Bifurcation analysis and chaos control in discrete chemical systems: Co-Dimension one and two insights from the Degn-Harrison model

Abstract

This research offers a well-developed analysis of the nonlinear dynamic behavior of a discrete-step Degn-Harrison system without diffusion. The continuous-time chemical reaction system is discretized by employing the explicit Euler scheme which results in a two-dimensional nonlinear map whose behavior depends on both the reaction parameters and the discretization step size. The presence and stability conditions of the interior equilibrium point are first established, and then the system’s bifurcation structure is explored in depth. Special attention is paid to codimension-one and codimension-two bifurcations, period-doubling, Neimark-Sacker and strong resonance cases (1:2, 1:3, and 1:4). To give a comprehensive characterization of these phenomena, we get the normal forms of these phenomena and check the non-degeneracy conditions of such phenomena, such that the validity of the theoretical results is ensured. Numerical simulations ensure the theoretical findings, such as bifurcation diagrams, phase portraits, and the maximum Lyapunov exponent, to give a picture of the switching from regular to chaotic dynamics. The outcomes of this study not only deepen the understanding of discretization-induced complexities in chemical kinetics as well as offering a structure for analyzing similar nonlinear systems within applied sciences.

Introduction

Nonlinear chemical oscillators have been used as paradigms of modeling complex dynamics in biological and chemical systems, since Turing first formulated his work on morphogenesis. ¹ The Degn-Harrison reaction model is one of the most insightful models. This reaction mechanism was first introduced in the context of oscillatory respiratory behavior in Klebsiella aerogenes, ² and this mechanism has been shown to embody the principles of nonlinear inhibitory dynamics in a biologically relevant system. The analysis of nonlinear dynamical systems has long been established as a pillar of the comprehension of interactions of complexity in biological and ecological phenomena. Over the past years, the discussion of population dynamics, predator-prey and epidemiology under both classical and fractional-order paradigms has received much attention. As an illustration, Berkal et al. ³ investigated a generalized conformable fractional SIR and showed that the stability results and chaotic transitions can be dramatically modified by the use of fractional dynamics. In another study, the Allee effect and Gompertz growth ecological complexities were also considered in predator-prey models, and the bifurcation diagrams and chaos control methods were developed in fractional and conformable derivative models by Almatrafi and Berkal.^4–6 Fractional dynamics have also been applied to other systems, where it has been observed to be useful in the study of activator inhibitor systems, where the combination of diffusion, reaction kinetics, and memory effects provides rich pattern formation and stability behavior, for example, the work by Berkal and Almatrafi. ⁷ These studies collectively illustrate the need to employ sophisticated tools of analysis, such as eigenvalues-based criteria of stability, bifurcation theory, and numerical solutions, to discover the complete spectrum of qualitative behavior of nonlinear biological models. Following this accumulating body of literature, the current paper seeks to further understand the dynamical characteristics of discrete ecological systems by examining equilibrium points, their stability and the bifurcation phenomenon thereof. For more detail see Refs.^8–11

In the Degn-Harrison model, ² the system is governed by a three-step reaction model:

$$\begin{array}{c}\multicolumn{1}{c}{A\to Y,B⟺X,X+Y\to P.}\end{array}$$

In this, X and Y signify the intermediate reactants that constitute the density of oxygen and nutrients. Components $A$ and $B$ are fixed inputs while P is the final output, whose concentration remains constant. Notably, the third step is repressed by an overabundance of oxygen, thus creating positive feedback, which is essential in the behavior of the model. The first and third reactions are considered irreversible, but the second one is reversible, which represents real biochemical reactions. Degn and Harrison suggested that the inhibition in the final step could be modeled using a nonlinear rate law of the form XY/(1+q $X^2$ ), where $q$ is the inhibition strength. When placed in a spatial domain with homogeneous Neumann boundary assumptions, the full reaction-diffusion system is expressed as:

$$\{\begin{array}{c}{X}_t-D_X\mathrm{\Delta }X=k_{2}B-k_{3}B-\frac{k_4\mathrm{XY}}{(1+{\mathrm{qX}}^2)},x\in \mathrm{\Omega },t>0,\hfill \\ \multicolumn{1}{c}{Y_t-D_Y\mathrm{\Delta }Y=k_{1}A-\frac{k_4\mathrm{XY}}{(1+{\mathrm{qX}}^2)},x\in \mathrm{\Omega },t>0,}\\ \multicolumn{1}{c}{\frac{\partial X}{\partial n}=\frac{\partial Y}{\partial n}=0,x\in \partial \mathrm{\Omega },t>0,}\end{array}$$

(1.1)

in which $A,B,X\mathrm{and}Y$ are non-dimensional concentrations, $k_i(\mathrm{for}i\in N)$ represents reaction rate constants, and $D_X$ and $D_Y$ are the respective diffusion coefficients of the intermediates. For more details, readers may refer to Refs.^12–19 The continuous-time dynamics of this system have been extensively studied using both analytical and numerical techniques, revealing a rich array of dynamical behaviors, including limit cycles, Hopf bifurcations, and excitability.^20,21 However, most of these studies have focused on the differential form of the model, whereas its discrete-time analog remains relatively underexplored.

This paper considers an analysis of a discrete dynamical system based on a reduced ordinary differential equation formulation of the model. This study is conducted through the application of a standard non-dimensionalization procedure ¹² with the following transformations:

$$\begin{array}{c}\multicolumn{1}{c}{\tau =k_{3}t,\mathcal{u}=\frac{X}{k_3},\mathcal{v}=\frac{Y}{k_3},\mathcal{a}=\frac{k_2}{k_3}B,\mathcal{b}=\frac{k_1}{k_3}A,}\\ \multicolumn{1}{c}{\mathcal{c}=k_4,1=k_3^{3}q,{\mathcal{d}}_1=\frac{1}{k_3}{D}_X,{\mathcal{d}}_2=\frac{1}{k_3}{D}_Y.}\end{array}$$

For simplicity, we continue to denote time by $t$ , system (1.1) yielding the following non-linear dimensionless PDE system:

$$\{\begin{array}{c}{\mathcal{u}}_t-{\mathcal{d}}_1\mathrm{\Delta }\mathcal{u}=\mathcal{a}-\mathcal{u}-\frac{\mathcal{c}\mathcal{uv}}{1+{\mathcal{u}}^2},x\in \mathrm{\Omega },t>0,\hfill \\ \multicolumn{1}{c}{{\mathcal{v}}_t-{\mathcal{d}}_2\mathrm{\Delta }\mathcal{v}=\mathcal{b}-\frac{\mathcal{c}\mathcal{uv}}{1+{\mathcal{u}}^2},x\in \mathrm{\Omega },t>0,}\\ \multicolumn{1}{c}{\frac{\partial \mathcal{u}}{\partial n}=\frac{\partial \mathcal{v}}{\partial n}=0,x\in \partial \mathrm{\Omega },t>0,}\\ \multicolumn{1}{c}{\mathcal{u}(x,0)={\mathcal{u}}_0\left(x\right)>0,\mathcal{v}(x,0)={\mathcal{v}}_0\left(x\right)>0,x\in \mathrm{\Omega }.}\end{array}$$

(1.2)

Focusing on a spatially homogeneous case, we simplify the system (1.2) to a two-dimensional nonlinear ODE model.

$$\begin{array}{c}{\mathcal{u}}_t=\mathcal{a}-\mathcal{u}-\frac{\mathcal{c}\mathcal{uv}}{1+{\mathcal{u}}^2},\hfill \\ \multicolumn{1}{c}{{\mathcal{v}}_t=\mathcal{b}-\frac{\mathcal{c}\mathcal{uv}}{1+{\mathcal{u}}^2},}\end{array}$$

(1.3)

which forms the basis for our discrete-time analysis. By employing the explicit Euler scheme for time discretization, we derive the subsequent discrete dynamical model:

$$\begin{array}{c}{\mathcal{u}}_{n+1}={\mathcal{u}}_n+\mathcal{h}(\mathcal{a}-{\mathcal{u}}_n-\frac{\mathcal{c}{\mathcal{u}}_n{\mathcal{v}}_n}{1+{\mathcal{u}}_n^2}),\hfill \\ \multicolumn{1}{c}{{\mathcal{v}}_{n+1}={\mathcal{v}}_n+\mathcal{h}(\mathcal{b}-\frac{\mathcal{c}{\mathcal{u}}_n{\mathcal{v}}_n}{1+{\mathcal{u}}_n^2}),}\end{array}$$

(1.4)

where $\mathcal{h}$ > 0 represents the time step size. Here, $\mathcal{a}$ and $\mathcal{b}$ are dimensionless inflow concentration of oxygen and nutrients, respectively and $\mathcal{c}$ is the dimensionless reaction rate constant for the autocatalytic or inhibitory step.

This discrete system has mainly been analyzed in this paper to study nonlinear dynamics with emphasis on fixed points, local stability, and bifurcation as parameters vary. In particular, we study codimension one bifurcations, for example, period-doubling bifurcation, Neimark-Sacker bifurcation and codimension two bifurcation, for example, strong resonances, 1:2, 1:3, and 1:4. This analysis not only closes a gap in the existing research but also provides information about the discrete versions of biologically motivated chemical systems.

The structure of the manuscript is outlined as follows.

In Section 2, we investigate the equilibrium points of the system (1.4) and analyze their local stability properties. In the third section, we examine the case of a bifurcation, codimension-one and codimension-two bifurcations, which are capable of rich dynamical transitions as a result of a parameter change or two. Section 4 presents the numerical simulations to justify and explain our analytical results.

Existence and stability of equilibrium points

In this section, we discussed the equilibrium points of model (1.4) and also assessed their stability, giving some insight into what happens in the system in the long term. Equilibrium equations derived from the model are solved to obtain the fixed points. Moreover, we are able to identify strong stability conditions of any fixed point by using the Jacobian matrix and eigenvalue analysis to better understand the behavior of the system. Model (1.4) has the following system with the equilibrium solutions:

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{u}=\mathcal{u}+\mathcal{h}(\mathcal{a}-\mathcal{u}-\frac{\mathcal{cuv}}{1+{\mathcal{u}}^2}),}\\ \multicolumn{1}{c}{\mathcal{v}=\mathcal{v}+\mathcal{h}(\mathcal{b}-\frac{\mathcal{cuv}}{1+{\mathcal{u}}^2}).}\end{array}$$

(2.1)

Upon simplification, a non-trivial interior fixed point is identified as:

$$({\mathcal{u}}^{*},{\mathcal{v}}^{*})=(\mathcal{a}-\mathcal{b},\frac{(1+{(\mathcal{a}-\mathcal{b})}^2)\mathcal{b}}{(\mathcal{a}-\mathcal{b})\mathcal{c}}),$$

which exists when the condition $\mathcal{a}>\mathcal{b}$ is met. For assessing the stability of the equilibria $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ , we calculate the variational matrix $\mathfrak{J}({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ for system (1.4) computed at this equilibrium. The Jacobian is expressed as follows:

$$\begin{array}{c}\multicolumn{1}{c}{\mathfrak{J}(\mathcal{a}-\mathcal{b},\frac{(1+{(\mathcal{a}-\mathcal{b})}^2)\mathcal{b}}{(\mathcal{a}-\mathcal{b})\mathcal{c}})=}\\ \multicolumn{1}{c}{\left(\begin{array}{cc}1+(\frac{2(\mathcal{a}-\mathcal{b})\mathcal{b}}{1+{(\mathcal{a}-\mathcal{b})}^2}-\frac{\mathcal{b}}{\mathcal{a}-\mathcal{b}}-1)\mathcal{h}\hfill & \frac{(\mathcal{b}-\mathcal{a})\mathcal{ch}}{1+{(\mathcal{a}-\mathcal{b})}^2}\hfill \\ \multicolumn{1}{c}{\frac{({(\mathcal{a}-\mathcal{b})}^2-1)\mathcal{bh}}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}}& 1+\frac{(\mathcal{b}-\mathcal{a})\mathcal{ch}}{1+{(\mathcal{a}-\mathcal{b})}^2}\hfill \end{array}\right).}\end{array}$$

Moreover, the characteristic polynomial related to $\mathfrak{J}(\mathcal{a}-\mathcal{b},\frac{(1+{(\mathcal{a}-\mathcal{b})}^2)\mathcal{b}}{(\mathcal{a}-\mathcal{b})\mathcal{c}})$ is expressed as:

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{P}\left(\lambda \right)={\lambda }^2-\mathcal{T}(\mathcal{a}-\mathcal{b},\frac{(1+{(\mathcal{a}-\mathcal{b})}^2)\mathcal{b}}{(\mathcal{a}-\mathcal{b})\mathcal{c}})\lambda }\\ \multicolumn{1}{c}{+\frac{\mathcal{a}+{\mathcal{a}}^3-\mathcal{b}-3{\mathcal{a}}^2\mathcal{b}+3\mathcal{a}{\mathcal{b}}^2-{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}{\mathcal{h}}^2}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}}\\ \multicolumn{1}{c}{-\frac{(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c})\mathcal{h}}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}.}\end{array}$$

(2.2)

Furthermore, from (2.2) it follows that

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{T}(\mathcal{a}-\mathcal{b},\frac{(1+{(\mathcal{a}-\mathcal{b})}^2)\mathcal{b}}{(\mathcal{a}-\mathcal{b})\mathcal{c}})=2}\\ \multicolumn{1}{c}{-\frac{(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c})\mathcal{h}}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})},}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{P}(-1)=4+}\\ \multicolumn{1}{c}{\frac{\mathcal{h}(-2(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3)-2{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}+{(\mathcal{a}-\mathcal{b})}^2\mathcal{ch})}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})},}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{P}\left(0\right)=\frac{\mathcal{a}+{\mathcal{a}}^3-\mathcal{b}-3{\mathcal{a}}^2\mathcal{b}+3\mathcal{a}{\mathcal{b}}^2-{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}{\mathcal{h}}^2}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}}\\ \multicolumn{1}{c}{-\frac{(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c})\mathcal{h}}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})},}\end{array}$$

and

$$\mathcal{P}\left(1\right)=\frac{(\mathcal{a}-\mathcal{b})\mathcal{c}{\mathcal{h}}^2}{1+{(\mathcal{a}-\mathcal{b})}^2}>0.$$

(2.3)

Furthermore, by utilizing the established connections between the coefficients and roots of the quadratic polynomials, the following Lemma 2.1 is formulated.^22,23

• Lemma 2.1 Consider the quadratic polynomial $\mathrm{\Psi }\left(\xi \right)={\xi }^2-\Gamma \xi +\mathrm{\Omega },$ satisfying $\mathrm{\Psi }\left(1\right)>0$ such that ${\xi }_1$ , ${\xi }_2$ are the roots of ${\xi }^2-\Gamma \xi +\mathrm{\Omega }=0.$ Then the following statements hold true:

•  (a)  $\left|{\xi }_1\right|<1$ & $\left|{\xi }_2\right|<1⟺\mathrm{\Psi }(-1)>0$ and $\mathrm{\Psi }\left(0\right)<1$ .

•  (b)  $\left|{\xi }_1\right|>1$ & $\left|{\xi }_2\right|>1⟺\mathrm{\Psi }(-1)>0$ and $\mathrm{\Psi }\left(0\right)>1$ .

•  (c)  $\left|{\xi }_1\right|<1$ & $\left|{\xi }_2\right|>1$ or ( $\left|{\xi }_1\right|>1$ & $\left|{\xi }_2\right|<1$ ) $⟺\mathrm{\Psi }(-1)<0$ .

•  (d)  ${\xi }_1=-1$ & $\left|{\xi }_2\right|\ne 1⟺\mathrm{\Psi }(-1)=0$ and $\mathrm{\Psi }\left(0\right)\ne \pm 1$ .

•  (e)  ${\xi }_1$ and ${\xi }_2$ are complex conjugate with $\left|{\xi }_1\right|=\left|{\xi }_2\right|=1⟺{\Gamma }^2-4\mathrm{\Omega }<0$ and $\mathrm{\Omega }=1$ .

Next, we assume that ${\xi }_1$ and ${\xi }_2$ are the roots of (2.2), then the interior equilibrium $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ of system (1.4) can be classified according to the absolute value of ${\xi }_1$ and ${\xi }_2$ as follows.

The equilibrium point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ is locally asymptotically stable (sink) if $\left|{\xi }_{1,2}\right|<1$ . $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ is called a repeller (source) if and only if $\left|{\xi }_{1,2}\right|>1$ and thus is always unstable. $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ is known as a saddle point if $\left|{\xi }_1\right|<1$ and $\left|{\xi }_2\right|>1$ or ( $\left|{\xi }_1\right|>1$ & $\left|{\xi }_2\right|<1$ ), and thus a saddle point is also unstable. $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ is known as non-hyperbolic if either $\left|{\xi }_1\right|=1$ or $\left|{\xi }_2\right|=1$ .

Consequently, from (2.3), it follows that $\mathcal{P}\left(1\right)=\frac{(\mathcal{a}-\mathcal{b})\mathcal{c}{\mathcal{h}}^2}{1+{(\mathcal{a}-\mathcal{b})}^2}>0$ . Therefore, we can apply Lemma 2.1 to prove the following stability results.

• Lemma 2.2. Assume that $\mathcal{a}>\mathcal{b}$ , then the interior equilibrium point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})=(\mathcal{a}-\mathcal{b},\frac{(1+{(\mathcal{a}-\mathcal{b})}^2)\mathcal{b}}{(\mathcal{a}-\mathcal{b})\mathcal{c}})$ of system (1.4) exists, and the following results hold.

•  i.  $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ is a stable spiral if and only if ${\mathcal{T}}^2-4\mathcal{P}\left(0\right)<0$ ,

$$\begin{array}{c}\multicolumn{1}{c}{4(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{b}-\mathcal{a})<\mathcal{h}(-2(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3)}\\ \multicolumn{1}{c}{-2{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}+{(\mathcal{a}-\mathcal{b})}^2\mathcal{ch}),}\end{array}$$

and

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\mathcal{a}+{\mathcal{a}}^3-\mathcal{b}-3{\mathcal{a}}^2\mathcal{b}+3\mathcal{a}{\mathcal{b}}^2-{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}{\mathcal{h}}^2}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}<1}\\ \multicolumn{1}{c}{+\frac{(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c})\mathcal{h}}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}.}\end{array}$$

•  ii.  $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ is a stable node if and only if ${\mathcal{T}}^2-4\mathcal{P}\left(0\right)>0$ ,

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\mathcal{a}+{\mathcal{a}}^3-\mathcal{b}-3{\mathcal{a}}^2\mathcal{b}+3\mathcal{a}{\mathcal{b}}^2-{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}{\mathcal{h}}^2}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}}\\ \multicolumn{1}{c}{<1+\frac{(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c})\mathcal{h}}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})},}\end{array}$$

and

$$\begin{array}{c}\multicolumn{1}{c}{4(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{b}-\mathcal{a})<\mathcal{h}(-2(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3)}\\ \multicolumn{1}{c}{-2{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}+{(\mathcal{a}-\mathcal{b})}^2\mathcal{ch}}\end{array}).$$

•  iii. $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ is an unstable spiral if and only if ${\mathcal{T}}^2-4\mathcal{P}\left(0\right)<0$ ,

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{h}(-2(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3)-2{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}+{(\mathcal{a}-\mathcal{b})}^2\mathcal{ch})>4}\\ \multicolumn{1}{c}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{b}-\mathcal{a}),}\end{array}$$

and

