Stability and bifurcation analysis of a discrete-time plant-herbivore model with harvesting effect

1.   Introduction

A plant-herbivore system is an ecological interaction in which plants function as main producers and herbivores consume them as their major source of living. The creation of this connection is vital for the efficient operation of the ecosystem and the conservation of biodiversity. Plants use photosynthesis to capture solar energy and convert it into organic matter. They serve as the fundamental basis of food webs, providing herbivores and other trophic levels with the energy and nutrition they need to survive. Herbivores are living things that eat plant tissues such as leaves, stems, fruits, and roots. Plant-herbivore interactions have a profound influence on ecosystem processes and resources. It appears that variations in the number of plants may have a substantial impact on the cycling of nutrients, the structure of the soil, and the availability of habitat, which in turn can have an effect on the abundance and distribution of other organisms within the ecosystem. Herbivores, in addition to being prey for predators, have the ability to impact the interactions that occur between predators and prey, as well as the trophic cascades that occur between different groupings of animals. The interaction between plants and herbivores is a fascinating and nuanced element of ecology, characterized by complex, frequently nonlinear relationships that may result in threshold effects. Minor changes on small scales may have significant implications on bigger scales ^[1,2,3,4]. Plant-herbivore interactions are important, but their dynamics may be challenging to understand due to complex ecological relationships ^{[5,6,7,8,9,10,11]}.

Plant-herbivore interactions are often described through modified versions of predator-prey models due to their similar dynamics and ecological principles. Both types of models aim to represent the relationships between predators (also known as herbivores) and their prey (also known as plants) at different trophic levels within an ecosystem. Although plant-herbivore interactions and predator-prey interactions have certain similarities, such as population control and feeding relationships, there are key differences that arise from the unique behaviors of each kind of interaction. In the context of plants and herbivores, there is no occurrence of a predator (herbivore) chasing its victim (plant) ^[12]. Contrary to predators, herbivores do not engage in the aggressive pursuit or killing of plants. Herbivory, on the other hand, refers to consuming leaves, stems, fruits, or roots using specialized mouthparts or feeding structures to obtain nutrients from plant tissues.

Various studies have explored the qualitative behavior of plant-herbivore models, examining phenomena such as bistability, bifurcations, and chaos control. These investigations have provided valuable insights into the dynamic complexities of plant-herbivore interactions, contributing to our understanding of ecosystem dynamics. Researchers have investigated plant-herbivore interactions using both continuous-time differential equations and discrete-time difference equations. Kartal ^[13] investigated the dynamical behavior of a plant-herbivore system, including both differential and difference equations. Beso et al. ^[14] investigated stability and various types of bifurcations in a plat-herbivore system with a strong Allee effect. Din ^[15] investigated the global behavior of a discrete-time plant-herbivore system. Khan et al. ^[16] conducted bifurcation analysis of a discrete-time plant-herbivore system. Hamada ^[17] investigated stability, bifurcation, and chaos in a discrete-time plant-herbivore system obtained from a continuous-time plant-herbivore system by applying the piecewise constant argument method. Similarly, for some other discussions related to the population dynamics of plant-herbivore systems, we refer the interested reader to ^{[18,19,20,21,22,23]} and references therein.

The apple twig borer (ATB), a pest infesting grapevines, serves as a case study illustrating the ecological dynamics of plant-herbivore systems. Previous studies ^[24,25] proposed fundamental frameworks based on the life cycle of adult ATBs, laying the groundwork for understanding the dynamics of plant-herbivore interactions. Understanding these dynamics is crucial for land management, environmental protection, and animal husbandry, prompting extensive research into dynamical system models to elucidate population and ecosystem behaviors. In ^[11], Din et al. investigated the stability, Neimark-Sacker bifurcation, and chaos in the following discrete-time plant-herbivore system with Holling type-Ⅱ functional response:

where $x_n$ and $y_n$, respectively, represent the population densities of grapevine and ATB. Moreover, $r, \kappa{, }$ and $\sigma$ are positive constants.

In population dynamics, harvesting refers to the intentional removal of individuals from a population for human use or consumption. This practice is common in various natural resource management contexts, including fisheries ^[26,27], forestry ^[28,29], and wildlife management ^[30,31,32]. Harvesting may be detrimental to the health of ecosystems and biodiversity if it is carried out in an excessive or indiscriminate manner. Understanding the consequences of such efforts is necessary for the implementation of sustainable management systems. Within the framework of plant-herbivore interactions, the act of harvesting may directly influence the populations of both plants and herbivores, resulting in changes to the composition, abundance, and distribution of these populations. The harvesting of plants has the potential to reduce population densities, disrupt reproductive cycles, and alter the structure of communities by favoring some species over others. There is also the possibility that these effects might occur simultaneously. Similarly, herbivores may experience changes in population dynamics, resource availability, and foraging behavior as a result of the pressure exerted by harvesting. Furthermore, the removal of individuals from both the plant and herbivore populations has the potential to upset the delicate balance that exists between them, which may have a domino effect on higher trophic levels as well as the general functioning of the ecosystem. In the study conducted by Virtala ^[33], the optimum harvesting in a plant-herbivore system was investigated. The study conducted by Asfaw et al. ^[34] investigated a plant-herbivore system that included herbivore harvesting effects. They provided evidence that the harvest rate of the herbivore population has a significant influence on the dynamics of herbivores.

