Synchronized Cycles of Generalized Nicholson-Bailey Model

Tahmineh Azizi¹, Gabriel Kerr²

¹Institute of Computational Comparative Medicine, Department of Anatomy and Physiology, Department of Mathematics, Kansas State University, Manhattan, Kansas.
²Department of Mathematics, Kansas State University, Manhattan, Kansas.

DOI: 10.4236/ajcm.2020.101009                

Abstract

In this paper, we study a drive-response discrete-time dynamical system which has been coupled using convex functions and we introduce a synchronization threshold which is crucial for the synchronizing procedure. We provide one application of this type of coupling in synchronized cycles of a generalized Nicholson-Bailey model. This model demonstrates a rich cascade of complex dynamics from stable fixed point to periodic orbits, quasi periodic orbits and chaos. We explain how this way of coupling makes these two chaotic systems starting from very different initial conditions, quickly get synchronized. We investigate the qualitative behavior of GNB model and its synchronized model using time series analysis and its long time dynamics by the help of bifurcation diagram.

Keywords

Generalized Nicholson Bailey Model, Synchronized Cycles, Synchronization Threshold, Complete Synchronization

1. Introduction

Traditionally, the study of long time behavior of a dynamical systems was based on the examples of ordinary differential equations with regular solutions and those solutions which remained in a bounded region of the phase space could be divided into two different types based on their local behavior: first, a stable equilibrium point and second, a periodic (or quasi-periodic) oscillation. Edward Lorenz in 1961, by working on a simplified version of atmospheric transfer model which was consisting of three nonlinear ordinary differential equations, numerically observed that a very small changing in the initial conditions of the system equations makes a huge difference on the long term behavior of their solutions [1]. Indeed, his finding was due to one of the major properties of chaotic dynamical systems which later called as sensitive dependence on initial conditions or butterfly effect.

Chaos is a complex nonlinear phenomenon that has been increasingly studied in the last three decades. During those years, many fields of science and engineering have been affected by chaos studies. One of the most important achievements in nonlinear and complex dynamics is the discovery of synchronized chaos. Synchronization happens when two events take place in synchrony at the same time and when time approaches infinity, the error between solutions of the first system and its synchronized one vanishes and approaches to zero. The synchronization between two dynamical systems is a well-known phenomena occurring in Physics, Biology or Engineering and refers to a phenomenon that may occur when two or more oscillators are coupled. For the first time, Christiaan Huygens discovered synchronization of coupled pendulum clocks in 1665 [2]. Occurrence of synchronization in coupled chaotic system composed of identical chaotic oscillators has been detected for the first time by Fujisaka and Yamada [3] and [4] and after that it has been reported by Pecora and Carroll [5].

The dynamics of coupled chaotic systems show properties which we cannot detect in the behavior of the individual elements [6]. Someone can find the same spatial synchronized fluctuations in biology, ecology and epidemiology [7] [8] [9] [10] [11]. Synchronization of complex population oscillations in natural systems has been examined widely by some researchers [12]. Bernd Blasius and Lewi Stone worked on a chaotic UPCA foodweb model and they claimed that the spatio-temporal structures associated with phase synchronization have important implications for conservation ecology. They proposed that even though perturbation of a local patch population can bring them to the brink of extinction, the periodicity of spatial phase synchronization can help to buffer the endangered population by colonizers. They also asserted that unlike this thought that population synchronization can cause global population extinction [13], however, phase synchronization can be useful for maintaining species persistence. Their finding indicated that synchronization can shape the distribution and abundance of species even in continental scale.

As we have already discussed, there are many varieties of synchronization. In this research, instead of exploring all of these different types of synchronization which have been proposed for different purposes and with various applications, we simply focus on the most fundamental case, and we will develop our new approach based on a basic mathematical concept. Indeed, the purpose of this paper is that after defining and setting the fundamental concepts which we need to establish the basis of our study, to demonstrate such configurations under a suitable coupling method is possible. Moreover, we explore the dynamical and ecological effects of synchronization of a host-parasitoid model which is a generalization of Nicholson-Bailey model (GNB). The GNB model demonstrates regular and irregular or chaotic oscillations. We define a lift function which is technically a convex function and maps the orbits of the drive system into the orbits of the response system. Using this convex function, we drive the response system which inherits all the complex qualitative dynamics of GNB model and mimics that certain properties of the motion which is shared between them. We investigate numerically that this method of coupling, synchronizes completely the stable and periodic cycles and even the chaotic motions of GNB model and to do that, we need to adjust the synchronization constant to be closer to zero. We demonstrate the complex dynamics of GNB model and its coupled system by conducting some time series and bifurcation analysis.

2. Drive-Response System Derivation

In this section, we derive the coupled system corresponding to the drive system by defining a convex link function. We consider the following drive system:

$X_{n+1}=g\left(X_n\right)\mathrm{,}$ (2.1)

where $X_n\in {ℝ}^n$ is the state vector of a general discrete-time drive system and $g\mathrm{<mo>\; :\; </mo>}{ℝ}^n\to {ℝ}^n$ is continuous. To find an appropriate response system, we provide the following definition [14]:

Definition 2.1. Assume $X\mathrm{,}Y\in {ℝ}^n$ are the state vectors of two non-linear discrete-time dynamical systems and a constant $0<s\le 1$. Then, a continuous function $h\left(X\mathrm{,}Y\mathrm{,}s\right)\mathrm{<mo>\; :\; </mo>}=\left(1-s\right)X+sY$ where $h\mathrm{<mo>\; :\; </mo>}{ℝ}^{2n+1}\to {ℝ}^n$ is called a link function which maps the orbits of first system keeping the same qualitative dynamics to the orbits of second system.

Using the definition (2.1), we develop a new system which inherits the qualitative features of the system (2.1) and has the following form

$Y_{n+1}=g\left(\left(1-s\right)X_n+s{Y}_n\right)$ (2.2)

and it is called response system. Now, for the following non-linear discrete-time dynamical system, we are going to develop a theorem which helps us to find the synchronization threshold. We consider the following drive-response system:

$\{\begin{array}{l}{X}_{n+1}=g\left(X_n\right),\\ Y_{n+1}=g\left(\left(1-s\right)X_n+s{Y}_n\right).\end{array}$ (2.3)

where $X_n\mathrm{,}{Y}_n\in {ℝ}^n$ are the state vectors of drive system (2.1) and response system (2.2) respectively, g is a mapping from ${ℝ}^n$ to itself and a constant $0<s\le 1$. The Jacobian matrix for drive-response system (2.3) has the following form:

$J\mathrm{<mo>\; :\; </mo>}=\left[\begin{array}{cc}\frac{\partial g\left(X\right)}{\partial X}& 0\\ \frac{\partial g\left(h\right)}{\partial X}& \frac{\partial g\left(h\right)}{\partial Y}\end{array}\right]$ (2.4)

Definition 2.2. We say that the drive system (2.1) and response system (2.2) are in complete synchronization if

$\underset{n\to \infty }{lim}\Vert Y_n-X_n\Vert =0$ (2.5)

Here, we imply to a general definition in synchronization theory which is crucial for proposed coupling method.

Definition 2.3. Let E be a Banach space. Then, the map $g\mathrm{<mo>\; :\; </mo>}E\to E$ is called a contraction mapping if there exists a constant $0\le \alpha <1$ such that for every pair of points $X\mathrm{,}Y\in E$, we have $\Vert g\left(X\right)-g\left(Y\right)\Vert \le \alpha \Vert X-Y\Vert$, where $\alpha$ is called a contraction constant of g on E [15].

Now, we have the following result for drive-response system (2.3).

Proposition 2.4. For drive-response system (2.3) if g is a contraction function, then the solutions eventually evolve identically in time and

$\underset{n\to \infty }{lim}\Vert e\left(n+1\right)\Vert =0$

where $e\left(n+1\right)=\Vert Y_{n+1}-X_{n+1}\Vert$ is the error between the solutions of the system (2.3).

Proof. For drive-response system (2.3), we have:

$e\left(n+1\right)=Y_{n+1}-X_{n+1}=g\left(\left(1-s\right)X_n+s{Y}_n\right)-g(Xn)$

We can easily see that for $0<s\le 1$ :

$\Vert \left(\left(1-s\right)X_n+s{Y}_n\right)-X_n\Vert \le s\Vert Y_n-X_n\Vert$

for $0\le \alpha <1$ and since g is a contraction mapping, we can write:

$\Vert e\left(n+1\right)\Vert =\Vert g\left(\left(1-s\right)X_n+s{Y}_n\right)-g\left(X_n\right)\Vert \le \alpha \Vert Y_n-X_n\Vert =\alpha \Vert e\left(n\right)\Vert$

That is to say

$\Vert e\left(n+1\right)\Vert \le \alpha \Vert e\left(n\right)\Vert$

It is obvious that

$\underset{n\to \infty }{\lim }\Vert e\left(n\right)\Vert =\underset{n\to \infty }{\lim }\Vert Y_n-X_n\Vert =0$

Therefore;

$\underset{n\to \infty }{lim}\Vert e\left(n+1\right)\Vert =0$

To obtain our rigorous results for complete synchronization, we need to find the normal form for drive-response system (2.3), we need to perform a few linear coordinate transformations that will put (2.3) into a form which is easier to work with [16]. First we transform the fixed point $\left(X^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{Y}^{\mathrm{<mo>\; *\; </mo>}}\right)$ of the system (2.3) to the origin by the translation $X=X^{\mathrm{<mo>\; *\; </mo>}}+\bar{X}$ and $Y=Y^{\mathrm{<mo>\; *\; </mo>}}+\bar{Y}$ under which drive-response system (2.3) becomes

$F:=\{\begin{array}{l}{\bar{X}}_{n+1}=g\left({\bar{X}}_n+X^{*}\right)-X^{*}\equiv G\left({\bar{X}}_n\right),\\ {\bar{Y}}_{n+1}=g\left(\left(1-s\right)\left({\bar{X}}_n+X^{*}\right)+s\left({\bar{Y}}_n+Y^{*}\right)\right)-Y^{*}\equiv G\left({\bar{X}}_n,{\bar{Y}}_n\right).\end{array}$ (2.6)

The Jacobian matrix for drive-response system (2.6) has the following form:

$\tilde{J}\mathrm{<mo>\; :\; </mo>}=\left[\begin{array}{cc}\frac{\partial G\left(\bar{X}\right)}{\partial \bar{X}}& 0\\ \frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{X}}& \frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{Y}}\end{array}\right]$ (2.7)

Then, we split off the linear part of the system (2.6) and write

$\begin{array}{c}\left[\begin{array}{c}\bar{X}\\ \bar{Y}\end{array}\right]=\left[\begin{array}{cc}\frac{\partial G\left(\bar{X}\right)}{\partial \bar{X}}& 0\\ \frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{X}}& \frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{Y}}\end{array}\right]\left[\begin{array}{c}\bar{X}\\ \bar{Y}\end{array}\right]+\left[\begin{array}{c}\tilde{G}\left(\bar{X}\right)\\ \tilde{G}\left(\bar{X}\mathrm{,}\bar{Y}\right)\end{array}\right]\\ =\left[\begin{array}{c}\frac{\partial G\left(\bar{X}\right)}{\partial \bar{X}}\bar{X}+\tilde{G}\left(\bar{X}\right)\\ \frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{X}}\bar{X}+\frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{Y}}\bar{Y}+\tilde{G}\left(\bar{X}\mathrm{,}\bar{Y}\right)\end{array}\right]\end{array}$ (2.8)

where

$\left[\begin{array}{c}\tilde{G}\left(\bar{X}\right)\\ \tilde{G}\left(\bar{X},\bar{Y}\right)\end{array}\right]=\left[\begin{array}{c}G\left(\bar{X}\right)\\ G\left(\bar{X},\bar{Y}\right)\end{array}\right]-\left[\begin{array}{cc}\frac{\partial G\left(\bar{X}\right)}{\partial \bar{X}}& 0\\ \frac{\partial G\left(\bar{X},\bar{Y}\right)}{\partial \bar{X}}& \frac{\partial G\left(\bar{X},\bar{Y}\right)}{\partial \bar{Y}}\end{array}\right]\left[\begin{array}{c}\bar{X}\\ \bar{Y}\end{array}\right]$ (2.9)

Let Q be the matrix that transforms the matrix $\tilde{J}$ into (real) Jordan canonical form which has the following form