$$\begin{array}{c}\multicolumn{1}{c}{1+\frac{(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c})\mathcal{h}}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}}\\ \multicolumn{1}{c}{<\frac{\mathcal{a}+{\mathcal{a}}^3-\mathcal{b}-3{\mathcal{a}}^2\mathcal{b}+3\mathcal{a}{\mathcal{b}}^2-{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}{\mathcal{h}}^2}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}.}\end{array}$$

•  iv.  $({\mathcal{u}}^{*},{\mathcal{v}}^{*})\mathrm{is}\mathrm{an}\mathrm{unstable}\mathrm{node}\mathrm{if}\mathrm{and}\mathrm{only}\mathrm{if}{\mathcal{T}}^2-4\mathcal{P}\left(0\right)>0,$

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{h}(-2(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3)}\\ \multicolumn{1}{c}{-2{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}+{(\mathcal{a}-\mathcal{b})}^2\mathcal{ch})>4(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{b}-\mathcal{a}),}\end{array}$$

and

$$\begin{array}{c}\multicolumn{1}{c}{1+\frac{(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c})\mathcal{h}}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}}\\ \multicolumn{1}{c}{<\frac{\mathcal{a}+{\mathcal{a}}^3-\mathcal{b}-3{\mathcal{a}}^2\mathcal{b}+3\mathcal{a}{\mathcal{b}}^2-{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}{\mathcal{h}}^2}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}.}\end{array}$$

•  v. $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ is a saddle point if and only if

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{h}(-2(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3)-2{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}+{(\mathcal{a}-\mathcal{b})}^2\mathcal{ch})}\\ \multicolumn{1}{c}{<4(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{b}-\mathcal{a}).}\end{array}$$

•  vi. Assume that ${\xi }_1$ and ${\xi }_2$ are the roots of (2.2), then ${\xi }_1=-1$ and $\left|{\xi }_2\right|\ne 1⟺$

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{c}=\frac{1}{{(\mathcal{a}-\mathcal{b})}^2\mathcal{h}}(4(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}\\ \multicolumn{1}{c}{+(\mathcal{a}+\mathcal{a}{(\mathcal{a}-\mathcal{b})}^2-2\mathcal{b})\mathcal{h}-(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3){\mathcal{h}}^2,}\end{array}$$

and

$$\begin{array}{c}\multicolumn{1}{c}{\frac{\mathcal{a}+{\mathcal{a}}^3-\mathcal{b}-3{\mathcal{a}}^2\mathcal{b}+3\mathcal{a}{\mathcal{b}}^2-{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}{\mathcal{h}}^2}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}}\\ \multicolumn{1}{c}{-\frac{(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3+{(\mathcal{a}-\mathcal{b})}^2\mathcal{c})\mathcal{h}}{(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}\ne \pm 1.}\end{array}$$

• vii. The roots of (2.2) are complex conjugate with $\left|{\xi }_1\right|=\left|{\xi }_2\right|=1$ if and only if

$$\mathcal{c}=\frac{\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3}{{(\mathcal{a}-\mathcal{b})}^2(-1+\mathcal{h})},$$

$$\begin{array}{l}({(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3)}^2-2(\mathcal{a}+\mathcal{a}{(\mathcal{a}-\mathcal{b})}^2-2\mathcal{b})\\ \frac{{(\mathcal{a}-\mathcal{b})}^2\mathcal{c}+{(\mathcal{a}-\mathcal{b})}^4{\mathcal{c}}^2){\mathcal{h}}^2}{{(1+{(\mathcal{a}-\mathcal{b})}^2)}^2{(\mathcal{a}-\mathcal{b})}^2<0.}\end{array}$$

The topological classification of the positive equilibrium point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ is depicted in Figure 1.

Figure 1.

Topology-based characterization of the non-negative equilibrium point $(\mathcal{a}-\mathcal{b},\frac{(1+{(\mathcal{a}-\mathcal{b})}^2)\mathcal{b}}{(\mathcal{a}-\mathcal{b})\mathcal{c}})$ for $\mathcal{a}=26,\mathcal{b}=24,$ with $0.001\le \mathcal{h}\le 1$ and $0.001\le \mathcal{c}\le 50$ .

Bifurcation analysis

Here, we examine the bifurcation of system (1.4) by varying the important parameters of the system. This model, when undiffused, exhibits numerous dynamic transitions which give an indication of qualitative changes in system behavior. These are a period-doubling bifurcation where a steady equilibrium point turns into an unsteady orbit and a period-two orbit is generated, and a Neimark-Sacker bifurcation, whereby closed growing curves are generated. Moreover, system (1.4) has higher-order resonances, that is, 1:2, 1:3, and 1:4 bifurcations, that are important in the dynamics of the system around particular critical parameter values. Our analysis begins with the study of the development and evolution of the period-doubling bifurcation.

Period-doubling bifurcation at $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$

Model (1.4) undergoes a period-doubling bifurcation when the Jacobian matrix evaluated at the fixed point (u*, v*) has an eigenvalue ${\lambda }_1=-1$ , in contrast, the second eigenvalue ${\lambda }_2\ne \pm 1$ .

We define the parameter set associated with this bifurcation as

$${\mathrm{\Psi }}_{\mathrm{PDB}}=\left\{\begin{array}{c}\multicolumn{1}{c}{(\mathcal{a},\mathcal{b},\mathcal{c},\mathcal{h})\in {\mathbb{R}}_{+}^4:\mathcal{c}={\mathcal{c}}_0}\end{array}\right\},$$

in the case that the critical parameters’ value ${\mathcal{c}}_0$ is given as

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{c}}_0=\frac{1}{{(\mathcal{a}-\mathcal{b})}^2\mathcal{h}}(4(1+{(\mathcal{a}-\mathcal{b})}^2)(\mathcal{a}-\mathcal{b})}\\ \multicolumn{1}{c}{+(\mathcal{a}+\mathcal{a}{(\mathcal{a}-\mathcal{b})}^2-2\mathcal{b})\mathcal{h}}\\ \multicolumn{1}{c}{-(\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3){\mathcal{h}}^2.}\end{array}$$

We now consider a small perturbation $\tilde{\mathcal{c}}$ of $\mathcal{c}$ , where $\left|\tilde{\mathcal{c}}\right|<<1,$ and study the following perturbed system:

$$\left[\begin{array}{c}\mathcal{u}\hfill \\ \multicolumn{1}{c}{\mathcal{v}}\end{array}\right]\to \left[\begin{array}{c}\mathcal{u}+\mathcal{h}(\mathcal{a}-\mathcal{u}-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_0)\mathcal{uv}}{1+{\mathcal{u}}^2}),\hfill \\ \multicolumn{1}{c}{\mathcal{v}+\mathcal{h}(\mathcal{b}-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_0)\mathcal{uv}}{1+{\mathcal{u}}^2}).}\end{array}\right].$$

(3.1)

To analyze system (3.1) near its fixed point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ , we shift toward the origin by defining $\mathcal{p}=\mathcal{u}-{\mathcal{u}}^{*}\mathrm{and}\mathcal{q}=\mathcal{v}-{\mathcal{v}}^{*}$ . Expanding the Taylor series around the origin yields the following expression for the map:

$$\left[\begin{array}{c}\mathcal{p}\hfill \\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right]\to \left[\begin{array}{cc}{\varrho }_1\hfill & {\varrho }_2\hfill \\ \multicolumn{1}{c}{{\sigma }_1}& {\sigma }_2\hfill \end{array}\right]\left[\begin{array}{c}\mathcal{p}\hfill \\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right]+\left[\begin{array}{c}{\mathcal{f}}_0(\mathcal{p},\mathcal{q},\tilde{\mathcal{c}})\hfill \\ \multicolumn{1}{c}{{\mathcal{g}}_0(\mathcal{p},\mathcal{q},\tilde{\mathcal{c}})}\end{array}\right],$$

(3.2)

where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{f}}_0(\mathcal{p},\mathcal{q},\tilde{\mathcal{c}})={\varrho }_{11}{\mathcal{p}}^2+{\varrho }_{12}\mathcal{pq}+{\varrho }_{13}\mathcal{p}\tilde{\mathcal{c}}+{\varrho }_{22}{\mathcal{q}}^2+{\varrho }_{23}\mathcal{q}\tilde{\mathcal{c}}+{\varrho }_{33}{\tilde{\mathcal{c}}}^2}\\ \multicolumn{1}{c}{+{\varrho }_{111}{\mathcal{p}}^3+{\varrho }_{112}{\mathcal{p}}^2\mathcal{q}+{\varrho }_{113}{\mathcal{p}}^2\tilde{\mathcal{c}}+{\varrho }_{122}\mathcal{p}{\mathcal{q}}^2+{\varrho }_{123}\mathcal{pq}\tilde{\mathcal{c}}}\\ \multicolumn{1}{c}{+{\varrho }_{222}{\mathcal{q}}^3+{\varrho }_{223}{\mathcal{q}}^2\tilde{\mathcal{c}}+{\varrho }_{333}{\tilde{\mathcal{c}}}^3+\mathcal{O}\left({(\left|\mathcal{p}\right|+\left|\mathcal{q}\right|+\left|\tilde{\mathcal{c}}\right|)}^4\right),}\\ \multicolumn{1}{c}{{\mathcal{g}}_0(\mathcal{p},\mathcal{q},\tilde{\mathcal{c}})={\sigma }_{11}{\mathcal{p}}^2+{\sigma }_{12}\mathcal{pq}+{\sigma }_{13}\mathcal{p}\tilde{\mathcal{c}}+{\sigma }_{22}{\mathcal{q}}^2+{\sigma }_{23}\mathcal{q}\tilde{\mathcal{c}}+{\sigma }_{33}{\tilde{\mathcal{c}}}^2}\\ \multicolumn{1}{c}{+{\sigma }_{111}{\mathcal{p}}^3+{\sigma }_{112}{\mathcal{p}}^2\mathcal{q}+{\sigma }_{113}{\mathcal{p}}^2\tilde{\mathcal{c}}+{\sigma }_{122}\mathcal{p}{\mathcal{q}}^2+{\sigma }_{123}\mathcal{pq}\tilde{\mathcal{c}}}\\ \multicolumn{1}{c}{+{\sigma }_{222}{\mathcal{q}}^3+{\sigma }_{223}{\mathcal{q}}^2\tilde{\mathcal{c}}+{\sigma }_{333}{\tilde{\mathcal{c}}}^3+\mathcal{O}\left({(\left|\mathcal{p}\right|+\left|\mathcal{q}\right|+\left|\tilde{\mathcal{c}}\right|)}^4\right),}\end{array}$$

(3.3)

and the coefficients are

