Bifurcation Analysis of a Discrete Basener–Ross Population Model: Exploring Multiple Scenarios

Abstract

The Basener and Ross mathematical model is widely recognized for its ability to characterize the interaction between the population dynamics and resource utilization of Easter Island. In this study, we develop and investigate a discrete-time version of the Basener and Ross model. First, the existence and the stability of fixed points for the present model are investigated. Next, we investigated various bifurcation scenarios by establishing criteria for the occurrence of different types of codimension-one bifurcations, including flip and Neimark–Sacker bifurcations. These criteria are derived using the center manifold theorem and bifurcation theory. Furthermore, we demonstrated the existence of codimension-two bifurcations characterized by 1:2, 1:3, and 1:4 resonances, emphasizing the model’s complex dynamical structure. Numerical simulations are employed to validate and illustrate the theoretical predictions. Finally, through bifurcation diagrams, maximal Lyapunov exponents, and phase portraits, we uncover a diversity of dynamical characteristics, including limit cycles, periodic solutions, and chaotic attractors.

1. Introduction

Recently, discrete-time models have played an effective and important role in the study of population dynamics, particularly in species that have one generation per year, such as annual plants, and species with pulsed reproductions, such as many wildlife species in seasonal environments [1,2,3]. Moreover, they exhibit more intricate dynamic patterns than those of continuous-time models, such as bifurcations and chaos [4,5,6,7].

The discussion surrounding societal (civilization) collapse as a potential future trajectory for the world has recently gained significant relevance [8]. Tonnelier employed bifurcation analysis to investigate the effects of the nature depletion rate and societal inequality on the dynamics of the HANDY model [9]. To emphasize the importance of social cohesion, Sargentis et al. employed a mathematical model based on Hurst–Kolmogorov dynamics. They suggested that social structural changes are more likely to be responsible for society’s ultimate decline [10]. Akhavan and York studied population collapse in Elite-Dominated societies [11]. They used a new Lyapunov function theorem to analyze a class of ordinary differential equation models characterized by the perpetual collapse of the worker population. Subsequently, they employed qualitative criteria, indicating that these conditions necessitate population collapse. In their study, Manzoor et al. investigated the consequences of aggregation through densely connected communities on a networked system of natural-resource consumption [12]. To investigate the association between drought scenarios and population dynamics, Kuil et al. employed a hydrology–demography model that had been calibrated to mimic realistic feedbacks for the Ancient Maya population drop in Central America [13]. Robert et al. employed a Bayesian model-driven approach to deduce the timing of events in relation to arguments about Rapa Nui collapse [14]. Michel introduced a class of models designed to explore the dynamics of population and resources, where birth decisions are guided by a rational, forward-looking evaluation of the anticipated future state of resources [15]. Roman et al. presented a dynamical model that incorporates land state, population dynamics, and workforce engaged in both swidden and intensive agriculture as well as monument construction to examine whether the interplay between societal dynamics and natural resource depletion sheds light on the decline of the Classic Maya civilization [16]. Several scholars have argued for various types of models that represented the relationship between consumers and their resources; regarding this, see [17] and the many references cited therein.

The mathematical modeling of the relationship between a population and natural resources consumption is based on the Lotka–Volterra predator–prey model, in which the population is the predator and resources are the prey. An isolated island, with its regulated environment and lack of external disruptions like migration, serves as a natural laboratory for examining ecosystem evolution. This isolation enables researchers to examine ecological dynamics and evolutionary processes devoid of external influences, providing significant insights into the development and adaptation of ecosystems across time. In 1998, Brander and Taylor developed an economic model to explore the environmental and societal collapse of the Polynesian civilization on Easter Island [18]. Theirs is a predator–prey-driven model, which is governed by the following equations:

$$\begin{array}{cc}\hfill {\displaystyle \frac{dP}{dt}}& =G(H,R)P,\hfill \\ \hfill {\displaystyle \frac{dR}{dt}}& =\rho \left(R\right)-H(R,P),\hfill \end{array}$$

(1)

where $P\left(t\right)$ and $R\left(t\right)$ are the human population and the stock of renewable resources, respectively; $G(H,R)$ represents the population growth rate; $\rho \left(R\right)$ denotes the renewal rate of resource; and $H=H(R,P)$ is the rate of resource harvesting. They have described resource dynamics through a logistic regeneration and a linear correlation between fertility and resource consumption as follows:

$$\begin{array}{cc}\hfill {\displaystyle \frac{dP}{dt}}& =(b-d+\varphi R)P,\hfill \\ \hfill {\displaystyle \frac{dR}{dt}}& =cR\left(1-{\displaystyle \frac{R}{K}}\right)-hPR.\hfill \end{array}$$

(2)

Here, b and d are birth and death rates, respectively; $\varphi$ is the conservation rate per capita resource consumption; c is the intrinsic regeneration rate of the resource; K is carrying capacity; and h is the rate of resource harvesting. Basener and Ross [19] employed the logistic growth to describe population dynamics as follows:

$$\begin{array}{cc}\hfill {\displaystyle \frac{dP}{dt}}& =aP\left(1-{\displaystyle \frac{P}{R}}\right),\hfill \\ \hfill {\displaystyle \frac{dR}{dt}}& =cR\left(1-{\displaystyle \frac{R}{K}}\right)-hP,\hfill \end{array}$$

(3)

where a is the population natural growth rate.

In this work, using the forward Euler approach, we obtain the discrete-time type of Model (3) as follows:

$$\begin{array}{cc}\hfill P_{n+1}& =P_n+a{P}_n\left(1-{\displaystyle \frac{P_n}{R_n}}\right),\hfill \\ \hfill R_{n+1}& =R_n+c{R}_n\left(1-{\displaystyle \frac{R_n}{K}}\right)-h{P}_n.\hfill \end{array}$$

(4)

Here, we consider the integral step size to be equal to 1. Our research conducts a comprehensive study of bifurcation analysis for the discrete-time Basener–Ross model (Equation (4)). According to Schaffer and Kot [20], it is critical to understand the periodic or chaotic dynamics that arise in ecological models. Their findings suggest that, far from being chaotic and disorderly, the chaotic trajectory structure could actually include crucial information about an ecosystem’s dynamics. The bifurcation of equilibria is a well-known cause of an ecosystem’s complicated dynamics. In [21], the authors numerically explored the bifurcation behaviors of codimension one for Model (4). Here, we examine analytically the codimension-one bifurcations, including flip and Neimark–Sacker bifurcations, and extend our study to codimension-two bifurcations characterized by $1:2$, $1:3$, and $1:4$ resonances, emphasizing the model’s complex dynamical structure. For analyzing codimension-one bifurcations, we use the center manifold theorem and bifurcation theory. This method is most effective for low-dimensional models. However, as the dimension of the model increases, the complexity of constructing and analyzing the center manifold grows significantly. Additionally, its utility is limited for higher-codimension bifurcations. Next, we examine the bifurcation behaviors of codimension two for the current model using the normal form approach and bifurcation theory. This approach does not need a transition into Jordan form and the computation of the model center manifold approximation. It is sufficient to compute the critical non-degeneracy coefficients to verify the existence of various bifurcation forms. Numerous studies have focused on bifurcation and chaotic behaviors in both discrete-time and continuous-time models. While numerous codimension-one bifurcation literature srouces have been considered, as shown in [4,5,6,7], only a small number of codimension-two bifurcation literature options are theoretically feasible [22,23,24,25]. To the best of our knowledge, there is very little literature on the topic of the bifurcation behaviors of the discrete-time Basener–Ross model as a function of two parameters that use both theoretical and numerical approaches, including continuation, invariant manifolds, maximal Lyapunov exponents, and normal forms. Not only that, numerical simulations are used to verify our theoretical results and describe other model behaviors like bifurcations of higher iterations (such as the third and fourth iterations).

The structure of this paper is outlined as follows. Section 2 focuses on examining the existence and stability of fixed points in Model (4). Section 3 presents the derivation of sufficient conditions for the occurrence of flip bifurcation and Neimark–Sacker bifurcation, utilizing the center manifold theorem and bifurcation theory. In Section 4, we discuss the codimension-two bifurcation. In Section 5, numerical simulations are conducted to validate the theoretical findings. Finally, concluding remarks are provided in Section 6.

2. Exploring the Existence and Stability of Fixed Points

This section explores the existence and stability of Model (4)’s fixed points. Model (4) exhibits two non-negative fixed points, which can be stated as follows:

• The semi-trivial fixed point $E_1=(0,K)$;
• The unique coexistence fixed point $E_2=(P^{*},R^{*})=(K-{\displaystyle \frac{hK}{c}},K-{\displaystyle \frac{hK}{c}})$, which exists when $c>h$.

In order to investigate the local stability of Model (4)’s fixed points, we examine the Jacobian matrix of Model (4) and compute the eigenvalues associated with each fixed point. At an arbitrary fixed point $(P,R)$, the Jacobian matrix is given by

$$\begin{array}{c}\hfill J(P,R)=\left(\begin{array}{cc}1+a-{\displaystyle \frac{2aP}{R}}& {\displaystyle \frac{a{P}^2}{R^2}}\\ -h& 1+c-{\displaystyle \frac{2cR}{K}}\end{array}\right).\end{array}$$

(5)

Now, we state the following lemmas that help us to understand the qualitative properties of the obtained fixed points.

Lemma 1 

([26]). Let $F\left(\lambda \right)={\lambda }^2+Tr\lambda +Det$. Suppose that $F\left(1\right)>0$ and that ${\lambda }_1$ and ${\lambda }_2$ are two roots of $F\left(\lambda \right)=0$. Then,

• $\left|{\lambda }_1\right|<1$ and $\left|{\lambda }_2\right|<1$ if and only if $F(-1)>0$, $Det<1$;
• $\left|{\lambda }_1\right|<1$ and $\left|{\lambda }_2\right|>1$ (or $\left|{\lambda }_1\right|>1$ and $\left|{\lambda }_2\right|<1$) if and only if $F(-1)<0$;
• $\left|{\lambda }_1\right|>1$ and $\left|{\lambda }_2\right|>1$ if and only if $F(-1)>0$ and $Det>1$;
• ${\lambda }_1=-1$ and $\left|{\lambda }_2\right|\ne 1$ if and only if $F(-1)=0$ and $Tr\ne 0,2$;
• ${\lambda }_1$ and ${\lambda }_2$ are complex and $\left|{\lambda }_1\right|=1$ and $\left|{\lambda }_2\right|=1$ if and only if ${\left(Tr\right)}^2-4Det<0$ and $Det=1$.

Lemma 2 

([26]). Assume that ${\lambda }_1$ and ${\lambda }_2$ are eigenvalues a fixed point $(P^{*},R^{*})$; then, $(P^{*},R^{*})$

• is referred to as a sink if $\left|{\lambda }_1\right|<1$ and $\left|{\lambda }_2\right|<1$, so the sink is deemed to be locally asymptotically stable;
• is referred to as a source if $\left|{\lambda }_1\right|>1$ and $\left|{\lambda }_2\right|>1$, so the source is deemed locally unstable;
• is referred to as a saddle if $\left|{\lambda }_1\right|>1$ and $\left|{\lambda }_2\right|<1$ (or $\left|{\lambda }_1\right|<1$ and $\left|{\lambda }_2\right|>1)$;
• is referred to as non-hyperbolic if either $\left|{\lambda }_1\right|=1$ or $\left|{\lambda }_2\right|=1$.

Based on the previous lemmas, we now present the following theorems.

Theorem 1. 

The semi-trivial fixed point $E_1(0,K)$ possesses the following properties:

(i) 

If $0<c<2$, then $E_1$ is a saddle point where $\left|{\lambda }_1\right|>1$ and $\left|{\lambda }_2\right|<1$;

(ii) 

If $c>2$, then $\left|{\lambda }_{1,2}\right|>1$ and $E_1$ is a source;

(iii) 

If $c=2$, then $E_1$ is a non-hyperbolic point where $\left|{\lambda }_1\right|=1+a$ and ${\lambda }_2=-1$.

Proof. 