$Q=\left[\begin{array}{cc}\frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{Y}}-\frac{\partial G\left(\bar{X}\right)}{\partial \bar{X}}& 0\\ \frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{X}}& 1\end{array}\right],Q^{-1}=\left[\begin{array}{cc}\frac{1}{\frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{Y}}-\frac{\partial G\left(\bar{X}\right)}{\partial \bar{X}}}& 0\\ \frac{\frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{X}}}{\frac{\partial G\left(\bar{X}\right)}{\partial \bar{X}}-\frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{Y}}}& 1\end{array}\right]$

Then, under the transformation

$\left[\begin{array}{c}\bar{X}\\ \bar{Y}\end{array}\right]=\left[\begin{array}{cc}\frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{Y}}-\frac{\partial G\left(\bar{X}\right)}{\partial \bar{X}}& 0\\ \frac{\partial G\left(\bar{X}\mathrm{,}\bar{Y}\right)}{\partial \bar{X}}& 1\end{array}\right]\left[\begin{array}{c}U\\ V\end{array}\right]$ (2.10)

(2.8) becomes

$\left[\begin{array}{c}U\\ V\end{array}\right]=\hat{J}\left[\begin{array}{c}U\\ V\end{array}\right]+Q^{-1}\left[\begin{array}{c}\tilde{G}\left(U\right)\\ \tilde{G}\left(U\mathrm{,}V\right)\end{array}\right]$ (2.11)

where $\hat{J}=Q^{-1}\tilde{J}Q$. We remark that the transformation (2.10) has simplified the linear part of (2.8) as much as possible. One can continue the task of simplifying the nonlinear part. However, for our purpose, we only need to focus on the linear part of the system. The schematic representation of the procedure of complete synchronization in a general discrete-time drive-response dynamical system has been demonstrated in Figure 1.

Hartman (1960) and Grobman (1959) proved that the orbit structure near a hyperbolic fixed point has the same qualitative structure of associated linearized system [17] [18] [19]. According to Hartman-Grobman theorem, the dynamical

Figure 1. The schematic diagram for complete synchronization in a discrete-time drive-response dynamical system.

systems behave similar to their linearization part around the fixed point. However, this theorem needs the linearization part without eigenvalues with real part zero for continuous-time system and for discrete-time systems, it needs the absolute value of the eigenvalues of linearized part not become one.

Remark 1. Because of the nature of the contraction mapping theorem and therefore the proposed coupling method, the Hartman-Grobman theorem can be applied directly into our problem.

Lemma 2.5. Using Hartman-Grobman theorem, to find the condition under which the drive-response system (2.11) achieves complete synchronization, we only need to look at the linear part of (2.11) which is

$L\left(U\mathrm{,}V\right)\equiv Q^{-1}\tilde{J}Q$ (2.12)

Theorem 2.6. Given Jacobian matrix $\hat{J}=Q^{-1}\tilde{J}Q$, for which the following inequality holds:

$\Vert \hat{J}\Vert :={\rho }_{\hat{J}}=\max \left|{\lambda }_i\right|\le \alpha <1$ (2.13)

where $\alpha$ is the contraction constant, ${\rho }_{\hat{J}}$ is the spectral radius of $\hat{J}$ and ${\lambda }_i$ for $i=\mathrm{1,}\cdots \mathrm{,}n$ are the eigenvalues of Jacobian matrix $\hat{J}$. Then, the mapping G is a contraction and the drive-response system (2.11) satisfies the complete synchronization properties, i.e.

$\underset{n\to \infty }{\lim }\Vert e\left(n+1\right)\Vert =\underset{n\to \infty }{\lim }\Vert {\bar{Y}}_{n+1}-{\bar{X}}_{n+1}\Vert =0$

Proof. Using the contraction mapping theorem [20].

In dynamical system point of view, it is possible for a point arbitrarily close to fixed point $Z^{\mathrm{<mo>\; *\; </mo>}}=\left(U^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{V}^{\mathrm{<mo>\; *\; </mo>}}\right)$ of system (2.11) to generate an orbit which stays arbitrarily close to $Z^{\mathrm{<mo>\; *\; </mo>}}$. An orbit which could circle around the equilibrium $Z^{\mathrm{<mo>\; *\; </mo>}}$ staying within the proposed bounded region, for initial conditions sufficiently close to $Z^{\mathrm{<mo>\; *\; </mo>}}$, could eventually approach $Z^{\mathrm{<mo>\; *\; </mo>}}$. In this case, $Z^{\mathrm{<mo>\; *\; </mo>}}$ and all invariant set of points which demonstrating the same attractive property called attractor. Using this statement, we define an attractor of drive-response system (2.11):

Definition 2.7. Let $\Lambda \subset {ℝ}^n$ such that $\Lambda$ is invariant under the function G; i.e., $G\left(\Lambda \right)\subseteq \Lambda$. We define the distance between $\Lambda$ and a point Z, as $d\left(\Lambda ,Z\right)={\min }_{z\in \Lambda }\Vert Z-z\Vert$. If there exists an $ϵ>0$ such that $d\left(\Lambda \mathrm{,}Z\right)<ϵ$ implies $li{m}_{\mathcal{l}\to \infty }d\left(\Lambda \mathrm{,}{G}^{\mathcal{l}}\left(Z\right)\right)=0$, then $\Lambda$ is called an attractor for drive-response system (2.11).

For a stable fixed point, it is of special interest to determine the set of initial conditions whose subsequent orbits end up at this fixed point and we call this set as the basin of attraction of stable fixed point which achieves by the following definition:

Definition 2.8. Given $G\mathrm{<mo>\; :\; </mo>}{ℝ}^{2n+1}\to {ℝ}^{2n}$ a continuously differentiable map, then, the compact and invariant set $\Gamma$ is called the attracting set of drive-response system (2.11) if there exists an open neighborhood B of $\Gamma$ such that $G\left(B\right)\subset B$ and ${\displaystyle {\cap }_{\mathcal{l}=0}^{\infty }{G}^{\mathcal{l}}\left(B\right)}=\Gamma$. The largest such B is called the basin of attraction for system.

Stable Threshold for Synchronization in Discrete-Time Dynamical Systems

We continue this section by defining a new concept in synchronization theory of discrete-time dynamical systems. The concept of stability in this study is similar to the one that we have in contraction mapping theorem which is different from the relative stability of equilibrium point and some nominal motion. We say that a system is stable if the final state of the system is independent on initial conditions and we call a system is attracting if the orbits of that system get pulled in or converge towards each other [21]. In general, stability can be interpreted as a property of solutions of a dynamical system which means all solutions converge towards each other [22].