When it comes to effectively managing plant-herbivore systems, there are a few different harvesting strategies to choose from. These include proportional, constant, nonselective, and selective harvesting ^{[35,36,37,38]}. Each approach has its own set of repercussions, both for the coexistence of organisms and for the dynamics of the ecosystem. The term selective harvesting refers to a technique that focuses on certain individuals or species, taking into consideration specific characteristics such as size, age, or quality. This strategy makes it possible to exert conservation and management activities in a focused manner. In contrast, nonselective harvesting entails the extraction of individuals without taking into account their distinct features or characteristics. This might potentially have unanticipated consequences for the structure of populations and the operation of ecosystems. Constant harvesting is the consistent and regular extraction of a set quantity from a population, which consistently affects the population dynamics independent of variations in population density. In contrast, proportionate harvesting maintains a harvest rate that is directly proportional to the size of the present population, ensuring that the intensity of harvesting increases or decreases with the abundance of the population. This adaptable strategy enables flexible administration by shifting population dynamics and fostering sustainable use of resources and stability in the population. Proportionate harvesting may provide more adaptability to changes in the environment and variations in population size, especially in ecosystems that are dynamic or changing.

Therefore, inspired by the preceding discussion, we are interested to inquire: What are the consequences for the dynamic properties when a harvesting impact is applied to the plant population in system (1.1)? Therefore, we modified system (1.1) by integrating the harvesting impact into the plant population. The result is the following modified system:

where $h > 0$ is the harvesting rate. It represents the intensity of harvesting. Ensuring nonnegative solutions is of utmost importance in ecological models, such as the plant-herbivore model, to maintain biological relevance. Negative solutions are not ecologically meaningful since populations cannot have negative densities. Ensuring solution positivity is crucial for aligning model predictions with ecological reality, which in turn facilitates the correct understanding and management of ecosystems. We assume that the initial values $x_0$ and $y_0$ are positive. Clearly, from the second equation of system (1.2), it is evident that $y_{n+1} > 0$ for all $n$. In the first equation, the term $-h x_n$ can produce negative values of $x_{n+1}$. One should select $h$ sufficiently small so that $x_{n+1}$ is positive. It can be checked that if $h$ is sufficiently small enough and $h < \frac{1}{r(1+y_0)+\kappa x_0}$, then $x_{n+1} > 0$ for all $n$. Another approach, given in Section 3.6 of ^[39], is that $x_{n+1}$ can be redefined as follows:

Moreover, in all our simulations performed here, solutions remained nonnegative. Through a combination of analytical techniques and numerical simulations, we seek to address the following research questions:

● What are the effects of plant harvesting on the fixed points of the plant-herbivore system (1.2)?

● How does plant harvesting influence the stability properties of fixed points, and under what conditions do bifurcation phenomena occur?

● What insights can be gained from bifurcation analysis regarding the long-term dynamics of the modified plant-herbivore system?

For the detailed analysis of stability and bifurcation in discrete-time systems, we refer the readers to ^{[40,41,42,43,44]} and references therein. The subsequent sections of the paper are arranged as follows: Section 2 is dedicated to investigating the existence and stability of fixed points. Section 3 examines the study of bifurcations that include period-doubling (PD), transcritical (TC), and Neimark-Sacker (NS) occurring at the positive fixed point (FP). In Section 4, numerical simulation results are presented to support the theoretical analysis and display the new and rich dynamic behavior. Moreover, the influence of harvesting on system dynamics is presented in Section 5. Finally, a brief conclusion is presented in Section 6.

2.   Existence and stability of FPs

Exploring the stability of FPs is very important in plant-herbivore environments. These FPs are states of equilibrium where both plant and herbivore populations are balanced. By examining their stability, we can predict long-term behavior in ecological systems, improving our understanding of the multiple components that affect ecosystem dynamics.

2.1.   Existence of FPs

The FPs of system (1.2) can be determined by solving the subsequent nonlinear equations:

for $x$ and $y$. The system (1.2) possesses three FPs: $E_0 = (0, 0), E_1 = (\frac{1-r-hr}{\kappa+\kappa h}, 0)$, and

$E_0 = (0, 0)$ is the trivial FP that exists always. $E_1 = (\frac{1-r-hr}{\kappa+\kappa h}, 0)$ is a boundary FP that exists if $r+rh < 1$. Moreover, $E_2 = \bigg(\frac{1-\sigma}{\sigma}, \frac{1}{r}\bigg(\kappa-\frac{\kappa}{\sigma}+\frac{1}{1+h}-r\bigg)\bigg)$ is coexistence FP, which exists if $r < 1, \ \frac{\kappa}{1+\kappa-r} < \sigma < 1$, and $h < \frac{\sigma}{\kappa-\kappa \sigma+\sigma r}-1$.