If Condition (iii) is met, then one of $E_1$’s eigenvalues is $-1$, while the other eigenvalue is different from 1. This condition suggests the occurrence of a flip bifurcation at the fixed point $E_1(0,K)$.    □

Next, we analyze the local stability of the positive fixed point $E_2(K-{\displaystyle \frac{hK}{c}},K-{\displaystyle \frac{hK}{c}})$. The characteristic equation of the positive fixed point $E_2$ is expressed as follows:

$$F\left(\lambda \right):={\lambda }^2+(a-2-A)\lambda +1+A+a(B-1)=0,$$

where $A=2h-c$ and $B=c-h$. Note that $F\left(1\right)=aB>0$ and $F(-1)=4+2A+a(B-2)$.

Applying Lemmas (1) and (2), we draw the following conclusions.

Theorem 2. 

Assume that $B>2$. The dynamic behaviors of Model (4) at the positive fixed point $E_2$ can be described as follows:

(i) 

When a is greater than ${\displaystyle \frac{4+2A}{2-B}}$ and less than ${\displaystyle \frac{A}{1-B}}$, the fixed point $E_2$ becomes a sink and is locally asymptotically stable;

(ii) 

When a is greater than both ${\displaystyle \frac{4+2A}{2-B}}$ and ${\displaystyle \frac{A}{1-B}}$, the fixed point $E_2$ becomes a source and is locally unstable;

(iii) 

When a is less than ${\displaystyle \frac{4+2A}{2-B}}$, the fixed point $E_2$ becomes a saddle;

(iv) 

When a equals ${\displaystyle \frac{4+2A}{2-B}}$ but is not equal to $A+2$ or $A+4$, the fixed point $E_2$ becomes non-hyperbolic, and Model (4) may undergo a flip bifurcation (period-doubling bifurcation);

(v) 

When ${(a-A)}^2<4aB$ and a equals ${\displaystyle \frac{A}{1-B}}$, the fixed point $E_2$ becomes non-hyperbolic, and System (4) may undergo a Neimark–Sacker bifurcation.

3. Codimension-One Bifurcation Analysis

Based on the analysis presented in Section 2, we establish the occurrence of flip bifurcation and Neimark–Sacker bifurcation in Model (4) at the positive fixed point $E_2$. This conclusion is supported by the application of the center manifold theorem and bifurcation theory, as discussed in [27,28].

3.1. Flip Bifurcation

This subsection clarifies the necessity for a flip bifurcation to occur at the positive fixed point $E_2$.

Theorem 3. 

When the parameter a assumes the value $a={\displaystyle \frac{4+2A}{2-B}}$, Model (4) exhibits a flip bifurcation at the positive fixed point $E_2$.

Proof. 

The eigenvalues of Model (4) at the positive fixed point $E_2=(P^{*},R^{*})$ are ${\lambda }_1=-1$ and $|{\lambda }_2| \ne 1$ when $a={\displaystyle \frac{4+2A}{2-B}}$. We need to apply a coordinate transformation to shift the fixed point $E_2$ into the origin. Assume that $u=P-P^{*}$, $v=R-R^{*}$, and substitute a with $a+a^{*}$ in Model (4). As a consequence, Model (4) can be expressed in the following form:

$$\left(\begin{array}{c}\hfill u\\ \hfill v\end{array}\right)\to \left\{\begin{array}{cc}\hfill a_{1}u& +a_{2}v+a_{13}u{a}^{*}+a_{23}v{a}^{*}+a_{11}{u}^2+a_{12}uv+a_{22}{v}^2+a_{113}{u}^2{a}^{*}+a_{123}uv{a}^{*}\hfill \\ \hfill & +a_{223}{v}^2{a}^{*}+O\left(\right(\left|u\right|+\left|v\right|+|a^{*}{\left|\right)}^3),\hfill \\ \hfill b_{1}u& +b_{2}v+b_{11}{u}^2+b_{12}uv+b_{22}{v}^2,\hfill \end{array}\right.$$

(6)

where

$$\left\{\begin{array}{cc}\hfill a_1& =1-a,a_2=a,a_{13}=-1,a_{23}=1,a_{11}=-{\displaystyle \frac{a}{R^{*}}},a_{22}=-{\displaystyle \frac{a}{R^{*}}},a_{12}={\displaystyle \frac{2a}{R^{*}}},a_{113}=-{\displaystyle \frac{1}{R^{*}}},\hfill \\ \hfill a_{223}& =-{\displaystyle \frac{1}{R^{*}}},a_{123}={\displaystyle \frac{2}{R^{*}}},b_1=-h,b_2=1-c+2h,b_{11}=0,b_{12}=0,b_{22}=-{\displaystyle \frac{c}{K}}.\hfill \end{array}\right.$$

(7)

Construct the following invertible matrix:

$$T=\left(\begin{array}{cc}{a}_2& a_2\\ -1-a_1& {\lambda }_2-a_1\end{array}\right).$$

Applying the transformation

$$\begin{array}{c}\hfill \left(\begin{array}{c}u\\ v\end{array}\right)=T\left(\begin{array}{c}\xi \\ \eta \end{array}\right),\end{array}$$

Model (6) undergoes the following form:

$$\begin{array}{c}\hfill \left(\begin{array}{c}\xi \\ \eta \end{array}\right)\to \left(\begin{array}{cc}-1& 0\\ 0& {\lambda }_2\end{array}\right)\left(\begin{array}{c}\xi \\ \eta \end{array}\right)+\left(\begin{array}{c}{f}_1(\xi ,\eta ,a^{*})\\ f_2(\xi ,\eta ,a^{*})\end{array}\right),\end{array}$$

(8)

where

$$\left\{\begin{array}{cc}\hfill f_1(\xi ,\eta ,a^{*})& ={\displaystyle \frac{({\lambda }_2-a_1)a_{13}}{a_2(1+{\lambda }_2)}}u{a}^{*}+{\displaystyle \frac{({\lambda }_2-a_1)a_{23}}{a_2(1+{\lambda }_2)}}v{a}^{*}+{\displaystyle \frac{({\lambda }_2-a_1)a_{11}-a_2{b}_{11}}{a_2(1+{\lambda }_2)}}{u}^2+{\displaystyle \frac{({\lambda }_2-a_1)a_{12}-a_2{b}_{12}}{a_2(1+{\lambda }_2)}}uv\hfill \\ \hfill & +{\displaystyle \frac{({\lambda }_2-a_1)a_{22}-a_2{b}_{22}}{a_2(1+{\lambda }_2)}}{v}^2+{\displaystyle \frac{({\lambda }_2-a_1)a_{113}}{a_2(1+{\lambda }_2)}}{u}^2{a}^{*}+{\displaystyle \frac{({\lambda }_2-a_1)a_{123}}{a_2(1+{\lambda }_2)}}uv{a}^{*}+{\displaystyle \frac{({\lambda }_2-a_1)a_{223}}{a_2(1+{\lambda }_2)}}{v}^2{a}^{*}\hfill \\ \hfill & +O\left(\right(\left|u\right|+\left|v\right|+|a^{*}{\left|\right)}^3),\hfill \\ \hfill f_2(\xi ,\eta ,a^{*})& ={\displaystyle \frac{(1+a_1)a_{13}}{a_2(1+{\lambda }_2)}}u{a}^{*}+{\displaystyle \frac{(1+a_1)a_{23}}{a_2(1+{\lambda }_2)}}v{a}^{*}+{\displaystyle \frac{(1+a_1)a_{11}+a_2{b}_{11}}{a_2(1+{\lambda }_2)}}{u}^2+{\displaystyle \frac{(1+a_1)a_{12}+a_2{b}_{12}}{a_2(1+{\lambda }_2)}}uv\hfill \\ \hfill & +{\displaystyle \frac{(1+a_1)a_{22}+a_2{b}_{22}}{a_2(1+{\lambda }_2)}}{v}^2+{\displaystyle \frac{(1+a_1)a_{113}}{a_2(1+{\lambda }_2)}}{u}^2{a}^{*}+{\displaystyle \frac{(1+a_1)a_{123}}{a_2(1+{\lambda }_2)}}uv{a}^{*}+{\displaystyle \frac{(1+a_1)a_{223}}{a_2(1+{\lambda }_2)}}{v}^2{a}^{*}\hfill \\ \hfill & +O\left(\right(\left|u\right|+\left|v\right|+|a^{*}{\left|\right)}^3),\hfill \end{array}\right.$$

and

$$\left\{\begin{array}{cc}\hfill u& =a_2(\xi +\eta ),v=-(1+a_1)\xi +({\lambda }_2-a_1)\eta ,u^2={a_2}^2({\xi }^2+\xi \eta +{\eta }^2),\hfill \\ \hfill uv& =a_2[-(1+a_1){\xi }^2+({\lambda }_2-2a_1-1)\xi \eta +({\lambda }_2-a_1){\eta }^2],\hfill \\ \hfill v^2& ={(1+a_1)}^2{\xi }^2-2(1+a_1)({\lambda }_2-a_1)\xi \eta +{({\lambda }_2-a_1)}^2{\eta }^2.\hfill \end{array}\right.$$

Using the center manifold theory [27,28], we can determine a center manifold $W^c(0,0,0)$ of Equation (8) at the fixed point $(0,0)$ within a small neighborhood of $a^{*}$. This center manifold can be expressed as follows:

$$W^c(0,0,0)=\{(\xi ,\eta ,a^{*})\in R^3,\eta =f(\xi ,a^{*}),f(0,0)=0,Df(0,0)=0\}.$$

Assume that the center manifold can be represented by the following expression:

$$\begin{array}{c}\hfill f(\xi ,a^{*})=c_1{\xi }^2+c_2\xi a^{*}+c_3{a^{*}}^2+O\left(\right(\left|\xi \right|+|a^{*}{\left|\right)}^3).\end{array}$$

(9)

The center manifold is required to meet

$$\begin{array}{c}\hfill f(-\xi +f_1(\xi ,f(\xi ,a^{*}),a^{*}))={\lambda }_{2}f(\xi ,a^{*})+f_2(\xi ,f(\xi ,a^{*}),a^{*}).\end{array}$$

(10)

By substituting from Equation (9) into Equation (10) and comparing the coefficients of same powers, we obtain

$$\begin{array}{cc}\hfill c_1& ={\displaystyle \frac{1}{a_2(1-{\lambda }_2^2)}}[(1+a_1)a_2^2{a}_{11}+a_2^3{b}_{11}-{(1+a_1)}^2{a}_2{a}_{12}-(1+a_1)a_2^2{b}_{12}+{(1+a_1)}^3{a}_{22}+{(1+a_1)}^2{a}_2{b}_{22}],\hfill \\ \hfill c_2& ={\displaystyle \frac{1}{a_2{(1+{\lambda }_2)}^2}}[{(1+a_1)}^2{a}_{23}-(1+a_1)a_2{a}_{13}],c_3=0.\hfill \end{array}$$

So, Model (8) is restricted to the following center manifold:

$$\begin{array}{c}\hfill G:\xi \to -\xi +{\delta }_{11}{\xi }^2+{\delta }_{12}\xi a^{*}+{\delta }_{111}{\xi }^3+{\delta }_{112}{\xi }^2{a}^{*}+{\delta }_{122}\xi {a^{*}}^2+O\left(\right(\left|\xi \right|+|a^{*}{\left|\right)}^4),\end{array}$$

(11)

where

$$\begin{array}{cc}\hfill {\delta }_{11}& ={\displaystyle \frac{1}{a_2(1+{\lambda }_2)}}[({\lambda }_2-a_1)a_2^2{a}_{11}-a_2^3{b}_{11}+({\lambda }_2-a_1){(1+a_1)}^2{a}_{22}-{(1+a_1)}^2{a}_2{b}_{22}\hfill \\ \hfill & -(1+a_1)({\lambda }_2-a_1)a_2{a}_{12}+(1+a_1)a_2^2{b}_{12}],\hfill \\ \hfill {\delta }_{12}& ={\displaystyle \frac{1}{a_2({\lambda }_2+1)}}[({\lambda }_2-1)a_2{a}_{13}-({\lambda }_2-a_1)(1+a_1)a_{23}],\hfill \\ \hfill {\delta }_{111}& ={\displaystyle \frac{c_1}{a_2(1+{\lambda }_2)}}[2({\lambda }_2-a_1)a_2^2{a}_{11}-2a_2^3{b}_{11}-2{({\lambda }_2-a_1)}^2(1+a_1)a_{22}+2({\lambda }_2-a_1)(1+a_1)a_2{b}_{22}\hfill \\ \hfill & +({\lambda }_2-a_1)({\lambda }_2-2a_1-1)a_2{a}_{12}-({\lambda }_2-2a_1-1)a_2^2{b}_{12}],\hfill \\ \hfill {\delta }_{112}& ={\displaystyle \frac{1}{a_2({\lambda }_2-a_1)}}[2({\lambda }_2-a_1)c_2{a}_2^2{a}_{11}-2c_2{a}_2^3{b}_{11}+({\lambda }_2-a_1)a_2^2{a}_{113}+({\lambda }_2-a_1){(1+a_1)}^2{a}_{223}\hfill \\ \hfill & -2{({\lambda }_2-a_1)}^2(1+a_1)c_2{a}_{11}+2({\lambda }_2-a_1)(1+a_1)c_2{a}_2{b}_{11}+{({\lambda }_2-a_1)}^2{c}_1{a}_{223}-({\lambda }_2-a_1)(1+a_1)a_2{a}_{123}\hfill \\ \hfill & +({\lambda }_2-a_1)a_2{c}_1{a}_{13}+({\lambda }_2-a_1)({\lambda }_2-2a_1-1)a_2{c}_2{a}_{12}+({\lambda }_2-2a_1-1)a_2^2{c}_2{b}_{12}],\hfill \\ \hfill {\delta }_{122}& ={\displaystyle \frac{1}{a_2(1+{\lambda }_2)}}[({\lambda }_2-a_1)a_2{c}_2{a}_{13}+{({\lambda }_2-a_1)}^2{c}_2{a}_{23}].\hfill \end{array}$$

For a flip bifurcation to occur, it is necessary that the following discriminant quantities:

$$\begin{array}{ccc}\hfill {\alpha }_1& =& {\left(2{\displaystyle \frac{{\partial }^{2}G}{\partial \xi \partial a^{*}}}+{\displaystyle \frac{\partial G}{\partial a^{*}}}{\displaystyle \frac{\partial G}{\partial \xi }}\right)}_{(0,0)}=2{\delta }_{12},\hfill \\ \hfill {\alpha }_2& =& {\left({\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{{\partial }^{2}G}{\partial {\xi }^2}}\right)}^2+{\displaystyle \frac{1}{6}}\left({\displaystyle \frac{{\partial }^{3}G}{\partial {\xi }^3}}\right)\right)}_{(0,0)}={\delta }_{111}+{{\delta }_{11}}^2,\hfill \end{array}$$

are not zero. Thus, as a result of the above analysis, Model (4) undergoes a flip bifurcation at $E_2$ when $a={\displaystyle \frac{4+A}{2-B}}$ and ${\alpha }_1\ne 0$, ${\alpha }_2\ne 0$.    □

Remark 1. 

Flip bifurcation as a biological phenomenon that happens when the population size fluctuates with periods 2, 4, 8, …, until it becomes completely chaotic.

3.2. Neimark–Sacker Bifurcation

In this subsection, we investigate the Neimark–Sacker bifurcation for Model (4) at the positive fixed point $E_2$, in particular when the parameter a is set to $a={\displaystyle \frac{A}{1-B}}$. Applying the transformations $u=P-P^{*}$, $v=R-R^{*}$ and $a=a+\overline{a}$, we may shift $E_2$ in System (4) to the origin. Consequently, we obtain the following expression:

$$\begin{array}{c}\hfill \left(\begin{array}{c}\hfill u\\ \hfill v\end{array}\right)\to \left\{\begin{array}{cc}& a_{1}u+a_{2}v+a_{11}{u}^2+a_{12}uv+a_{22}{v}^2+O\left(\right(\left|u\right|+{\left|v\right|)}^3),\hfill \\ & b_{1}u+b_{2}v+b_{11}{u}^2+b_{22}{v}^2,\hfill \end{array}\right.\end{array}$$

(12)

The values of $a_1,a_2,a_{11},a_{12},a_{22},b_1,b_2$, and $b_{22}$ are defined in (7) with $a=a+\overline{a}$. The linearized system in (12) at the point $(u,v)=(0,0)$ yields the following characteristic equation:

$${\lambda }^2+m\left(\overline{a}\right)\lambda +n\left(\overline{a}\right)=0,$$

where

$$m\left(\overline{a}\right)=a+\overline{a}-A-2,n\left(\overline{a}\right)=1+A+(a+\overline{a})(B-1).$$

Based on Theorem (2), the eigenvalues of Model (12) are a complex conjugate pair, denoted as ${\lambda }_1$ and ${\lambda }_2$, having modulus 1, where

$${\lambda }_{1,2}=-{\displaystyle \frac{m\left(\overline{a}\right)}{2}}\pm {\displaystyle \frac{i}{2}}\sqrt{4n\left(\overline{a}\right)-{\left(m\left(\overline{a}\right)\right)}^2}.$$

Moreover,

$$|{\lambda }_{1,2}{|}_{\overline{a}=0}=\sqrt{n\left(0\right)}=1,{\displaystyle \frac{d|{\lambda }_{1,2}|}{d\overline{a}}}{|}_{\overline{a}=0}={\displaystyle \frac{1}{2}}(B-1)\ne 0.$$

To satisfy the conditions ${\lambda }_{1,2}^j\ne 1$ for $j=1,2,3,4$, which implies that $m\left(0\right)\ne -2,0,1,2$, we only need to require that

$$m\left(0\right)\ne 0,1.$$

(13)

Assume that $a-A-2\ne 0,1.$ In this case, the eigenvalues ${\lambda }_{1,2}$ of the fixed point $(0,0)$ in System (12) do not fall within the intersection of the unit circle, and the conditions in (13) are satisfied. In order to obtain the normal form of Model (12), we set

$$\mu =-{\displaystyle \frac{m\left(0\right)}{2}},\omega ={\displaystyle \frac{1}{2}}\sqrt{4n\left(0\right)-{\left(m\left(0\right)\right)}^2},$$

at $\overline{a}=0$. Assume an invertible matrix as follows:

$$\begin{array}{c}\hfill T=\left(\begin{array}{cc}{a}_2& 0\\ \mu -a_1& -\omega \end{array}\right),\end{array}$$

Using the following transformation:

$$\begin{array}{c}\hfill \left(\begin{array}{c}u\\ v\end{array}\right)=T\left(\begin{array}{c}\xi \\ \eta \end{array}\right),\end{array}$$

Model (12) is transformed into the following form:

$$\begin{array}{c}\hfill \left(\begin{array}{c}\xi \\ \eta \end{array}\right)\to \left(\begin{array}{cc}\mu & -\omega \\ \omega & \mu \end{array}\right)\left(\begin{array}{c}\xi \\ \eta \end{array}\right)+\left(\begin{array}{c}\overline{f_1}(\xi ,\eta )\\ \overline{f_2}(\xi ,\eta )\end{array}\right),\end{array}$$

(14)

where

$$\begin{array}{cc}\hfill \overline{f_1}(\xi ,\eta )& ={\displaystyle \frac{1}{a_2}}[a_{11}{u}^2+a_{12}uv+a_{22}{v}^2]+O\left(\right(\left|\xi \right|+{\left|\eta \right|)}^3),\hfill \\ \hfill \overline{f_2}(\xi ,\eta )& =\left[{\displaystyle \frac{a_{11}(\mu -a_1)-a_2{b}_{11}}{a_2\omega }}\right]u^2+\left[{\displaystyle \frac{a_{12}(\mu -a_1)}{a_2\omega }}\right]uv+\left[{\displaystyle \frac{a_{22}(\mu -a_1)-a_2{b}_{22}}{a_2\omega }}\right]v^2,\hfill \end{array}$$

and

$$\begin{array}{cc}\hfill u& =a_2\xi ,v=(\mu -a_1)\xi -\omega \eta ,u^2={a_2}^2{\xi }^2,uv=a_2\xi [(\mu -a_1)\xi -\omega \eta ]\hfill \\ \hfill v^2& ={(\mu -a_1)}^2{\xi }^2-2(\mu -a_1)\omega \xi \eta +{\omega }^2{\eta }^2.\hfill \end{array}$$

For System (14) to undergo a Neimark–Sacker bifurcation, it is necessary that the following discriminatory quantity is nonzero:

$$\theta =-Re\left({\displaystyle \frac{(1-2\lambda ){\overline{\lambda }}^2}{1-\lambda }}{L}_{11}{L}_{12}\right)-{\displaystyle \frac{1}{2}}|L_{11}{|}^2-{\left|L_{21}\right|}^2+Re\left(\overline{\lambda }{L}_{22}\right),$$

where

$$\begin{array}{ccc}\hfill L_{11}& =& {\displaystyle \frac{1}{4}}(({\overline{f}}_{1\xi \xi }+{\overline{f}}_{1\eta \eta })+i({\overline{f}}_{2\xi \xi }+{\overline{f}}_{2\eta \eta })),\hfill \\ \hfill L_{12}& =& {\displaystyle \frac{1}{8}}(({\overline{f}}_{1\xi \xi }-{\overline{f}}_{1\eta \eta }+2{\overline{f}}_{2\xi \eta })+i({\overline{f}}_{2\xi \xi }-{\overline{f}}_{2\eta \eta }-2{\overline{f}}_{1\xi \eta })),\hfill \\ \hfill L_{21}& =& {\displaystyle \frac{1}{8}}(({\overline{f}}_{1\xi \xi }-{\overline{f}}_{1\eta \eta }-2{\overline{f}}_{2\xi \eta })+i({\overline{f}}_{2\xi \eta }-{\overline{f}}_{2\eta \eta }+2{\overline{f}}_{1\xi \eta })),\hfill \\ \hfill L_{22}& =& {\displaystyle \frac{1}{16}}(({\overline{f}}_{1\xi \xi \xi }+{\overline{f}}_{1\xi \eta \eta }+{\overline{f}}_{2\xi \xi \eta }+{\overline{f}}_{2\eta \eta \eta })+i({\overline{f}}_{2\xi \xi \xi }+{\overline{f}}_{2\xi \eta \eta }-{\overline{f}}_{1\xi \xi \eta }-{\overline{f}}_{1\eta \eta \eta })),\hfill \end{array}$$

Based on the analysis presented above, we can conclude the following result.

Theorem 4. 

If Condition (13) is satisfied and θ is nonzero, System (4) undergoes a Neimark–Sacker bifurcation at the fixed point $E_2(P^{*},R^{*})$. Additionally, if θ is negative (resp., positive), a closed invariant curve emerges at $a={\displaystyle \frac{A}{1-B}}$, which is subcritical (resp., supercritical) and asymptotically stable (resp., unstable).

Remark 2. 

From a sustainability perspective, a stable invariant curve arises from the coexistence fixed point $E_2$ once a surpasses the critical value ${\displaystyle \frac{A}{1-B}}$. This indicates a stable and reproducing cohabitation between the human population and resources. On the other hand, ecological imbalance will result from human population and their resource instability will occur if the invariant curve bifurcated from $E_2$ is unstable as a approaches the critical value ${\displaystyle \frac{A}{1-B}}$.

4. Codimension-Two Bifurcation Analysis

In this section, at the positive fixed point $E_2(P^{*},R^{*})$, we investigate the strong resonance phenomena for Model (4), as mentioned in [22,23,24,25]. Consider the characteristic equation of Model (4):

$${\lambda }^2+p(P^{*},R^{*})\lambda +q(P^{*},R^{*})=0,$$

(15)

where $p=a-A-2$, $q=1+A+a(B-1)$, and A, B are previously defined. Referring to Equation (15), we obtain the following expression:

$${\lambda }_{1,2}={\displaystyle \frac{-p(P^{*},R^{*})\pm \sqrt{p^2(P^{*},R^{*})-4q(P^{*},R^{*})}}{2}}.$$

When the values of $p(P^{*},R^{*})$ and $q(P^{*},R^{*})$ are such that $p(P^{*},R^{*})=2$ and $q(P^{*},R^{*})=1$, the eigenvalues ${\lambda }_1$ and ${\lambda }_2$ become equal, and both are equal to $-1$. This signifies that Model (4) exhibits $1:2$ resonance. If $p(P^{*},R^{*})=q(P^{*},R^{*})=1$, then ${\lambda }_{1,2}={\displaystyle \frac{(-1)\pm i\sqrt{3}}{2}}$. This signifies that Model (4) exhibits $1:3$ resonance. If $p(P^{*},R^{*})=0$ and $q(P^{*},R^{*})=1$, then ${\lambda }_{1,2}=\pm i$. This signifies that Model (4) exhibits $1:4$ resonance. Thus, strong resonance phenomena are characterized by the conditions specified by the following sets, $F_{1j}$, defined as follows:

$$F_{1j}=\left\{(a,c,h,K)\in {\mathbb{R}}_{+}^4\mid p(P^{*},R^{*})=4-j,q(P^{*},R^{*})=1\right\},j=2,3,4.$$

4.1. The 1:2 Resonance of Model (4) at the Positive Fixed Point $E_2(P^{*},R^{*})$

In this subsection, we explore the conditions for the occurrence of 1:2 resonance at $E_2$ by applying bifurcation theory and normal form theory, as discussed in references [27,28,29]. We randomly select parameters $(\tilde{a},c,\tilde{h},K)$ from the set $F_{12}$. Consequently, Model (4) can be expressed as follows:

$$\begin{array}{cc}\hfill P_{n+1}& =P_n+\tilde{a}{P}_n(1-{\displaystyle \frac{P_n}{R_n}}),\hfill \\ \hfill R_{n+1}& =R_n+c{R}_n(1-{\displaystyle \frac{R_n}{K}})-\tilde{h}{P}_n.\hfill \end{array}$$

(16)

Model (16) has two eigenvalues ${\lambda }_1={\lambda }_2=-1$ at the fixed point $E_2(P^{*},R^{*})$. To investigate the effects of bifurcation parameters $\tilde{a}$ and $\tilde{h}$, we introduce a small perturbation into Model (16). This perturbed model can be described as follows:

$$\begin{array}{cc}\hfill P_{n+1}& =P_n+(\tilde{a}+a^{*})P_n(1-{\displaystyle \frac{P_n}{R_n}}),\hfill \\ \hfill R_{n+1}& =R_n+c{R}_n(1-{\displaystyle \frac{R_n}{K}})-(\tilde{h}+h^{*})P_n,\hfill \end{array}$$

(17)

where $\left|a^{*}\right|,\left|h^{*}\right|\ll 1$ are small perturbation parameters. Assume that $\tilde{u}\left(n\right)=P\left(n\right)-P^{*}$, $\tilde{v}\left(n\right)=R\left(n\right)-R^{*}$, $a=\tilde{a}+a^{*}$, and $h=\tilde{h}+h^{*}$. We perform a coordinate transformation to shift the fixed point $E_2$ into the origin. This transformation allows us to simplify Model (17), resulting in the following expression:

$$\begin{array}{cc}\hfill \tilde{u}(n+1)& =\left(1+a-{\displaystyle \frac{2a{P}^{*}}{R^{*}}}\right)\tilde{u}\left(n\right)+{\displaystyle \frac{a{P^{*}}^2}{{R^{*}}^2}}\tilde{v}\left(n\right)-{\displaystyle \frac{a}{R^{*}}}\widetilde{u^2}\left(n\right)+{\displaystyle \frac{2a{P}^{*}}{{R^{*}}^2}}\tilde{u}\left(n\right)\tilde{v}\left(n\right)\hfill \\ \hfill & -{\displaystyle \frac{a{P^{*}}^2}{{R^{*}}^3}}\widetilde{v^2}\left(n\right)+O\left(\right(|\tilde{u}\left(n\right)|+|\tilde{v}{\left(n\right)\left|\right)}^3),\hfill \\ \hfill \tilde{v}(n+1)& =-h\tilde{u}\left(n\right)+\left(1+c-{\displaystyle \frac{-2c{R}^{*}}{K}}\right)\tilde{v}\left(n\right)-{\displaystyle \frac{c}{K}}\widetilde{v^2}\left(n\right).\hfill \end{array}$$

(18)

Let T be an invertible matrix is given by

$$\begin{array}{c}\hfill T=\left(\begin{array}{cc}-{\displaystyle \frac{a_2}{1+a_1}}& -{\displaystyle \frac{a_2}{{(1+a_1)}^2}}\\ 1& 0\end{array}\right),\end{array}$$

and use the translation

$$\begin{array}{c}\hfill \left(\begin{array}{c}\tilde{u}\\ \tilde{v}\end{array}\right)=T\left(\begin{array}{c}\tilde{P}\\ \tilde{R}\end{array}\right).\end{array}$$

Then, Map (18) is transformed into

$$\begin{array}{c}\hfill \left(\begin{array}{c}\tilde{P}\\ \tilde{R}\end{array}\right)\mapsto \left(\begin{array}{cc}-1+a_{10}\left(a,h\right)& 1+a_{01}\left(a,h\right)\\ b_{10}(a,h)& -1+b_{01}\left(a,h\right)\end{array}\right)\left(\begin{array}{c}\tilde{P}\\ \tilde{R}\end{array}\right)+\left(\begin{array}{c}\tilde{f}\left(\tilde{P},\tilde{R},a,h\right)\\ \tilde{g}\left(\tilde{P},\tilde{R},a,h\right)\end{array}\right),\end{array}$$

(19)

where

$$\tilde{f}\left(\tilde{P},\tilde{R},a,h\right)=\sum _{2\le j+k\le 3}{\tilde{f}}_{j,k}\left(a,h\right){\tilde{P}}^j{\tilde{R}}^k,$$

$$\tilde{g}\left(\tilde{P},\tilde{R},a,h\right)=\sum _{2\le j+k\le 3}{\tilde{g}}_{j,k}\left(a,h\right){\tilde{P}}^j{\tilde{R}}^k,$$

and

$$\left\{\begin{array}{cc}\hfill a_{10}& =2+c+{\displaystyle \frac{a_{2}h}{1+a_1}}-{\displaystyle \frac{2c{R}^{*}}{K}},a_{01}=-1+{\displaystyle \frac{a_{2}h}{{(1+a_1)}^2}},\hfill \\ \hfill b_{10}& =-{\displaystyle \frac{1}{a_{2}K{R^{*}}^2}}\left\{\left[a_2^{2}hk+c(1+a_1)(K-2R^{*})\right]{R^{*}}^2+aK{(1+a_1)}^2{P^{*}}^2+2a_2{P}^{*}{R}^{*}-a_2{R^{*}}^2\right\},\hfill \\ \hfill b_{01}& =1-{\displaystyle \frac{1}{R^{*}(1+a_1)}}\left\{a(1+a_1)(2P^{*}-R^{*})-R^{*}(1+a_1-a_{2}h)\right\},{\tilde{f}}_{20}=-{\displaystyle \frac{c}{K}},\hfill \\ \hfill {\tilde{f}}_{11}& ={\tilde{f}}_{02}={\tilde{f}}_{30}={\tilde{f}}_{12}={\tilde{f}}_{21}={\tilde{f}}_{03}=0,\hfill \\ \hfill {\tilde{g}}_{20}& ={\displaystyle \frac{1}{a_{2}K{R^{*}}^3}}\left\{a_{2}c(1+a_1){R^{*}}^3+aK{(P^{*}+a_1{P}^{*}+a_2{R}^{*})}^2\right\},\hfill \\ \hfill {\tilde{g}}_{11}& ={\displaystyle \frac{1}{(1+a_1){R^{*}}^2}}\left\{2a(P^{*}+a_1{P}^{*}+a_2{R}^{*})\right\},{\tilde{g}}_{02}={\displaystyle \frac{a{a}_2}{{(1+a_1)}^2{R}^{*}}},\hfill \\ \hfill {\tilde{g}}_{30}& ={\tilde{g}}_{12}={\tilde{g}}_{21}={\tilde{g}}_{03}=0.\hfill \end{array}\right.$$

Now, we utilize a nonsingular linear coordinate transformation, which is provided by

$$\begin{array}{c}\hfill \left(\begin{array}{c}\tilde{P}\\ \tilde{R}\end{array}\right)=\left(\begin{array}{cc}1+a_{01}\left(a,h\right)& 0\\ -a_{10}\left(a,h\right)& 1\end{array}\right)\left(\begin{array}{c}\widehat{P}\\ \widehat{R}\end{array}\right).\end{array}$$

As a result, Equation (19) is transformed into

$$\begin{array}{c}\hfill \left(\begin{array}{c}\widehat{P}\\ \widehat{R}\end{array}\right)\mapsto \left(\begin{array}{cc}-1& 1\\ {\theta }_1\left(a,h\right)& -1+{\theta }_2\left(a,h\right)\end{array}\right)\left(\begin{array}{c}\widehat{P}\\ \widehat{R}\end{array}\right)+\left(\begin{array}{c}\widehat{f}\left(\widehat{P},\widehat{R},a,h\right)\\ \widehat{g}\left(\widehat{P},\widehat{R},a,h\right)\end{array}\right),\end{array}$$

(20)

where

$$\widehat{f}\left(\widehat{P},\widehat{R},a,h\right)=\sum _{2\le j+k\le 3}{\widehat{f}}_{j,k}\left(a,h\right){\widehat{P}}^j{\widehat{R}}^k,$$

$$\widehat{g}\left(\widehat{P},\tilde{R},a,h\right)=\sum _{2\le j+k\le 3}{\widehat{g}}_{j,k}\left(a,h\right){\widehat{P}}^j{\widehat{R}}^k,$$

and

$$\left\{\begin{array}{cc}\hfill {\theta }_1& ={\displaystyle \frac{-1}{K{a}_2{R^{*}}^2}}[a_2{R^{*}}^2(K(1+2c+a_{1}c+a_{2}h)-2c{R}^{*}(2+a_1))+a(2c{a}_2{R^{*}}^2(2P^{*}-R^{*})\hfill \\ \hfill & +K((1+2a_1+{a_1}^2+a_{2}h){P^{*}}^2+2a_2(a_1-c)P^{*}{R}^{*}+a_2(-a_1+c){R^{*}}^2)\left)\right],\hfill \\ \hfill {\theta }_2& =a+c-{\displaystyle \frac{2a{P}^{*}}{R^{*}}}-{\displaystyle \frac{2c{R}^{*}}{K}},{\widehat{f}}_{20}=-{\displaystyle \frac{c(1+2a_1+{a_1}^2+a_{2}h)}{K{(1+a_1)}^2}},\hfill \\ \hfill {\widehat{f}}_{11}& ={\widehat{f}}_{20}={\widehat{f}}_{30}={\widehat{f}}_{12}={\widehat{f}}_{21}={\widehat{f}}_{03},\hfill \\ \hfill {\widehat{g}}_{20}& ={\displaystyle \frac{1}{a_2{K}^2{R^{*}}^3{(1+a_1)}^2}}\{a_{2}c{R^{*}}^3(1+2a_1+{a_1}^2+a_{2}h)(a_{1}K-cK+2c{R}^{*})+a(2c{a}_2{R^{*}}^2\hfill \\ \hfill & +K((1+2a_1+{a_1}^2+a_{2}h)P^{*}+a_2(a_1-c)R^{*}){)}^2\},\hfill \\ \hfill {\widehat{g}}_{11}& ={\displaystyle \frac{1}{K{R^{*}}^2{(1+a_1)}^2}}\left\{2a(2c{a}_2{R^{*}}^2+K((1+2a_1+{a_1}^2+a_{2}h)P^{*}+a_2(a_1-c)R^{*}))\right\},\hfill \\ \hfill {\widehat{g}}_{02}& ={\displaystyle \frac{a{a}_2}{R^{*}{(1+a_1)}^2}},{\widehat{g}}_{30}={\widehat{g}}_{12}={\widehat{g}}_{21}={\widehat{g}}_{03}=0.\hfill \end{array}\right.$$

Finally, we apply the following transformation:

$$\begin{array}{c}\hfill \widehat{P}=\xi +\sum _{2\le j+k\le 3}{\varphi }_{jk}\left(\tilde{a},\tilde{h}\right){\xi }^j{\eta }^k,\\ \hfill \widehat{R}=\eta +\sum _{2\le j+k\le 3}{\psi }_{jk}\left(\tilde{a},\tilde{h}\right){\xi }^j{\eta }^k,\end{array}$$

(21)

The coefficients ${\varphi }_{jk}$ and ${\psi }_{jk}$ can be obtained through the following relations:

$$\left\{\begin{array}{cc}\hfill {\varphi }_{20}& ={\displaystyle \frac{1}{2}}{\widehat{f}}_{20}+{\displaystyle \frac{1}{4}}{\widehat{g}}_{20},{\varphi }_{11}={\displaystyle \frac{1}{2}}{\widehat{f}}_{20}+{\displaystyle \frac{1}{2}}{\widehat{f}}_{11}+{\displaystyle \frac{1}{2}}{\widehat{g}}_{20}+{\displaystyle \frac{1}{4}}{\widehat{g}}_{11},\hfill \\ \hfill {\varphi }_{02}& ={\displaystyle \frac{1}{4}}{\widehat{f}}_{11}+{\displaystyle \frac{1}{2}}{\widehat{f}}_{02}+{\displaystyle \frac{1}{8}}{\widehat{g}}_{20}+{\displaystyle \frac{1}{4}}{\widehat{g}}_{11}+{\displaystyle \frac{1}{4}}{\widehat{g}}_{02},\hfill \\ \hfill {\psi }_{20}& ={\displaystyle \frac{1}{2}}{\widehat{g}}_{20},{\psi }_{11}={\displaystyle \frac{1}{2}}{\widehat{g}}_{20}+{\displaystyle \frac{1}{2}}{\widehat{g}}_{11},{\psi }_{02}={\displaystyle \frac{1}{4}}{\widehat{g}}_{11}+{\displaystyle \frac{1}{2}}{\widehat{g}}_{02}\hfill \end{array}\right.$$

By applying Transformation (21) and its inverse transformation, we obtain the following expression:

$$\begin{array}{c}\hfill \left(\begin{array}{c}\xi \\ \eta \end{array}\right)\mapsto \left(\begin{array}{cc}-1& 1\\ {\theta }_1\left(a,h\right)& 1+{\theta }_2\left(a,h\right)\end{array}\right)\left(\begin{array}{c}\xi \\ \eta \end{array}\right)+\left(\begin{array}{c}\mathsf{\Gamma }\left(\widehat{P},\widehat{R},a,h\right)\\ \mathsf{\Sigma }\left(\widehat{P},\widehat{R},a,h\right)\end{array}\right),\end{array}$$

(22)

where

$$\begin{array}{c}\hfill \mathsf{\Gamma }\left(\widehat{P},\widehat{R},a,h\right)=\sum _{2\le j+k\le 3}{\gamma }_{jk}{\xi }^j{\eta }^k+O\left(\left(\right|\xi |+{\left|\eta \right|)}^4\right),\\ \hfill \mathsf{\Sigma }\left(\widehat{P},\widehat{R},a,h\right)=\sum _{2\le j+k\le 3}{\sigma }_{jk}{\xi }^j{\eta }^k+O\left(\left(\right|\xi |+{\left|\eta \right|)}^4\right),\end{array}$$

where ${\gamma }_{jk}$ and ${\sigma }_{jk}$ are given below:

$$\left\{\begin{array}{cc}\hfill {\gamma }_{20}(a,h)& =\widehat{f_{20}}+{\psi }_{20}-2{\varphi }_{20}-{\varphi }_{02}{{\theta }_1}^2+{\varphi }_{11}{\theta }_1,\hfill \\ \hfill {\gamma }_{11}(a,h)& =\widehat{f_{11}}+{\psi }_{11}-2{\varphi }_{02}{\theta }_1(1+{\theta }_2)+{\varphi }_{11}({\theta }_2-{\theta }_1)+2{\varphi }_{20},\hfill \\ \hfill {\gamma }_{02}(a,h)& =\widehat{f_{02}}+{\psi }_{02}-{\varphi }_{02}(1+{(1+{\theta }_2)}^2)+{\varphi }_{11}({\theta }_2+1)-{\varphi }_{20},\hfill \\ \hfill {\sigma }_{20}(a,h)& =\widehat{g_{20}}-{\psi }_{02}{{\theta }_1}^2+{\psi }_{11}{\theta }_1+{\psi }_{20}{\theta }_2+{\psi }_{20}{\theta }_1,\hfill \\ \hfill {\sigma }_{11}(a,h)& =\widehat{g_{11}}-2{\psi }_{02}{\theta }_1({\theta }_2+1)+{\psi }_{11}(2-{\theta }_1+2{\theta }_2)+2{\psi }_{20}+{\varphi }_{11}{\theta }_1,\hfill \\ \hfill {\sigma }_{02}(a,h)& =\widehat{g_{02}}-{\psi }_{02}{\theta }_2(1+{\theta }_2)-{\psi }_{11}(1+{\theta }_2)-{\psi }_{20}+{\varphi }_{02}{\theta }_1.\hfill \end{array}\right.$$

By setting ${\gamma }_{jk}={\sigma }_{jk}=0$ with $j+k=2$, we can eliminate all quadratic terms. Also, all cubic terms are annihilated except those resonant terms by setting ${\gamma }_{jk}={\sigma }_{jk}=0$ with $j+k=3$. For $j+k=3$, this results in a system for ${\varphi }_{jk}$ and ${\psi }_{jk}$, from which we may obtain ${\varphi }_{jk}$ and ${\psi }_{jk}$. From the above transformations, Map (22) is transformed to the $1:2$ resonance normal form as follows:

$$\begin{array}{c}\hfill \left(\begin{array}{c}\xi \\ \eta \end{array}\right)\mapsto \left(\begin{array}{c}-\xi +\eta \\ {\theta }_1\left(a,h\right)\xi +\left(-1+{\theta }_2\left(a,h\right)\right)\eta +C\left(a,h\right){\xi }^3+D\left(a,h\right){\xi }^2\eta \end{array}\right)+\begin{array}{c}O{\left(\left(\right|\xi |+|\eta \left|\right)\right)}^4,\end{array}\end{array}$$

(23)

where $C\left(a,h\right)={\sigma }_{30}$ and $D\left(a,h\right)={\sigma }_{21}$. When $\left(a,h\right)=\left(\tilde{a},\tilde{h}\right)$, we obtain ${\theta }_1\left(a,h\right)={\theta }_2\left(a,h\right)=0,$

$$C\left(\tilde{a},\tilde{h}\right)={\widehat{g}}_{30}+{\widehat{f}}_{20}{\widehat{g}}_{20}+{\displaystyle \frac{1}{2}}{{\widehat{g}}_{20}}^2+{\displaystyle \frac{1}{2}}{\widehat{g}}_{20}{\widehat{g}}_{11},$$

$$\begin{array}{ccc}\hfill D\left(\tilde{a},\tilde{h}\right)& =& {\widehat{g}}_{21}+3{\widehat{f}}_{30}+{\displaystyle \frac{1}{2}}{\widehat{f}}_{20}{\widehat{g}}_{11}+{\displaystyle \frac{5}{4}}{\widehat{g}}_{20}{\widehat{g}}_{11}\hfill \\ & +& {\widehat{g}}_{20}{\widehat{g}}_{02}+3{\widehat{f}}_{20}+{\displaystyle \frac{5}{2}}{\widehat{f}}_{20}{\widehat{g}}_{20}+{\displaystyle \frac{5}{2}}{\widehat{f}}_{11}{\widehat{g}}_{20}+{{\widehat{g}}_{20}}^2+{\displaystyle \frac{1}{2}}{{\widehat{g}}_{11}}^2.\hfill \end{array}$$

By utilizing Lemma 9.10 and Theorem 9.3 from reference [28], we can establish the following theorem.

Theorem 5. 

Suppose that $C\left(\tilde{a},\tilde{h}\right)\ne 0$ and $D\left(\tilde{a},\tilde{h}\right)+3C\left(\tilde{a},\tilde{h}\right)\ne 0$; then, Map (4) contains the subsequent dynamical behaviors:

• A flip bifurcation curve is present at $F=\{\left({\theta }_1,{\theta }_2\right):{\theta }_1=0\},$ and there exists a non-trivial fixed point for ${\theta }_1<0;$
• A non-degenerate Neimark–Sacker bifurcation curve can be identified:

  $$H=\{\left({\theta }_1,{\theta }_2\right):{\theta }_1=-{\theta }_2+O\left(\left(\right|{\theta }_1|+|{\theta }_2{\left|\right)}^2\right),{\theta }_1<0\};$$
• A heteroclinic bifurcation curve exists:

  $$HL=\{\left({\theta }_1,{\theta }_2\right):{\theta }_1=-{\displaystyle \frac{5}{3}}{\theta }_2+O\left(\left(\right|{\theta }_1|+|{\theta }_2{\left|\right)}^2\right),{\theta }_1<0\}.$$

Remark 3. 

The existence of a 1:2 strong resonance signifies that Model (4) is acutely responsive to variations in bifurcation parameters, influencing its intricate dynamics. The non-degenerate Neimark–Sacker bifurcation has important biological consequences, causing periodic or quasi-periodic fluctuations in the population–resource system as the bifurcation parameters (a,h) move around the ($\tilde{a}$,$\tilde{h}$) region. These fluctuations can lead to long-period variations, significant population surges, and even chaotic behavior in the population–resource system. These arise from periodic oscillations with periods of 2, 4, and 8, or due to the presence of a homoclinic structure.

4.2. The 1:3 Resonance of Model (4) at the Positive Fixed Point $E_2(P^{*},R^{*})$

In this subsection, we explain the dynamical behavior that occurs near 1:3 resonance for Model (4). Taking both $\left(\widehat{a},\widehat{h}\right)$ as bifurcation parameters from $F_{13}$, Model (4) with parameters $\left(\widehat{a},c,\widehat{h},K\right)$ is formed as

$$\begin{array}{cc}\hfill P_{n+1}& =P_n+\widehat{a}{P}_n(1-{\displaystyle \frac{P_n}{R_n}}),\hfill \\ \hfill R_{n+1}& =R_n+c{R}_n(1-{\displaystyle \frac{R_n}{K}})-\widehat{h}{P}_n,\hfill \end{array}$$

(24)

Model (24) at the positive fixed point $E_2$ has two eigenvalues: ${\lambda }_{1,2}={\displaystyle \frac{-1\pm \sqrt{3}i}{2}}.$ Let $\widehat{u}\left(n\right)=P\left(n\right)-P^{*}$, $\widehat{v}\left(n\right)=R\left(n\right)-R^{*}$. Then, Model (24) is transformed to

$$\begin{array}{cc}\hfill \widehat{u}(n+1)& =\left(1+\widehat{a}-{\displaystyle \frac{2\widehat{a}{P}^{*}}{R^{*}}}\right)\widehat{u}\left(n\right)+{\displaystyle \frac{\widehat{a}{P^{*}}^2}{{R^{*}}^2}}\widehat{v}\left(n\right)-{\displaystyle \frac{\widehat{a}}{R^{*}}}\widehat{u^2}\left(n\right)+{\displaystyle \frac{2\widehat{a}{P}^{*}}{{R^{*}}^2}}\widehat{u}\left(n\right)\widehat{v}\left(n\right)-{\displaystyle \frac{\widehat{a}{P^{*}}^2}{{R^{*}}^3}}\widehat{v^2}\left(n\right)+O\left(\right(|\widehat{u}\left(n\right)|+|\widehat{v}{\left(n\right)\left|\right)}^3),\hfill \\ \hfill \widehat{v}(n+1)& =-\widehat{h}\widehat{u}\left(n\right)+\left(1+c-{\displaystyle \frac{2c{R}^{*}}{K}}\right)\widehat{v}\left(n\right)-{\displaystyle \frac{c}{K}}\widehat{v^2}\left(n\right).\hfill \end{array}$$

(25)

The Jacobian matrix of Model (24) at $E_2$ is

$$\begin{array}{c}\hfill A(\widehat{a},\widehat{h})=\left(\begin{array}{cc}1+\widehat{a}-{\displaystyle \frac{2\widehat{a}{P}^{*}}{R^{*}}}& {\displaystyle \frac{\widehat{a}{P^{*}}^2}{{R^{*}}^2}}\\ -\widehat{h}& 1+c-{\displaystyle \frac{2c{R}^{*}}{K}}\end{array}\right).\end{array}$$

Also, we compute a pair of adjoint eigenvectors $q(\widehat{a},\widehat{h}),p(\widehat{a},\widehat{h})\in {\mathbb{C}}^2$ to achieve the following relations:

$$\begin{array}{c}\hfill A(\widehat{a},\widehat{h})q(\widehat{a},\widehat{h})={\displaystyle \frac{-1+\sqrt{3}i}{2}}q(\widehat{a},\widehat{h}),A^T(\widehat{a},\widehat{h})p(\widehat{a},\widehat{h})=-{\displaystyle \frac{1+\sqrt{3}i}{2}}p(\widehat{a},\widehat{h}),\end{array}$$

$<p(\widehat{a},\widehat{h}),q(\widehat{a},\widehat{h})>=1,$ where $<.,.>$ means the standard scalar product in ${\mathbb{C}}^2$:$<p,q>=\overline{p_1}{q}_1+\overline{p_2}{q}_2.$

After some calculations, we can choose

$$\begin{array}{c}\hfill q(\widehat{a},\widehat{h})=\left(\begin{array}{c}-{\displaystyle \frac{a{P^{*}}^2}{{R^{*}}^2}}\\ a-{\displaystyle \frac{2a{P}^{*}}{R^{*}}}+{\displaystyle \frac{3-\sqrt{3}i}{2}}\end{array}\right),\end{array}$$

$$\begin{array}{c}\hfill p(\widehat{a},\widehat{h})=\left(\begin{array}{c}-{\displaystyle \frac{{R^{*}}^2}{6a{P^{*}}^2}}(\sqrt{3}i(2a+3)-3)\\ {\displaystyle \frac{-\sqrt{3}i}{3}}\end{array}\right),\end{array}$$

that, for any vector $x\in {\mathbb{R}}^2$, can give the expression

$$x=zq(\widehat{a},\widehat{h})+\overline{z}\overline{q(\widehat{a},\widehat{h})},z\in \mathbb{C}.$$

Map (25) can be transformed into the complex form

$$\begin{array}{c}\hfill z\to {\displaystyle \frac{\sqrt{3}i-1}{2}}z+\sum _{2\le l+k\le 3}{\displaystyle \frac{1}{k!l!}}{g}_{kl}{z}^k{\overline{z}}^l,\end{array}$$

(26)

where

$$\left\{\begin{array}{c}{g}_{20}(a,h)={\displaystyle \frac{1}{12{R^{*}}^5}}[{\displaystyle \frac{\sqrt{3}ic{R}^{*}{[(-3+\sqrt{3}i){R^{*}}^2+2a({P^{*}}^2-{R^{*}}^2)]}^2}{K}}+(3i-3\sqrt{3}i+2\sqrt{3}a)(3i+3\sqrt{3}){R^{*}}^4\hfill \\ +2a{R^{*}}^2(3i+3\sqrt{3})(-{P^{*}}^2+P^{*}{R}^{*}+{R^{*}}^2)+2i{a}^2{(-{P^{*}}^2+P^{*}{R}^{*}+{R^{*}}^2)}^2],\hfill \\ g_{11}(a,h)={\displaystyle \frac{1}{6{R^{*}}^5}}[{\displaystyle \frac{\sqrt{3}ic{R}^{*}{[(-3+\sqrt{3}i){R^{*}}^2+2a({P^{*}}^2-{R^{*}}^2)]}^2}{K}}+(3i-3\sqrt{3}i+2\sqrt{3}a)(3i+3\sqrt{3}){R^{*}}^4\hfill \\ +2a{R^{*}}^2(3i+\sqrt{3})(-{P^{*}}^2+P^{*}{R}^{*}+{R^{*}}^2)+2i{a}^2{(-{P^{*}}^2+P^{*}{R}^{*}+{R^{*}}^2)}^2],\hfill \\ g_{02}(a,h)={\displaystyle \frac{1}{6{R^{*}}^5}}[{\displaystyle \frac{\sqrt{3}ic{R}^{*}{[(-3+\sqrt{3}i){R^{*}}^2+2a({P^{*}}^2-{R^{*}}^2)]}^2}{K}}+(3i-3\sqrt{3}i+2\sqrt{3}a)(3i+3\sqrt{3}){R^{*}}^4\hfill \\ +2a{R^{*}}^2(3i+\sqrt{3})(-{P^{*}}^2+P^{*}{R}^{*}+{R^{*}}^2)+2i{a}^2{(-{P^{*}}^2+P^{*}{R}^{*}+{R^{*}}^2)}^2]\hfill \\ g_{30}(a,h)=g_{12}(a,h)=g_{21}(a,h)=g_{03}(a,h)=0.\hfill \end{array}\right.$$

Now, to remove all quadratic terms, introduce the following transformation:

$$\begin{array}{c}\hfill z=\omega +{\displaystyle \frac{1}{2}}{h}_{20}{\omega }^2+h_{11}\omega \overline{\omega }+{\displaystyle \frac{1}{2}}{h}_{02}{\overline{\omega }}^2,\end{array}$$

(27)

by using Transformation (27) and its inverse transformation. Map (26) changed into

$$\begin{array}{c}\hfill \omega \to {\displaystyle \frac{\sqrt{3}i-1}{2}}\omega +\sum _{2\le l+k\le 3}{\displaystyle \frac{1}{k!l!}}{\varrho }_{kl}{\omega }^k{\overline{\omega }}^l+{O\left(\right|\omega |}^4),\end{array}$$

(28)

where

$${\varrho }_{20}=g_{20}+\sqrt{3}i{h}_{20},{\varrho }_{11}=2g_{11}+(\sqrt{3}i-3)h_{11}and{\varrho }_{02}=g_{02}.$$

By setting

$$h_{20}={\displaystyle \frac{\sqrt{3}i}{3}}{g}_{20},h_{11}={\displaystyle \frac{3+\sqrt{3}i}{6}}{g}_{11}and{h}_{02}=0.$$

according to the previous values of $\varrho$ and h, we deduce that ${\varrho }_{20}={\varrho }_{11}=0$ and ${\varrho }_{02}=g_{02}$ and ${\varrho }_{30},{\varrho }_{21},{\varrho }_{12},{\varrho }_{03}$ can be simplified as follows:

$$\left\{\begin{array}{c}{\varrho }_{30}={\displaystyle \frac{3-\sqrt{3}i}{2}}{g}_{11}{\overline{g}}_{02}+\sqrt{3}i{g}_{20}^2+g_{30},\hfill \\ {\varrho }_{21}={\displaystyle \frac{3+2\sqrt{3}i}{2}}{g}_{11}{g}_{20}+{\displaystyle \frac{3-\sqrt{3}i}{3}}{|g_{11}|}^2+g_{21},\hfill \\ {\varrho }_{12}={\displaystyle \frac{3+\sqrt{3}i}{6}}{g}_{20}{g}_{02}+{\displaystyle \frac{3-\sqrt{3}i}{3}}{\overline{g}}_{11}{g}_{02}+{\displaystyle \frac{3+\sqrt{3}i}{3}}{g}_{11}^2-{\displaystyle \frac{\sqrt{3}i}{3}}{\overline{g}}_{20}{g}_{11},\hfill \\ {\varrho }_{03}=\sqrt{3}i{g}_{11}{g}_{02}-\sqrt{3}i{g}_{02}{\overline{g}}_{20}+g_{03}.\hfill \end{array}\right.$$

To simplify Equation (28), we introduce the following transformation:

$$\begin{array}{c}\hfill \omega =\xi +{\displaystyle \frac{1}{6}}{h}_{30}{\xi }^3+{\displaystyle \frac{1}{2}}{h}_{21}{\xi }^2\overline{\xi }+{\displaystyle \frac{1}{2}}{h}_{12}\xi {\overline{\xi }}^2+{\displaystyle \frac{1}{6}}{h}_{03}{\overline{\xi }}^3.\end{array}$$

(29)

After using Equation (29) and its inverse, Map (28) is formed as

$$\begin{array}{c}\hfill \xi ={\displaystyle \frac{\sqrt{3}i-1}{2}}\xi +{\displaystyle \frac{1}{2}}{g}_{02}{\overline{\xi }}^2+\sum _{l+k=3}{\displaystyle \frac{1}{k!l!}}{\tilde{\varrho }}_{kl}{\xi }^k{\overline{\xi }}^l+{O\left(\right|\xi |}^4),\end{array}$$

(30)

where

$$\begin{array}{cc}\hfill {\tilde{\varrho }}_{30}& ={\varrho }_{30}+{\displaystyle \frac{\sqrt{3}i-3}{2}}{h}_{30},{\tilde{\varrho }}_{21}={\varrho }_{21},\hfill \\ \hfill {\tilde{\varrho }}_{12}& ={\varrho }_{12}+\sqrt{3}i{h}_{12}and{\tilde{\varrho }}_{03}={\varrho }_{03}+{\displaystyle \frac{\sqrt{3}i-3}{2}}{h}_{03}.\hfill \end{array}$$

By putting

$$\begin{array}{c}\hfill h_{30}={\displaystyle \frac{3+\sqrt{3}i}{6}}{\varrho }_{30},h_{12}={\displaystyle \frac{\sqrt{3}i}{3}}{\varrho }_{12},h_{03}={\displaystyle \frac{3+\sqrt{3}i}{6}}{\varrho }_{03}and{h}_{21}=0,\end{array}$$

we then have ${\tilde{\varrho }}_{30}={\tilde{\varrho }}_{21}={\tilde{\varrho }}_{03}=0.$

Finally, we introduce the normal form of the bifurcation with 1:3 resonance as follows:

$$\begin{array}{c}\hfill \xi \to {\displaystyle \frac{\sqrt{3}i-1}{2}}\xi +B(\widehat{a},\widehat{h})\overline{\xi }+C(\widehat{a},\widehat{h})\xi {\left|\xi \right|}^2+{O\left(\right|\xi |}^4),\end{array}$$

(31)

where

$$\begin{array}{c}\hfill B(\widehat{a},\widehat{h})={\displaystyle \frac{1}{2}}{g}_{02}(a,h),C(\widehat{a},\widehat{h})={\displaystyle \frac{(3+\sqrt{3}i)g_{20}{g}_{11}}{6}}+{\displaystyle \frac{(3-\sqrt{3}i){|g_{11}|}^2}{6}}+{\displaystyle \frac{1}{2}}{g}_{21}.\end{array}$$

If

$$\begin{array}{cc}\hfill B_1(\widehat{a},\widehat{h})& ={\displaystyle \frac{-3}{2}}(\sqrt{3}i+1)B(\widehat{a},\widehat{h}),\hfill \\ \hfill C_1(\widehat{a},\widehat{h})& =-3{\left|B\right(\widehat{a},\widehat{h}\left)\right|}^2-{\displaystyle \frac{3}{2}}(\sqrt{3}i+1)C(\widehat{a},\widehat{h}).\hfill \end{array}$$

By Lemma (9.11) and Lemma (9.13) in [28], we can obtain the following theorem.

Theorem 6. 

If $B_1(\widehat{a},\widehat{h})\ne 0$ and $Re\left(C_1(\widehat{a},\widehat{h})\right)\ne 0$, then Map (4) has the the following complex dynamical behaviors at the fixed point $E^{*}$:

(a) 

There is a Neimark–Sacker bifurcation at the trivial fixed point $E_1$ of map (31);

(b) 

There is a saddle cycle of period three, corresponding to the saddle fixed point $E_k(k=1,2,3)$ of map (31);

(c) 

There is a homoclinic structure formed by the stable and unstable invariant manifolds of the period-three cycle, intersecting transversally in an exponentially narrow parameter region.

Remark 4. 

In the 1:3 resonance scenario, the meeting of stable and unstable manifolds creates an infinite number of orbits with a period of three, leading to a homoclinic tangency. This reveals that a period-three cycle can cause chaos. Biologically, this means that the population–resource system may experience periodic or quasi-periodic fluctuations due to the non-degenerate Neimark–Sacker bifurcation. The period-three fluctuations, resulting from a saddle cycle of period three, can generate chaotic sets. These chaotic sets contribute to long-term fluctuations, population explosions, and overall chaos, all due to the presence of the homoclinic structure.

4.3. The 1:4 Resonance of Model (4) at the Positive Fixed Point $E_2(P^{*},R^{*})$

Here, we clear that the 1:4 resonance of Model (4) at the positive fixed point $E_3(P^{*},R^{*})$ when parameters a and h vary within a small neighborhood of $F_{14}$. Also, we choose parameters $\left(\tilde{a},c,\tilde{h},K\right)$ arbitrarily from $F_{14}$, taking $\tilde{a}$ and $\tilde{h}$ as bifurcation parameters.

Model (4) with $\left(\tilde{a},c,\tilde{h},K\right)$ is given by

$$\begin{array}{cc}\hfill P_{n+1}& =P_n+\tilde{a}{P}_n(1-{\displaystyle \frac{P_n}{R_n}}),\hfill \\ \hfill R_{n+1}& =R_n+c{R}_n(1-{\displaystyle \frac{R_n}{K}})-\tilde{h}{P}_n,\hfill \end{array}$$

(32)

where $|a-\tilde{a}|$, $|h-\tilde{h}|$$\ll 1$, and the eigenvalues of Model (32) at the positive point are ${\lambda }_{1,2}=\pm i$.

Let $\tilde{u}=P-P^{*}$ and $\tilde{v}=R-R^{*}$. Then, we turn the point $E_3(P^{*},R^{*})$ into the origin point.

Map (32) becomes

$$\begin{array}{cc}\hfill \tilde{u}(n+1)& =\left(1+\tilde{a}-{\displaystyle \frac{2\tilde{a}{P}^{*}}{R^{*}}}\right)\tilde{u}\left(n\right)+{\displaystyle \frac{\tilde{a}{P^{*}}^2}{{R^{*}}^2}}\tilde{v}\left(n\right)-{\displaystyle \frac{\tilde{a}}{R^{*}}}\widetilde{u^2}+{\displaystyle \frac{2\tilde{a}{P}^{*}}{{R^{*}}^2}}\tilde{u}\left(n\right)\tilde{v}\left(n\right)-{\displaystyle \frac{\tilde{a}{P^{*}}^2}{{R^{*}}^3}}\widetilde{v^2}\left(n\right)+O\left(\right(|\tilde{u}\left(n\right)|+|\tilde{v}{\left(n\right)\left|\right)}^3),\hfill \\ \hfill \tilde{v}(n+1)& =-\tilde{h}\tilde{u}\left(n\right)+\left(1+c-{\displaystyle \frac{2c{R}^{*}}{K}}\right)\tilde{v}\left(n\right)-{\displaystyle \frac{c}{K}}\widetilde{v^2}\left(n\right).\hfill \end{array}$$

(33)

At $E_2(\tilde{a},\tilde{h})$, the Jacobian matrix of Model (32) is given by

$$\begin{array}{c}\hfill A(\tilde{a},\tilde{h})=\left(\begin{array}{cc}1+\tilde{a}-{\displaystyle \frac{2\tilde{a}{P}^{*}}{R^{*}}}& {\displaystyle \frac{\tilde{a}{P^{*}}^2}{{R^{*}}^2}}\\ -\tilde{h}& 1+c-{\displaystyle \frac{2c{R}^{*}}{K}}\end{array}\right).\end{array}$$

Furthermore, we can obtain a pair of adjoint eigenvectors $q(\tilde{a},\tilde{h}),p(\tilde{a},\tilde{h})\in {\mathbb{C}}^2$, such that $A(\tilde{a},\tilde{h})q(\tilde{a},\tilde{h})=iq(\tilde{a},\tilde{h})$, $A^T(\tilde{a},\tilde{h})p(\tilde{a},\tilde{h})=-ip(\tilde{a},\tilde{h})$, $<p(\tilde{a},\tilde{h}),q(\tilde{a},\tilde{h})>=1.$ According to the previous relations, we have

$$\begin{array}{c}\hfill q(\tilde{a},\tilde{h})=\left(\begin{array}{c}-{\displaystyle \frac{a{P^{*}}^2}{{R^{*}}^2}}\\ 1+a-{\displaystyle \frac{2a{P}^{*}}{R^{*}}}-i\end{array}\right),\end{array}$$

$$\begin{array}{c}\hfill p(\tilde{a},\tilde{h})=\left(\begin{array}{c}{\displaystyle \frac{{R^{*}}^2}{2a{P^{*}}^2}}(i(1+a)-1)\\ -{\displaystyle \frac{i}{2}}\end{array}\right),\end{array}$$

For any vector ${\mathrm{Y}}_{\mathrm{n}}={(\tilde{P},\tilde{R})}^T\in {\mathbb{R}}^2$, this can be written as

$${\mathrm{Y}}_{\mathrm{n}}=zq(\tilde{a},\tilde{h})+\overline{z}\overline{q}(\tilde{a},\tilde{h}),z\in \mathbb{C}.$$

Then, Model (33) transforms to the following complex relation:

$$\begin{array}{c}\hfill z⟼iz+\sum _{2\le k+l\le 3}{\displaystyle \frac{1}{k!l!}}{g}_{kl}{z}^k{\overline{z}}^l,\end{array}$$

(34)

where

$$\left\{\begin{array}{c}{g}_{20}={\displaystyle \frac{-i}{2{R^{*}}^3}}[{\displaystyle \frac{-c{R}^{*}{(2+a(P^{*}-R^{*})-(1-i)R^{*})}^2}{K}}+((1+i)+a){(-2+((1-i)+a)R^{*})}^2],\hfill \\ g_{11}={\displaystyle \frac{-i}{{R^{*}}^3}}[{\displaystyle \frac{-c{R}^{*}{(2+a(P^{*}-R^{*})-(1-i)R^{*})}^2}{K}}+((1+i)+a){(-2+((1-i)+a)R^{*})}^2],\hfill \\ g_{02}={\displaystyle \frac{-i}{2{R^{*}}^3}}[{\displaystyle \frac{-c{R}^{*}{(2+a(P^{*}-R^{*})-(1-i)R^{*})}^2}{K}}+((1+i)+a){(-2+((1-i)+a)R^{*})}^2],\hfill \\ g_{30}=g_{12}=g_{21}=g_{03}=0.\hfill \end{array}\right.$$

To remove the quadratic terms, we introduce the following equation:

$$\begin{array}{c}\hfill z=\omega +{\displaystyle \frac{1}{2}}{h}_{20}{\omega }^2+h_{11}\omega \overline{\omega }+{\displaystyle \frac{1}{2}}{h}_{02}{\overline{\omega }}^2,\end{array}$$

(35)

where coefficients $h_{kl}$ with $k+l=2$ are presented in the following relation.

By using Transformation (35) and its inverse transformation, Map (34) can be turned into the form

$$\begin{array}{c}\hfill \omega ⟼i\omega +\sum _{2\le +l\le 3}{\displaystyle \frac{1}{k!l!}}{\varrho }_{kl}{\omega }^k{\overline{\omega }}^l+O\left(\right(\left|\omega \right|+|\overline{\omega }{\left|\right)}^4),\end{array}$$

(36)

where

$$\left\{\begin{array}{c}{\varrho }_{20}={\displaystyle \frac{1}{2}}(g_{20}+(\lambda -{\lambda }^2)h_{20}),{\varrho }_{11}=(g_{11}{+(\lambda -|\lambda |}^2\left)h_{11}\right),\hfill \\ {\varrho }_{20}={\displaystyle \frac{1}{2}}(g_{02}+(\lambda -{\overline{\lambda }}^2)h_{02}),{\varrho }_{30}={\displaystyle \frac{1}{6}}(g_{30}+(\lambda -{\lambda }^2)h_{30}),\hfill \\ {\varrho }_{21}={\displaystyle \frac{1}{2}}(g_{21}+(\lambda -\lambda {\left|\lambda \right|}^2)h_{21}),{\varrho }_{12}={\displaystyle \frac{1}{2}}(g_{12}+(\lambda -\overline{\lambda }{\left|\lambda \right|}^2)h_{12}),\hfill \\ {\varrho }_{03}={\displaystyle \frac{1}{6}}(g_{03}+(\lambda -{\overline{\lambda }}^3)h_{03}).\hfill \end{array}\right.$$

By setting

$$h_{20}={\displaystyle \frac{1}{2}}{g}_{20}(i-1),h_{11}={\displaystyle \frac{1}{2}}{g}_{11}(i+1),h_{02}={\displaystyle \frac{1}{2}}{g}_{02}(i+1),$$

then, we obtain ${\varrho }_{20}={\varrho }_{11}={\varrho }_{02}=0.$

To eliminate some cubic terms, we introduce the following transformation:

$$\begin{array}{c}\hfill \omega =\xi +{\displaystyle \frac{1}{6}}{h}_{30}{\xi }^3+{\displaystyle \frac{1}{2}}{h}_{21}{\xi }^2\overline{\xi }+{\displaystyle \frac{1}{2}}{h}_{12}\xi {\overline{\xi }}^2+{\displaystyle \frac{1}{3}}{h}_{03}{\overline{\xi }}^3.\end{array}$$

(37)

By using Equation (37) and its inverse transformation, we obtain

$$\begin{array}{c}\hfill \omega =i\xi +\sum _{n+l=3}{\displaystyle \frac{1}{k!l!}}{\overline{\varrho }}_{kl}{\xi }^k{\overline{\xi }}^l+O{\left(\right|\omega \left|\right)}^4,\end{array}$$

(38)

where

$${\overline{\varrho }}_{30}={\varrho }_{30}+2i{h}_{30},{\overline{\varrho }}_{21}={\varrho }_{21},$$

$${\overline{\varrho }}_{12}={\varrho }_{12}+2i{h}_{12},{\overline{\varrho }}_{03}={\varrho }_{03}.$$

Let

$$h_{30}={\displaystyle \frac{i}{2}}{\varrho }_{30},h_{12}={\displaystyle \frac{i}{2}}{\varrho }_{12},h_{03}=h_{21}=0,$$

from the previous relation, we obtain ${\overline{\varrho }}_{30}={\overline{\varrho }}_{12}=0.$

Finally, near 1:4 resonance, Model (38) can be changed into

$$\begin{array}{c}\hfill \xi \to i\xi +C(\overline{a},\overline{h}){\xi \left|\xi \right|}^2+C(\overline{a},\overline{h})|\overline{\xi }{|}^3+O{\left(\right|\xi \left|\right)}^4,\end{array}$$

(39)

where

$$\begin{array}{cc}\hfill C(\overline{a},\overline{h})& =i{g}_{11}^2-{\displaystyle \frac{1}{2}}{g}_{11}\overline{g_{20}}(1+i)+{\overline{g}}_{11}{g}_{02}+g_{02}{g}_{11}(i-1)-{\displaystyle \frac{1}{2}}{g}_{11}{g}_{20}(1-2i),\hfill \\ \hfill D(\overline{a},\overline{h})& ={\displaystyle \frac{i-1}{4}}{g}_{02}{g}_{11}-{\displaystyle \frac{i+1}{4}}{g}_{20}{g}_{11}+{\displaystyle \frac{1}{6}}{g}_{03}.\hfill \end{array}$$

Let $C_1(\overline{a},\overline{h})=-4iC(\overline{a},\overline{h})$, $D_1(\overline{a},\overline{h})=-4iD(\overline{a},\overline{h})$. At $D_1\ne 0$, we conclude that $A(\overline{a},\overline{h})={\displaystyle \frac{C_1(\overline{a},\overline{h})}{|D_1(\overline{a},\overline{h})|}}.$

Theorem 7. 

Assume that $(\overline{a},\overline{h})\in F_{14}$. If $ReA(\overline{a},\overline{h})\ne 0$ and $ImA(\overline{a},\overline{h})\ne 0$, then Model (4) has many dynamical behaviors (see [28]):

• There is a Neimark–Sacker bifurcation curve at trivial fixed point $E_1$ in (39). Also, if $\lambda =-i$, there is an invariant circle; at $\lambda =i$, the invariant circle disappears,
• If $\left|A\right|>1$, System (39) has eight non-trivial equilibrium points $S_k,E_k,=1,2,3,4$. The eight non-trivial fixed points appear or disappear as the fold bifurcation at the corresponding parameter values;
• There are Neimark–Sacker bifurcations at fixed points $E_k,=1,2,3,4$.

Remark 5. 

The presence of 1:4 resonance signifies the existence of a nondegenerate Neimark–Sacker bifurcation, allowing for the formation of an invariant cycle with a period of four in a specific parameter range. In biological systems, the nondegenerate Neimark–Sacker bifurcation can lead to periodic or quasiperiodic fluctuations in the population–resource system. Moreover, the presence of an invariant cycle with a period of four indicates that a stable state of the population–resource system would transition into a state that repeats (almost) after every four time intervals.

5. Numerical Simulations

In this section, we conduct numerical simulations to validate the theoretical analysis proposed earlier for System (4). We utilize various techniques, including bifurcation diagrams, phase portraits, and maximal Lyapunov exponents (MLEs) with specific parameter values. Moreover, we investigate novel and complex dynamics within the system.

Fix $K=5$, $c=1.3$, and $h=0.5$; then, vary the value of parameter a. We choose the initial value of Model (4) to be $(P_0,R_0)=(0.5,0.5)$ in all our numerical simulations. Based on the theoretical analysis presented in Section 3, a flip bifurcation occurs when $a=2.833$, as depicted in Figure 1a. To further investigate the dynamics, we calculate the maximal Lyapunov exponents (MLEs) corresponding to Figure 1a and plot them in Figure 1b. It can be observed that MLEs are negative for $a<3.33$, suggesting stability within the system. However, a small portion of MLEs are positive for $a>3.33$, indicating a cascade of period-doubling bifurcations in System (4). Let us consider the parameter values $a=0.3$, $c=3$, and $h=1$ in System (4). In this case, the bifurcation diagram in Figure 2a reveals that the system starts with a unique point, which then undergoes a period-two cycle bifurcation. This cycle loses stability and transforms into a period-four cycle. As the bifurcation progresses, a period doubling occurs, leading to the emergence of chaotic behavior. The corresponding maximal Lyapunov exponents in Figure 2b provide further insight into the system’s dynamics.

We consider a new set of parameters: $a=1$, $h=1$, and different values of c. In Figure 3a, Figure 4a, Figure 5a and Figure 6a, we observe the Neimark–Sacker bifurcation occurring in the $(a,P)$ plane. This indicates that System (4) loses its stability through Neimark–Sacker bifurcation. To gain further insights into the dynamics, we examine the corresponding maximal Lyapunov exponents in Figure 3b, Figure 4b, Figure 5b and Figure 6b. We select a new parameter set with $a=1$, $h=1$, and various values of c. In Figure 3a, Figure 4a, Figure 5a and Figure 6a, the Neimark–Sacker bifurcation is observed in the $(a,P)$ plane, indicating the loss of stability in System (4). To analyze the system’s behavior further, we examine the corresponding maximal Lyapunov exponents in Figure 3b, Figure 4b, Figure 5b and Figure 6b.

From Figure 1, Figure 2, Figure 3, Figure 4, Figure 5 and Figure 6, when the rate at which resources grow, represented as $c,$ is lower than the rate at which the population grows, represented as $a,$ increasing the harvesting rate will cause both the resource and population to vanish together, as the equilibrium point moves towards the origin. Subcritical flip bifurcation occurs when the rate of resource growth reduces to the rate of growth in the population. Also, Neimark–Sacker bifurcation occurs when the rate of resource growth exceeds the rate of growth in population, and the equilibrium point is shifted away from the bifurcation as the harvesting rate is increased.

Figure 7 displays the phase portraits of the system (4) near the Neimark-Sacker with a = 0.0045, c = 3 and different values of h, which corresponds to the scenario depicted in Figure 1a. Figure 8 displays the chaotic attractors of the system (4) with a = 0.3, c = 3, and varying values of h, which corresponds to the scenario depicted in Figure 2a. When we set a = 1, c = 3:5, and different values of h, Figure 9a–e represent the attractors for system (4). These figures are associated with the parameter configuration shown in Figure 6a. Additionally, Figure 9f illustrates the strange attractor observed when h = 1:5. From Figure 7, Figure 8 and Figure 9, while increasing the parameter “h”, a thought-provoking inquiry emerges concerning the closed invariant curve’s vanishing point. As “h” continues to increase, gradual cusps appear along the closed curve. Additionally, a secondary Neimark-Sacker bifurcation occurs within the confines of this curve, encompassing a periodic orbit. In excess of a particular threshold value of the parameter, numerical investigations indicate the emergence of chaotic behavior. As we further explore the dynamics, a secondary Neimark-Sacker bifurcation occurs, followed by the emergence of chaotic behavior beyond a specific threshold value of the parameter, as indicated by numerical studies. The intricate nature of this process can be elucidated through an examination of codimension two bifurcations, revealing additional details about the system’s behavior.

Numerical Continuation for Strong Resonances

In this subsection, numerical continuation simulations are provided to illustrate the analytical findings mentioned above and for analyzing strong resonances using the MATLAB package MATCONTM [29]. This analysis utilizes continuation methods to trace the solution manifolds of fixed points while varying specific parameters of the map. By employing this approach, we observe that, in the two-parameter configuration, the boundaries of the stabilizing domains for the cycles are typically represented by bifurcation curves obtained through the MATCONTM software. These bifurcation curves need to be computed using numerical continuation methods, as they cannot be obtained analytically.

By a continuation of the fixed point $E_2=(6.60846,6.60846)$ and putting the values of parameters as $a=1.5,c=2.2$, $K=10$, and h as free and changing from 0 to 1, we note that the positive fixed point $E_2$ is stable when $0<h<0.04$. It loses stability via a subcritical period doubling point (detected as PD) and the corresponding normal form coefficient is $-1.419573\times {10}^{-2}<0$. Furthermore, the fixed point $E_2$ is denoted as a subcritical Neimark–Sacker bifurcation; it is detected as NS when $h=0.8$, and the corresponding normal-form coefficient is $-1.999130\times {10}^{-2}<0$. Figure 10 visually depicts. The MATCONTM report is given by

label=PD, x=(9.818182 9.818182 0.040000)

normal form coefficient of NS=-1.999130e-02

label=NS, x=(6.363636 6.363636 0.800000)

normal form coefficient of NS=-1.999130e-02

For $h>0.04$, stable 4-cycle and 8-cycle versions are presented in Figure 11.

Furthermore, we calculate the Neimark–Sacker curve by varying the free parameters a and h, starting from the Neimark–Sacker (NS) point. The MATCONTM report for this computation is presented below:

label=R3, x=(7.890409 7.890409 1.728220 0.464110 -0.500000)

Normal form coefficient of R3: Re(c_1)=-6.470793e-01

label=R2, x=(9.929222 9.929222 1.831142 0.015571 -1.000000)

Normal form coefficient of R2:[c, d]=3.105740e-02,-3.219911e-01

label=R4, x=(6.227983 6.227983 1.459688 0.829844 0.000000)

Normal form coefficient of R4: A =-9.464392e-01-5.573169e-01i

Based on the previous report generated by MATCONTM, we have identified several codimension-2 bifurcations occurring along the Neimark–Sacker curve. These include 1:2 resonance (R2), 1:3 resonance (R3), and 1:4 resonance (R4). Figure 12 visually depicts these findings.

Here, we present the codimension-two bifurcation analysis conducted from the previously identified PD point to determine the period-doubling curve. The analysis is performed at specific parameter values, namely, $c=2.2$ and $K=10$, while considering a and h as free parameters. Figure 13 visually depicts. The MATCONTM report for this analysis is provided below:

label=R2, x=(9.929222 9.929222 1.831142 0.015571 -1.000000)

Normal form coefficient of R2:[c, d]=3.105740e-02,-3.219911e-01

label=GPD, x=(9.761393 9.761393 1.288254 0.052494)

Normal form coefficient of GPD =1.969966e-03

label=LPPD, x=(9.545455 9.545455 -0.000000 0.100000)

Normal form coefficient for LPPD : [a/e, be]=1.00801e-7, 8.985363e-11,

First Lyapunov coefficient for second iterate =8.985363e-11,

Based on the MATCONTM report, our observations indicate that the period-doubling curve exhibits various bifurcation phenomena, including 1:2 resonance (R2), generalized period doubling (GPD), and a fold-flip (LPPD). These findings are illustrated in Figure 14.

6. Conclusions

In this study, we have conducted a comprehensive analysis of the dynamical behavior of a discrete population model for human population and its resource, represented by Model (4). We have performed an elaborate analysis of the existence and stability of model fixed points. For the positive fixed point, we have examined how it undergoes several codimension-one bifurcations, including flip bifurcation, Neimark–Sacker bifurcation, and the emergence of chaotic attractors. Moreover, we have provided a thorough examination of codimension-two bifurcations related to $1:2$, $1:3$, and $1:4$ resonance and their characteristics. To achieve this objective, a successful technique was implemented. Specifically, we have analyzed the dynamics of Model (4) using the normal-form approach. The process of normalizing a model involves distilling it to its essential components. Utilizing the normal form, we identified the criteria that govern the occurrence of subcritical or supercritical bifurcations. To further support the complexity of Model (4) dynamics, we have numerically computed the maximal Lyapunov exponents. Finally, numerical simulations were conducted to validate our analytical findings, ensuring their consistency with the actual behavior of the system.

The authors in [19] contend that the renewable rate of resources could serve as either a stabilizing or destabilizing factor in Model (3). However, our theoretical analysis indicates that the impact of the renewable rate may vary, potentially stabilizing or destabilizing Model (4), contingent upon the values of other model parameters. Even with substantial increases in the renewable rate of resource, the human population may not gain advantages, as this would result in a corresponding rise in population due to the availability of harvested resources. As the human population surges, resource depletion may lead to extinction, which would have dire consequences for the human population. Moreover, according to the Poincaré–Bendixon theorem [30], two-dimensional BR continuous-time Model (3) exhibits either stable coexistence or oscillations. The continuous-time BR system in (3) exhibits no additional complicated dynamics or multistability. Our discrete-time Model (4) exhibited complex dynamical behaviors, including periodicity, quasiperiodicity, and chaos. The bifurcation diagram demonstrated the presence of periodic bubbling and periodic windows, leading to chaotic behavior. The maximal Lyapunov exponents validated the existence of non-periodic dynamics in the system. These behaviors highlight that, when population growth and resource growth coexist, they give rise to highly intricate patterns.

Fractional-order predator–prey models incorporate memory effects and hereditary features, making them more suitable for simulating real-world phenomena where previous states affect current dynamics. For a fractional-order version of Model (4), fractional derivative memory effects can significantly affect dynamics, bifurcation thresholds, stability features, and oscillatory behaviors. It would be interesting to tackle such kinds of modeling in future work. Contemporary research continues to advance our understanding of these dynamics by employing more sophisticated mathematical models and advanced data analysis techniques. Factors such as ecological changes and climatic variations on the ecosystem are also being considered in these investigations. The implications of these findings extend to various fields, including economics, biology, and ecology, where they can find practical applications.