Definition 2.9. A variable response system to the dynamical system ${\bar{X}}_0\mathrm{,}{\bar{X}}_1\mathrm{,}\cdots \mathrm{,}{\bar{X}}_n$, with the map $\Phi \mathrm{<mo>\; :\; </mo>}D\subseteq {ℝ}^n\to D\subseteq {ℝ}^n$ with respect to the sequence $s=\left(s_0\mathrm{,}{s}_1\mathrm{,}\cdots \mathrm{,}{s}_n\right)$ where $0<s_i\le 1$ is the sequence ${\tilde{Y}}_0=\left({\bar{Y}}_1,\cdots ,{\bar{Y}}_n\right),\cdots ,{\tilde{Y}}_{n+1}=\Phi \left(Con{t}_{{\bar{X}}_n}^{s_n}\left({\bar{Y}}_n\right)\right)$ where the map $Con{t}_u^s\mathrm{<mo>\; :\; </mo>}{ℝ}^n\to {ℝ}^n$ can be represented via $Con{t}_u^s\left(v\right)=\left(1-s\right)\stackrel{\to }{u}+s\stackrel{\to }{v}$.

Definition 2.10. The threshold for the coupled dynamical system with drive system ${\bar{X}}_0\mathrm{,}{\bar{X}}_1\mathrm{,}\cdots \mathrm{,}{\bar{X}}_n$, with the map $\Phi \mathrm{<mo>\; :\; </mo>}D\subseteq {ℝ}^n\to D\subseteq {ℝ}^n$ is s if given any $\tilde{s}<s$ and any ${\bar{Y}}_0$, there exists a sequence $s=\left(s_0\mathrm{,}{s}_1\mathrm{,}\cdots \mathrm{,}{s}_n\right)$ with $0<s_i\le 1$ and $s_i=\tilde{s}$ for $i\gg 0$ such that

$\underset{n\to \infty }{lim}\Vert {\bar{Y}}_n-{\bar{X}}_n\Vert =0$ (2.14)

Using the definitions (2.9) and (2.10), we state the main results of this section.

Theorem 2.11. Consider a linear discrete-time dynamical system (drive-response system) as following form:

$X_{n+1}=A{X}_n\mathrm{.}$ (2.15)

$Y_{n+1}=A\left(\left(1-s\right)X_n+s{Y}_n\right)\mathrm{.}$ (2.16)

where matrix A is similar to a diagonal matrix. For the values $s<{\tilde{s}}_1=\frac{1}{\rho \left(A\right)}$ where ${\tilde{s}}_1$ represents the synchronization threshold, we have

$\underset{n\to \infty }{lim}\Vert Y_n-X_n\Vert =0$

and passing this threshold decreases the stability of synchronization and consequently, the drive-response system (2.15) and (2.16) lose the complete synchronization properties.

Proof. Considering the drive system (2.15) and (2.16), we have

$\Vert Y_{n+1}-X_{n+1}\Vert \le sA\Vert Y_n-X_n\Vert \mathrm{.}$ (2.17)

To have a complete synchronization, at first we need to clarify the concept of the norm of a matrix. To find the behavior of the sequence of $\left\{A\right\}$, we need to look at the modulus of the largest eigenvalue of A. But, for the initial value $X_0=\bar{0}$, we have $X_0=X_1=X_2=\cdots =X_n=\bar{0}$. So, (2.16) can be written as $Y_{n+1}=A\left(s{Y}_n\right)$. Thus, for ${ℝ}^n$ to be the direct sum of ${ℝ}^n=V_{{\lambda }_1}\oplus V_{{\lambda }_2}\oplus \cdots \oplus V_{{\lambda }_n}$ where ${\lambda }_i$ are the eigenvalues of the matrix A, we have $Y_{n+1}=sA{Y}_n$. For the sequences $Y_1\mathrm{,}{Y}_2\mathrm{,}\cdots \mathrm{,}{Y}_n$, we have:

$\begin{array}{l}{Y}_1=sA{Y}_0\\ Y_2=sA{Y}_1=sA\left(sA\right)Y_0={\left(sA\right)}^2{Y}_0\\ ⋮\\ Y_n={\left(sA\right)}^n{Y}_0\end{array}$

we need to find minimal k where $Y_0\in V_{{\lambda }_{i_1}}\oplus V_{{\lambda }_{i_2}}\oplus \cdots V_{{\lambda }_{i_k}}$. So, for $Y_0={\displaystyle {\sum }_{j=1}^k}{v}_j$ and for all $v_j\in V_{{\lambda }_{i_j}}$, we have $Y_n={\displaystyle {\sum }_{j=1}^k}{\left(s{\lambda }_{i_j}\right)}^n{v}_j$. So, for $A{v}_i={\lambda }_{i_j}{v}_i$, we can write, $li{m}_{n\to \infty }{Y}_n=0$ if and only if $li{m}_{n\to \infty }{\left(s{\lambda }_{i_j}\right)}^n{v}_j=0$ which gives us $\left|s{\lambda }_{i_j}\right|<1$

if and only if $s\le \frac{1}{\left|{\lambda }_{i_j}\right|}$. If $Y_0$ is generic, we get $s\le \frac{1}{\left|\lambda \right|}$ for all $\lambda$ which gives $s<\frac{1}{\rho \left(A\right)}$, where, $\rho \left(A\right)$ is the spectral radius of A and can be written as

$\rho \left(A\left(n\right)\right)=max\left|{\lambda }_i\right|$. Therefore, the threshold for (2.15) and (2.16) to be completely

synchronized is $s<\frac{1}{\rho \left(A\right)}={\tilde{s}}_1$. When we pass this threshold, two systems (2.15) and (2.16) cannot be completely synchronized anymore.

Now, for the non-linear discrete-time dynamical system, we are going to develop a theorem which helps us to find the synchronization threshold.

Theorem 2.12. Given the non-linear coupled dynamical system (2.3), where the map $g\mathrm{<mo>\; :\; </mo>}D\subseteq {ℝ}^n\to D\subseteq {ℝ}^n$, for the values $s<{\tilde{s}}_2=\frac{1}{{\rho }_A+\alpha }$, we get