$$\begin{array}{c}\multicolumn{1}{c}{{\varrho }_1=1-\mathcal{h}+\frac{2(-2+\mathcal{h})(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{(-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}},}\\ \multicolumn{1}{c}{{\varrho }_{13}=\frac{\mathcal{h}(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{{(1+{\mathcal{u}}^{*2})}^2},}\\ \multicolumn{1}{c}{{\varrho }_2=-\frac{2(-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})}{(-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}},}\\ \multicolumn{1}{c}{{\varrho }_{22}=0,{\varrho }_{23}=-\frac{\mathcal{h}{\mathcal{u}}^{*}}{1+{\mathcal{u}}^{*2}},}\\ \multicolumn{1}{c}{{\varrho }_{11}=-\frac{2(-2+\mathcal{h}){\mathcal{u}}^{*}(-3+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{(-2+\mathcal{h}){\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^2+2(-1+{\mathcal{u}}^{*4}){\mathcal{v}}^{*}},}\\ \multicolumn{1}{c}{{\varrho }_{113}=-\frac{\mathcal{h}{\mathcal{u}}^{*}(-3+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{{(1+{\mathcal{u}}^{*2})}^3},}\\ \multicolumn{1}{c}{{\varrho }_{12}=\frac{2(-2+\mathcal{h})(-1+{\mathcal{u}}^{*2})}{(-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}},}\\ \multicolumn{1}{c}{{\varrho }_{123}=\frac{\mathcal{h}(-1+{\mathcal{u}}^{*2})}{(1+{\mathcal{u}}^{*2}{)}^2},}\\ \multicolumn{1}{c}{{\varrho }_{111}=\frac{2(-2+\mathcal{h})(1-6{\mathcal{u}}^{*2}+{\mathcal{u}}^{*4}){\mathcal{v}}^{*}}{(1+{\mathcal{u}}^{*2}{)}^2((-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*})},}\\ \multicolumn{1}{c}{{\varrho }_{122}=0,{\varrho }_{33}=0,}\\ \multicolumn{1}{c}{{\varrho }_{112}=-\frac{2(-2+\mathcal{h}){\mathcal{u}}^{*}(-3+{\mathcal{u}}^{*2})}{(-2+\mathcal{h}){\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^2+2(-1+{\mathcal{u}}^{*4}){\mathcal{v}}^{*}},}\\ \multicolumn{1}{c}{{\varrho }_{222}=0,{\varrho }_{223}=0,{\varrho }_{333}=0,}\\ \multicolumn{1}{c}{{\sigma }_1=\frac{2(-2+\mathcal{h})(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{(-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}},}\\ \multicolumn{1}{c}{{\sigma }_{13}=\frac{\mathcal{h}(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{{(1+{\mathcal{u}}^{*2})}^2},}\\ \multicolumn{1}{c}{{\sigma }_2=1-\frac{2(-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})}{(-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}},}\\ \multicolumn{1}{c}{{\sigma }_{22}=0,{\sigma }_{23}=(-\frac{\mathcal{h}{\mathcal{u}}^{*}}{1+{\mathcal{u}}^{*2}}),}\\ \multicolumn{1}{c}{{\sigma }_{11}=-\frac{2(-2+\mathcal{h}){\mathcal{u}}^{*}(-3+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{(-2+\mathcal{h}){\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^2+2(-1+{\mathcal{u}}^{*4}){\mathcal{v}}^{*}},}\\ \multicolumn{1}{c}{{\sigma }_{113}=-\frac{\mathcal{h}{\mathcal{u}}^{*}(-3+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{{(1+{\mathcal{u}}^{*2})}^3},}\\ \multicolumn{1}{c}{{\sigma }_{12}=\frac{2(-2+\mathcal{h})(-1+{\mathcal{u}}^{*2})}{(-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}},}\\ \multicolumn{1}{c}{{\sigma }_{33}=0,{\sigma }_{123}=\frac{\mathcal{h}(-1+{\mathcal{u}}^{*2})}{(1+{\mathcal{u}}^{*2}{)}^2},}\\ \multicolumn{1}{c}{{\sigma }_{111}=\frac{2(-2+\mathcal{h})(1-6{\mathcal{u}}^{*2}+{\mathcal{u}}^{*4}){\mathcal{v}}^{*}}{(1+{\mathcal{u}}^{*2}{)}^2((-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*})},}\\ \multicolumn{1}{c}{{\sigma }_{122}=0,}\\ \multicolumn{1}{c}{{\sigma }_{112}=-\frac{2(-2+\mathcal{h}){\mathcal{u}}^{*}(-3+{\mathcal{u}}^{*2})}{(-2+\mathcal{h}){\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^2+2(-1+{\mathcal{u}}^{*4}){\mathcal{v}}^{*}},}\\ \multicolumn{1}{c}{{\sigma }_{222}=0,{\sigma }_{223}=0,{\sigma }_{333}=0.}\end{array}$$

To transform the coefficient matrix $\left[\begin{array}{cc}{\varrho }_1\hfill & {\varrho }_2\hfill \\ \multicolumn{1}{c}{{\sigma }_1}& {\sigma }_2\hfill \end{array}\right]$ from the map (3.2) into its normal form, we apply the following coordinate transformation:

$$\left[\begin{array}{c}\multicolumn{1}{c}{\mathcal{p}}\\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right]=\mathrm{T}\left[\begin{array}{c}\multicolumn{1}{c}{\mathcal{s}}\\ \multicolumn{1}{c}{\mathcal{t}}\end{array}\right],$$

(3.4)

where $\mathrm{T}$ is an invertible transformation matrix formulated as

$$\mathrm{T}=\left[\begin{array}{cc}{\varrho }_2\hfill & {\varrho }_2\hfill \\ \multicolumn{1}{c}{{\lambda }_1-{\varrho }_1}& {\lambda }_2-{\varrho }_1\hfill \end{array}\right],$$

and ${\lambda }_1$ and ${\lambda }_2$ are the eigenvalues of the matrix $\left[\begin{array}{cc}{\varrho }_1\hfill & {\varrho }_2\hfill \\ \multicolumn{1}{c}{{\sigma }_1}& {\sigma }_2\hfill \end{array}\right]$ . Applying this transformation, the map (3.2) is rewritten in the new coordinates $(\mathcal{s},\mathcal{t})$ as

$$\begin{array}{c}\multicolumn{1}{c}{\left[\begin{array}{c}\mathcal{s}\hfill \\ \multicolumn{1}{c}{\mathcal{t}}\end{array}\right]\to \left[\begin{array}{cc}-1\hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \frac{-((-2+\mathcal{h})(-1+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}))+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{(-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}\hfill \end{array}\right]}\\ \multicolumn{1}{c}{\left[\begin{array}{c}\mathcal{s}\hfill \\ \multicolumn{1}{c}{\mathcal{t}}\end{array}\right]+\left[\begin{array}{c}{\mathcal{f}}_1(\mathcal{s},\mathcal{t},\tilde{\mathcal{c}})\hfill \\ \multicolumn{1}{c}{{\mathcal{g}}_1(\mathcal{s},\mathcal{t},\tilde{\mathcal{c}})}\end{array}\right],}\end{array}$$

(3.5)

where

$$\{\begin{array}{c}{\mathcal{f}}_1(\mathcal{s},\mathcal{t},\tilde{\mathcal{c}})={\Gamma }_1{\mathcal{f}}_0(\mathcal{p},\mathcal{q},\tilde{\mathcal{c}})+{\Gamma }_2{\mathcal{g}}_0(\mathcal{p},\mathcal{q},\tilde{\mathcal{c}}),\hfill \\ \multicolumn{1}{c}{}\\ \multicolumn{1}{c}{{\mathcal{g}}_1(\mathcal{s},\mathcal{t},\tilde{\mathcal{c}})={\Gamma }_3{\mathcal{f}}_0(\mathcal{p},\mathcal{q},\tilde{\mathcal{c}})+{\Gamma }_4{\mathcal{g}}_0(\mathcal{p},\mathcal{q},\tilde{\mathcal{c}}),}\end{array}$$

(3.6)

and the coefficients are

$$\begin{array}{c}\multicolumn{1}{c}{{\Gamma }_1=\frac{2{\mathcal{v}}^{*}((-2+\mathcal{h}){\mathcal{u}}^{*}(-1+{\mathcal{u}}^{*4})+2(-1+{\mathcal{u}}^{*2}{)}^2{\mathcal{v}}^{*})}{(-2+\mathcal{h}){\mathcal{u}}^{*}((-2+{\mathcal{h})}^2{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^2-4(-1+{\mathcal{u}}^{*4}){\mathcal{v}}^{*})},}\\ \multicolumn{1}{c}{{\Gamma }_2=\frac{(-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{(-2+{\mathcal{h})}^2{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})-4(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}},}\\ \multicolumn{1}{c}{{\Gamma }_3=\frac{(-2+\mathcal{h})(-((-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}))-2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*})}{2(-2+{\mathcal{h})}^2{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})-8(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}},}\\ \multicolumn{1}{c}{{\Gamma }_4=-\frac{(-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{{(-2+\mathcal{h})}^2{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})-4(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}.}\end{array}$$

Let ${\mathcal{W}}^c(0,0,0)$ denote the center manifold of the transformed structure (3.5), computed at the origin $(0,0)$ and within a slight vicinity around $\tilde{\mathcal{c}}=0$ . The given topological space can be locally approximated by

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{W}}^c(0,0,0)=\{(\mathcal{s},\mathcal{t},\tilde{\mathcal{c}})\in {\mathbb{R}}^3:\mathcal{t}}\\ \multicolumn{1}{c}{={\mathcal{m}}_1{\mathcal{s}}^2+{\mathcal{m}}_2\mathcal{s}\tilde{\mathcal{c}}+{\mathcal{m}}_3{\tilde{\mathcal{c}}}^2}\end{array}+\mathcal{O}\left({(\left|\mathcal{s}\right|+\left|\tilde{\mathcal{c}}\right|)}^3\right)\}.$$

(3.7)

Let

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{t}=\mathcal{h}(\mathcal{s},\tilde{\mathcal{c}})={\mathcal{m}}_1{\mathcal{s}}^2+{\mathcal{m}}_2\mathcal{s}\tilde{\mathcal{c}}+{\mathcal{m}}_3{\tilde{\mathcal{c}}}^2.}\end{array}$$

(3.8)

Then, the function $\mathcal{h}(\mathcal{s},\tilde{\mathcal{c}})$ satisfies the following equation

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{P}\left[\mathcal{h}(\mathcal{s},\tilde{\mathcal{c}})\right]=\mathcal{h}(-\mathcal{s}+{\mathcal{f}}_1(\mathcal{s},\mathcal{h}(\mathcal{s},\tilde{\mathcal{c}}),\tilde{\mathcal{c}}),\tilde{\mathcal{c}})}\\ \multicolumn{1}{c}{-(1+\frac{1}{2}\mathcal{a}(-\mathcal{b}+\mathcal{d})\mathcal{c}{\mathcal{p}}^{*}{\mathcal{q}}^{*})\mathcal{h}(\mathcal{s},\tilde{\mathcal{c}})}\\ \multicolumn{1}{c}{-{\mathcal{g}}_1(\mathcal{s},\mathcal{h}(\mathcal{s},\tilde{\mathcal{c}}),\tilde{\mathcal{c}})=0.}\end{array}$$

(3.9)

By comparing the coefficients, we derive the following values for ${\mathcal{m}}_1$ , ${\mathcal{m}}_2$ and ${\mathcal{m}}_3$ ,

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{m}}_1=\frac{2(-2+{\mathcal{h})}^2{\mathcal{u}}^{*}((-2+\mathcal{h})(-1+{\mathcal{u}}^{*4})+2{\mathcal{u}}^{*}(-3+{\mathcal{u}}^{*2}){\mathcal{v}}^{*})}{(-2+{\mathcal{h})}^3({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}{)}^2+2(-4+\mathcal{h})(-2+\mathcal{h}){\mathcal{u}}^{*}(-1+{\mathcal{u}}^{*4}){\mathcal{v}}^{*}-8(-1+{\mathcal{u}}^{*2}{)}^2{\mathcal{v}}^{*2}},}\\ \multicolumn{1}{c}{{\mathcal{m}}_2=\frac{(-2+\mathcal{h}){\mathcal{u}}^{*}((-2+\mathcal{h})\mathcal{h}{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2\mathcal{h}(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}{)}^2}{2(1+{\mathcal{u}}^{*2})((-2+{\mathcal{h})}^2{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})-4(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}{)}^2},}\\ \multicolumn{1}{c}{{\mathcal{m}}_3=0.}\end{array}$$

As a result, the map restricted to the center manifold ${\mathcal{W}}^c(0,0,0),$ reduces to a one-dimensional form involving only the variable $\mathcal{s}$ , given by

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{F}:\mathcal{s}\to -\mathcal{s}+{\tau }_{11}{\mathcal{s}}^2+{\tau }_{12}\mathcal{s}\tilde{\mathcal{c}}+{\tau }_{112}{\mathcal{s}}^2\tilde{\mathcal{c}}+{\tau }_{122}\mathcal{s}{\tilde{\mathcal{c}}}^2}\\ \multicolumn{1}{c}{+{\tau }_{111}{\mathcal{s}}^3+\mathcal{O}\left({(\left|\mathcal{s}\right|+\left|\tilde{\mathcal{c}}\right|)}^4\right),}\end{array}$$

(3.10)

where the coefficients are

$$\{\begin{array}{l}{\tau }_{11}=\frac{4{(-2+\mathcal{h})}^2{\mathcal{u}}^{\ast }(-((-2+\mathcal{h})(-1+{\mathcal{u}}^{\ast 4}))-2{\mathcal{u}}^{\ast }(-3+{\mathcal{u}}^{\ast 2}){\mathcal{v}}^{\ast })}{{(-2+\mathcal{h})}^3{({\mathcal{u}}^{\ast }+{\mathcal{u}}^{\ast 3})}^2+2(-4+\mathcal{h})(-2+\mathcal{h}){\mathcal{u}}^{\ast }(-1+{\mathcal{u}}^{\ast 4}){\mathcal{v}}^{\ast }-8{(-1+{\mathcal{u}}^{\ast 2})}^2{\mathcal{v}}^{\ast 2}},\\ {\tau }_{12}=-\frac{\mathcal{h}{\left(\left(-2+\mathcal{h}\right){\mathcal{u}}^{\ast }\left(1+{\mathcal{u}}^{\ast 2}\right)+2\left(-1+{\mathcal{u}}^{\ast 2}\right){\mathcal{u}}^{\ast }\right)}^2}{{\left(1+{\mathcal{u}}^{\ast 2}\right)}^2\left({\left(-2+\mathcal{u}\right)}^2{\mathcal{u}}^{\ast }\left(1+{\mathcal{u}}^{\ast 2}\right)-4\left(-1+{\mathcal{u}}^{\ast 2}\right){\mathcal{v}}^{\ast }\right)},\\ {\tau }_{112}=\frac{1}{({\left(1+{\mathcal{u}}^{\ast 2}\right)}^2{\left({\left(-2+\mathcal{h}\right)}^2{\mathcal{u}}^{\ast }\left(1+{\mathcal{u}}^{\ast 2}\right)-4\left(-1+{\mathcal{u}}^{\ast 2}\right){\mathcal{v}}^{\ast }\right)}^3))}(-((4{\left(-2+\mathcal{h}\right)}^2\mathcal{h}{\mathcal{u}}^{\ast }\\ ({(-2+\mathcal{h})}^3(-1+\mathcal{h}){\mathcal{u}}^{\ast 2}(-1+{\mathcal{u}}^{\ast 2}){(1+{\mathcal{u}}^{\ast 2})}^3+(-2+\mathcal{h}){\mathcal{u}}^{\ast }{(1+{\mathcal{u}}^{\ast 2})}^2(8+2(-3+\mathcal{h})\\ \mathcal{h}+(-28+(36-13\mathcal{h})\mathcal{h}){\mathcal{u}}^{\ast 2}+(12+\mathcal{h}(-14+5\mathcal{h})){\mathcal{u}}^{\ast 4}){\mathcal{v}}^{\ast }+2(-1+{\mathcal{u}}^{\ast 4})\\ (4+{\mathcal{u}}^{\ast 2}(\mathcal{h}(-10+3\mathcal{h})(-3+{\mathcal{u}}^{\ast 2})+4(-8+3{\mathcal{u}}^{\ast 2}))){\mathcal{v}}^{\ast 2}-8{\mathcal{u}}^{\ast }(-3+{\mathcal{u}}^{\ast 2}){\left(-1+{\mathcal{u}}^{\ast 2}\right)}^2{\mathcal{v}}^{\ast 3})),\\ {\tau }_{122}=-\frac{{\mathcal{h}}^4{\mathcal{u}}^{\ast }\left(-1+{\mathcal{u}}^{\ast 2}\right){\mathcal{v}}^{\ast }{\left(\left(-2+\mathcal{h}\right){\mathcal{u}}^{\ast }\left(1+{\mathcal{u}}^{\ast 2}\right)+2\left(-1+{\mathcal{u}}^{\ast 2}\right){\mathcal{v}}^{\ast }\right)}^3}{{\left({\left(-2+\mathcal{h}\right)}^2{\mathcal{u}}^{\ast }{\left(1+{\mathcal{u}}^{\ast 2}\right)}^2-4\left(-1+{\mathcal{u}}^{\ast 4}\right){\mathcal{v}}^{\ast }\right)}^3},\\ {\tau }_{111}=\frac{-1}{{({(-2+\mathcal{h})}^3{({\mathcal{u}}^{\ast }+{\mathcal{u}}^{\ast 3})}^2+2(-4+\mathcal{h})(-2+\mathcal{h}){\mathcal{u}}^{\ast }(-1+{\mathcal{u}}^{\ast 4}){\mathcal{v}}^{\ast }-8{(-1+{\mathcal{u}}^{\ast 2})}^2{\mathcal{v}}^{\ast 2})}^2)}\\ ((8{(-2+\mathcal{h})}^3{\mathcal{u}}^{\ast }({(-2+\mathcal{h})}^3{\mathcal{u}}^{\ast }{(1+{\mathcal{u}}^{\ast 2})}^2(1-5{\mathcal{u}}^{\ast 2}+2{\mathcal{u}}^{\ast 4})+4(-2+\mathcal{h})(1+{\mathcal{u}}^{\ast 2})\\ (-1+5(-2+\mathcal{h}){\mathcal{u}}^{\ast 2}+(19-9\mathcal{h}){\mathcal{u}}^{\ast 4}+2(-2+\mathcal{h}){\mathcal{u}}^{\ast 6}){\mathcal{v}}^{\ast }+8{\mathcal{u}}^{\ast }(-2+9(-2+\mathcal{h})\\ {\mathcal{u}}^{\ast 2}+(14-6\mathcal{h}){\mathcal{u}}^{\ast 4}+(-2+\mathcal{h}){\mathcal{u}}^{\ast 6}){\mathcal{v}}^{\ast 2})).\end{array}$$

To rigorously establish the emergence of a 2-cycle bifurcation, we evaluate the following conditions,

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{P}}_1=\frac{\partial \mathcal{F}}{\partial \mathcal{s}\partial \tilde{\mathcal{c}}}+\frac{1}{2}\frac{\partial \mathcal{F}}{\partial \tilde{\mathcal{c}}}\frac{{\partial }^2\mathcal{F}}{\partial {\mathcal{s}}^2}}\\ \multicolumn{1}{c}{=-\frac{\mathcal{h}{((-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*})}^2}{{(1+{\mathcal{u}}^{*2})}^2({(-2+\mathcal{h})}^2{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})-4(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*})}\ne 0,}\end{array}$$

and

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{P}}_2=\frac{1}{6}\frac{{\partial }^3\mathcal{F}}{\partial {\mathcal{s}}^3}+\frac{1}{2}\frac{{\partial }^2\mathcal{F}}{\partial {\mathcal{s}}^2}=\frac{8(-2+{\mathcal{h})}^3{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})((-2+\mathcal{h})(1+{\mathcal{u}}^{*2})+4{\mathcal{u}}^{*}{\mathcal{v}}^{*})}{((-2+{\mathcal{h})}^2{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})-4(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*})((-2+\mathcal{h}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}{)}^2}\ne 0}\end{array}.$$

These conditions ensure the non-degeneracy and transversality required for a period-doubling bifurcation. We consolidate our findings in the succeeding theorem.

• Theorem 3.1. System (1.4) exhibits a period-doubling bifurcation at $(\mathcal{a}-\mathcal{b},\frac{(1+{(\mathcal{a}-\mathcal{b})}^2)\mathcal{b}}{(\mathcal{a}-\mathcal{b})\mathcal{c}}),$ as $\mathcal{c}$ varies near the critical value ${\mathcal{c}}_0\in {\psi }_{\mathrm{PDB}}$ , provided that the bifurcation conditions ${\mathcal{P}}_1\ne 0\mathrm{and}{\mathcal{P}}_2\ne 0$ are satisfied. Furthermore, if ${\mathcal{P}}_2>0$ , a supercritical period-doubling bifurcation happens at $\mathcal{c}={\mathcal{c}}_0,$ leading to the emergence of an attracting 2-cycle for $\mathcal{c}>{\mathcal{c}}_0$ . In contrast, when ${\mathcal{P}}_2<0,$ a subcritical bifurcation occurs, and a repelling 2-cycle exists for $\mathcal{c}<{\mathcal{c}}_0$ .

Neimark-Sacker bifurcation at $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$

The Neimark-Sacker bifurcation (NSB) is common in discrete dynamical systems. It happens when an equilibrium point becomes unstable: two complex-conjugate eigenvalues of the Jacobian matrix of the equilibrium point cross the unit circle in the complex plane. It is the transition that leads to the appearance of closed invariant curves, which is the case of quasiperiodic dynamics.

The bifurcation parameter associated with NSB is identified as the critical value

$${\mathcal{c}}_1=\frac{\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3}{{(\mathcal{a}-\mathcal{b})}^2(-1+\mathcal{h})}.$$

Accordingly, we define the critical parameter set associated with the NSB as

$${\mathrm{\Psi }}_{\mathrm{NSB}}=\{(\mathcal{a},\mathcal{b},\mathcal{c},\mathcal{h})\in {\mathbb{R}}_{+}^4:\mathcal{c}={\mathcal{c}}_1\}.$$

To investigate the local dynamics near this bifurcation point, we introduce a small perturbation $\tilde{\mathcal{c}}$ of the bifurcation parameter $\mathcal{c}$ , such that $\left|\tilde{\mathcal{c}}\right|<<1,$ and study the following perturbed version of the model

$$\left[\begin{array}{c}\mathcal{u}\hfill \\ \multicolumn{1}{c}{\mathcal{v}}\end{array}\right]\to \left[\begin{array}{c}\mathcal{u}+\mathcal{h}(\mathcal{a}-\mathcal{u}-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{uv}}{1+{\mathcal{u}}^2}),\hfill \\ \multicolumn{1}{c}{\mathcal{v}+\mathcal{h}(\mathcal{b}-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{uv}}{1+{\mathcal{u}}^2}).}\end{array}\right].$$

(3.11)

For the case of shifting the fixed point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ of map (3.11) to the origin, we define the new parameters $\mathcal{p}=\mathcal{u}-{\mathcal{u}}^{*}\mathrm{and}\mathcal{q}=\mathcal{v}-{\mathcal{v}}^{*}$ . Expanding the map (3.11) using a Taylor series around the origin yields

$$\left[\begin{array}{c}\mathcal{p}\hfill \\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right]\to \left[\begin{array}{cc}{\varrho }_1\hfill & {\varrho }_2\hfill \\ \multicolumn{1}{c}{{\sigma }_1}& {\sigma }_2\hfill \end{array}\right]\left[\begin{array}{c}\mathcal{p}\hfill \\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right]+\left[\begin{array}{c}{\mathcal{f}}_0(\mathcal{p},\mathcal{q})\hfill \\ \multicolumn{1}{c}{{\mathcal{g}}_0(\mathcal{p},\mathcal{q})}\end{array}\right],$$

(3.12)

where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{f}}_0(\mathcal{p},\mathcal{q})={\varrho }_{11}{\mathcal{p}}^2+{\varrho }_{12}\mathcal{pq}+{\varrho }_{22}{\mathcal{q}}^2++{\varrho }_{111}{\mathcal{p}}^3+{\varrho }_{112}{\mathcal{p}}^2\mathcal{q}+{\varrho }_{122}\mathcal{p}{\mathcal{q}}^2+{\varrho }_{222}{\mathcal{q}}^3+\mathcal{O}\left({(\left|\mathcal{p}\right|+\left|\mathcal{q}\right|)}^4\right),}\\ \multicolumn{1}{c}{{\mathcal{g}}_0(\mathcal{p},\mathcal{q})={\sigma }_{11}{\mathcal{p}}^2+{\sigma }_{12}\mathcal{pq}+{\sigma }_{22}{\mathcal{q}}^2+{\sigma }_{111}{\mathcal{p}}^3+{\sigma }_{112}{\mathcal{p}}^2\mathcal{q}+{\sigma }_{122}\mathcal{p}{\mathcal{q}}^2+{\sigma }_{222}{\mathcal{q}}^3+\mathcal{O}\left({(\left|\mathcal{p}\right|+\left|\mathcal{q}\right|)}^4\right),}\end{array}$$

(3.13)

and the coefficients are

$$\{\begin{array}{c}\multicolumn{1}{c}{{\varrho }_1=1+\mathcal{h}(-1+\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{{(1+{\mathcal{u}}^{*2})}^2}),{\varrho }_2=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}{\mathcal{u}}^{*}}{1+{\mathcal{u}}^{*2}},{\varrho }_{12}=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}(-1+{\mathcal{u}}^{*2})}{{(1+{\mathcal{u}}^{*2})}^2},}\\ \multicolumn{1}{c}{{\varrho }_{11}=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}{\mathcal{u}}^{*}(-3+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{(1+{\mathcal{u}}^{*2}{)}^3},{\varrho }_{22}=0,{\varrho }_{111}=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}(1-6{\mathcal{u}}^{*2}+{\mathcal{u}}^{*4}){\mathcal{v}}^{*}}{(1+{\mathcal{u}}^{*2}{)}^4},{\varrho }_{122}=0,}\\ \multicolumn{1}{c}{{\varrho }_{222}=0,{\varrho }_{112}=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}{\mathcal{u}}^{*}(-3+{\mathcal{u}}^{*2})}{{(1+{\mathcal{u}}^{*2})}^3},{\sigma }_1=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{{(1+{\mathcal{u}}^{*2})}^2},{\sigma }_{22}=0,}\\ \multicolumn{1}{c}{{\sigma }_2=1-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}{\mathcal{u}}^{*}}{1+{\mathcal{u}}^{*2}}{\sigma }_{11}=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}{\mathcal{u}}^{*}(-3+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{{(1+{\mathcal{u}}^{*2})}^3},{\sigma }_{12}=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}(-1+{\mathcal{u}}^{*2})}{{(1+{\mathcal{u}}^{*2})}^2},}\\ \multicolumn{1}{c}{{\sigma }_{111}=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}(1-6{\mathcal{u}}^{*2}+{\mathcal{u}}^{*4}){\mathcal{v}}^{*}}{(1+{\mathcal{u}}^{*2}{)}^4},{\sigma }_{112}=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}{\mathcal{u}}^{*}(-3+{\mathcal{u}}^{*2})}{(1+{\mathcal{u}}^{*2}{)}^3},{\sigma }_{122}=0,{\sigma }_{222}=0.}\end{array}$$

The characteristic polynomial for the normalized model corresponding to (3.12), calculated at the equilibrium point (0, 0), is defined as

$$\mathcal{F}\left(\lambda \right)={\lambda }^2-\mathcal{T}\left(\tilde{\mathcal{c}}\right)\lambda +\mathcal{D}\left(\tilde{\mathcal{c}}\right),$$

(3.14)

where