2.2.   Stability of FPs

The Jacobian matrix of the system

is the matrix given below:

Thus, the Jacobian matrix $J(x, y)$ of the system (1.2) evaluated at any FP $(x, y)$ is as follows:

The eigenvalues $\xi_{1, 2}$ of the Jacobian matrix $J$ are helpful in determining the stability of FPs. The FP $(x, y)$ is called a sink if $|\xi_{1, 2}| < 1$ and a source if $|\xi_1| > 1$ and $|\xi_2| > 1$. Furthermore, the FP $(x, y)$ is classified as a saddle point (SP) if $|\xi_1| > 1$ and $|\xi_2| < 1$ (or $|\xi_1| < 1$ and $|\xi_2| > 1$). In the case of a non-hyperbolic point (NHBP) $(x, y)$, either $|\xi_1| = 1$ or $|\xi_2| = 1$. However, if the eigenvalues are in complicated form, making direct use of this definition challenging, the following result proves to be very useful in determining the stability of the FP.

Lemma 2.1. ^[45] Let $\Lambda(\xi) = \xi^2+K_1\xi+K_0$. Assume that $\Lambda(1) > 0$. If $\xi_1$ and $\xi_2$ are solutions of $\Lambda(\xi) = 0$, then

(1) $|\xi_1| < 1$ together with $|\xi_2| < 1$ if $\Lambda(-1) > 0 \wedge K_0 < 1$,

(2) $|\xi_1| < 1 \wedge |\xi_2| > 1$ (or $|\xi_1| > 1 \wedge |\xi_2| < 1$) if $\Lambda(-1) < 0$,

(3) $|\xi_{1, 2}| > 1$ if $\Lambda(-1) > 0 \wedge K_0 > 1$,

(4) $|\xi_2|\neq 1 \wedge \xi_1 = -1$ if $\Lambda(-1) = 0 \wedge K_1 \neq 0, 2$,

(5) $\xi_1, \xi_2 \in \mathbb{C}$ along with $|\xi_{1, 2}| = 1$ if $K_1^2-4K_0 < 0 \wedge K_0 = 1$.

A sink reflects the long-term coexistence of herbivores and plants, suggesting a balanced environment in which populations maintain sustainable levels throughout time. An SP suggests intermittent stability, with the system alternating between plenty and shortage, similar to natural population fluctuations. An unstable source suggests vulnerability to unexpected changes, stressing the potential for population collapses or rapid expansion that disrupts ecological equilibrium. Furthermore, non-hyperbolic points may indicate intricate and nuanced interactions, demonstrating the complexity of ecological dynamics and the difficulty of predicting population behavior. Through computations, it is obtained that

One can easily obtain the following result for the topological classification of $E_0$:

Theorem 2.2. $E_0$ is a

(1) sink if $\sigma < 1$ and $\frac{1-r}{r} < h < \frac{1+r}{r}$,

(2) source if one of the below-listed parametric conditions is fulfilled:

(i) $\sigma > 1, r < 1$, and $h < \frac{1-r}{r}$,

(ii) $\sigma > 1$ and $h > \frac{1+r}{r}$,

(3) SP if one of the listed below parametric conditions is fulfilled:

(i) $\sigma < 1, r < 1$, and $h < \frac{1-r}{r}$,

(ii) $\sigma < 1$ and $h > \frac{1+r}{r}$,

(iii) $\sigma > 1$ and $\frac{1-r}{r} < h < \frac{1+r}{r}$,

(4) NHBP if one of the listed below parametric conditions is fulfilled:

(i) $\sigma = 1$,

(ii) $h = \frac{1+r}{r}$,

(iii) $r < 1$, $h = \frac{1-r}{r}$.

The eigenvectors corresponding to eigenvalues $\xi_1 = -h+\frac{1}{r}$ and $\xi_2 = \sigma$ of $J(E_0)$ are $v_1 = < 1, 0 >$ and $v_2 = < 0, 1 >$. The stable $W^s$ and unstable $W^u$ manifolds in the saddle case (3-ⅲ) are given by

Remark 2.3. At NHBP, the system experiences different types of bifurcations depending on the nature of the eigenvalues of the Jacobian matrix. If (4-ⅰ) is true, it follows that one of the eigenvalues of the matrix $J(E_0)$ is $1$. Consequently, a TC bifurcation occurs at $E_0$. Similarly, a TC bifurcation occurs at $E_0$ if (4-ⅲ) is true. Moreover, if (4-ⅱ) is true, it follows that one of the eigenvalues of the matrix $J(E_0)$ is $-1$. Consequently, a PD bifurcation occurs at $E_0$.

Next, we obtain

Since $E_1$ exists if $r < \frac{1}{1+h}$. Thus, we obtain the following result:

Theorem 2.4. $E_1$ is a

(1) sink if $\frac{h-1}{(h+1)^2} < r < \frac{1}{h+1}$ and $0 < \sigma < \frac{\kappa+\kappa h}{\kappa+\kappa h-r-hr+1}$,

(2) source if $r < \frac{h-1}{(h+1)^2}, h > 1$ and $\sigma > \frac{\kappa+\kappa h}{\kappa+\kappa h-r-hr+1}$,

(3) SP if any of the following criteria is met:

(i) $r < \frac{h-1}{(h+1)^2}, h > 1$, and $0 < \sigma < \frac{\kappa+\kappa h}{\kappa+\kappa h-r-hr+1}$,

(ii) $\sigma > \frac{\kappa+\kappa h}{\kappa+\kappa h-r-hr+1}$ and $\frac{h-1}{(h+1)^2} < r < \frac{1}{h+1}$,

(4) NHBP if one of the listed below criteria is met:

(i) $r = \frac{h-1}{(h+1)^2}$,

(ii) $\sigma = \frac{\kappa+\kappa h}{\kappa+\kappa h-r-hr+1}$.

The eigenvectors corresponding to eigenvalues $\xi_1 = r+h^2 r+h (-1+2 r)$ and $\xi_2 = \sigma -\frac{\sigma (-1+r+h r)}{\kappa (1+h)}$ are $v_1 = < 1, 0 >$ and

The stable $W^s$ and unstable $W^u$ manifolds in the saddle case (3-ⅱ) are given by

Remark 2.5. If (4-i) is true, then it follows that one of the eigenvalues of the matrix $J(E_1)$ is $-1$. Consequently, a PD bifurcation occurs at $E_1$. Moreover, if (4-ii) is true, it follows that one of the eigenvalues of the matrix $J(E_1)$ is $1$. Consequently, a TC bifurcation occurs at $E_1$.

The Jacobian matrix at positive FP is

The characteristic polynomial of $J(E_2)$ is

where

Through computations, it is obtained that

After applying Lemma 2.1, we obtain the following result for the classification of the positive FP, $E_2$, of system (1.2).

Theorem 2.6. The positive FP

(1) $E_2$ is a sink if $r < \frac{1}{\sigma (1+h)}(\sigma +\kappa \sigma +\kappa \sigma h-2\kappa-2h)$ and

(2) $E_2$ is a source if $r > \frac{1}{\sigma (1+h)}(\sigma +\kappa \sigma +\kappa \sigma h-2\kappa-2h)$ and

(3) $E_2$ is an SP if $r < \frac{1}{(\sigma -1)(1+h)^2} \bigg(-5+\sigma -h+\sigma h+\frac{\kappa(\sigma -3)(\sigma -1)(1+h)^2}{\sigma }\bigg)$,

(4) $E_2$ is NHBP and experiences PD bifurcation if $\sigma \neq \frac{\kappa(1+h)^2}{2+\kappa(1+h)^2}, \frac{\kappa(1+h)^2}{4+\kappa(1+h)^2}$ and

(5) $E_2$ is NHBP and experiences NS bifurcation if $0 < \kappa < \frac{4}{(1-\sigma)(1+h)^2}$ and

3.   Bifurcation analysis

This section aims to thoroughly investigate bifurcation phenomena, including TC bifurcation at $E_1$ as well as PD and NS bifurcations at $E_2$. For a comprehensive understanding of bifurcation analysis, we suggest referring to the sources ^{[46,47,48,49]}. Bifurcation is a crucial factor in determining the dynamics of the system, revealing situations where even little changes in parameters may lead to significant changes in plant-herbivore interactions. Moreover, exploring bifurcation processes improves our comprehension of ecosystem dynamics, facilitating the development of effective strategies for preserving and controlling plant and herbivore populations in a sustainable manner.

3.1.   TC bifurcation at $E_1$

We start our examination of the TC bifurcation at $E_1$ by considering condition (4-ⅱ) outlined in Theorem 2.4. We can rewrite it as $r = \frac{-\kappa -h \kappa +\sigma +\kappa \sigma +h \kappa \sigma }{(1+h) \sigma }$. Introducing a sufficiently small perturbation $\varepsilon$ into the bifurcation parameter $r$ around $R_{10} = \frac{-\kappa -h \kappa +\sigma +\kappa \sigma +h \kappa \sigma }{(1+h) \sigma }$, the system (1.2) takes the subsequent form:

We shift $E_1$ to $(0, 0)$ by taking $u_n = x_n-\frac{1-(R_{10}+\varepsilon)-h(R_{10}+\varepsilon)}{\kappa+\kappa h}, \ v_n = y_n$. Consequently, the system (3.1) becomes

where

Next, the system (3.2) is diagonalized through the consideration of the following transformation:

Upon applying the mapping (3.3), the system (3.2) undergoes the alteration as follows:

where

Now, we use the center manifold theory to obtain the center manifold $M^C$ of (3.2) at $(0, 0)$ in a close neighborhood of $\varepsilon = 0$. The center manifold $M^C$ can be obtained using