$\underset{n\to \infty }{lim}\Vert Y_n-X_n\Vert =0$

means that passing the synchronization threshold ${\tilde{s}}_2$ makes the drive-response system (2.3) lose the complete synchronization properties.

Proof. Suppose the following $C^r$ maps which have a fixed point at the origin:

$X_{n+1}=A{X}_n+\hat{F}\left(X_n\right)\mathrm{,}$ (2.18)

$Y_{n+1}=A\left(\left(1-s\right)X_n+s{Y}_n\right)+\hat{F}\left(\left(1-s\right)X_n+s{Y}_n\right)\mathrm{,}$ (2.19)

and the $\hat{F}\left(X_n\right)={\hat{F}}_2\left(X_n\right)+\cdots +{\hat{F}}_{r-1}\left(X_n\right)+O\left({\left|\left(X_n\right)\right|}^r\right)$ which including the vector-valued homogeneous polynomials of degree $2,\cdots \mathrm{,}r$. Consider the following equation for the error:

$\begin{array}{c}e\left(n+1\right)=Y_{n+1}-X_{n+1}\\ =\left(1-s\right)A{X}_n+sA{Y}_n-A{X}_n+\hat{F}\left(\left(1-s\right)X_n+s{Y}_n\right)-\hat{F}\left(X_n\right)\\ =sA\left(Y_n-X_n\right)+\hat{F}\left(\left(1-s\right)X_n+s{Y}_n\right)-\hat{F}(Xn)\end{array}$

By triangular inequality we can write:

$\begin{array}{c}\Vert e\left(n+1\right)\Vert =\Vert sA\left(Y_n-X_n\right)+\hat{F}\left(\left(1-s\right)X_n+s{Y}_n\right)-\hat{F}\left(X_n\right)\Vert \\ \le s{\rho }_A\Vert Y_n-X_n\Vert +s\alpha \Vert Y_n-X_n\Vert \\ =s{\rho }_A\Vert e\left(n\right)\Vert +s\alpha \Vert e\left(n\right)\Vert \end{array}$

where, ${\rho }_A$ is the spectral radius of A which is equal to ${\rho }_A=max\left|{\lambda }_i\right|$ where $\lambda$ is the root of characteristic polynomial or eigenvalue for A and $0\le \alpha <1$ and $0<s_i\le 1$ (to find the behavior of the sequence of $\left\{A\right\}$, we need to look at the modulus of the largest eigenvalue of A). So,

$\Vert e\left(n+1\right)\Vert \le s{\rho }_A\Vert e\left(n\right)\Vert +s\alpha \Vert e\left(n\right)\Vert$

We know that

$\underset{n\to \infty }{\lim }\Vert e\left(n\right)\Vert =\underset{n\to \infty }{\lim }\Vert Y_n-X_n\Vert =0$

Therefore,

$\Vert e\left(n+1\right)\Vert \le s\left({\rho }_A+\alpha \right)\Vert e\left(n\right)\Vert$

Thus, for $s\left({\rho }_A+\alpha \right)<1$ we have

$\underset{n\to \infty }{lim}\Vert Y_n-X_n\Vert =0$

for which, $s<\frac{1}{{\rho }_A+\alpha }={\tilde{s}}_2$. Here, ${\tilde{s}}_2=\beta$, which we discussed in the beginning

of this section. After passing ${\tilde{s}}_2$, we lose the complete synchronization in system (2.3).

Lemma 2.13. If drive system (2.15) becomes periodic. Then, for the values $s<{\left(\frac{1}{{\rho }_A+\alpha }\right)}^{\frac{1}{k}}:={\tilde{s}}_3$, where ${\tilde{s}}_3$ implies to the synchronization threshold, the

non-linear coupled dynamical system (2.3) becomes completely synchronized. In other word,

$\underset{n\to \infty }{lim}\Vert Y_n-X_n\Vert =0$

3. Synchronized Cycles in Generalized Nicholson Bailey (GNB) Model: Description of the Model

Generalized Nicholson Bailey (GNB) model is a generalization of the work presented by Beddington, Free and Lawton in 1975 [23]. They have investigated the complex dynamics of a host-parasitoid model which was an extension work of Nicholson-Bailey model in 1935 [24]. This model depends on three biological parameters a, k and r and has the following form

$H\left(n+1\right)=H\left(n\right){\text{e}}^{r\left(1-\frac{H\left(n\right)}{k}\right)-aP\left(n\right)},P\left(n+1\right)=H\left(n\right)\left(1-{\text{e}}^{-aP\left(n\right)}\right)$ (3.1)

where $H\left(n\right)$ presents the host population after being attacked by the parasitoid and $P\left(n\right)$ implies to the parasitoids population before they die because of biological reasons like shortage of food and or some other natural biological reasons at the end of the season n. k is the carrying capacity and shows maximum population size that can be supported by availability of all the potentially limiting resources. It is usually limited by the intensity of light and space. The fractions of hosts not parasitized is $exp\left(-aP\left(n\right)\right)$ where a is called the searching efficiency which is the probability that a given parasitoid confronts a host whole of the lifetime.

Here, we focus on the above model including a new parameter b which has the following form

$H\left(n+1\right)=H\left(n\right){\text{e}}^{r\left(1-\frac{H\left(n\right)}{k}\right)-bP\left(n\right)},P\left(n+1\right)=H\left(n\right)\left(1-{\text{e}}^{-aP\left(n\right)}\right)$ (3.2)

Without loss of generality, we assume $b=r$ and follow the same model assumptions as Asheghi in 2014 [25]. Therefore, (3.2) can be written as the form:

$H\left(n+1\right)=H\left(n\right){\text{e}}^{r\left(1-\frac{H\left(n\right)}{k}\right)-rP\left(n\right)},P\left(n+1\right)=H\left(n\right)\left(1-{\text{e}}^{-aP\left(n\right)}\right)$ (3.3)

The local dynamics of system (3.3), have been studied by different authors numerically and analytically [25] [26]. We replace $H\mathrm{,}P$ in (3.3) with $x\mathrm{,}y$ respectively and re wright (3.3) as the following form

$x\left(n+1\right)=x\left(n\right){\text{e}}^{r\left(1-\frac{x\left(n\right)}{k}\right)-ry\left(n\right)},y\left(n+1\right)=x\left(n\right)\left(1-{\text{e}}^{-ay\left(n\right)}\right)$ (3.4)