$$\mathcal{T}\left(\tilde{\mathcal{c}}\right)=2+\frac{-\mathcal{h}(1+{\mathcal{u}}^{*2})(1+{\mathcal{u}}^{*}(\tilde{\mathcal{c}}+{\mathcal{c}}_1+{\mathcal{u}}^{*}))+(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{\begin{array}{c}\multicolumn{1}{c}{(1+{\mathcal{u}}^{*2}{)}^2}\end{array}},$$

(3.15)

and

$$\mathcal{D}\left(\tilde{\mathcal{c}}\right)=\frac{(-1+\mathcal{h})(-1+(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}{\mathcal{u}}^{*}-{\mathcal{u}}^{*2})(1+{\mathcal{u}}^{*2})+(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{\begin{array}{c}\multicolumn{1}{c}{(1+{\mathcal{u}}^{*2}{)}^2}\end{array}}.$$

(3.16)

Let $(\mathcal{a},\mathcal{b},\mathcal{c},\mathcal{h})\in {\mathrm{\Psi }}_{\mathrm{NSB}}$ be given. The characteristic equation (3.14) has solutions (3.14) ${\lambda }_1$ and ${\lambda }_{2,}$ which are complex conjugates satisfying $\left|{\lambda }_{1,2}\right|=1$ . These roots are expressed as

$${\lambda }_{1,2}=\frac{\mathcal{T}\left(\tilde{\mathcal{c}}\right)}{2}\pm i\sqrt{\frac{4\mathcal{D}\left(\tilde{\mathcal{c}}\right)-{\mathcal{T}}^2\left(\tilde{\mathcal{c}}\right)}{4}}.$$

As a result, we derive

$$\left|{\lambda }_{1,2}\right|=\sqrt{\mathcal{D}\left(\tilde{\mathcal{c}}\right)}.$$

Differentiating the magnitudes of these eigenvalues with respect to $\tilde{\mathcal{c}}$ at $\tilde{\mathcal{c}}=0$ yields

$${\frac{d\left|{\lambda }_{1,2}\right|}{d\tilde{\mathcal{c}}}|}_{\tilde{\mathcal{c}}=0}=\frac{(-1+\mathcal{h})\mathcal{h}{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+\mathcal{h}(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{2\sqrt{\mathcal{D}\left(0\right)}}>0.$$

Additionally, we assume that $\mathcal{T}\left(0\right)\ne \{0,-1,2,-2\}$ , which ensures that ${\lambda }_1^m,{\lambda }_2^m\ne 1$ for all $\mathcal{m}=1,2,3,4$ at $\tilde{\mathcal{c}}=0$ , that is, the eigenvalues do not lie at the resonance points on the unit circle. This requirement is formally stated

$$\begin{array}{c}\multicolumn{1}{c}{2+\frac{-\mathcal{h}(1+{\mathcal{u}}^{*2})(1+{\mathcal{u}}^{*}(\tilde{\mathcal{c}}+{\mathcal{c}}_1+{\mathcal{u}}^{*}))+(\tilde{\mathcal{c}}+{\mathcal{c}}_1)\mathcal{h}(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}}{(1+{\mathcal{u}}^{*2}{)}^2}\ne 0,-1.}\end{array}$$

(3.17)

To derive the simplified form of a structure (3.12) on $\tilde{\mathcal{c}}=0$ , we set $\alpha =\frac{\mathcal{T}\left(0\right)}{2}$ and $\beta =\frac{1}{2}\sqrt{4\mathcal{D}\left(0\right)-{\mathcal{T}}^2\left(0\right)}$ , and apply the coordinate transformation

$$\left[\begin{array}{c}\mathcal{p}\hfill \\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right]=\mathrm{T}\left[\begin{array}{c}\mathcal{s}\hfill \\ \multicolumn{1}{c}{\mathcal{t}}\end{array}\right],$$

(3.18)

where $\mathrm{T}=\left[\begin{array}{cc}{\varrho }_2\hfill & 0\hfill \\ \multicolumn{1}{c}{\alpha -{\varrho }_1}& -\beta \hfill \end{array}\right]$ is an invertible matrix. By applying the transformation in (3.18), the map described in (3.12) can be re-expressed in a new form

$$\begin{array}{c}\multicolumn{1}{c}{\left[\genfrac{}{}{0ex}{}{\mathcal{s}}{\mathcal{t}}\right]\to \left[\begin{array}{cc}\alpha \hfill & -\beta \hfill \\ \multicolumn{1}{c}{\beta }& \alpha \hfill \end{array}\right]\left[\begin{array}{c}\mathcal{s}\hfill \\ \multicolumn{1}{c}{\mathcal{t}}\end{array}\right]+\left[\begin{array}{c}{\mathcal{f}}_1(\mathcal{s},\mathcal{t})\hfill \\ \multicolumn{1}{c}{{\mathcal{g}}_1(\mathcal{s},\mathcal{t})}\end{array}\right],}\end{array}$$

(3.19)

where

$$\{\begin{array}{c}\multicolumn{1}{c}{{\mathcal{f}}_1(\mathcal{s},\mathcal{t})=(-\frac{1+{\mathcal{u}}^{*2}}{{\mathcal{c}}_1\mathcal{h}{\mathcal{u}}^{*}}){\mathcal{f}}_0(\mathcal{p},\mathcal{q}),}\\ \multicolumn{1}{c}{{\mathcal{g}}_1(\mathcal{s},\mathcal{t})=\mathrm{K}{\mathcal{f}}_0(\mathcal{p},\mathcal{q})+\mathrm{E}{\mathcal{g}}_0(\mathcal{p},\mathcal{q}),}\end{array}$$

(3.20)

and the coefficients are

$$\begin{array}{c}\multicolumn{1}{c}{\mathrm{K}=\frac{-(1+{\mathcal{u}}^{*2}{)}^2+{\mathcal{c}}_1({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}-{\mathcal{v}}^{*}+{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}{{\mathcal{c}}_1{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})\sqrt{{\mathcal{h}}^2(-1-\frac{{\mathcal{c}}_1^2{({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+{\mathcal{v}}^{*}-{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}^2}{{(1+{\mathcal{u}}^{*2})}^4}+\frac{2{\mathcal{c}}_1({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}-{\mathcal{v}}^{*}+{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}{{(1+{\mathcal{u}}^{*2})}^2})}},}\\ \multicolumn{1}{c}{\mathrm{E}=-\frac{2}{\sqrt{{\mathcal{h}}^2(-1-\frac{{\mathcal{c}}_1^2{({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+{\mathcal{v}}^{*}-{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}^2}{{(1+{\mathcal{u}}^{*2})}^4}+\frac{2{\mathcal{c}}_1({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}-{\mathcal{v}}^{*}+{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}{{(1+{\mathcal{u}}^{*2})}^2})}},}\\ \multicolumn{1}{c}{\mathcal{p}=-\frac{{\mathcal{c}}_1\mathcal{hs}{\mathcal{u}}^{*}}{1+{\mathcal{u}}^{*2}},\mathcal{q}=\frac{1}{2}\mathcal{hs}(1-\frac{{\mathcal{c}}_1({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}-{\mathcal{v}}^{*}+{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}{{(1+{\mathcal{u}}^{*2})}^2})-\frac{1}{2}}\\ \multicolumn{1}{c}{\mathcal{t}\sqrt{{\mathcal{h}}^2(-1-\frac{{\mathcal{c}}_1^2{({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+{\mathcal{v}}^{*}-{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}^2}{{(1+{\mathcal{u}}^{*2})}^4}+\frac{2{\mathcal{c}}_1({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}-{\mathcal{v}}^{*}+{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}{{(1+{\mathcal{u}}^{*2})}^2})}.}\end{array}$$

Now we define the subsequent non-zero real number:

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{L}={\left[(-\mathrm{Re}\left(\frac{(1-2{\lambda }_1){\lambda }_2^2}{1-{\lambda }_1}{\mathfrak{X}}_{20}{\mathfrak{X}}_{11}\right)-\frac{1}{2}{\left|{\mathfrak{X}}_{11}\right|}^2-{\left|{\mathfrak{X}}_{02}\right|}^2+\mathrm{Re}\left({\lambda }_2{\mathfrak{X}}_{21}\right))\right]}_{\tilde{\mathcal{c}}=0}\ne 0,}\end{array}$$

(3.21)

where the coefficients are defined as:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathfrak{X}}_{20}=\frac{1}{8}[{\mathcal{f}}_{\mathcal{ss}}-{\mathcal{f}}_{\mathcal{tt}}+2{\mathcal{g}}_{\mathcal{st}}+i({\mathcal{g}}_{\mathcal{ss}}-{\mathcal{g}}_{\mathcal{tt}}+2{\mathcal{f}}_{\mathcal{st}})],}\\ \multicolumn{1}{c}{{\mathfrak{X}}_{11}=\frac{1}{4}[{\mathcal{f}}_{\mathcal{ss}}+{\mathcal{f}}_{\mathcal{tt}}+i({\mathcal{g}}_{\mathcal{ss}}+{\mathcal{g}}_{\mathcal{tt}})],}\\ \multicolumn{1}{c}{{\mathfrak{X}}_{02}=\frac{1}{8}[{\mathcal{f}}_{\mathcal{ss}}-{\mathcal{f}}_{\mathcal{tt}}-2{\mathcal{g}}_{\mathcal{st}}+i({\mathcal{g}}_{\mathcal{ss}}-{\mathcal{g}}_{\mathcal{tt}}+2{\mathcal{f}}_{\mathcal{st}})],}\\ \multicolumn{1}{c}{{\mathfrak{X}}_{21}=\frac{1}{16}[{\mathcal{f}}_{\mathcal{sss}}+{\mathcal{f}}_{\mathcal{stt}}+{\mathcal{g}}_{\mathcal{sst}}+{\mathcal{g}}_{\mathcal{ttt}}+i({\mathcal{g}}_{\mathcal{sss}}+{\mathcal{g}}_{\mathcal{stt}}-{\mathcal{f}}_{\mathcal{sst}}-{\mathcal{f}}_{\mathcal{ttt}})].}\end{array}$$

Based on the above analysis, we formulate the subsequent outcomes:

• Theorem 3.2. Suppose that $\mathcal{L}\ne 0$ and the non-resonance condition (3.17) are satisfied, then system (1.4) undergoes Neimark-Sacker bifurcation at $(\mathcal{a}-\mathcal{b},\frac{(1+{(\mathcal{a}-\mathcal{b})}^2)\mathcal{b}}{(\mathcal{a}-\mathcal{b})\mathcal{c}}),$ as the bifurcation parameter $\mathcal{c}$ fluctuates in the neighborhood of ${\mathcal{c}}_1=\frac{\mathcal{a}+{\mathcal{a}}^3-4{\mathcal{a}}^2\mathcal{b}+5\mathcal{a}{\mathcal{b}}^2-2{\mathcal{b}}^3}{{(\mathcal{a}-\mathcal{b})}^2(-1+\mathcal{h})}.$ Furthermore, if $\mathcal{L}<0$ , an attracting invariant closed curve bifurcates from the fixed point for $\mathcal{c}>{\mathcal{c}}_1$ . Conversely, if $\mathcal{L}>0$ , a repelling invariant closed curve bifurcates from the fixed point for $\mathcal{c}<{\mathcal{c}}_1$ .

Strong resonance 1:2

In the following segment, we apply bifurcation theory outlined in Ref. ²⁴ to examine the codimension-2 bifurcation related to 1:2 strong resonance at the interior fixed point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ of model (1.4). We select c and h as the bifurcation parameters. Under this setup, the parameter set corresponding to the 1:2 resonance is given by

$${\mathrm{\Psi }}_{1:2}=\left\{\begin{array}{c}\multicolumn{1}{c}{(\mathcal{a},\mathcal{b},\mathcal{c},\mathcal{h})\in {\mathbb{R}}_{+}^4:\mathcal{h}={\mathcal{h}}_2,\mathcal{c}={\mathcal{c}}_2}\end{array}\right\},$$

where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{h}}_2=\frac{2({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3+\eta )}{{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3},}\\ \multicolumn{1}{c}{{\mathcal{c}}_2=\frac{(1+{\mathcal{u}}^{*2}{)}^4}{-2(1+{\mathcal{u}}^{*2}{)}^2{\mathcal{v}}^{*}+2\eta +(1+{\mathcal{u}}^{*2}{)}^3({\mathcal{u}}^{*}+{\mathcal{v}}^{*})},}\end{array}$$

with $\eta =\sqrt{{\mathcal{u}}^{*}(-1+{\mathcal{u}}^{*2}){(1+{\mathcal{u}}^{*2})}^5{\mathcal{v}}^{*}}$ .

To explore the local bifurcation structure near this resonance, we introduce small perturbations $\tilde{\mathcal{h}}$ and $\tilde{\mathcal{c}}$ to the critical parameter values ${\mathcal{h}}_2$ and ${\mathcal{c}}_2$ , respectively, such that $\mathcal{h}={\mathcal{h}}_2+\tilde{\mathcal{h}}$ and $\mathcal{c}={\mathcal{c}}_2+\tilde{\mathcal{c}}$ . Substituting these into system (1.4), the resulting perturbed system is given by:

$$\left[\begin{array}{c}\mathcal{u}\hfill \\ \multicolumn{1}{c}{\mathcal{v}}\end{array}\right]\to \left[\begin{array}{c}\mathcal{u}+({\mathcal{h}}_2+\tilde{\mathcal{h}})(\mathcal{a}-\mathcal{u}-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)\mathcal{uv}}{1+{\mathcal{u}}^2}),\hfill \\ \multicolumn{1}{c}{\mathcal{v}+({\mathcal{h}}_2+\tilde{\mathcal{h}})(\mathcal{b}-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)\mathcal{uv}}{1+{\mathcal{u}}^2}).}\end{array}\right].$$

(3.22)

To facilitate a local analysis near the interior fixed point, the fixed point $({\mathcal{u}}^{**},{\mathcal{v}}^{**})=(\mathcal{a}-\mathcal{b},\frac{(1+{(\mathcal{a}-\mathcal{b})}^2)\mathcal{b}}{(\mathcal{a}-\mathcal{b})({\mathcal{c}}_0+\tilde{\mathcal{c}})})$ transforms to the origin by using the transformation $\mathcal{p}=\mathcal{u}-{\mathcal{u}}^{**}\mathrm{and}\mathcal{q}=\mathcal{v}-{\mathcal{v}}^{**}$ . Expanding system (3.22) in a Taylor series about the origin yields the following form:

$$\left[\begin{array}{c}\mathcal{p}\hfill \\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right]\to \left[\begin{array}{cc}{\varrho }_1\hfill & {\varrho }_2\hfill \\ \multicolumn{1}{c}{{\sigma }_1}& {\sigma }_2\hfill \end{array}\right]\left[\begin{array}{c}\mathcal{p}\hfill \\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right]+\left[\begin{array}{c}{\mathcal{f}}_0(\mathcal{p},\mathcal{q})\hfill \\ \multicolumn{1}{c}{{\mathcal{g}}_0(\mathcal{p},\mathcal{q})}\end{array}\right],$$

(3.23)

where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{f}}_0(\mathcal{p},\mathcal{q})={\vartheta }_{20}{\mathcal{p}}^2+{\vartheta }_{11}\mathcal{pq}+{\vartheta }_{02}{\mathcal{q}}^2+{\vartheta }_{30}{\mathcal{p}}^3+{\vartheta }_{21}{\mathcal{p}}^2\mathcal{q}}\\ \multicolumn{1}{c}{+{\vartheta }_{12}\mathcal{p}{\mathcal{q}}^2+{\vartheta }_{03}{\mathcal{q}}^3+\mathcal{O}\left({(\left|\mathcal{p}\right|+\left|\mathcal{q}\right|)}^4\right),}\\ \multicolumn{1}{c}{{\mathcal{g}}_0(\mathcal{p},\mathcal{q})={\sigma }_{20}{\mathcal{p}}^2+{\sigma }_{11}\mathcal{pq}+{\sigma }_{02}{\mathcal{q}}^2+{\sigma }_{30}{\mathcal{p}}^3+{\sigma }_{21}{\mathcal{p}}^2\mathcal{q}}\\ \multicolumn{1}{c}{+{\sigma }_{12}\mathcal{p}{\mathcal{q}}^2{\sigma }_{03}{\mathcal{q}}^3+\mathcal{O}\left({(\left|\mathcal{p}\right|+\left|\mathcal{q}\right|)}^4\right),}\end{array}$$

(3.24)

and the coefficients are

$$\begin{array}{c}\multicolumn{1}{c}{{\varrho }_1=1+(\tilde{\mathcal{h}}+{\mathcal{h}}_2)(-1+\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^2}),}\\ \multicolumn{1}{c}{{\varrho }_2=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}}{1+{\mathcal{u}}^{**2}},}\\ \multicolumn{1}{c}{{\sigma }_1=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2)(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{(1+{\mathcal{u}}^{**2}{)}^2},}\\ \multicolumn{1}{c}{{\sigma }_2=1-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}}{1+{\mathcal{u}}^{**2}},}\\ \multicolumn{1}{c}{{\vartheta }_{20}=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{(1+{\mathcal{u}}^{**2}{)}^3},}\\ \multicolumn{1}{c}{{\vartheta }_{11}=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2)(-1+{\mathcal{u}}^{**2})}{(1+{\mathcal{u}}^{**2}{)}^2},}\\ \multicolumn{1}{c}{{\vartheta }_{02}=0,{\vartheta }_{30}=(\frac{\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2)(1-6{\mathcal{u}}^{**2}+{\mathcal{u}}^{**4}){\mathcal{v}}^{**}}{(1+{\mathcal{u}}^{**2}{)}^4},}\\ \multicolumn{1}{c}{{\vartheta }_{12}=0,{\vartheta }_{03}=0,{\vartheta }_{21}=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2})}{{(1+{\mathcal{u}}^{**2})}^3},}\\ {\sigma }_{11}=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2)(-1+{\mathcal{u}}^{**2})}{{(1+{\mathcal{u}}^{**2})}^2},{\sigma }_{02}=0,{\sigma }_{30}=\frac{(\mathfrak{b}-\mathcal{d})(1+{\mathfrak{b}}^2+\mathcal{d}(3+\mathcal{d})-\mathfrak{b}(3+2\mathcal{d})){({\mathcal{c}}_2+\tilde{\mathcal{c}})}^3{\mathcal{v}}^{**}}{6},\\ \multicolumn{1}{c}{{\sigma }_{20}=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^3},}\\ \multicolumn{1}{c}{{\sigma }_{12}=0,{\sigma }_{03}=0,}\\ \multicolumn{1}{c}{{\sigma }_{21}=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2})}{{(1+{\mathcal{u}}^{**2})}^3}.}\end{array}$$

From the linearized system in (3.23), the Jacobian matrix $\mathcal{A}$ of the relation on translated equilibrium (0, 0) is expressed as

$$\mathcal{A}=\left(\begin{array}{cc}1+(\tilde{\mathcal{h}}+{\mathcal{h}}_2)(-1+\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^2})\hfill & -\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}}{1+{\mathcal{u}}^{**2}}\hfill \\ \multicolumn{1}{c}{\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2)(-1+{\mathcal{u}}^{**}{)}^2){\mathcal{v}}^{**}}{(1+{\mathcal{u}}^{**2}{)}^2}}& 1-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}}{1+{\mathcal{u}}^{**2}}\hfill \end{array}\right).$$

(3.25)

Evaluating this matrix at $\tilde{\mathcal{h}}=\tilde{\mathcal{c}}=0$ , we obtain the unperturbed Jacobian matrix

$${\mathcal{A}}_0=\left(\begin{array}{cc}{\theta }_{11}\hfill & {\theta }_{12}\hfill \\ \multicolumn{1}{c}{{\theta }_{21}}& {\theta }_{22}\hfill \end{array}\right),$$

(3.26)

where the coefficients are given by