Through calculations, we obtain $p_1 = \frac{d_1}{1-\xi }, p_2 = \frac{d_4}{1-\xi }, p_3 = 0$. Consequently, the system (3.2) is restricted to $M^C$ as follows:

Through simple computations, we obtain

Thus, the system (1.2) experiences TC bifurcation at $E_1$. The next result gives the conditions for the existence and direction of TC bifurcation in the system (1.2) at $E_1$.

Theorem 3.1. If condition (4-ii) in Theorem 2.4 is fulfilled, then (1.2) undergoes TC bifurcation at $E_1$ when $r$ differs in a close neighborhood of $R_{10} = \frac{-\kappa -h \kappa +\sigma +\kappa \sigma +h \kappa \sigma }{(1+h) \sigma }$.

The boundary FP is unstable for $r < R_{10}$. It collides with positive FP at $r = R_{10}$, and subsequently they exchange stability. The result shows that small variations in environmental or biological parameters may cause major modifications in plant-herbivore population dynamics, culminating in a TC bifurcation at the border fixed point. At the transition, essential ecosystem dynamics shift, causing the system to shift from a state of stable coexistence between plants and herbivores to one in which plants can live on their own.

3.2.   PD bifurcation at $E_2$

Next, we investigate the PD bifurcation at $E_2$ by considering condition (4) presented in Theorem 2.6. Biologically, PD bifurcations indicate significant alterations in plant-herbivore dynamics. Consider a situation in which plant and herbivore populations vary regularly over time. During a PD bifurcation, these variations become more intense, and the space between peak and trough densities widens, perhaps doubling. This may have a severe environmental effect. If environmental circumstances unexpectedly favor plant growth, herbivore populations may increase, thereby depleting plant resources. Plant populations rebound when herbivore numbers drop. However, when plant numbers recover, herbivore populations rebound, resulting in increased grazing pressure and a consequent drop in plant populations. Introducing a sufficiently small perturbation $\varepsilon$ into the bifurcation parameter $r$ around $R_1 = \frac{1}{(\sigma -1)(1+h)^2} \bigg(-5+\sigma -h+\sigma h+\frac{\kappa(\sigma -3)(\sigma -1)(1+h)^2}{\sigma }\bigg)$, the system (1.2) takes the subsequent form:

To translate the point $E_2$ at $(0, 0)$, we use the change of variables $u_n = x_n-\frac{1-\sigma }{\sigma }, \ v_n = y_n-\frac{1}{(R_1+\varepsilon)}\bigg(\kappa-\frac{\kappa}{\sigma }+\frac{1}{1+h}-(R_1+\varepsilon)\bigg)$. Consequently, the system (3.6) is expressed as below:

where

Subsequently, the system (3.7) goes through diagonalization after a consideration of the following transformation:

After implementing the map (3.8), the system (3.7) takes the subsequent form:

where

Now, we use the center manifold theory to obtain the center manifold $M^C$ of (3.9) at $(0, 0)$ in a close neighborhood of $\varepsilon = 0$. Its approximation is as follows:

where

Consequently, the system (3.9) is restricted to $M^C$ as follows:

The mapping (3.10) undergoes PD bifurcation if the following are nonzero:

Thus, we obtain the following result:

Theorem 3.2. If condition (4) of Theorem 2.6 is fulfilled, then (1.2) undergoes PD bifurcation at $E_2$ if $l_1, l_2$ provided in (3.11) and (3.12) are nonzero and $r$ fluctuates in a small neighborhood of

Moreover, if $\ l_2 > 0$ (respectively, $l_2 < 0$), then an orbit of period-2 emanates from $E_2$, which is stable (respectively, unstable).

The result illustrates how minor fluctuations in parameters can lead to substantial shifts in the system's dynamics, specifically causing a doubling of the population oscillation periods. This happens when the herbivore population increases and feeds more on the plants, thus leading to a reduction in the plant population. With fewer plants available, the herbivore population begins to decrease due to the lack of food. When the herbivore numbers decrease, the plant population is given the chance to increase once again, and the cycle begins again. This cycle reflects the delicate balance and complex interactions that exist within the environment.

3.3.   NS bifurcation at $E_2$

Next, we investigate NS bifurcation at $E_2$ by considering condition (5) presented in Theorem 3.3. Biologically, NS bifurcations reveal complex dynamics within a plant-herbivore system, often marked by unexpected changes and unpredictable patterns. Imagine a situation in which the changes in the populations of plants and herbivores cannot be accurately predicted. NS bifurcation may cause population numbers to fluctuate significantly due to even minor alterations in environmental circumstances or plant-herbivore interactions, making it difficult to predict these changes. The instability has a dramatic effect on the dynamics of ecosystems, which presents issues for ecologists attempting to forecast how populations will respond to changes in the environment. For example, sudden changes in the environment may cause unpredictable and disorderly variations in populations of plants and herbivores, leading to disruptions in the interactions between different trophic levels and unpredictable changes in how the ecosystem functions. Introducing a sufficiently small perturbation $\varepsilon$ into the bifurcation parameter $r$ around $R_2 = \frac{1}{\sigma (1+h)}(\sigma +\kappa \sigma +\kappa \sigma h-2\kappa-2\kappa h)$, the system (1.2) takes the subsequent form:

We use the translation $u_n = x_n-\frac{1-\sigma }{\sigma }, \ v_n = y_n-\frac{1}{(R_2+\varepsilon)}\bigg(\kappa-\frac{\kappa}{\sigma }+\frac{1}{1+h}-(R_2+\varepsilon)\bigg)$ to translate the point $E_2$ to origin. Consequently, the system (3.13) becomes

where

The characteristic equation of the linearized part of (3.14) computed at $(0, 0)$ is

where

The solutions of (3.15) are given by

We then obtain $|\xi_{1, 2}| = \sqrt{\alpha_2(\varepsilon)}$. Since $\sigma < 1$, then

Moreover, it is necessary that $\xi_{1, 2}^i \neq 1$ when $\varepsilon = 0$ for $i = 1, 2, 3, 4,$ which is equivalent to $\alpha_1(0) \neq -2, 0, 1, 2$. Through simple computations, we obtain

Clearly, $\alpha_1(0)\neq 2$. Thus, we only require that

To determine the canonical form of (3.14) at $\varepsilon = 0$, we employ the map

After applying the map (3.18), the system (3.14) converts to

where

In order to investigate the direction of the NS bifurcation, we consider the first Lyapunov exponent, which is derived as follows:

where

Based on the aforementioned study, the following conclusion can be drawn:

Theorem 3.3. Assume that condition (5) of Theorem 2.6 is fulfilled. If condition (3.17) is satisfied and $L$ presented in (3.20) is nonzero, then (1.2) undergoes NS bifurcation at $E_2$ when the parameter $r$ varies in a small neighborhood of $R_2 = \frac{1}{\sigma (1+h)}(\sigma +\kappa \sigma +\kappa \sigma h-2\kappa-2\kappa h)$. Furthermore, if $L < 0$ (or $L > 0$), the NS bifurcation at $E_2$ is classified as supercritical (or subcritical), and a single closed invariant curve bifurcates from $E_2$, which is attracting (or repelling).

The result shows that small changes in environmental or biological conditions may generate significant changes in plant-herbivore population dynamics, resulting in an NS bifurcation. This indicates that the system shifts from steady-state oscillations to quasiperiodic oscillations, in which populations no longer follow simple cycles but instead show increasingly complicated and interconnected patterns.

4.   Numerical examples

The purpose of this section is to verify our theoretical results through numerical simulations. Calculations are done using MATHEMATICA, while MATLAB is used for plotting graphs.

4.1.   Stability and bifurcation analysis varying $r$

Assume that $h = 0.15, \kappa = 0.01, \sigma = 0.04$ and $x_0 = 23.4, y_0 = 0.4$. For these parametric values, the boundary FP, $E_1$, is obtained as $(24, 0)$. The system (1.2) experiences TC bifurcation at $E_1$ when $r \approx 0.629565$. The eigenvalues of $J(E_1)$ are $\xi_1 = 1$ and $\xi_2 = 0.6826$. Moreover, some careful calculations give

It verifies Theorem 3.1. We present phase portraits in Figure 1 to graphically verify the TC bifurcation at $E_1$. Figure 1a illustrates that $E_1$ is unstable while $E_2$ is stable. In Figure 1b, it is shown that $E_1$ and $E_2$ collide at the critical value $r = 0.629565$ and, subsequently, they exchange stability. By choosing $r = 0.65$, Figure 1c demonstrates that $E_1$ becomes stable while $E_2$ turns unstable. This phenomenon is known as TC bifurcation.

Our modified system (1.2) is a generalization of the system (1.1) studied in ^[11]. We can obtain the original system (1.1) by substituting $h = 0$ into the modified system (1.2). In this analysis, we consider the same parameter values $\kappa = 0.01, \sigma = 0.04$ and initial conditions $x_0 = 23.4, y_0 = 0.4$, which were used in ^[11]. Moreover, we assume that $h = 0.15$. We vary $r$ in a bigger interval $[0.2, 0.9]$ instead of $[0.2, 0.7]$ to address additional important dynamics that are not reported in the numerical analysis presented in ^[11]. Our investigation shows that harvesting is helpful in stabilizing the plant-herbivore system. The system (1.2) undergoes NS bifurcation for $r \approx 0.379565$ at FP $E_2 = (24, 0.658648)$. The eigenvalues of $J(E_2)$ are $\xi_{1, 2} = 0.8413\pm 0.540569 i$ with $|\xi_{1, 2}| = 1$. Moreover, some careful calculations give