Now, we apply this coupling method on the system (3.4). Consider the following drive-response system:

$x_1\left(n+1\right)=x_1\left(n\right){\text{e}}^{r\left(1-\frac{x_1\left(n\right)}{k}\right)-r{y}_1\left(n\right)}$ (3.5)

$y_1\left(n+1\right)=x_1\left(n\right)\left(1-{\text{e}}^{-a{y}_1\left(n\right)}\right)$ (3.6)

$x_2\left(n+1\right)=p{\text{e}}^{r\left(1-\frac{p}{k}\right)-rq}$ (3.7)

$y_2\left(n+1\right)=p\left(1-{\text{e}}^{-aq}\right)$ (3.8)

where

$p=\left(1-s\right)x_1\left(n\right)+s{x}_2\left(n\right)$ (3.9)

$q=\left(1-s\right)y_1\left(n\right)+s{y}_2\left(n\right)$ (3.10)

Here, $p\mathrm{<mo>\; :\; </mo>}{ℝ}^2\to ℝ$ and $q\mathrm{<mo>\; :\; </mo>}{ℝ}^2\to ℝ$ are two continuous functions. So, if we consider the drive system $G_1\left(x_1\mathrm{,}{y}_1\right)\mathrm{<mo>\; :\; </mo>}{ℝ}^2\to {ℝ}^2$, the synchronized system would be $G_2\left(p\mathrm{,}q\right)\mathrm{<mo>\; :\; </mo>}{ℝ}^4\to {ℝ}^2$, where $p=\left(x_1\mathrm{,}{x}_2\right)$ and $q=\left(y_1\mathrm{,}{y}_2\right)$. The local stability results for the drive system (3.5) - (3.6) and (3.7) - (3.8) are the same. To investigate the qualitative dynamics of the solutions of system (3.5) - (3.8), we use several dynamical systems tools.

Synchronized Cycles in GNB Model without Parasitoid

Here, we show that in system (3.4), when the parasitoid populations go extinct (because of severe intraspecific competition or due to external factors), the dynamics of (3.4) inherits all different behaviors of Ricker curves from stable fixed point to cascade of period doubling bifurcations and then chaos [27] [28]. We rewrite the drive-response system (3.5) - (3.8) in one dimension as the following form:

$x_1\left(n+1\right)=x_1\left(n\right){\text{e}}^{r\left(1-\frac{x_1\left(n\right)}{k}\right)}$ (3.11)

$x_2\left(n+1\right)=p{\text{e}}^{r\left(1-\frac{p}{k}\right)}$ (3.12)

where

$p=\left(1-s\right)x_1\left(n\right)+s{x}_2\left(n\right)$ (3.13)

Here, $p\mathrm{<mo>\; :\; </mo>}{ℝ}^2\to ℝ$ is a continuous function. So, if we consider the drive system $R_1\left(x_1\right)\mathrm{<mo>\; :\; </mo>}ℝ\to ℝ$, the synchronized system would be $R_2\left(p\right)\mathrm{<mo>\; :\; </mo>}{ℝ}^2\to ℝ$, where p is a function of $x_1$ and $x_2$.

We have demonstrated the time-series corresponding to the solution of the drive-response system (3.11) - (3.12) in Figure 2. As we can see, with increasing the growth rate r, the behavior of system changes from stable equilibrium point to periodic behavior and then to irregular and chaotic dynamics which was expected since the system (3.11) - (3.12) has the same form and so dynamics of Ricker model.

Also, we performed a one co-dimensional bifurcation analysis for system (3.11) - (3.12) with respect to growth rate r in Figure 3(a) to discover the long term behavior of the system and we have compared the solution of drive system (3.11) with the response system (3.12) by showing the error between the solutions in Figure 3(b). As we will discuss also in Section 4, when synchronization constant s is larger, the drive-response system cannot be synchronized completely. Moreover, we have shown the Lyapunov Exponent corresponding the drive-response system (3.11) - (3.12) in Figure 3(c) which is the best place to

Figure 2. Evolution of host population $x_1$ and its coupled $x_2$ in time for for drive-response system (3.11) - (3.12) when $s=0.95$, $k=10$.

Figure 3. (a): bifurcation diagram for drive-response system (3.11) - (3.12) when $s=0.95$, $k=10$, red (drive system) and black (response system); (b): the error between the solutions of drive system and response system receptively; (c): the Lyapunov Exponent corresponding to drive-response system (3.11) - (3.12) when $s=0.95$, $k=10$, red (drive system) and black (response system).

investigate the stability or chaotic behavior of system (3.11) - (3.12). As we know, the negative Lyapunov Exponent implies to stable behavior and when it is positive, we expect to see chaotic behavior.

4. Numerical Simulations

In chaotic systems, it seems to not being possible to reproduce exactly the same initial conditions and parameters and force the orbits converge. However, in this section, we will numerically show that by using a sufficiently strong coupling method, we can change the track of the orbits to converge. Therefore, there exists a possible way to get a complete synchronization in chaotic systems whereas they have been coupled by a suitable coupling method.

In this section, we demonstrate some numerical simulation to describe the qualitative behavior of drive-response system (3.5) - (3.8). The orbits of the system (3.5) - (3.8) in chaotic regime can be considered as chaotic oscillations. Now we want to study the evolution of the dynamic variables $x_1\mathrm{,}{x}_2$ corresponding to the host population of drive system and response system, and $y_1\mathrm{,}{y}_2$ which are corresponding to the parasitoid population of drive system and response system respectively. All analysis and numerical simulations which have been conducted in this section, are expected to reveal the type of attractor from equilibrium point, periodic and quasi-periodic orbits, and chaotic attractors for which the dynamics will eventually settle down and remain forever.

In Figure 4, we can see the evolution of the attractors of drive system (red color) and response system (black color) when we change the growth rate r. As we see, as long as we are increasing r, the dynamics of the system change from the stable equilibrium point which loses stability and arises to a limit cycle. This figure demonstrates that the drive and response system with different initial conditions, and for the smaller value of the synchronization constant s, become completely synchronized.

Figure 4. Attractors for drive-response system (3.5) - (3.8) when $s=0.5$, $a=40$, $k=10$, from up left side $r=\mathrm{1.08,1.087,1.095}$, and from down left side $r=1.099,1.1,1.15$.