$$\begin{array}{c}\multicolumn{1}{c}{{\theta }_{11}=\frac{-{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3+3(-1+{\mathcal{u}}^{*2})(1+{\mathcal{u}}^{*2}{)}^2{\mathcal{v}}^{*}-2\eta }{(1+{\mathcal{u}}^{*2}{)}^2({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+{\mathcal{v}}^{*}-{\mathcal{u}}^{*2}{\mathcal{v}}^{*})},}\\ \multicolumn{1}{c}{{\theta }_{12}=-\frac{2({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3-\eta )}{{(1+{\mathcal{u}}^{*2})}^2({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+{\mathcal{v}}^{*}-{\mathcal{u}}^{*2}{\mathcal{v}}^{*})},}\\ \multicolumn{1}{c}{{\theta }_{21}=\frac{2(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}({\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3+\eta )}{{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})(-2(1+{\mathcal{u}}^{*2}{)}^2{\mathcal{v}}^{*}+2\eta +(1+{\mathcal{u}}^{*2}{)}^3({\mathcal{u}}^{*}+{\mathcal{v}}^{*}))},}\\ \multicolumn{1}{c}{{\theta }_{22}=\frac{(1+{\mathcal{u}}^{*2}{)}^2({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+{\mathcal{v}}^{*}-{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}{{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3+(-1+{\mathcal{u}}^{*2})(1+{\mathcal{u}}^{*2}{)}^2{\mathcal{v}}^{*}+2\eta }.}\end{array}$$

The characteristic roots of ${\mathcal{A}}_0$ are ${\lambda }_{1,2}=-1$ . The eigenvector and generalized eigenvectors of ${\mathcal{A}}_0$ corresponding to eigenvalue -1 are ${\mathcal{E}}_1=\left(\begin{array}{c}-\frac{{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3}{\eta }\hfill \\ \multicolumn{1}{c}{1}\end{array}\right)$ and ${\mathcal{E}}_2=\left(\begin{array}{c}\frac{{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3}{2\eta }\hfill \\ \multicolumn{1}{c}{0}\end{array}\right)$ . Additionally, ${\mathcal{G}}_1=\left(\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{1}\end{array}\right)$ and ${\mathcal{G}}_2=\left(\begin{array}{c}\frac{2\eta }{{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3}\hfill \\ \multicolumn{1}{c}{2}\end{array}\right)$ are eigenvectors and generalized eigenvectors for ${\mathcal{A}}_0^T$ associated with eigenvalue $-1$ , where ${\mathcal{E}}_i\mathrm{and}{\mathcal{G}}_i(i=0,1)$ satisfy the following relations,

$$\{\begin{array}{c}\multicolumn{1}{c}{{\mathcal{A}}_0{\mathcal{E}}_1=-{\mathcal{E}}_1,}\\ \multicolumn{1}{c}{{\mathcal{A}}_0{\mathcal{E}}_2=-{\mathcal{E}}_2+{\mathcal{E}}_1,}\\ \multicolumn{1}{c}{{\mathcal{A}}_0^T{\mathcal{G}}_2=-{\mathcal{G}}_2,}\\ \multicolumn{1}{c}{{\mathcal{A}}_0^T{\mathcal{G}}_1=-{\mathcal{G}}_1+{\mathcal{G}}_2,}\\ \multicolumn{1}{c}{〈{\mathcal{E}}_1,{\mathcal{G}}_1〉=〈{\mathcal{E}}_2,{\mathcal{G}}_2〉=1,}\\ \multicolumn{1}{c}{〈{\mathcal{E}}_1,{\mathcal{G}}_2〉=〈{\mathcal{E}}_2,{\mathcal{G}}_1〉=0.}\end{array}$$

(3.27)

Now, if

$$\left[\begin{array}{c}\mathcal{p}\hfill \\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right]=\mathcal{m}{\mathcal{E}}_1+\mathcal{n}{\mathcal{E}}_2=\left[\begin{array}{cc}0\hfill & 1\hfill \\ \multicolumn{1}{c}{\frac{2\eta }{{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3}}& 2\hfill \end{array}\right]\left[\begin{array}{c}\mathcal{m}\hfill \\ \multicolumn{1}{c}{\mathcal{n}}\end{array}\right],$$

(3.28)

with $\mathcal{x}=\left(\begin{array}{c}\mathcal{p}\hfill \\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right)$ , then by a straightforward calculation, the coefficients $\mathcal{m}$ and $\mathcal{n}$ are given by,

$$\{\begin{array}{c}\multicolumn{1}{c}{\mathcal{m}=〈{\mathcal{x},\mathcal{G}}_1〉=\mathcal{q},}\\ \multicolumn{1}{c}{\mathcal{n}=〈{\mathcal{x},\mathcal{G}}_2〉=\frac{2(\eta \mathcal{p}+{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3\mathcal{q})}{{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3}.}\end{array}$$

Furthermore, in a more general setting, (3.23) can be rewritten as,

$$\left[\begin{array}{c}\mathcal{m}\hfill \\ \multicolumn{1}{c}{\mathcal{n}}\end{array}\right]\to \left[\begin{array}{cc}-1+{\kappa }_{10}\hfill & 1+{\kappa }_{01}\hfill \\ \multicolumn{1}{c}{{\tau }_{10}}& -1+{\tau }_{01}\hfill \end{array}\right]\left[\begin{array}{c}\mathcal{m}\hfill \\ \multicolumn{1}{c}{\mathcal{n}}\end{array}\right]+\left[\begin{array}{c}{\mathcal{f}}_1(\mathcal{m},\mathcal{n})\hfill \\ \multicolumn{1}{c}{{\mathcal{g}}_1(\mathcal{m},\mathcal{n})}\end{array}\right],$$

(3.29)

where

$$\{\begin{array}{c}\multicolumn{1}{c}{{\mathcal{f}}_1(\mathcal{m},\mathcal{n})={\mathcal{g}}_0(-\frac{(2\mathcal{m}-\mathcal{n}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3}{2\eta },\mathcal{m}),}\\ \multicolumn{1}{c}{{\mathcal{g}}_1(\mathcal{m},\mathcal{n})=\frac{2\eta }{{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3}{\mathcal{f}}_0(-\frac{(2\mathcal{m}-\mathcal{n}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3}{2\eta },\mathcal{m})}\\ \multicolumn{1}{c}{+2{\mathcal{g}}_0(-\frac{(2\mathcal{m}-\mathcal{n}){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3}{2\eta },\mathcal{m}),}\end{array}$$

(3.30)

and

$$\begin{array}{c}\multicolumn{1}{c}{{\kappa }_{10}=2-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}}{1+{\mathcal{u}}^{**2}}}\\ \multicolumn{1}{c}{-\frac{1}{{(1+{\mathcal{u}}^{**2})}^2\eta }((\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){(1+{\mathcal{u}}^{*2})}^3}\\ \multicolumn{1}{c}{{\mathcal{u}}^{*}(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}),}\\ \multicolumn{1}{c}{{\kappa }_{01}=-1+\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{2(1+{\mathcal{u}}^{**2}{)}^2\eta },}\\ \multicolumn{1}{c}{{\tau }_{10}=2(\tilde{\mathcal{h}}+{\mathcal{h}}_2)(1-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2){\mathcal{u}}^{**}}{1+{\mathcal{u}}^{**2}}-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2){\mathcal{u}}^{**}\eta }{{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3(1+{\mathcal{u}}^{**2})}}\\ \multicolumn{1}{c}{-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^2}-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2){\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^2\eta },}\\ \multicolumn{1}{c}{{\tau }_{01}=2+\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{(1+{\mathcal{u}}^{**2}{)}^2\eta }}\\ \multicolumn{1}{c}{+(\tilde{\mathcal{h}}+{\mathcal{h}}_2)(-1+\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^2}).}\end{array}$$

(3.31)

Moreover, the calculation shows that at $(\tilde{\mathcal{h}}=\tilde{\mathcal{c}}=0)$ , ${\kappa }_{10}={\kappa }_{01}={\tau }_{10}={\tau }_{01}=0$ . Furthermore, when defining

$$\mathcal{B}=\left(\begin{array}{cc}-1+{\kappa }_{01}\hfill & 0\hfill \\ \multicolumn{1}{c}{-{\kappa }_{10}}& 1\hfill \end{array}\right),$$

(3.32)

with

$$\left[\begin{array}{c}\multicolumn{1}{c}{\mathcal{m}}\\ \multicolumn{1}{c}{\mathcal{n}}\end{array}\right]=\mathcal{B}\left[\begin{array}{c}\multicolumn{1}{c}{\mathcal{s}}\\ \multicolumn{1}{c}{\mathcal{t}}\end{array}\right].$$

(3.33)

Map (3.29) can be written as

$$\left[\begin{array}{c}\mathcal{s}\hfill \\ \multicolumn{1}{c}{\mathcal{t}}\end{array}\right]\to \left[\begin{array}{cc}-1\hfill & 1\hfill \\ \multicolumn{1}{c}{\vartheta }& -1+\phi \hfill \end{array}\right]\left[\begin{array}{c}\mathcal{s}\hfill \\ \multicolumn{1}{c}{\mathcal{t}}\end{array}\right]+\left[\begin{array}{c}{\mathcal{f}}_2(\mathcal{s},\mathcal{t})\hfill \\ \multicolumn{1}{c}{{\mathcal{g}}_2(\mathcal{s},\mathcal{t})}\end{array}\right],$$

(3.34)

where

$$\begin{array}{c}\multicolumn{1}{c}{\vartheta =-\frac{1}{{(1+{\mathcal{u}}^{**2})}^2}((-2+\tilde{\mathcal{h}}+{\mathcal{h}}_2)}\\ \multicolumn{1}{c}{(-2+(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}-2{\mathcal{u}}^{**2})}\\ \multicolumn{1}{c}{(1+{\mathcal{u}}^{**2})+2(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2)(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}),}\\ \multicolumn{1}{c}{\phi =\frac{-1}{(1+{\mathcal{u}}^{**2}{)}^2}((1+{\mathcal{u}}^{**2})(-3+\tilde{\mathcal{h}}+{\mathcal{h}}_2}\\ \multicolumn{1}{c}{+(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}+(-3+\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**2}))}\\ \multicolumn{1}{c}{+(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2)(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**})+1,}\end{array}$$

(3.35)

and

$$\{\begin{array}{c}\multicolumn{1}{c}{{\mathcal{f}}_2(\mathcal{s},\mathcal{t})=\zeta {\mathcal{f}}_1(\mathcal{m},\mathcal{n})}\\ \multicolumn{1}{c}{{\mathcal{g}}_2(\mathcal{s},\mathcal{t})=\xi {\mathcal{f}}_1(\mathcal{m},\mathcal{n})+{\mathcal{g}}_1(\mathcal{m},\mathcal{n}),}\end{array}$$

(3.36)

where

$$\begin{array}{c}\multicolumn{1}{c}{\zeta =\frac{2{(1+{\mathcal{u}}^{**2})}^2\eta }{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}},}\\ \multicolumn{1}{c}{\xi =-2-\frac{2(-2+(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}-2{\mathcal{u}}^{**2})(1+{\mathcal{u}}^{**2})\eta }{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}},}\\ \multicolumn{1}{c}{m=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2)s{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{2(1+{\mathcal{u}}^{**2}{)}^2\eta },}\\ \multicolumn{1}{c}{\mathcal{n}=t+s(-2+\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){\mathcal{u}}^{**}}{1+{\mathcal{u}}^{**2}}+\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_2)(\tilde{\mathcal{h}}+{\mathcal{h}}_2){(1+{\mathcal{u}}^{*2})}^3(-1+{\mathcal{u}}^{**2}){\mathcal{u}}^{*}{\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^2\eta }).}\end{array}$$

Taking into account $\vartheta$ and $\phi$ , we define the following matrix

$$\det (\mathcal{C}({\tilde{\mathcal{h}},\tilde{\mathcal{c}})|}_{(0,0)})=\left(\begin{array}{cc}\frac{\partial \vartheta }{\partial \tilde{\mathcal{h}}}\hfill & \frac{\partial \vartheta }{\partial \tilde{\mathcal{c}}}\hfill \\ \multicolumn{1}{c}{\frac{\partial \phi }{\partial \tilde{\mathcal{h}}}}& \frac{\partial \phi }{\partial \tilde{\mathcal{c}}}\hfill \end{array}\right).$$

Then, by simple computation, $\det (\mathcal{C}({\tilde{\mathcal{h}},\tilde{\mathcal{c}})|}_{(0,0)})$ is obtain as

$$\begin{array}{c}\multicolumn{1}{c}{\det (\mathcal{C}({\tilde{\mathcal{h}},\tilde{\mathcal{c}})|}_{(0,0)})=}\\ \multicolumn{1}{c}{\frac{{\mathcal{h}}_2^2{\mathcal{u}}^{**}({(1+{\mathcal{u}}^{**2})}^2-{\mathcal{c}}_2({\mathcal{u}}^{**}+{\mathcal{u}}^{**3}+{\mathcal{v}}^{**}-{\mathcal{u}}^{**2}{\mathcal{v}}^{**}))}{(1+{\mathcal{u}}^{**2}{)}^3}\ne 0.}\end{array}$$

(3.37)

Condition (3.37) is referred to as the transversality assumption, and we assume it holds. Furthermore, if

$$\beta =\left(\begin{array}{c}{\beta }_1\hfill \\ \multicolumn{1}{c}{{\beta }_2}\end{array}\right)=\left(\begin{array}{c}\vartheta \hfill \\ \multicolumn{1}{c}{\phi }\end{array}\right),$$

(3.38)

then ${\beta }_1={\beta }_2=0$ , at $({\mathcal{h}}_2,{\mathcal{c}}_2)$ and $\tilde{\mathcal{c}}$ and $\tilde{\mathcal{h}}$ can be expressed as a function of ${\beta }_1$ and ${\beta }_2$ as follows:

$$\tilde{\mathcal{h}}=\frac{1}{2}(4-2{\mathcal{h}}_2-{\beta }_2-\frac{\Xi }{\sqrt{{\mathcal{u}}^{**}}\Upsilon }),$$

$$\begin{array}{c}\multicolumn{1}{c}{\tilde{\mathcal{c}}=\frac{1}{\left(2{({\mathcal{u}}^{**}+{\mathcal{u}}^{**3}+{\mathcal{v}}^{**}-{\mathcal{u}}^{**2}{\mathcal{v}}^{**})}^2\Lambda \right)}}\\ \multicolumn{1}{c}{(-2{\mathcal{u}}^{**6}({\mathcal{c}}_2-{\mathcal{v}}^{**})\Lambda -2{\mathcal{v}}^{**}(1+{\mathcal{c}}_2{\mathcal{v}}^{**})\Lambda }\\ \multicolumn{1}{c}{-2{\mathcal{u}}^{**4}(-{\mathcal{v}}^{**}+{\mathcal{c}}_2(2+{\mathcal{v}}^{**2}))\Lambda +2{\mathcal{u}}^{**2}}\\ \multicolumn{1}{c}{(-{\mathcal{v}}^{**}+{\mathcal{c}}_2(-1+2{\mathcal{v}}^{**2}))\Lambda -3{\mathcal{u}}^{**3}}\\ \multicolumn{1}{c}{(8+2{\beta }_1+(-4+{\beta }_2){\beta }_2)-{\mathcal{u}}^{**7}(8+2{\beta }_1+(-4+{\beta }_2){\beta }_2)}\\ \multicolumn{1}{c}{-{\mathcal{u}}^{**}(8+2{\beta }_1+(-4+{\beta }_2){\beta }_2+4{\mathcal{c}}_2{\mathcal{v}}^{**}\Lambda )}\\ \multicolumn{1}{c}{+{\mathcal{u}}^{**5}(4{\mathcal{c}}_2{\mathcal{v}}^{**}\Lambda -3(8+2{\beta }_1+(-4+{\beta }_2){\beta }_2))+\sqrt{{\mathcal{u}}^{**}}}\\ \multicolumn{1}{c}{\Upsilon (-4+{\beta }_2)\Xi 2{\mathcal{u}}^{**5/2}\Upsilon (-4+{\beta }_2)\Xi +{\mathcal{u}}^{**9/2}\Upsilon (-4+{\beta }_2)\Xi ),}\end{array}$$

(3.39)

where $\Upsilon =\sqrt{1+{\mathcal{u}}^{**2}}$ , $\Lambda =(-4+{\beta }_1+2{\beta }_2)$ and

$\Xi =\sqrt{-4(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}\Lambda +{\mathcal{u}}^{**}(1+{\mathcal{u}}^{**2})(4{\beta }_1+{\beta }_2^2)}$ . Using (3.38) in (3.34), we have the following mapping:

$$\left[\begin{array}{c}\multicolumn{1}{c}{\mathcal{s}}\\ \multicolumn{1}{c}{\mathcal{t}}\end{array}\right]\to \left[\begin{array}{cc}-1\hfill & 1\hfill \\ \multicolumn{1}{c}{{\beta }_1}& -1+{\beta }_2\hfill \end{array}\right]\left[\begin{array}{c}\multicolumn{1}{c}{\mathcal{s}}\\ \multicolumn{1}{c}{\mathcal{t}}\end{array}\right]+\left[\begin{array}{c}{\mathcal{f}}_3(\mathcal{s},\mathcal{t})\hfill \\ \multicolumn{1}{c}{{\mathcal{g}}_3(\mathcal{s},\mathcal{t})}\end{array}\right],$$

(3.40)

where

$$\{\begin{array}{c}\multicolumn{1}{c}{{\mathcal{f}}_3(\mathcal{s},\mathcal{t})={\mathcal{g}}_{20}{\mathcal{s}}^2+{\mathcal{g}}_{11}\mathcal{st}+{\mathcal{g}}_{02}{\mathcal{t}}^2+{\mathcal{g}}_{30}{\mathcal{s}}^3+{\mathcal{g}}_{21}{\mathcal{s}}^2\mathcal{t}}\\ \multicolumn{1}{c}{+{\mathcal{g}}_{12}\mathcal{s}{\mathcal{t}}^2+{\mathcal{g}}_{03}{\mathcal{t}}^3+\mathcal{O}\left({(\left|\mathcal{s}\right|+\left|\mathcal{t}\right|)}^4\right),}\\ \multicolumn{1}{c}{{\mathcal{g}}_3(\mathcal{s},\mathcal{t})={\mathcal{h}}_{20}{\mathcal{s}}^2+{\mathcal{h}}_{11}\mathcal{st}+{\mathcal{h}}_{02}{\mathcal{t}}^2+{\mathcal{h}}_{30}{\mathcal{s}}^3+{\mathcal{h}}_{21}{\mathcal{s}}^2\mathcal{t}}\\ \multicolumn{1}{c}{+{\mathcal{h}}_{12}\mathcal{s}{\mathcal{t}}^2+{\mathcal{h}}_{03}{\mathcal{t}}^3+\mathcal{O}\left({(\left|\mathcal{s}\right|+\left|\mathcal{t}\right|)}^4\right),}\end{array}$$

(3.41)

and the coefficients at $\beta =0$ are

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{g}}_{20}=\frac{1}{\left(\sqrt{{\mathcal{u}}^{**}}(-1+{\mathcal{u}}^{**4})\eta {({\mathcal{u}}^{**}+{\mathcal{u}}^{**3}+{\mathcal{v}}^{**}-{\mathcal{u}}^{**2}{\mathcal{v}}^{**})}^2\right)}}\\ \multicolumn{1}{c}{(2{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3((-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}\\ \multicolumn{1}{c}{+\sqrt{{\mathcal{u}}^{**}}\mathrm{\Upsilon }\mathrm{\Theta }(\Upsilon \Theta +\sqrt{{\mathcal{u}}^{**}}(1-{\mathcal{u}}^{**2}-{\mathcal{u}}^{**4}+{\mathcal{u}}^{**6}}\\ \multicolumn{1}{c}{-3{\mathcal{u}}^{**}{\mathcal{v}}^{**}+4{\mathcal{u}}^{**3}{\mathcal{v}}^{**}-{\mathcal{u}}^{**5}{\mathcal{v}}^{**}+{\mathcal{u}}^{**3/2}\mathrm{\Upsilon }\mathrm{\Theta }))),}\\ \multicolumn{1}{c}{{\mathcal{g}}_{02}=-\frac{{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3{\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2})}{2(-1+{\mathcal{u}}^{**4})\eta },}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{g}}_{11}=\frac{1}{\sqrt{{\mathcal{u}}^{**}}(-1+{\mathcal{u}}^{**4})\eta ({\mathcal{u}}^{**}+{\mathcal{u}}^{**3}+{\mathcal{v}}^{**}-{\mathcal{u}}^{**2}{\mathcal{v}}^{**})}}\\ \multicolumn{1}{c}{({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3(\Upsilon \Theta +\sqrt{{\mathcal{u}}^{**}}(1-{\mathcal{u}}^{**2}}\\ \multicolumn{1}{c}{-{\mathcal{u}}^{**4}+{\mathcal{u}}^{**6}-6{\mathcal{u}}^{**}{\mathcal{v}}^{**}+8{\mathcal{u}}^{**3}{\mathcal{v}}^{**}-2{\mathcal{u}}^{**5}{\mathcal{v}}^{**})}\\ \multicolumn{1}{c}{(+4{\mathcal{u}}^{**3/2}\mathrm{\Upsilon }\mathrm{\Theta }-{\mathcal{u}}^{**7/2}\mathrm{\Upsilon }\mathrm{\Theta })),}\\ \multicolumn{1}{c}{{\mathcal{g}}_{21}=\frac{-1}{(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}{({\mathcal{u}}^{**}+2{\mathcal{u}}^{**3}+{\mathcal{u}}^{**5}+{\mathcal{v}}^{**}-{\mathcal{u}}^{**4}{\mathcal{v}}^{**})}^2}}\\ \multicolumn{1}{c}{{(\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})+(9{\mathcal{u}}^{**}{\mathcal{v}}^{**}+11{\mathcal{u}}^{**3}}\\ \multicolumn{1}{c}{{\mathcal{v}}^{**}+3{\mathcal{u}}^{**5}{\mathcal{v}}^{**}{\mathcal{u}}^{**7})+3{\mathcal{v}}^{**2}-21{\mathcal{u}}^{**2}{\mathcal{v}}^{**2}+21{\mathcal{u}}^{**4}{\mathcal{v}}^{**2}}\\ \multicolumn{1}{c}{-3{\mathcal{u}}^{**6}{\mathcal{v}}^{**2}-4{\mathcal{u}}^{**7/2}-6{\mathcal{u}}^{**3/2}\mathrm{\Upsilon }\mathrm{\Theta }}\\ \multicolumn{1}{c}{+2{\mathcal{u}}^{**11/2}\mathrm{\Upsilon }\mathrm{\Theta }+28{\mathcal{u}}^{**5/2}\Upsilon {\mathcal{v}}^{**}\Theta -4{\mathcal{u}}^{**9/2}\Upsilon {\mathcal{v}}^{**}\Theta ),}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{g}}_{12}=\frac{-1}{2(-1+{\mathcal{u}}^{*2})(-1+{\mathcal{u}}^{**2}){(1+{\mathcal{u}}^{**2})}^2{\mathcal{v}}^{*}({\mathcal{u}}^{**}+{\mathcal{u}}^{**3}+{\mathcal{v}}^{**}-{\mathcal{u}}^{**2}{\mathcal{v}}^{**})}}\\ \multicolumn{1}{c}{({\mathcal{u}}^{**}(3-{\mathcal{u}}^{**2}-({\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}3{\mathcal{u}}^{**4}+{\mathcal{u}}^{**6})}\\ \multicolumn{1}{c}{-3(-1+7{\mathcal{u}}^{**2}-7{\mathcal{u}}^{**4}+{\mathcal{u}}^{**6}){\mathcal{v}}^{**}-14{\mathcal{u}}^{**5/2}\mathrm{\Upsilon }\mathrm{\Theta }-2{\mathcal{u}}^{**9/2}\mathrm{\Upsilon }\mathrm{\Theta }),}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{h}}_{20}=\frac{1}{\left({\mathcal{u}}^{**}{(1+{\mathcal{u}}^{**2})}^{5/2}\eta {({\mathcal{u}}^{**}+{\mathcal{u}}^{**3}+{\mathcal{v}}^{**}-{\mathcal{u}}^{**2}{\mathcal{v}}^{**})}^2\right)}}\\ \multicolumn{1}{c}{(-((4{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3((-1+{\mathcal{u}}^{**2})\Upsilon {\mathcal{v}}^{**}}\\ \multicolumn{1}{c}{+\sqrt{{\mathcal{u}}^{**}}(1+{\mathcal{u}}^{**2})\Theta )(\sqrt{{\mathcal{u}}^{**}}\Upsilon (-1+{\mathcal{u}}^{**4})\Theta }\\ \multicolumn{1}{c}{+{\mathcal{v}}^{**}(1+2{\mathcal{u}}^{**2}+{\mathcal{u}}^{**4}+3{\mathcal{u}}^{**3/2}\mathrm{\Upsilon }\mathrm{\Theta }-{\mathcal{u}}^{**7/2}\mathrm{\Upsilon }\mathrm{\Theta }))),}\\ \multicolumn{1}{c}{{\mathcal{h}}_{11}=\frac{1}{{\mathcal{u}}^{**}{(1+{\mathcal{u}}^{**2})}^2\eta ({\mathcal{u}}^{**}+{\mathcal{u}}^{**3}+{\mathcal{v}}^{**}-{\mathcal{u}}^{**2}{\mathcal{v}}^{**})}}\\ \multicolumn{1}{c}{(2{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3(-\sqrt{{\mathcal{u}}^{**}}\Upsilon (-1+{\mathcal{u}}^{**4})\Theta +}\\ \multicolumn{1}{c}{{\mathcal{v}}^{**}(-1-5{\mathcal{u}}^{**2}-3{\mathcal{u}}^{**4}+{\mathcal{u}}^{**6}-6{\mathcal{u}}^{**3/2}\mathrm{\Upsilon }\mathrm{\Theta }+2{\mathcal{u}}^{**7/2}\mathrm{\Upsilon }\mathrm{\Theta })),}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{h}}_{02}=\frac{{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3\sqrt{{\mathcal{u}}^{**}}(-3+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^{3/2}\eta \Theta },}\\ \multicolumn{1}{c}{{\mathcal{h}}_{21}=\frac{1}{\left((-1+{\mathcal{u}}^{*2})\sqrt{{\mathcal{u}}^{**}}{(1+{\mathcal{u}}^{**2})}^{7/2}{\mathcal{v}}^{*}{({\mathcal{u}}^{*}{\mathcal{u}}^{**}+{\mathcal{v}}^{**}-{\mathcal{u}}^{**2}{\mathcal{v}}^{**})}^2\right)}}\\ \multicolumn{1}{c}{(2{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})((-1+{\mathcal{u}}^{**2})\mathrm{\Upsilon }}\\ \multicolumn{1}{c}{{\mathcal{v}}^{**}+\sqrt{{\mathcal{u}}^{**}}(1+{\mathcal{u}}^{**2})\mathrm{\Theta })(3\sqrt{{\mathcal{u}}^{**}}{\mathcal{v}}^{**}+13{\mathcal{u}}^{**5/2}{\mathcal{v}}^{**}}\\ \multicolumn{1}{c}{+9{\mathcal{u}}^{**9/2}{\mathcal{v}}^{**}-{\mathcal{u}}^{**13/2}{\mathcal{v}}^{**}-6{\mathcal{u}}^{**}\mathrm{\Upsilon }\mathrm{\Theta }-4{\mathcal{u}}^{**3}}\\ \multicolumn{1}{c}{\Upsilon \Theta +2{\mathcal{u}}^{**5}\mathrm{\Upsilon }\mathrm{\Theta }-3\Upsilon {\mathcal{v}}^{**}\Theta +18{\mathcal{u}}^{**2}\Upsilon {\mathcal{v}}^{**}\Theta -3{\mathcal{u}}^{**4}\Upsilon {\mathcal{v}}^{**}\Theta )),}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{h}}_{12}=\frac{-1}{\left((-1+{\mathcal{u}}^{*2})\sqrt{{\mathcal{u}}^{**}}{(1+{\mathcal{u}}^{**2})}^{5/2}{\mathcal{v}}^{*}({\mathcal{u}}^{**}+{\mathcal{u}}^{**3}+{\mathcal{v}}^{**}-{\mathcal{u}}^{**2}{\mathcal{v}}^{**})\right)}}\\ \multicolumn{1}{c}{({\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2})(14-4{\mathcal{u}}^{**5/2}\Upsilon {\mathcal{v}}^{**}-2{\mathcal{u}}^{**9/2}\Upsilon {\mathcal{v}}^{**}-3{\mathcal{u}}^{**}\Theta }\\ \multicolumn{1}{c}{-2{\mathcal{u}}^{**3}\Theta +{\mathcal{u}}^{**5}\Theta -3{\mathcal{v}}^{**}\Theta +18{\mathcal{u}}^{**2}{\mathcal{v}}^{**}\Theta -3{\mathcal{u}}^{**4}{\mathcal{v}}^{**}\Theta )).}\end{array}$$

where $\Theta =\sqrt{(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}$ . In accordance with Lemma 9.9 of Ref., ²⁴ there exists an identity transformation under which system (3.40) can be equivalently reformulated as follows:

$$\left[\genfrac{}{}{0ex}{}{\mathcal{w}}{\mathcal{z}}\right]\to \left[\begin{array}{cc}-1\hfill & 1\hfill \\ \multicolumn{1}{c}{{\beta }_1}& {\beta }_2-1\hfill \end{array}\right]\left[\genfrac{}{}{0ex}{}{\mathcal{w}}{\mathcal{z}}\right]+\left[\begin{array}{c}0\hfill \\ \multicolumn{1}{c}{\mathcal{M}\left(\beta \right){\mathcal{w}}^3+\mathcal{N}(\beta ){\mathcal{w}}^2\mathcal{z}}\end{array}\right],$$

(3.42)

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{M}({\beta }_1,{\beta }_2)={\mathcal{g}}_{20}({\beta }_1,{\beta }_2){\mathcal{h}}_{20}({\beta }_1,{\beta }_2)}\\ \multicolumn{1}{c}{+\frac{1}{2}{\mathcal{h}}_{20}^2({\beta }_1,{\beta }_2)+\frac{1}{2}{\mathcal{h}}_{20}({\beta }_1,{\beta }_2){\mathcal{h}}_{11}({\beta }_1,{\beta }_2),}\\ \multicolumn{1}{c}{\mathcal{N}({\beta }_1,{\beta }_2)=\frac{1}{2}{\mathcal{g}}_{20}({\beta }_1,{\beta }_2){\mathcal{h}}_{11}({\beta }_1,{\beta }_2)}\\ \multicolumn{1}{c}{+\frac{5}{4}{\mathcal{h}}_{20}({\beta }_1,{\beta }_2){\mathcal{h}}_{11}({\beta }_1,{\beta }_2)+{\mathcal{h}}_{20}^2({\beta }_1,{\beta }_2)}\\ \multicolumn{1}{c}{+\frac{1}{2}{\mathcal{h}}_{11}^2({\beta }_1,{\beta }_2)+{\mathcal{h}}_{20}({\beta }_1,{\beta }_2){\mathcal{h}}_{02}({\beta }_1,{\beta }_2)}\\ \multicolumn{1}{c}{+3{\mathcal{g}}_{20}^2({\beta }_1,{\beta }_2)+\frac{5}{2}{\mathcal{g}}_{20}({\beta }_1,{\beta }_2){\mathcal{h}}_{20}({\beta }_1,{\beta }_2)}\\ \multicolumn{1}{c}{+\frac{5}{2}{\mathcal{g}}_{11}({\beta }_1,{\beta }_2){\mathcal{h}}_{20}({\beta }_1,{\beta }_2).}\end{array}$$

(3.43)

In summary the aforementioned analysis in the resulting theorem.

• Theorem 3.3. Assume the conditions $\mathcal{M}(0,0)\ne 0$ , $\mathcal{N}(0,0)+3\mathcal{M}(0,0)\ne 0$ and $\det (\mathcal{C}({\tilde{\mathcal{h}},\tilde{\mathcal{c}})|}_{(0,0)})\ne 0$ , are satisfied. Then, the discrete model (1.4) shows the following bifurcation phenomena near the trivial fixed point (0, 0) as:

•  (i) A non-degenerate flip bifurcation takes place along the curve $\{({\epsilon }_1,{\epsilon }_2)\in {\mathbb{R}}^2:{\epsilon }_2\ne 0,{\epsilon }_1=0\}$ .

•  (ii) A non-degenerate Neimark-Sacker bifurcation happens along the curve $\{({\epsilon }_1,{\epsilon }_2)\in {\mathbb{R}}^2:{\epsilon }_1+{\epsilon }_2=0,{\epsilon }_1<0\}$ , and this bifurcation is supercritical if $\mathcal{N}(0,0)+3\mathcal{M}(0,0)>0$ .

•  (iii) Moreover, if $\mathcal{M}(0,0)<0$ , then there occurs a bifurcation loop $\{({\epsilon }_1,{\epsilon }_2)\in {\mathbb{R}}^2:{\epsilon }_2=(\frac{\mathcal{N}(0,0)}{\mathcal{M}(0,0)}+2){\epsilon }_1+O\left({\epsilon }_1^2\right),{\epsilon }_1>0\}$ alongside the 2-period orbits experiences a Neimark-Sacker bifurcation, if $\mathcal{N}(0,0)+3\mathcal{M}(0,0)<0$ then bifurcation is said to be supercritical.

Strong resonance 1:4

In the subsequent subsection, we explore the codimension-2 bifurcation related to $1:4$ strong resonance in structure (1.4) at fixed point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ . We select $\mathcal{c}$ and $\mathcal{h}$ as the bifurcation parameters. Under this formulation, the parameter set corresponding to the 1:2 resonances is given by

$${\mathrm{\Psi }}_{1:2}=\left\{\begin{array}{c}\multicolumn{1}{c}{(\mathcal{a},\mathcal{b},\mathcal{c},\mathcal{h})\in {\mathbb{R}}_{+}^4:\mathcal{h}={\mathcal{h}}_4,\mathcal{c}={\mathcal{c}}_4}\end{array}\right\},$$

where ${\mathcal{h}}_4=\frac{{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3-\mathfrak{A}}{{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3},{\mathcal{c}}_4=\frac{(-1+{\mathcal{u}}^{*2})(1+{\mathcal{u}}^{*2}{)}^2{\mathcal{v}}^{*}+\mathfrak{A}}{{({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+{\mathcal{v}}^{*}-{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}^2},$

with $\mathfrak{A}=\sqrt{-{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^5({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+2{\mathcal{v}}^{*}-2{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}$ .

To explore the local bifurcation structure near this resonance, we introduce small perturbations $\tilde{\mathcal{h}}$ and $\tilde{\mathcal{c}}$ to the critical parameter values ${\mathcal{h}}_4$ and ${\mathcal{c}}_4$ , respectively, such that $\mathcal{h}={\mathcal{h}}_4+\tilde{\mathcal{h}}$ and $\mathcal{c}={\mathcal{c}}_4+\tilde{\mathcal{c}}$ . Substituting these into system (1.4), the resulting perturbed system is given by

$$\left[\begin{array}{c}\mathcal{u}\hfill \\ \multicolumn{1}{c}{\mathcal{v}}\end{array}\right]\to \left[\begin{array}{c}\mathcal{u}+({\mathcal{h}}_4+\tilde{\mathcal{h}})(\mathcal{a}-\mathcal{u}-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)\mathcal{uv}}{1+{\mathcal{u}}^2}),\hfill \\ \multicolumn{1}{c}{\mathcal{v}+({\mathcal{h}}_4+\tilde{\mathcal{h}})(\mathcal{b}-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)\mathcal{uv}}{1+{\mathcal{u}}^2}).}\end{array}\right].$$

(3.44)

To facilitate a local analysis near the interior fixed point, we transform the equilibrium point $({\mathcal{u}}^{**},{\mathcal{v}}^{**})=(\mathcal{a}-\mathcal{b},\frac{(1+{(\mathcal{a}-\mathcal{b})}^2)\mathcal{b}}{(\mathcal{a}-\mathcal{b})({\mathcal{c}}_0+\tilde{\mathcal{c}})})$ to the origin by using the following transformation $\mathcal{p}=\mathcal{u}-{\mathcal{u}}^{**}\mathrm{and}\mathcal{q}=\mathcal{v}-{\mathcal{v}}^{**}$ . Expanding the system (3.44) in a Taylor series about the origin yields the following form:

$$\left[\begin{array}{c}\mathcal{p}\hfill \\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right]\to \left[\begin{array}{cc}{\varrho }_1\hfill & {\varrho }_2\hfill \\ \multicolumn{1}{c}{{\sigma }_1}& {\sigma }_2\hfill \end{array}\right]\left[\begin{array}{c}\mathcal{p}\hfill \\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right]+\left[\begin{array}{c}{\mathcal{f}}_0(\mathcal{p},\mathcal{q})\hfill \\ \multicolumn{1}{c}{{\mathcal{g}}_0(\mathcal{p},\mathcal{q})}\end{array}\right],$$

(3.45)

where

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{f}}_0(\mathcal{p},\mathcal{q})={\vartheta }_{20}{\mathcal{p}}^2+{\vartheta }_{11}\mathcal{pq}+{\vartheta }_{02}{\mathcal{q}}^2+{\vartheta }_{30}{\mathcal{p}}^3}\\ \multicolumn{1}{c}{+{\vartheta }_{21}{\mathcal{p}}^2\mathcal{q}+{\vartheta }_{12}\mathcal{p}{\mathcal{q}}^2+{\vartheta }_{03}{\mathcal{q}}^3+\mathcal{O}\left({(\left|\mathcal{p}\right|+\left|\mathcal{q}\right|)}^4\right),}\\ \multicolumn{1}{c}{{\mathcal{g}}_0(\mathcal{p},\mathcal{q})={\sigma }_{20}{\mathcal{p}}^2+{\sigma }_{11}\mathcal{pq}+{\sigma }_{02}{\mathcal{q}}^2+{\sigma }_{30}{\mathcal{p}}^3}\\ \multicolumn{1}{c}{+{\sigma }_{21}{\mathcal{p}}^2\mathcal{q}+{\sigma }_{12}\mathcal{p}{\mathcal{q}}^2{\sigma }_{03}{\mathcal{q}}^3+\mathcal{O}\left({(\left|\mathcal{p}\right|+\left|\mathcal{q}\right|)}^4\right),}\end{array}$$

(3.46)

$$\begin{array}{c}\multicolumn{1}{c}{{\varrho }_1=1+(\tilde{\mathcal{h}}+{\mathcal{h}}_4)(-1+\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^2}),}\\ {\varrho }_2=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4){\mathcal{u}}^{**}}{1+{\mathcal{u}}^{**2}},\hfill \\ {\sigma }_1=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4)(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{(1+{\mathcal{u}}^{**2}{)}^2},\hfill \\ {\sigma }_2=1-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4){\mathcal{u}}^{**}}{1+{\mathcal{u}}^{**2}},\hfill \\ {\vartheta }_{20}=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4){\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^3},\hfill \\ {\vartheta }_{11}=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4)(-1+{\mathcal{u}}^{**2})}{{(1+{\mathcal{u}}^{**2})}^2},\hfill \\ \multicolumn{1}{c}{{\vartheta }_{02}=0,{\vartheta }_{30}=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4)(1-6{\mathcal{u}}^{**2}+{\mathcal{u}}^{**4}){\mathcal{v}}^{**}}{(1+{\mathcal{u}}^{**2}{)}^4},}\\ \multicolumn{1}{c}{{\vartheta }_{12}=0,{\vartheta }_{03}=0,}\\ \multicolumn{1}{c}{{\vartheta }_{21}=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4){\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2})}{{(1+{\mathcal{u}}^{**2})}^3},}\\ \multicolumn{1}{c}{{\sigma }_{11}=\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4)(-1+{\mathcal{u}}^{**2})}{{(1+{\mathcal{u}}^{**2})}^2},}\\ \multicolumn{1}{c}{{\sigma }_{02}=0,}\\ \multicolumn{1}{c}{{\sigma }_{30}=\frac{(\mathfrak{b}-\mathcal{d})(1+{\mathfrak{b}}^2+\mathcal{d}(3+\mathcal{d})-\mathfrak{b}(3+2\mathcal{d})){({\mathcal{c}}_4+\tilde{\mathcal{c}})}^3{\mathcal{v}}^{**}}{6},}\\ \multicolumn{1}{c}{{\sigma }_{20}=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4){\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^3},}\\ \multicolumn{1}{c}{{\sigma }_{12}=0,{\sigma }_{03}=0,}\\ {\sigma }_{21}=-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4){\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2})}{{(1+{\mathcal{u}}^{**2})}^3}.\hfill \end{array}$$