Thus, it is obtained that $L = -0.016401 < 0$, which verifies Theorem 3.3. Thus, NS bifurcation is supercritical, and a closed invariant attracting curve will emerge from $E_2$. Bifurcation diagrams are depicted in Figures 2a and 2b. These illustrate that $r$ stabilizes the system (1.2). The corresponding maximum Lyapunov exponent (MLE) graph is given in Figure 2c.

Next, Figure 3 shows phase portraits of (1.2) for different values of $r$. It can be seen that $E_2$ is a sink for $r > 0.379565$, yet it becomes unstable at $r \approx 0.379565$ due to the occurrence of NS bifurcation. For $r \leq 0.379565$, an invariant curve emanates from $E_2$, with its radius expanding as $r$ is decreasing.

In ^[11], for the same parametric values and initial conditions, the positive FP is a sink for $0.51 < r < 0.76$. Their investigation shows that the system experiences NS bifurcation at $r \approx 0.51$. Although it is not highlighted in their paper, positive FP vanishes at $r \approx 0.76$ and a boundary FP appears due to TC bifurcation at the boundary FP. Later, this boundary FP also disappears for $r \geq 1$ and the system has only a trivial FP, $E_0$, which is then stable for $r > 1$. The same phenomenon is observed in our modified system (1.2). By adding the harvesting effect, the positive FP of the modified system (1.2) is a sink for $0.379565 < r < 0.629565$. The system experiences NS bifurcation at positive FP, $E_2$, for $r \approx 0.379565$. The positive FP $E_2$ disappears for $r \geq 0.629565$. The system (1.2) is experiencing TC bifurcation at boundary FP $E_1$ for $r \approx 0.629565$. Before $r \approx 0.629565$, $E_2$ is the only nontrivial FP, which is a sink for $0.379565 < r < 0.629565$. After $r \approx 0.629565$, $E_1$ is the only nontrivial FP, which is a sink for $0.629565 < r < 0.869565$. Later, the boundary FP $E_1$ also disappears for $r \geq 0.869565$. The system (1.2) is experiencing TC bifurcation at trivial FP $E_0$ for $r \approx 0.869565$. For $0.629565 < r < 0.869565$, the system possesses one boundary FP $E_1$ and one trivial FP $E_0$. Here, $E_1$ is stable, but $E_0$ is unstable. For $r > 0.869565$, the boundary FP, $E_1$, disappears, but the trivial FP $E_0$ exists, and it also becomes stable.

4.2.   Stability and bifurcation analysis varying $h$

In the previous section, it was shown that the harvesting effect has a stabilizing impact on plant-herbivore systems. One natural question arises: What will happen if we harvest without any set upper and lower limits? Will it still stabilize the system or destabilize if we significantly increase or decrease the harvesting limit? Ecologically, threshold harvesting is beneficial for both plants and herbivores. In this numerical analysis, it is shown that a balanced amount of harvesting is advantageous for both plants and herbivores. For this purpose, assume that $r = 0.022, \kappa = 0.01, \sigma = 0.5$, and vary $h$. For these parametric values, the NS bifurcation value is $h_1 = 18.2308$ and the PD bifurcation value is $h_2 = 18.5673$. The FP $E_2$ is a sink if $18.2308 < h < 18.5673$. The FP $E_2$ becomes unstable for $h \leq 18.2308$ due to NS bifurcation. Moreover, it is unstable for $h \geq 18.5673$ due to NS bifurcation.

The positive FP is obtained as $E_2 = (1, 0.909091)$ for $h = h_1$. The eigenvalues of $J(E_2)$ are $\xi_{1, 2} = -0.849112\pm 0.528212i$ with $|\xi_{1, 2}| = 1$. Moreover, some careful calculations give

Thus, it is obtained that $L = -93.568 < 0$, which verifies Theorem 3.3. Thus, NS bifurcation is supercritical, and a closed invariant attracting curve will emerge from $E_2$. The bifurcation diagrams of system (1.2) are given in Figures 4a and 4b by using initial conditions $x_0 = 1.1, y_0 = 0.9$, and varying $h$ in $[18.20, 18.30]$. The corresponding MLE graph is given in Figure 4e. Next, the positive FP is obtained as $E_2 = (1, 0.86844)$ for $h = h_2$. The eigenvalues of $J(E_2)$ are $\xi_1 = -1$ and $\xi_2 = -0.828791$. Moreover, some careful calculations give $l_1 = -10.6168\neq 0$ and $l_2 = 6.90668 \neq 0$, which verifies Theorem 3.2. Since $l_2 > 0$, a stable orbit of period-2 emanates from $E_2$. The bifurcation diagrams are given in Figures 4c and 4d by using initial conditions $x_0 = 1.1, y_0 = 0.85$, and varying $h$ in $[18.50, 18.65]$. The corresponding MLE graph is given in Figure 4f. Figures depict phase portraits of (1.2) for different values of $r$. These illustrate that $E_2$ is a sink if $18.2308 < h < 18.5673$. The system undergoes NS bifurcation at $h \approx 18.2308$ and PD bifurcation at $h \approx 18.5673$.