In Figure 5, we increase the value of the synchronization constant s and we notice that two systems keep the same qualitative behavior from stable equilibrium point to limit cycle but they are not completely synchronized.

In Figure 6, we demonstrated the evolution of host population for drive-response system (3.5) - (3.8) (The interested parameters values have been selected from the given bifurcation diagram). For fixed parameter values $s=0.5$, $a=40$, $k=10$, and different growth rate r, we can easily see how the time series for $x_1\mathrm{,}{x}_2$ change from stable and periodic oscillations to chaotic motions. However, as we discussed before, for smaller synchronization constant s, both $x_1$ and its coupled $x_2$ are completely synchronized. With previous illustration about the chaotic synchronization, we conclude that using this method of coupling, we can expect to get a complete synchronization in chaotic regime for smaller synchronization constant s which this has been demonstrated for growth rate $r=3.8$ in Figure 6.

However, when we compare Figure 6 with Figure 7, the coupling method which has been explained in Section 2, is successful when the synchronization constant s has smaller values and is closer to zero. With increasing the synchronization constant s, we noticed that two convex functions (3.9) and (3.10) which we introduced in Section 2 cannot make a complete synchronization between the drive and response systems (3.5) - (3.8). For different growth rate r, bifurcation diagrams for drive system (3.5), (3.6) and response system (3.7), (3.8), have been demonstrated for fixed synchronization constant $s=0.5$ in Figure 8 and for fixed synchronization constant $s=0.95$ in Figure 9. The one-co-dimensional bifurcation diagram helps us to know about the dependence of the drive-response dynamical systems (3.5) - (3.8) on the certain parameter which here is the growth rate r. Here, we are expecting to get completely synchronization for $s=0.5$ (Figure 8) and we do not expect to get a complete synchronization phase for $s=0.95$ (Figure 9).

As it can be easily seen in Figure 8 and Figure 9, the dynamics of host population for drive-response system (3.5) - (3.8), for some range of parameter values

Figure 5. Attractors for drive-response system (3.5) - (3.8) when $s=0.99995$, $a=40$, $k=10$, from up left side $r=\mathrm{1.08,1.087,1.095}$, and from down left side $r=1.099,1.1,1.15$.

Figure 6. Evolution of host population $x_1$ and its coupled $x_2$ in time for drive-response system (3.5) - (3.8) when $s=0.5$, $a=40$, $k=10$.

Figure 7. Evolution of host population $x_1$ and its coupled $x_2$ in time for drive-response system (3.5) - (3.8) when $s=0.5$, $a=40$, $k=10$.

Figure 8. Up: bifurcation diagram for drive-response system (3.5) - (3.8) for $a=40$, $s=0.5$, $k=10$. Down: the error between the solutions of drive system and response system receptively.

Figure 9. Up: bifurcation diagram for drive-response system (3.5) - (3.8) for $a=40$, $s=0.95$, $k=10$, red (drive system) and blue (response system). Down: the error between the solutions of drive system and response system receptively.

r and when $\left(y_1\mathrm{,}{y}_2\right)=\left(\mathrm{0,0}\right)$ (without parasitoid), is similar to classical bifurcation diagram of Ricker model where the routes to chaos happen through the cascade of period-doubling bifurcations and crisis corresponding to the extinction of parasitoid for drive and response system. However, for smaller values of growth rate r, the host and parasitoid in both systems can coexist through the periodic and quasi-periodic cycles of the Neimark-Sacker bifurcation of interior steady states. We also can easily observe that the sudden changes of attractors (crisis) happen frequently when we increase the parameter values of r.

5. Conclusion

In this paper, we developed a drive-response system by defining a convex continuous link function which maps the orbits of the drive system into the orbits of its coupled system and keeps the same qualitative dynamics. We found an appropriate normal form for drive-response system and we obtained the conditions under which the solutions of drive and response system become completely synchronized. We provided a new concept in chaos synchronization, called, synchronization threshold, which means that the solutions of drive and response system diverge from each other and lose the complete synchronization properties when they pass the threshold. Also, we studied a coupled discrete-time two-dimensional host-parasitoid model which is a generalization of famous Nicholson-Bailey model. One of our objectives in this paper was to investigate the rich dynamics of drive-response system (3.5) - (3.8) around its equilibrium points and achieving the chaos synchronization. We developed a new drive-response system by defining a convex continuous link function which maps the orbits of the drive system keeping the same qualitative properties such as stability and periodicity into the orbits of its coupled system. We observed that this coupling method can be successful for drive-response system (3.5) - (3.8) to get a complete synchronization when the synchronization constant has smaller values, closer to zero. Moreover, numerical verification is performed to show the existence of wide range of dynamics of drive-response system (3.5) - (3.8) around the positive equilibrium point. We also changed the values of synchronization constant s in its range between (0, 1) and we observed that the response system (3.7) - (3.8) for smaller values of synchronization constant s is completely synchronized with its original drive system (3.5) - (3.6) and when we increased the values of synchronization constant s, we noticed that the qualitative behaviors of both systems remain the same, however, we do not get a complete synchronization between the solutions of drive and response system (3.5) - (3.8). In chaotic regime, for larger values of synchronization constant s, closer to one, we could not get a complete synchronization. But, for smaller synchronization constant s, closer to zero, we have shown that two systems are in complete synchronization when the dynamic is chaotic.

Acknowledgements

This work was supported by the Institute of Computational Comparative Medicine (ICCM) and Department of Mathematics of Kansas State University. With a special thanks to Dr. Majid Jaberi-Douraki for his full support.

Appendix

We have demonstrated different types of attractors of drive-response system (3.5) - (3.8) in Figures A1-A4 when we are changing the threshold s.

Figure A1. Attractors for drive-response system (3.5) - (3.8) when $s=0.94$, $a=5$, $k=1.5$, from up left side $r=\mathrm{2.0,2.2,2.3}$, and from down left side $r=\mathrm{2.5,2.7,2.8}$.

Figure A2. Attractors for drive-response system (3.5) - (3.8) when $s=0.95$, $a=5$, $k=1.5$, from up left side $r=\mathrm{2.0,2.2,2.3}$, and from down left side $r=\mathrm{2.5,2.7,2.8}$.

Figure A3. Attractors for drive-response system (3.5) - (3.8) when $s=0.96$, $a=5$, $k=1.5$, from up left side $r=\mathrm{2.0,2.2,2.3}$, and from down left side $r=\mathrm{2.5,2.7,2.8}$.