From the linearized system in (3.45), the Jacobian matrix $\mathcal{A}$ for a structure at the translated equilibrium point (0,0) takes the form

$$\mathcal{A}=\left(\begin{array}{cc}1+(\tilde{\mathcal{h}}+{\mathcal{h}}_4)(-1+\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^2})\hfill & -\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4){\mathcal{u}}^{**}}{1+{\mathcal{u}}^{**2}}\hfill \\ \multicolumn{1}{c}{\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4)(-1+{\mathcal{u}}^{**2}){\mathcal{v}}^{**}}{{(1+{\mathcal{u}}^{**2})}^2}}& 1-\frac{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4){\mathcal{u}}^{**}}{1+{\mathcal{u}}^{**2}}\hfill \end{array}\right).$$

(3.47)

By evaluating this matrix at $\tilde{\mathcal{h}}=\tilde{\mathcal{c}}=0$ , we obtain the unperturbed Jacobian matrix

$${\mathcal{A}}_2=\left(\begin{array}{cc}{\theta }_{11}\hfill & {\theta }_{12}\hfill \\ \multicolumn{1}{c}{{\theta }_{21}}& {\theta }_{22}\hfill \end{array}\right),$$

(3.48)

where the coefficients are given by

$$\begin{array}{c}\multicolumn{1}{c}{{\theta }_{11}=\frac{(-1+{\mathcal{u}}^{*2}){(1+{\mathcal{u}}^{*2})}^2{\mathcal{v}}^{*}+\mathfrak{A}}{(1+{\mathcal{u}}^{*2}{)}^2({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+{\mathcal{v}}^{*}-{\mathcal{u}}^{*2}{\mathcal{v}}^{*})},}\\ \multicolumn{1}{c}{{\theta }_{12}=-\frac{{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3+\mathfrak{A}}{{(1+{\mathcal{u}}^{*2})}^2({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+{\mathcal{v}}^{*}-{\mathcal{u}}^{*2}{\mathcal{v}}^{*})},}\\ \multicolumn{1}{c}{{\theta }_{21}=\frac{(-1+{\mathcal{u}}^{*2}){\mathcal{v}}^{*}({\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3+\mathfrak{A})}{{\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+{\mathcal{v}}^{*}-{\mathcal{u}}^{*2}{\mathcal{v}}^{*})},}\\ \multicolumn{1}{c}{{\theta }_{22}=\frac{(-1+{\mathcal{u}}^{*2}){(1+{\mathcal{u}}^{*2})}^2{\mathcal{v}}^{*}+\mathfrak{A}}{(1+{\mathcal{u}}^{*2}{)}^2{\mathcal{v}}^{*}({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+{\mathcal{v}}^{*}-{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}.}\end{array}$$

${\lambda }_{1,2}=\pm i$ are eigenvalues of ${\mathcal{A}}_0$ . The eigenvectors of ${\mathcal{A}}_2$ corresponding to eigenvalues are

$\mathcal{E}=\left(\begin{array}{c}\multicolumn{1}{c}{\frac{(1-i)({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3-\mathfrak{Z})}{(-1+{\mathcal{u}}^{*2}){(1+{\mathcal{u}}^{*2})}^2{\mathcal{v}}^{*}}}\\ \multicolumn{1}{c}{2}\end{array}\right)$ and $\mathcal{G}=\left(\begin{array}{c}\multicolumn{1}{c}{\frac{i({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3+\mathfrak{Z})}{4{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^2}}\\ \multicolumn{1}{c}{\frac{1-i}{4}}\end{array}\right)$ , with $\mathfrak{Z}=\sqrt{{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^5({\mathcal{u}}^{*}+{\mathcal{u}}^{*3}+2{\mathcal{v}}^{*}-2{\mathcal{u}}^{*2}{\mathcal{v}}^{*})}$ , satisfying the following relations:

$$\{\begin{array}{c}\multicolumn{1}{c}{{\mathcal{A}}_2\mathcal{E}=i\mathcal{E},}\\ \multicolumn{1}{c}{{\mathcal{A}}_2\overline{\mathcal{E}}=-i\overline{\mathcal{E}},}\\ \multicolumn{1}{c}{{\mathcal{A}}_2^T\mathcal{G}=-i\mathcal{G},}\\ \multicolumn{1}{c}{{\mathcal{A}}_2^T\overline{\mathcal{G}}=i\overline{\mathcal{G}},}\\ \multicolumn{1}{c}{〈\mathcal{E},\mathcal{G}〉=1.}\end{array}$$

(3.49)

Now, if $\mathcal{Z}=\left(\begin{array}{c}\mathcal{p}\hfill \\ \multicolumn{1}{c}{\mathcal{q}}\end{array}\right)\in {\mathbb{R}}^2$ where $\mathcal{Z}=\mathcal{zE}+\overline{\mathcal{z}}\overline{\mathcal{E}}$ , can represented it, then becomes

$$\mathcal{z}\to i\mathcal{z}+\mathcal{g}(\mathcal{z},\overline{\mathcal{z}},P_4),$$

(3.50)

where

$$\mathcal{g}(\mathcal{z},\overline{\mathcal{z}},P_4)=〈\mathcal{G},\mathcal{F}(\mathcal{zE}+\overline{\mathcal{z}}\overline{\mathcal{E}},P_4)〉=\sum _{2\le j+k\le 3}\frac{{\mathcal{l}}_{\mathrm{jk}}}{j!k!}{\mathcal{z}}^j{\overline{\mathcal{z}}}^k,$$

(3.51)

with

$$\begin{array}{c}\multicolumn{1}{c}{\mathcal{F}(\mathcal{zE}+\overline{\mathcal{z}}\overline{\mathcal{E}},P_4)=}\\ \multicolumn{1}{c}{\left[\begin{array}{c}{\mathcal{f}}_0(\mathcal{z}\frac{(1-i)({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3-\mathfrak{Z})}{(-1+{\mathcal{u}}^{*2}){(1+{\mathcal{u}}^{*2})}^2{\mathcal{v}}^{*}}+\overline{\mathcal{z}}\frac{(1+i)({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3-\mathfrak{Z})}{(-1+{\mathcal{u}}^{*2}){(1+{\mathcal{u}}^{*2})}^2{\mathcal{v}}^{*}},2(\mathcal{z}+\overline{\mathcal{z}}))\hfill \\ \multicolumn{1}{c}{{\mathcal{g}}_0(\mathcal{z}\frac{(1-i)({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3-\mathfrak{Z})}{(-1+{\mathcal{u}}^{*2}){(1+{\mathcal{u}}^{*2})}^2{\mathcal{v}}^{*}}+\overline{\mathcal{z}}\frac{(1+i)({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3-\mathfrak{Z})}{(-1+{\mathcal{u}}^{*2}){(1+{\mathcal{u}}^{*2})}^2{\mathcal{v}}^{*}},2(\mathcal{z}+\overline{\mathcal{z}}))}\end{array}\right],}\end{array}$$

(3.52)

and

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{l}}_{20}=\frac{1}{{\mathcal{u}}^{*}{(-1+{\mathcal{u}}^{*2})}^2{(1+{\mathcal{u}}^{*2})}^7{(1+{\mathcal{u}}^{**2})}^3{\mathcal{v}}^{*2}}}\\ \multicolumn{1}{c}{i(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4)({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3-\mathfrak{Z}){\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3}\\ \multicolumn{1}{c}{+\mathfrak{iZ}((-1-i)(-1+{\mathcal{u}}^{*2}){(1+{\mathcal{u}}^{*2})}^2(-1+{\mathcal{u}}^{**4}){\mathcal{v}}^{*}}\\ \multicolumn{1}{c}{-{\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2})(-{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3+\mathfrak{Z}){\mathcal{v}}^{**}),}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{l}}_{11}=\frac{1}{{\mathcal{u}}^{*}{(-1+{\mathcal{u}}^{*2})}^3{(1+{\mathcal{u}}^{*2})}^9{(1+{\mathcal{u}}^{**2})}^4{\mathcal{v}}^{*3}}}\\ \multicolumn{1}{c}{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4)({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3-\mathfrak{Z})({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3}\\ \multicolumn{1}{c}{+\mathfrak{iZ})((1-i)(-1+{\mathcal{u}}^{*2}{)}^2(1+{\mathcal{u}}^{*2}{)}^4(-1+{\mathcal{u}}^{**2})}\\ \multicolumn{1}{c}{(1+{\mathcal{u}}^{**2}{)}^2{\mathcal{v}}^{*2}-(4-8i)(-1+{\mathcal{u}}^{*2})(1+{\mathcal{u}}^{*2}{)}^2}\\ \multicolumn{1}{c}{{\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2})(1+{\mathcal{u}}^{**2}){\mathcal{v}}^{*}({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3-\mathfrak{Z})}\\ \multicolumn{1}{c}{+i(-1+{\mathcal{u}}^{*2})(1+{\mathcal{u}}^{*2}{)}^2{\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2})}\\ \multicolumn{1}{c}{(1+{\mathcal{u}}^{**2}){\mathcal{v}}^{*}({\mathcal{u}}^{*}(1+{\mathcal{u}}^{*2}{)}^3-\mathfrak{Z}){\mathcal{v}}^{**}}\\ \multicolumn{1}{c}{-6i(1-6{\mathcal{u}}^{**2}+{\mathcal{u}}^{**4}){(-{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3+\mathfrak{Z})}^2{\mathcal{v}}^{**}),}\end{array}$$

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{l}}_{02}=\frac{1}{{\mathcal{u}}^{*}{(-1+{\mathcal{u}}^{*2})}^2{(1+{\mathcal{u}}^{*2})}^7{(1+{\mathcal{u}}^{**2})}^3{\mathcal{v}}^{*2}}}\\ \multicolumn{1}{c}{(\tilde{\mathcal{c}}+{\mathcal{c}}_4)(\tilde{\mathcal{h}}+{\mathcal{h}}_4)({\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3-\mathfrak{Z})(-i{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3}\\ \multicolumn{1}{c}{+\mathfrak{Z})((-1+i)(-1+{\mathcal{u}}^{*2}){(1+{\mathcal{u}}^{*2})}^2(-1+{\mathcal{u}}^{**4}){\mathcal{v}}^{*}}\\ \multicolumn{1}{c}{-{\mathcal{u}}^{**}(-3+{\mathcal{u}}^{**2})(-{\mathcal{u}}^{*}{(1+{\mathcal{u}}^{*2})}^3+\mathfrak{Z}){\mathcal{v}}^{**}).}\end{array}$$

Next, as established in Lemma 9.13 of Ref., ²⁴ there is a smooth transformation of variables, depending on the parameter exists and alters the equation (3.50) into the following normal form

$$\omega \to i\omega +\mathcal{C}\left({\alpha }_2\right)\omega {\left|\omega \right|}^2+\mathcal{D}\left({\alpha }_2\right){\overline{\omega }}^3+\left(O{\left|\omega \right|}^4\right),$$

(3.53)

where

$$\mathcal{C}\left(P_4\right)=\left(\frac{1+3i}{4}\right){\mathcal{l}}_{20}{\mathcal{l}}_{11}+\left(\frac{1-i}{2}\right){\left|{\mathcal{l}}_{11}\right|}^2-\left(\frac{1+i}{4}\right){\left|{\mathcal{l}}_{02}\right|}^2,$$

$$\mathcal{D}\left(P_4\right)=\left(\frac{i-1}{4}\right){\mathcal{l}}_{02}{\mathcal{l}}_{11}-\left(\frac{1+i}{4}\right){\mathcal{l}}_{02}{\overline{\mathcal{l}}}_{20}.$$

Next, we introduce the following quantities

$${\mathcal{C}}_1\left(P_4\right)=-4i\mathcal{C}\left(P_4\right),$$

$${\mathcal{D}}_1\left(P_4\right)=-4i\mathcal{D}\left(P_4\right),$$

under the condition ${\mathcal{D}}_1\left(P_4\right)\ne 0$ , the Jacobian matrix can be expressed as

$$\mathfrak{J}\left(P_4\right)=\frac{{\mathcal{C}}_1\left(P_4\right)}{\left|{\mathcal{D}}_1\left(P_4\right)\right|}.$$

From the above investigation, the subsequent theorem is stated as follows.

• Theorem 3.4. Assume the subsequent no-degeneracy assumptions hold:

${\mathcal{C}}_1\left(P_4\right)\ne 0,\mathrm{Re}\mathfrak{J}\left(P_4\right)\ne 0$ and $\mathrm{Im}\mathfrak{J}\left(P_4\right)\ne 0$ .

Under these conditions, the local behavior of model (1.4) near its interior fixed point can be characterized by the standard form map (3.53), particularly in the neighborhood of the critical values ${\mathcal{h}}_4$ and ${\mathcal{c}}_4$ . This map exhibits several notable dynamical phenomena.

• The trivial equilibrium point of structure (3.53) experiences a Neimark-Sacker bifurcation, giving rise to a stable invariant circle when the linearized system has eigenvalues $\lambda =-i$ . This circle collapses as the eigenvalues pass through $\lambda =i$ .

• Provided that the spectral radius of the Jacobian is larger than unity, that is, $\left|\mathfrak{J}\left(P_4\right)\right|>1$ , the map (3.53) has eight distinct non-trivial equilibria, that change by fold bifurcations with variations in the parameters.

• Besides, the map (3.53) has four non-trivial isolated equilibrium points which experience NSB, resulting in the formation of four small invariant loops. These loops then fade away by homoclinic bifurcations.

Numerical simulations

In this section, we carry out a detailed numerical simulation of system (1.4) to further verify and validate the analytical findings derived in the aforementioned Sections. Our analysis of bifurcation diagrams and maximum Lyapunov exponents of a system by varying the key system parameter is used to identify and describe qualitative transitions in the system dynamics systematically, for example, the development of periodic oscillations, quasi-periodicity and strong resonances.

Period-doubling bifurcation

We define the values for system’s variables as $\mathcal{a}=0.9,\mathcal{b}=0.6,\mathrm{and}\mathcal{h}=0.5$ , with the bifurcation parameter $\mathcal{c}\in [5.0,12.0]$ . Numerical simulations reveal that structure (1.4) undertakes a period-doubling bifurcation on $\mathcal{c}=6.44444$ . The corresponding bifurcation figures and the MLE are illustrated in Figure 2(a) to (c). For the specific case where $\mathcal{c}=3.97531$ and the fixed point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})=(0.3,0.338276),$ the system (1.4) yields the following characteristic expression:

$${\lambda }^2+0.221713\lambda -0.778287=0.$$

(4.1)

Figure 2.

(a and b) demonstrate that the period-doubling bifurcation in model (1.4) for $\mathcal{a}=0.9,\mathcal{b}=0.6,\mathcal{h}=0.5,$ and $\mathcal{c}\in [5.0,12.0]$ and (c) represents the corresponding maximum Lyapunov exponents.

The eigenvalues are ${\lambda }_{1,2}=-1\mathrm{and}0.778287$ , satisfying the required criteria for period-doubling bifurcation near the interior fixed point. Moreover, the conditions stated in Theorem 3.1 are verified: ${\mathcal{P}}_1-0.26175\ne 0$ and ${\mathcal{P}}_2=4.50991\ne 0,$ with ${\mathcal{P}}_2>0.$ This establishes the emergence of a supercritical bifurcation with 2-cycle, whereby the emerging period-two orbit is locally attracting in a neighborhood of the bifurcating point.

Neimark-Sacker bifurcation

Choosing parameters $\mathcal{a}=14,\mathcal{b}=12,\mathcal{h}=0.1$ and allowing the bifurcation parameter $\mathcal{c}\in [6.6,7.4]$ , the system (1.4) exhibits a NSB at $\mathcal{c}=7.2222$ . Figure 3(a) and (b), present the bifurcation diagrams of ${\mathcal{u}}_{\mathcal{n}}\mathrm{and}{\mathcal{v}}_{\mathcal{n}}$ respectively, while the associated MLE is shown in Figure 3(c). On the equilibrium point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})=(2,4.15385),$ the characteristic equation is represented by

$${\lambda }^2-1.97111\lambda +1=0.$$

(4.2)

Figure 3.

(a and b) depict the Neimark-Sacker bifurcation within structure (1.4) for $\mathcal{a}=14,\mathcal{b}=12,\mathcal{h}=0.1$ and $\mathcal{c}\in [6.6,7.4],$ and (c) represents the corresponding maximum Lyapunov exponents.

The eigenvalues are ${\lambda }_{1,2}=0.985556\pm 0.169352i$ with $\left|{\lambda }_{1,2}\right|=1$ , confirming the NSB condition. Furthermore, the non-degeneracy condition $\mathcal{L}=-0.00179829\ne 0$ is satisfied, as required by Theorem 3.2. Since $\mathcal{L}<0$ , this confirms the occurrence of a supercritical NSB, whereby the system undergoes a bifurcation where a stable periodic orbit originates from the equilibrium point as the parameter $\mathcal{c}$ decreases below the critical value $\mathcal{c}=7.2222$ .

1:2 Strong resonance

To study the $1:2$ strong resonance in model (1.4), we set $\mathcal{a}=31.5,\mathfrak{b}=30.2,\mathcal{h}\in [0.57,0.59]$ and $\mathcal{c}\in [24.5015,24.5019]$ . At $\mathcal{c}=24.50168$ and $\mathcal{h}=0.58$ , the model (1.4) experiences a 1:2 resonance near the equilibrium point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})=(1.3,2.55047)$ . The bifurcation and MLE are depicted in Figure 4(a) to (c). The characteristic polynomial is

$${\lambda }^2+2\lambda +1=0,$$

(4.3)

Figure 4.

(a and b) depict 1:2 strong resonance in structure (1.4) for $\mathcal{a}=31.5,\mathfrak{b}=30.2,\mathcal{h}=0.58$ and $\mathcal{c}\in [24.5015,24.5019,]$ and (c) represents the corresponding maximum Lyapunov exponents.

yielding a double root at ${\lambda }_{1,2}=-1$ . Further, the non-degeneracy conditions $\det (\mathcal{C}({\tilde{\mathcal{h}},\tilde{\mathcal{c}})|}_{(0,0)})=-0.666689\ne 0$ , $\mathcal{M}(0,0)=2.23066\ne 0$ and $\mathcal{N}(0,0)+3\mathcal{M}(0,0)=12.8974\ne 0$ are satisfied, confirming the 1:2 resonance.

1:3 Strong resonance

For 1:3 strong resonance in model (1.4), we select $\mathcal{a}=22,\mathfrak{b}=20,26.6331\le \mathcal{c}\le 26.6336$ and $0.5302\le \mathcal{h}\le 0.5308$ . At $\mathcal{c}=26.6333$ and $\mathcal{h}=0.5307,$ the framework (1.4) exhibits a $1:3$ strong resonance around equilibrium point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})=(2,1.87735)$ . The bifurcation diagrams are displayed in Figure 5(a) to (c).

Figure 5.

(a and b) show the 1:3 strong resonance in system (1.4) for select $\mathcal{a}=22,\mathfrak{b}=20,\mathcal{h}=0.5307$ and $\mathcal{c}\in [26.6331,26.6336]$ and (c) represents the corresponding maximum Lyapunov exponents.

On this point the computed Jacobian takes the form

$$\mathcal{J}(2,1.87735)=\left(\begin{array}{cc}-0.06019\hfill & -5.65331\hfill \\ \multicolumn{1}{c}{3.18397}& -4.65331\hfill \end{array}\right),$$

(4.4)

with eigenvalues ${\lambda }_{1,2}=-\frac{1}{2}\pm \frac{\sqrt{3}}{2}i$ , fulfilling the 1:3 resonance condition.

1:4 Strong resonance

Finally, we investigate 1:4 strong resonance in model (1.4), with $\mathcal{a}=18,\mathfrak{b}=16,18.3\le \mathcal{c}\le 19.5$ and $0.5081\le \mathcal{h}\le 0.5089$ . At $\mathcal{c}=19.331$ and $\mathcal{h}=0.5086$ , the model (1.4) experiences a $1:4$ strong resonance near the non-negative equilibrium point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})=(2,2.06916)$ . The bifurcation diagrams are presented in Figure 6(a) and (b), and MLE are shown in Figure 6(c).

Figure 6.

(a and b) Show the 1:4 strong resonance in framework (1.4) when parameters are defined as $\mathcal{a}=18,\mathfrak{b}=16,\mathcal{h}=0.5089$ and $\mathcal{c}\in [18.3,19.5]$ , and (c) represents the corresponding MLE.

Evaluating the variation matrix at the equilibrium point yields

$$\mathcal{J}(2,2.06916)=\left(\begin{array}{cc}2.9326\hfill & -3.9326\hfill \\ \multicolumn{1}{c}{2.4411}& -2.9326\hfill \end{array}\right),$$

(4.5)

with purely complex eigenvalues ${\lambda }_{1,2}=\pm i$ , confirming the 1:4 resonance condition. Moreover, requirements for Theorem 3.4 are fulfilled at the critical point, as the relevant quantities are computed as follows:

${\mathcal{C}}_1\left(P_4\right)=-58.865-95.863i\ne 0,\mathrm{Re}\mathfrak{J}\left(P_4\right)=-7.40147\ne 0$ and $\mathrm{Im}\mathfrak{J}\left(P_4\right)=-12.0534\ne 0.$

Moreover, in Figure 7, phase portraits of system (1.4) are plotted for different values of bifurcation parameter $\mathcal{c}\in [6.6,7.4]$ while other parameters $(\mathcal{a},\mathcal{b},\mathcal{h})=(14,12,0.1)$ and initial conditions $(2,4.15385)$ are remaining same in each case.

Figure 7.

Phase portraits of system (1.4) for different values of bifurcation parameter $c$ . Phase portrait for $\mathcal{c}=7.0$ . Phase portrait for $\mathcal{c}=7.2222$ . Phase portrait for $\mathcal{c}=7.225$ . Phase portrait for $\mathcal{c}=7.23$ . Phase portrait for $\mathcal{c}=7.3$ . Phase portrait for $\mathcal{c}=7.4$ .

Chaos control

Optimization of dynamical structures to achieve specific performance goals with the lessening of chaotic behavior is an important objective in diverse scientific and engineering disciplines. Among the many methods that have been developed for this purpose, chaos control techniques have found great popularity in both the theoretical and applied fields of research. This section is devoted to one such method that has been very influential, the OGY control technique originally proposed by Ott, Grebogi and Yorke. To implement this control technique, system (1.4) is reformulated in the following form:

$$\begin{array}{c}{\mathcal{u}}_{n+1}={\mathcal{u}}_n+\mathcal{h}(\mathcal{a}-{\mathcal{u}}_n-\frac{\mathcal{c}{\mathcal{u}}_n{\mathcal{v}}_n}{1+{\mathcal{u}}_n^2})=\mathcal{f}({\mathcal{u}}_n,{\mathcal{v}}_n,\mathcal{c}),\hfill \\ \multicolumn{1}{c}{{\mathcal{v}}_{n+1}={\mathcal{v}}_n+\mathcal{h}(\mathcal{b}-\frac{\mathcal{c}{\mathcal{u}}_n{\mathcal{v}}_n}{1+{\mathcal{u}}_n^2})=\mathcal{g}({\mathcal{u}}_n,{\mathcal{v}}_n,\mathcal{c}).}\end{array}$$

(5.1)

In the mentioned framework, the parameter $c$ is the control variable, that is, it allows the regulation of the behavior of the system by means of minimal perturbations. To ensure adequate control, $\mathcal{c}$ is confined within a narrow interval $\mathcal{c}\in ({\mathcal{c}}_0-\zeta ,{\mathcal{c}}_0+\zeta ),$ where $\mathcal{c}>0$ and ${\mathcal{c}}_0$ represents a minor value situated within a chaotic region. A steadying feedback switch mechanism is then applied to direct the model’s routes toward the anticipated periodic trajectory. Consider the unstable fixed point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ of model (1.4), which emerges within the chaotic regime as a result of both Neimark-Sacker and period-doubling bifurcations. Near this unstable equilibrium point, the dynamics of model (1.5) can be approximated by a linear map:

$$\left[\begin{array}{c}{\mathcal{u}}_{n+1}-{\mathcal{u}}^{*}\hfill \\ \multicolumn{1}{c}{{\mathcal{v}}_{n+1}-{\mathcal{v}}^{*}}\end{array}\right]\approx \mathfrak{J}\left[\begin{array}{c}{\mathcal{u}}_n-{\mathcal{u}}^{*}\hfill \\ \multicolumn{1}{c}{{\mathcal{v}}_n-{\mathcal{v}}^{*}}\end{array}\right]+\mathcal{B}[\mathcal{c}-{\mathcal{c}}_0],$$

(5.2)

where

$$\mathfrak{J}=\left[\begin{array}{cc}\frac{\partial \mathcal{f}({\mathcal{u}}^{*},{\mathcal{v}}^{*},{\mathcal{c}}_0)}{\partial \mathcal{u}}\hfill & \frac{\partial \mathcal{f}({\mathcal{u}}^{*},{\mathcal{v}}^{*},{\mathcal{c}}_0)}{\partial \mathcal{v}}\hfill \\ \multicolumn{1}{c}{\frac{\partial \mathcal{g}({\mathcal{u}}^{*},{\mathcal{v}}^{*},{\mathcal{c}}_0)}{\partial \mathcal{u}}}& \frac{\partial \mathcal{g}({\mathcal{u}}^{*},{\mathcal{v}}^{*},{\mathcal{c}}_0)}{\partial \mathcal{v}}\hfill \end{array}\right],$$

(5.3)

$$\mathcal{B}=\left[\begin{array}{c}\frac{\partial \mathcal{f}({\mathcal{u}}^{*},{\mathcal{v}}^{*},{\mathcal{c}}_0)}{\partial \mathcal{c}}\hfill \\ \multicolumn{1}{c}{\frac{\partial \mathcal{g}({\mathcal{u}}^{*},{\mathcal{v}}^{*},{\mathcal{c}}_0)}{\partial \mathcal{c}}}\end{array}\right].$$

(5.4)

Framework (5.1) is said to be controllable if $[\mathcal{B}:\mathfrak{J}\mathcal{B}]$ be a matrix with rank 2. Next, by setting $[\mathcal{c}-{\mathcal{c}}_0]=-\mathfrak{R}\left[\begin{array}{c}{\mathcal{u}}_n-{\mathcal{u}}^{*}\hfill \\ \multicolumn{1}{c}{{\mathcal{v}}_n-{\mathcal{v}}^{*}}\end{array}\right]$ , where $\mathfrak{R}=\left[\begin{array}{cc}{s}_1& s_2\end{array}\right]$ , system (5.2) can be expressed as:

$$\left[\begin{array}{c}{\mathcal{u}}_{n+1}-{\mathcal{u}}^{*}\hfill \\ \multicolumn{1}{c}{{\mathcal{v}}_{n+1}-{\mathcal{v}}^{*}}\end{array}\right]\approx [\mathfrak{J}-\mathcal{B}\mathfrak{R}]\left[\begin{array}{c}{\mathcal{u}}_n-{\mathcal{u}}^{*}\hfill \\ \multicolumn{1}{c}{{\mathcal{v}}_n-{\mathcal{v}}^{*}}\end{array}\right].$$

(5.5)

The corresponding controlled version of (5.1) is given by:

$$\begin{array}{c}\multicolumn{1}{c}{{\mathcal{u}}_{n+1}={\mathcal{u}}_n+\mathcal{h}(\mathcal{a}-{\mathcal{u}}_n-\frac{(\mathcal{c}-s_1({\mathcal{u}}_n-{\mathcal{u}}^{*})-s_2({\mathcal{v}}_n-v^{*})){\mathcal{u}}_n{\mathcal{v}}_n}{1+{\mathcal{u}}_n^2}),}\\ \multicolumn{1}{c}{{\mathcal{v}}_{n+1}={\mathcal{v}}_n+\mathcal{h}(\mathcal{b}-\frac{(\mathcal{c}-s_1({\mathcal{u}}_n-{\mathcal{u}}^{*})-s_2({\mathcal{v}}_n-v^{*})){\mathcal{u}}_n{\mathcal{v}}_n}{1+{\mathcal{u}}_n^2}).}\end{array}$$

(5.6)

Fixed point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ is locally asymptotically stable if and only if ${\lambda }_1$ and ${\lambda }_2$ for the matrix $\mathfrak{J}-\mathcal{B}\mathfrak{R}$ exists inside the unit disk. Characteristic equation corresponding to $\mathfrak{J}-\mathcal{B}\mathfrak{R}$ is given as

$$K\left(\omega \right)={\omega }^2-\mathrm{T}/\omega +\mathfrak{D}.$$

(5.7)

Here, $\mathrm{T}/$ and $\mathfrak{D}$ denotes the determinant of the matrix $\mathfrak{J}-\mathcal{B}\mathfrak{R}$ . Boundaries for the marginal stability region are determined with evaluating the characteristic polynomial at specific points:

$$l_1=\mathcal{F}\left(0\right)=1,l_2=\mathcal{F}\left(1\right)=0\mathrm{and}{l}_3=\mathcal{F}(-1)=0.$$

• Theorem 5.1. A non-negative equilibrium point $({\mathcal{u}}^{*},{\mathcal{v}}^{*})$ of the controlled structure (5.6) is locally asymptotically stable if the eigenvalues ${\lambda }_1{\mathrm{and}\mathrm{\lambda }}_2$ are positioned inside the triangular area in the $s_1{s}_2$ -plane enclosed with the lines $l_1,l_2$ and $l_3$ for given parameter values (Figure 8).

Figure 8.

(a) depicts the controllable region associated with period-doubling bifurcation for parameters $\mathcal{a}=0.9,\mathfrak{b}=0.6,\mathcal{h}=0.5$ and $\mathcal{c}=6.4444$ with $-20\le s_1\le 15$ and $-20\le s_2\le 15$ and (b) represents the controllable region corresponding to the Neimark-Sacker bifurcation, with $\mathcal{a}=14,\mathfrak{b}=12,\mathcal{h}=0.1$ and $\mathcal{c}=7.2222$ under the same bounds: $-20\le s_1\le 15$ and $-20\le s_2\le 15$ .

Results and discussions

In this study, we have conducted a detailed bifurcation and stability investigation of a discretized model derived from the chemically inspired Degn-Harrison reaction model. By employing the forward Euler discretization scheme, the original continuous system was transformed into a planar nonlinear map, enabling the investigation of its dynamic properties in discrete time. This analysis demonstrates that even a fairly simple chemical reaction mechanism can have colorful and complex dynamics when discretized and that the mechanism can have many paths to complexity and chaos. Through the combination of strong theoretical tools and many numerical simulations, we were able to locate the parameter regions of the system where codimension one bifurcations, namely Neimark-Sacker and period-doubling, and codimension two strong resonances of type 1:2, 1:3, and 1:4 are encountered close to the interior equilibrium point. These results emphasize the importance of the discretization in the determination of the observational nature of nonlinear structure and exhibit the presence of phenomena that may not be evident for their continuous counterparts.

Furthermore, the forward discretization scheme of the classical Degn-Harrison model, based on equation (2.1) of the aforementioned work, when applied to this model, yields a discrete-time dynamical system with a much more abundant and complex set of behaviors than the continuous version. Though the system under analysis has been reported to exhibit Hopf bifurcation, resulting in the existence of stable limit cycles, in our discrete analysis system, we observe the appearance of Neimark-Sacker and period-doubling bifurcations, highly non-trivial strong resonances (co-dimension two), and routes to deterministic chaos, which are successfully stabilized by contemporary chaos control techniques. This extension of regular to chaotic dynamics and their reappearances as order has far-reaching implications on the possible behavior of the biochemical regulatory network of bio products in the face of discrete sampling, as well as innate discrete-time dynamics. The stability analysis of the isolated positive steady state, which resembles a homeostatic metabolic state in which oxygen and nutrient uptake are balanced, is based on this. This equilibrium becomes unstable as the main kinetic parameters are varied by changing the reaction rate $c$ . In the discrete model, this destabilization is realized through two main mechanisms: a pair of complex conjugate eigenvalues crossing over the unit circle, which creates a Neimark-Sacker bifurcation, and a real eigenvalue that goes through $-1$ , which creates a period-doubling bifurcation. The Neimark-Sacker bifurcation is biologically defined as the origin of quasi-periodic oscillations of the chemical concentrations. This is a multi-frequency respiration rate rhythmicity in bacteria that may arise due either to the nonlinearity of the interaction between the substrate inhibition kinetics and the discrete-time character of some cellular processes, for example, the feedback between gene expression and metabolic periodic sampling. To these, the period-doubling bifurcation is the beginning of a cascade of period-doubling oscillations, culminating in long-term periodic chaotic behavior. This movement is a pathological loss of normal metabolic oscillation into an irregular, erratic respiratory condition. This kind of chaos will mimic a malfunctioning regulatory circuit in which the feedback is inhibitory and is massively out of control, resulting in highly irregular rates of oxygen uptake that cannot maintain stable operation, similar to some of the arrhythmic or fibrillation states of biological systems. These bifurcation curves in the space of parameters form degenerate codimension two points, namely strong resonances of order 1:2, 1:3, and 1:4. It is at these critical points where the nonlinear interaction of the emerging frequencies results in the so-called mode-locking. In this case, the intrinsic oscillatory vulnerability of the system will be entrained to a subharmonic of the system, which will yield stable, tightly periodic orbits of low period. These resonances are dynamical regimes in a biochemical sense, which are highly structured and stable. They may indicate those states in which an internal metabolic oscillator gets strongly coupled to an external periodic driver, like a circadian rhythm or periodic nutrient flow, or in which two or more internal oscillatory modules become strongly coupled, with one in a perfect phase relationship. Such entrainment defines tightly concerned parameters to have these resonance tongues, which can be evolutionarily chosen to provide good timing and communication in cellular activities.

The chaotic attractors in the parameter space of the discrete model demonstrate a very important discovery: the simple, well-understood Degn-Harrison kinetics can exhibit deterministic chaos when placed in a discrete-time context. This commotion is a condition of extreme dynamical disease to the hypothetical bacterial culture, which is (apparently) characterized by random and unpredictable changes in both oxygen and nutrient levels that would be extremely inefficient and could be fatal. The fact that chaos control scheme, namely, the Ott-Grebogi-Yorke control scheme, have been successfully implemented, proves the fact that this complex, undesired state is not final. Through the application of minute, time sensitive perturbation to an accessible parameter of the system, given its present state, the chaotic trace can be brought to a suitable periodic course or even manipulated back to the steady state. This has a strong biological connotation. The theory indicates that a pathological chaotic metabolic network could be treated not by the enormous, long-term intervention, but by fine adjustments, which occur at the right time and location, akin to a dynamical pacemaker, which would restore rhythmicity by providing small but essential feedback.

There is strong dynamical evidence on the basis of numerical simulations which reveal this whole dynamical story. A graph of the landscape is bifurcation diagrams, which illustrate a smooth passage between stability to quasi-periodic through the Neimark-Sacker bifurcation and the period-doubling cascades which narrow into wide, chaotic bands. Phase portraits make a graphic visualization of the closed curves of invariant quasiperiodic motion, the highly detailed geometry of chaotic attractors rich in fractals, against a backdrop of the simple point equilibrium. The control laws are most effectively visualized by the action of projecting a trajectory that starts chaotically wandering in the chaotic attractor into a domesticated and bounded periodic cycle. Finally, our discrete-time study of the Degn-Harrison model continues its legacy that the nonlinearity of the degree of the model can, through its interactions with discrete-time evolution, coordinate all the complex behaviors of rhythmicity through to chaos. The co-dimension two phenomena analysis gives a fine-grained map of the complex resonant state of the system and the successful chaos control provides an optimistic principled way to theoretically restore the functionality of life, in terms of concentrations but in the language of nonlinear dynamical systems and their control.

Future directions

Building upon this work, several promising research avenues can be explored:

• Alternative Discretization Methods: Examining the effect of different implicit schemes on the bifurcation structure and stability characteristics.

• Inclusion of diffusion and spatial effects: Extending the model to reaction-diffusion systems to analyze pattern formation and spatiotemporal chaos in discrete environments.

• Stochastic Perturbations: Incorporating noise or random fluctuations to analyze the robustness of bifurcation scenarios under uncertainty.

• Control and Stabilization: Developing control strategies to suppress undesirable oscillations or chaos in discrete chemical systems.

• Applications to Real Chemical Networks: Using the theoretical framework on experimentally validated reaction mechanisms to connect theory with practice.

These directions aim to improve the theoretical understanding of discrete-time chemical models and broaden their applicability in real-world chemical and biological systems.