Figures 4 and 5 illustrate that the higher or lower harvesting effect destabilizes the system, resulting in more complicated dynamical behaviors. Thus, a balanced amount of harvesting is advantageous for both plants and herbivores. Ecologically, under moderate harvesting, there is a stable coexistence of both plants and herbivores, represented by the positive FP. As harvesting increases, herbivores struggle to find food, eventually leading to extinction. Continued harvesting eventually causes the plant population to decline as well, resulting in the extinction of plants. This indicates that moderate harvesting can sustain both populations, but excessive harvesting leads to the collapse of the entire ecosystem.

5.   Influence of harvesting effect

The aim of this section is to investigate the influence of the harvesting effect on the dynamics of system (1.2). Incorporating a harvesting effect into a plant-herbivore system yields insightful findings that underscore the intricate dynamics of species interactions in ecological systems. While the harvesting effect is directly applied to the plant population, its influence extends beyond that, indirectly impacting herbivore dynamics through the interaction between plant and herbivore populations. One can observe in Figure 4 that system (1.2) undergoes PD and NS bifurcations at the positive FP $E_2$ by varying the harvesting effect parameter $h$. These show that the harvesting effect can stabilize or destabilize the system. Notably, a moderate level of the harvesting effect appears to be advantageous for both plant and herbivore populations, as evidenced by the occurrence of these bifurcations. Furthermore, the positive FP $E_2 = (x_2, y_2)$ is dependent on the harvesting effect parameter $h$. Through simple calculations, one can obtain that $\frac{dx_2}{dh} = 0$ and $\frac{dy_2}{dh} = -\frac{1}{r(1+h)^2} < 0$. Ecologically, it represents that when the harvesting effect strengthens in the plant population, it might become slightly more difficult for herbivores to survive. A decrease in plant availability could lead to reduced reproduction rates or increased mortality among herbivores, ultimately resulting in a decrease in the population growth rate of herbivores.

6.   Conclusions

This study explores the implications of adding a harvesting impact to the dynamics of a plant-herbivore system. Our findings highlight the intricate interplay between species interactions in ecological systems, emphasizing the significance of the harvesting effect. The harvesting effect affects the plant population directly, but it also has an impact on herbivore dynamics because of the complex interactions that occur between herbivores and plant populations. The existence and stability of FPs are examined. Furthermore, a comprehensive analysis of local bifurcations at the positive FP is conducted. The investigation reveals that the system (1.2) undergoes PD, TC, and NS bifurcations. These bifurcations demonstrate the system's sensitivity to small parameter changes, which result in significant changes in population dynamics. The role of transcritical bifurcation is to emphasize the critical transition from coexistence to dominance of plant populations; PD bifurcations show the formation of complex oscillatory patterns. NS bifurcations, in turn, indicate the transition from steady-state to quasiperiodic oscillations, displaying the complexity and adaptability of ecological interactions.

The plant-herbivore system without harvesting effect was investigated in ^[11]. Our plant-herbivore system is more general, and we can obtain the results of ^[11] by just replacing $h = 0$ into our results. Although we obtain similar results as in ^[11], the following are the main differences observed in both plant-herbivore systems:

● With the same parametric values and initial conditions, the positive FP of system (1.1) is a sink inside the range $0.51 < r < 0.76$. However, when the harvesting effect is included, the positive FP of system (1.2) is a sink inside the range $0.379565 < r < 0.629565$ for $h = 0.15$. This illustrates that harvesting helps to stabilize the plant-herbivore system for lower values of $r$.

● In the absence of harvesting, the system (1.1) does not experience PD bifurcation. However, when the harvesting is added, the system (1.2) experiences PD bifurcation at the positive FP.

An interesting finding in our paper is that a moderate amount of the harvesting effect benefits both plant and herbivore populations, emphasizing the intricate impact of harvesting on ecosystem dynamics. These results have significant biological consequences, as they reveal the intricate relationship between human involvement and the capacity of ecosystems to recover. They also highlight the need to adopt sustainable harvesting methods to safeguard the long-term well-being and stability of ecological systems. In the future, one can try to investigate the global dynamics and codimension-two bifurcations for the plant-herbivore system (1.2).

Author contributions

Mohammed Alsubhi: Formal Analysis, Validation, Writing-original draft, Methodology; Rizwan Ahmed: Conceptualization, Investigation, Writing-review & editing, Software, Supervision; Ibrahim Alraddadi: Conceptualization, Validation, Writing-review & editing, Visualization, Supervision, Methodology; Faisal Alsharif: Formal Analysis, Software, Visualization, Methodology; Muhammad Imran: Formal Analysis, Investigation, Writing-original draft, Software. All authors have read and approved the final version of the manuscript for publication.

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 extend their appreciation to the Deputyship for Research & Innovation, Ministry of Education in Saudi Arabia for funding this research work through the project number (445-9-923).

Conflict of interest

All authors declare no conflicts of interest in this paper.