Figure A4. Attractors for drive-response system (3.5) - (3.8) when $s=0.99$, $a=5$, $k=1.5$, from up left side $r=\mathrm{2.0,2.2,2.3}$, and from down left side $r=\mathrm{2.5,2.7,2.8}$.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Lorenz, E.N. (1963) Deterministic Nonperiodic Flow. Journal of the Atmospheric Sciences, 20, 130-141. https://doi.org/10.1175/1520-0469(1963)020<0130:DNF>2.0.CO;2
[2]	Pikovsky, A., Kurths, J., Rosenblum, M. and Kurths, J. (2003) Synchronization: A Universal Concept in Nonlinear Sciences. Cambridge University Press, Cambridge, 12.
[3]	Fujisaka, H. and Yamada, T. (1983) Stability Theory of Synchronized Motion in Coupled-Oscillator Systems. Progress of Theoretical Physics, 69, 32-47. https://doi.org/10.1143/PTP.69.32
[4]	Yamada, T. and Fujisaka, H. (1983) Stability Theory of Synchronized Motion in Coupled-Oscillator Systems. II: The Mapping Approach. Progress of Theoretical Physics, 70, 1240-1248. https://doi.org/10.1143/PTP.70.1240
[5]	Pecora, L.M. and Carroll, T.L. (1990) Synchronization in Chaotic Systems. Physical Review Letters, 64, 821. https://doi.org/10.1103/PhysRevLett.64.821
[6]	Balmforth, N.J., Jacobson, A. and Provenzale, A. (1999) Synchronized Family Dynamics in Globally Coupled Maps. Chaos: An Interdisciplinary Journal of Nonlinear Science, 9, 738-754. https://doi.org/10.1063/1.166448
[7]	Moran, P.A.P. (1953) The Statistical Analysis of the Canadian Lynx Cycle. Australian Journal of Zoology, 1, 291-298. https://doi.org/10.1071/ZO9530291
[8]	Bulmer, M.G. (1974) A Statistical Analysis of the 10-Year Cycle in Canada. The Journal of Animal Ecology, 47, 701-718. https://doi.org/10.2307/3532
[9]	Korpimäki, E. and Krebs, C.J. (1996) Predation and Population Cycles of Small Mammals: A Reassessment of the Predation Hypothesis. BioScience, 46, 754-764. https://doi.org/10.2307/1312851
[10]	Ranta, E., Kaitala, V. and Lundberg, P. (1997) The Spatial Dimension in Population Fluctuations. Science, 278, 1621-1623. https://doi.org/10.1126/science.278.5343.1621
[11]	Gurney, W.S.C., Crowley, P.H. and Nisbet, R.M. (1992) Locking Life-Cycles onto Seasons: Circle-Map Models of Population Dynamics and Local Adaptation. Journal of Mathematical Biology, 30, 251. https://doi.org/10.1007/BF00176151
[12]	Blasius, B. and Stone, L. (2000) Chaos and Phase Synchronization in Ecological Systems. International Journal of Bifurcation and Chaos, 10, 2361-2380. https://doi.org/10.1142/S0218127400001511
[13]	Earn, D.J.D., Rohani, P. and Grenfell, B.T. (1998) Persistence, Chaos and Synchrony in Ecology and Epidemiology. Proceedings of the Royal Society of London. Series B: Biological Sciences, The Royal Society, 265, 7-10. https://doi.org/10.1098/rspb.1998.0256
[14]	Azizi, T. and Kerr, G. (2020) Chaos Synchronization in Discrete-Time Dynamical Systems with Application in Population Dynamics. Journal of Applied Mathematics and Physics, 8, 406-423. https://doi.org/10.4236/jamp.2020.83031
[15]	Kolmogorov, A.N. and Fomin, S.V. (1957) Elements of the Theory of Functions and Functional Analysis.
[16]	Wiggins, S. (2003) Introduction to Applied Nonlinear Dynamical Systems and Chaos. Princeton University Press, Princeton, NJ, 2.
[17]	Grobman, D.M. (1959) Homeomorphism of Systems of Differential Equations. Doklady Akademii Nauk SSSR, 128, 880-881.
[18]	Hartman, P. (1960) On Local Homeomorphisms of Euclidean Spaces. Boletín de la Sociedad Matemática Mexicana, 5, 220-241.
[19]	Hartman, P. (1960) A Lemma in the Theory of Structural Stability of Differential Equations. Proceedings of the American Mathematical Society, 11, 610-620. https://doi.org/10.1090/S0002-9939-1960-0121542-7
[20]	Conrad, K. (2014) The Contraction Mapping Theorem, II.
[21]	Russo, G. (2010) Analysis, Control and Synchronization of Nonlinear Systems and Networks via Contraction Theory: Theory and Applications. University of Naples Federico II, Italy.
[22]	Angeli, D. (2002) A Lyapunov Approach to Incremental Stability Properties. IEEE Transactions on Automatic Control, 47, 410-421. https://doi.org/10.1109/9.989067
[23]	Beddington, J.R., Free, C.A. and Lawton, J.H. (1975) Host-Parasitoid Dynamics of a Generalized Thompson Model. Nature, 255, 719-732.
[24]	Nicholson, A.J. and Bailey, V.A. (1935) The Balance of Animal Populations Part I. Proceedings of the Zoological Society of London, 105, 551-598. https://doi.org/10.1111/j.1096-3642.1935.tb01680.x
[25]	Asheghi, R. (2014) Bifurcations and Dynamics of a Discrete Predator: Prey System. Journal of Biological Dynamics, 8, 161-186. https://doi.org/10.1080/17513758.2014.927596
[26]	Kapçak, S., Ufuktepe, ü. and Elaydi, S. (2013) Stability and Invariant Manifolds of a Generalized Beddington Host-Parasitoid Model. Journal of Biological Dynamics, 148, 233-253. https://doi.org/10.1080/17513758.2013.849764
[27]	Ricker, W.E. (1954) Stock and Recruitment. Journal of the Fisheries Board of Canada, 11, 559-623. https://doi.org/10.1139/f54-039
[28]	Azizi, T., et al. (2015) Dynamics of a Discrete-Time Plant-Herbivore Model. Caspian Journal of Mathematical Sciences, 4, 241-256.