Next Article in Journal
Three-Dimensional Lorentz-Invariant Velocities
Next Article in Special Issue
Numerical Solution of External Boundary Conditions Inverse Multilayer Diffusion Problems
Previous Article in Journal
A Unified Hardware Design for Multiplication, Division, and Square Roots Using Binary Logarithms
Previous Article in Special Issue
Improving Influenza Epidemiological Models under Caputo Fractional-Order Calculus
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Invariant Sets, Global Dynamics, and the Neimark–Sacker Bifurcation in the Evolutionary Ricker Model

1
Department of Mathematics, University of Madeira, Campus Universitário da Penteada, 9020-105 Funchal, Portugal
2
Center for Mathematical Analysis, Geometry, and Dynamical Systems, Instituto Superior Tecnico, Technical University of Lisbon, Av. Rovisco Pais, 1049-001 Lisbon, Portugal
3
Department of Mathematics, California State University Bakersfield, Bakersfield, CA 93311-1022, USA
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Symmetry 2024, 16(9), 1139; https://doi.org/10.3390/sym16091139
Submission received: 5 August 2024 / Revised: 26 August 2024 / Accepted: 29 August 2024 / Published: 2 September 2024
(This article belongs to the Special Issue Symmetry in Mathematical Models)

Abstract

:
In this paper, we study the local, global, and bifurcation properties of a planar nonlinear asymmetric discrete model of Ricker type that is derived from a Darwinian evolution strategy based on evolutionary game theory. We make a change of variables to both reduce the number of parameters as well as bring symmetry to the isoclines of the mapping. With this new model, we demonstrate the existence of a forward invariant and globally attracting set where all the dynamics occur. In this set, the model possesses two symmetric fixed points: the origin, which is always a saddle fixed point, and an interior fixed point that may be globally asymptotically stable. Moreover, we observe the presence of a supercritical Neimark–Sacker bifurcation, a phenomenon that is not present in the original non-evolutionary model.

1. Introduction

The standard one-dimensional Ricker equation refers to the discrete population model x n + 1 = r x n e a x n with parameters r, a and an initial condition x 0 > 0 . This model is named after Ricker’s work fitting the parameters to reproductive curves for fish populations, including the Pacific Herring, Pink Salmon, Coho Salmon, and Sockeye Salmon, among others [1]. From a mathematical perspective, an increase in r leads to period doubling bifurcations and eventually chaos [2]. In a controlled laboratory, such bifurcations were observed in Tribolium flour beetles [3]. However, as noted in [4], ecologists have rarely observed chaos in animal populations in nature. This suggested that, in a real-world setting, one-dimensional models such as the Ricker model are missing important factors.
In a series of papers [5,6,7,8,9,10,11], Jim Cushing, along with some coauthors, established the foundations of discrete evolutionary models using the idea of Darwinian evolution strategy [12] based on evolutionary game theory [13]. These evolutionary game models differ in significant ways from classical theory, while classical games focus on strategies that optimize players’ payoffs, evolutionary games emphasize strategies that persist over time. Through births and deaths, players come and go, but their strategies are passed down from generation to generation [14].
The first appearance of evolutionary models in discrete dynamical systems was in [6], applied to the classic Beverton–Holt (discrete logistic) model. Using standard methods for modeling evolution, the one-dimensional model becomes a planar system of coupled difference equations, consisting of a Beverton–Holt type equation for the population dynamics and a difference equation for the dynamics of the mean phenotype trait. By studying the case when the trait equation uncouples from the population dynamic equation, Cushing was able to obtain criteria under which the evolutionary system has globally asymptotically stable equilibria or periodic solutions.
Recently, Ch-Chaoui and Moknni [15] studied the coupled case of the evolutionary Beverton–Holt model and discussed the existence of the positive fixed point, exploring its asymptotic stability. Under certain parametric conditions, they discussed the Neimark–Sacker bifurcation, a phenomenon that does not occur in the classic Beverton–Holt model.
This strategy of deriving an evolutionary model was extended to more general discrete models [7,9,11]. In particular, models with the Allee effect [10], models of Ricker type [8,16,17] and multi-species models [17,18,19].
In all these studies, the main mathematical goals are to find and classify the stability of equilibria, identify bifurcations, and when possible, describe the global dynamics of the models. In the case of a lone locally stable equilibrium, a great deal of studying has been performed on whether this fixed point is also globally attracting. While there exist many general theorems on global stability in one-dimensional mappings [20,21,22], this problem is more difficult in higher dimensions. Some successful examples of establishing global stability in higher dimensions include of those of monotone and mixed-monotone systems [19,23,24,25,26,27]. For non-monotone systems, singularity curves [28,29,30,31] have also been used to establish global stability if the mapping satisfies certain topological constraints [32,33,34,35]; these curves will also play an important role in the present work. However, in full generality, the question of global stability even in the planar case is far from understood. Some open problems and conjectures in two-dimensions may be found in [36].
In a recent paper [37], the author proposed a series of open problems and conjectures concerning the global stability of the following evolutionary discrete model of Ricker type for a single-species x (with a single mean trait u and individual trait v)
x ( t + 1 ) = x ( t ) e α ( t ) u 2 ( t ) 2 c 0 x ( t ) u ( t + 1 ) = ( 1 σ 2 ) u ( t ) c 1 σ 2 x ( t ) ,
where the parameters α ( t ) are the sequence of growth rates (which are functions of the trait), c 0 is the competition coefficient, σ 2 0 is the speed of evolution, and c 1 is the sensitivity of the competition. A complete derivation and explanation of these parameters may be found in [8].
Notice that if
α ( t ) = α , for all t = 0 , 1 , 2 , ,
then Model (1) is autonomous whereas if
α ( t + p ) = α ( t ) , for all t = 1 , 2 , and some integer p > 1 ,
then Model (1) is p -periodic.
We remark that the first Equation in (1) is known in the literature as the inherent fertility equation while the second equation is the trait equation. Other names for the trait equation include the canonical equation for evolution [38], Lande’s Equation [39,40] and Fisher’s Equation [41].
In the autonomous case, for the parameter values 1 < σ 2 < 2 , Model (1) is mixed-monotone, and the following result addresses the global stability of the model:
Theorem 1 
([19]). Let α ( t ) = α > 0 for all t = 0 , 1 , 2 , such that 0 < α < 1 , 1 < σ 2 < 2 , and c 1 0 . Then, the interior equilibrium point E = ( x , c 1 x ) of Model (1) is globally asymptotically stable if it is locally asymptotically stable, where x = c 0 + c 0 2 + 2 α c 1 2 / c 1 2 .
The conditions of local stability required in Theorem 1 are
c 0 ( 2 + σ 2 ) x < 2 ( ( α 1 ) σ 2 + 2 )
and
c 0 ( 1 + σ 2 ) x > ( 2 α 1 ) σ 2 ,
with α > 0 , σ 2 < 2 and c 1 0 as stated in Theorem 2.1 in [19]. In particular, Theorem 1 shows that local stability implies global stability in the parameter range 1 < σ 2 < 2 . In this paper, we address the case of 0 < σ 2 < 1 , where the mapping is no longer mixed-monotone and new techniques are needed.
In this paper, we solve one of the conjectures presented in [37]; specifically, we demonstrate the global stability of the interior fixed point of Model (1) when 0 < α < 1 and σ 2 < 1 , thereby completing Theorem 1 in the global stability approach. Furthermore, we go beyond the concept of global stability by characterizing the invariant manifolds of the fixed points and also the bifurcation of the interior fixed point. As is seen in populations in nature, we will prove that this model does not undergo a period doubling cascade to chaos.
In order to achieve these goals, in Section 2 we make a change in the variables in Model (1) and write a simplified model that will help us understand both the local and global properties of the original model.
Using the concept of critical curves [29], in Section 3 we identify a forward-invariant and globally attracting subset under the action of the system’s mapping. Subsequently, in Section 4 and Section 5, we analyze the stability and bifurcation of the fixed points within this invariant set. Finally, in the last two sections, we provide the conclusion and discuss future work.
We believe that the techniques followed in this paper will aid in the study of similar models, namely, the use of critical curves and isoclines for identifying the maximal invariant set where all the dynamics take place; the characterization of the invariant (stable and unstable) manifolds that attract all orbits in the invariant set; and the classification of the Neimark–Sacker bifurcation using a geometric argument, which, to our knowledge, is the first time that a supercritical case has been classified via such an argument.

2. Notation and Preliminaries

By using a change in variables in Model (1), it is possible to reduce the number of parameters involved, thereby simplifying the study.
Let us consider the autonomous system
w n + 1 = w n e α 1 2 u n 2 c 0 w n u n + 1 = ( 1 σ 2 ) u n c 1 σ 2 w n ,
equivalent to Model (1), for α > 0 , c 0 0 , c 1 0 , and 0 < σ 2 < 1 .
Let x n = c 0 w n and y n = c 0 c 1 u n . Then,
x n + 1 = c 0 w n + 1 = c 0 w n e α 1 2 u n 2 c 0 w n = x n e α 1 2 u n 2 x n = x n e α c 1 2 2 c 0 2 y n 2 x n
and
y n + 1 = c 0 c 1 u n + 1 = c 0 c 1 ( 1 σ 2 ) u n + σ 2 c 0 w n = ( 1 σ 2 ) y n + σ 2 x n .
Hence, the properties of Model (2) are the same as the following system
x n + 1 = x n e α c 1 2 2 c 0 2 y n 2 x n y n + 1 = ( 1 σ 2 ) y n + σ 2 x n .
Let F be the mapping representing Model (3) given by
F ( x , y ) = f ( x , y ) , g ( x , y ) = x e α x β y 2 , ( 1 σ 2 ) y + σ 2 x ,
where β = c 1 2 2 c 0 2 for x 0 .
In the following sections, we will study the dynamics of invariant sets under the mapping F representing Model (3), as well as some properties of the fixed points of F.

3. Critical Curves and Invariant Sets

We start by examining the preimages of F representing Model (3). Initially, we follow the work of [29]. They assume there is a mapping F that is a noninvertible planar mapping whose preimage at any point consists of a finite number of isolated points. They then divide the plane into regions Z k , where ( x , y ) Z k if and only if ( x , y ) has k preimages under F. Following their notation, the boundaries of these regions is denoted by L C , with preimage L C 1 = F 1 ( L C ) . For smooth mappings, they show that L C 1 can be determined by the set of points where the determinant of the Jacobian of F vanishes, in which case L C = F ( L C 1 ) .
For the map F representing Model (3), we will show that L C separates the first and fourth quadrant into two distinct regions Z 2 and Z 0 , as is illustrated in Figure 1. All points in the interior of the shaded region have two distinct preimages while points in the unshaded region have no preimage, as will be shown in Theorem 2.
The Jacobian of the mapping F is given by
D F = J ( x , y ) = ( 1 x ) e α x β y 2 2 β x y e α x β y 2 σ 2 1 σ 2 .
Then,
det ( J ( x , y ) ) = ( 1 x ) ( 1 σ 2 ) + 2 β σ 2 x y e α x β y 2
and this determinant vanishes along the graph of
y = ϕ ( x ) = ( x 1 ) ( 1 σ 2 ) 2 β σ 2 x
with x > 0 . Consequently, L C 1 = { ( x , y ) R 2 : y = ϕ ( x ) , x > 0 } and L C = F ( L C 1 ) . We now look at the geometry of L C 1 and L C in more detail.
Lemma 1. 
Let α , β > 0 , and let σ ( 0 , 1 ) . Then the graph of L C 1 is a monotone increasing and concave curve. It has a vertical asymptote at x = 0 and a horizontal asymptote at y = 1 σ 2 2 β σ 2 .
Proof. 
The first claim follows by calculating derivatives of ϕ . The asymptotes follow by taking the appropriate limits. □
Lemma 2. 
Let α , β > 0 , and let σ ( 0 , 1 ) . Then, the curve L C can be decomposed into curves γ 1 from ( 0 , ) to ( e α 1 , σ 2 ) and γ 2 from ( 0 , ) to ( e α 1 , σ 2 ) where γ 1 is strictly increasing and γ 2 is strictly decreasing. At the intersection point ( e α 1 , σ 2 ) , L C has a vertical tangent line.
Proof. 
The curve L C has parametrization
r ( t ) = f 1 ( t ) , g 1 ( t ) = t e α t β ϕ ( t ) 2 , ( 1 σ 2 ) ϕ ( t ) + σ 2 t
for t ( 0 , ) . We claim that γ 1 consists of the portion of r ( t ) for t ( 0 , 1 ] and that γ 2 consists of the part with t [ 1 , ) . We calculate
d f 1 ( t ) d t = e α t β ϕ ( t ) 2 1 t 2 β t ϕ ( t ) ϕ ( t )
and
d g 1 ( t ) d t = ( 1 σ 2 ) ϕ ( t ) + σ 2 .
It is clear that d g 1 d t > 0 for all t. Also, for t < 1 , both 1 t > 0 and ϕ ( t ) ϕ ( t ) > 0 , so d f 1 d t > 0 for t ( 0 , 1 ) . Similarly, for t > 1 , both 1 t < 0 and ϕ ( t ) ϕ ( t ) < 0 , so d f 1 d t < 0 for t > 1 . Since d f 1 ( t ) d t vanishes at t = 1 , the curve has a vertical tangent line at r ( 1 ) = ( e α 1 , σ 2 ) . □
We remark that the relative positions of the critical curves as well their properties, as described in Lemmas 1 and 2, are clearly illustrated in Figure 2. As we can observe, the critical curves do not intersect, and there are points on the left side of L C with the same y component that have two distinct preimages symmetric about a point on the critical curve L C 1 .
Before we proceed with the analysis of our model, let us take an overall look at a concrete example:
Example 1. 
Let α = 1 2 , β = 2 and σ = 3 4 . Then, the critical curve L C 1 = { ( x , y ) R 2 : y = 7 ( x 1 ) 36 x } and a parametrization of L C is written as
r ( t ) = f 1 ( t ) , g 1 ( t ) = t e 49 ( t 1 ) 2 648 t 2 t + 1 2 , 324 t 2 + 49 t 49 576 t , t ( 0 , ) .
Clearly, the conclusion of the precedent theorem is satisfied since g 1 ( t ) = 324 t 2 + 49 576 t 2 and
f 1 ( t ) = ( t 1 ) 324 t 2 + 49 324 t 2 e 49 ( t 1 ) 2 648 t 2 t + 1 2 .
A representation of the graphs of L C 1 and L C is shown in Figure 2. As one can observe the curve, L C divides the first and fourth quadrant into two distinct regions. Now, let A = 1 50 , 21 32 , B = 3 25 , 21 32 , C = 1 3 , 21 32 , and D = x ¯ , 21 32 , all taken with the same y-coordinate, where x ¯ is the intersection of the line y = 21 32 and the curve L C . To calculate x ¯ , we use the parametrization g 1 ( t ) and solve 324 t ¯ 2 + 49 t ¯ 49 576 t ¯ = 21 32 , yielding t ¯ = 329 648 + 7 3505 648 . Then, x ¯ = f 1 ( t ¯ ) 0.5998 .
A forward computation shows that the first three points have two preimages as follows
F 1 ( A ) { ( 0.25026 , 1.17824 ) , ( 2.12565 , 1.23297 ) }
F 1 ( B ) { ( 0.50788 , 0.84702 ) , ( 1.81746 , 0.83673 ) }
F 1 ( C ) { ( 0.75425 , 0.53025 ) , ( 1.55036 , 0.49333 ) }
while D has a single preimage on L C 1 given approximately by
F 1 ( D ) ( 1.14725 , 0.02496 ) .
All of the preimages fall on the line y = 3 2 9 7 x , marked as L in Figure 2.
The precedent example motivates the following result:
Theorem 2. 
Let α , β > 0 , and let σ ( 0 , 1 ) . Then, the curve L C divides the first and fourth quadrant into two regions. For the region adjacent to the y-axis, each point with positive x coordinate has two distinct preimages. The other region is not in the range of F.
Proof. 
The conclusion that L C divides the plane into two regions follows from Lemma 2. To prove the remaining claim, we consider the set of lines with slope σ 2 1 σ 2 , whose image we will show is a horizontal line. Let
L ( x , y ) = ( x ( 1 σ 2 ) x t , y + σ 2 x t ) : < t 1 1 σ 2 .
Since L C 1 is monotonically increasing, each line intersects L C 1 exactly once, so it suffices to consider the set of lines L ( x , ϕ ( x ) ) .
Noting that 2 β x ϕ ( x ) σ 2 = ( x 1 ) ( 1 σ 2 ) , we have
f L ( x , ϕ ( x ) ) = x 1 ( 1 σ 2 ) t e α x + ( 1 σ 2 ) x t β ( ϕ ( x ) ) 2 2 β σ 2 ϕ ( x ) x t β σ 4 x 2 t 2 = x e α x β ( ϕ ( x ) ) 2 1 ( 1 σ 2 ) t e ( 1 σ 2 ) x t 2 β σ 2 ϕ ( x ) x t β σ 4 x 2 t 2 = f ( x , ϕ ( x ) ) 1 ( 1 σ 2 ) t e ( 1 σ 2 ) x t ( x 1 ) ( 1 σ 2 ) t β σ 4 x 2 t 2 = f ( x , ϕ ( x ) ) 1 ( 1 σ 2 ) t e ( 1 σ 2 ) t β σ 4 x 2 t 2
and
g L ( x , ϕ ( x ) ) = ϕ ( x ) + σ 2 x t + σ 2 x ( 1 σ 2 ) x t ϕ ( x ) σ 2 x t = ϕ ( x ) + σ 2 x t + σ 2 x σ 2 x t + σ 4 x t σ 2 ϕ ( x ) σ 4 x t = ϕ ( x ) + σ 2 x σ 2 ϕ ( x ) = g ( x , ϕ ( x ) ) .
Thus, the image of the line L ( x , ϕ ( x ) ) is horizontal since the y coordinate is fixed at g ( x , ϕ ( x ) ) . To study what happens to the x coordinate, we look at the multiplicative factor
h ( t ) = 1 ( 1 σ 2 ) t e ( 1 σ 2 ) t β σ 4 x 2 t 2 .
When t = 0 , h ( t ) = 1 , so the x coordinate f L ( x , ϕ ( x ) ) is f ( x , ϕ ( x ) ) , i.e., it lies on L C . Since
h ( t ) = 2 β σ 4 x 2 t ( 1 ( 1 σ 2 ) t ) ( 1 σ 2 ) 2 t e ( 1 σ 2 ) t β σ 4 x 2 t 2 ,
h is increasing for t < 0 and decreasing for 0 < t 1 1 σ 2 . Consequently, the image of the line is always to the left of L C . Moreover, it traces the line from 0 , g ( x , ϕ ( x ) ) to f ( x , ϕ ( x ) ) , g ( x , ϕ ( x ) ) exactly once for < t < 0 , and the image traces the horizontal line in the reverse direction for 0 < t 1 1 σ 2 . □
To summarize, every point to the left of L C with positive x coordinate has exactly two preimages, one to the upper left of L C 1 and one to its lower right. The map’s orientation is preserving in the upper left and the orientation is reversing in the lower right, best seen by the images of the lines L ( x , ϕ ( x ) ) . From the symmetry of the image of the mapping about these line, we only need to consider the mapping F restricted to Z 2 L C . In fact, for sufficiently large n the y coordinate of F n is always between 0 and 1. This motivates the final theorem of this subsection.
Theorem 3. 
Let U denote the closed region bounded by the coordinate axes, the curve L C , and the line y = 1 . Let α ( 0 , 1 ) , β > 0 , and σ ( 0 , 1 ) . Then, U is forward invariant under F, F : U U is injective, and U is globally attracting.
Proof. 
First, we establish that U is invariant. Let ( x , y ) U . From the discussion above, f ( x , y ) always lies to the left of L C , so we must prove that ( x , y ) U implies g ( x , y ) = ( 1 σ 2 ) y + σ 2 x [ 0 , 1 ] . Since x lies to the left of L C , we must have 0 x < 1 . Then, we have
0 ( 1 σ 2 ) y + σ 2 x 1 σ 2 + σ 2 1
establishing that U is forward invariant.
The graph of L C 1 crosses the x-axis at the point ( 1 , 0 ) . Thus, the graph of L C 1 is always to the right of U. Consequently, at most one preimage of a point in U can lie in U since the other preimage must lie to the lower right of L C 1 . This shows that the restriction of F to U is injective.
Finally, we show U is globally attracting. First, observe that the first quadrant is forward invariant. Note that if y < 0 , then g ( x , y ) = y + σ 2 ( x y ) > y . If F n ( x , y ) never maps to the first quadrant, then the set of y coordinates of { F n ( x , y ) } n Z is a bounded and strictly monotone increasing sequence that must converge to a value y 0 . If y < 0 , this is a contradiction since there is no invariant set of points with fixed y negative coordinate by the inequality above in g. Consequently, either the sequence { F n ( x , y ) } n Z has its y coordinate approach 0, or it eventually maps to a positive value.
Thus, it suffices to consider ( x , y ) with y 0 . We will show there exists an integer k such that the y coordinate of F k ( x , y ) is less than or equal to 1. Since f ( x , y ) < 1 (since the x coordinate lies to the left of L C ), if g ( x , y ) > 1 , then we have
g f ( x , y ) , g ( x , y ) = g ( x , y ) + σ 2 f ( x , y ) g ( x , y ) < g ( x , y ) .
If the y value of F n ( x , y ) always exceeds 1, then similar to before, the y coordinates of { F n ( x , y ) } n Z + form a bounded and strictly monotone decreasing sequence that must converge to a value y 1 . As there is no invariant set with y > 1 , the theorem is proved. □
In Figure 3, we depict a prototype of the closed invariant region U mentioned in Theorem 3, where the mapping F is injective. As is clearly shown, only a portion of the critical curve L C on the boundary of this region will move strictly according to the values of the parameters, while the remaining three segments of the boundary are defined by conditions that do not depend on the parameters. This region U is crucial for the global analysis of the model, as all the dynamics of (3) occur within it.

4. Isoclines and Fixed Points

Let us denote the restriction of F to U, mentioned in the previous section, by F U . The isoclines of F ( x , y ) = f ( x , y ) , g ( x , y ) are given by the set of points that satisfy f ( x , y ) = x and the points that satisfy g ( x , y ) = y . For F U : U U with
F U ( x , y ) = x e α x β y 2 , ( 1 σ 2 ) y + σ 2 x
the isoclines can be broken into three distinct sets:
I 1 = { ( x , y ) U : y = x } , I 2 = { ( x , y ) U : x + β y 2 = α }
and
I 3 = { ( x , y ) U : x = 0 } .
From the isoclines, we easily obtain that there are two fixed points, i.e., points ( x , y ) U with F ( x , y ) = ( x , y ) .
Proposition 1. 
Let α ( 0 , 1 ) , β > 0 , and σ ( 0 , 1 ) . Then, the mapping F from model (3) has exactly two fixed points given by ( 0 , 0 ) and ( x , x ) with x = 1 + 4 α β 1 2 β = 2 α 1 + 1 + 4 α β .
Proof. 
To be a fixed point, one either needs ( x , y ) I 1 I 2 or I 1 I 3 . In the first case, the fixed point satisfies x = α β ( x ) 2 with x 0 . The positive solution to this quadratic yields x = 1 + 4 α β 1 2 β = 2 α 1 + 1 + 4 α β . The second case I 1 I 3 immediately yields the origin as the other fixed point. □
Note that the isoclines I 1 and I 2 each break U into two regions, which we denote as follows:
I 1 = { ( x , y ) U : y x } I 1 + = { ( x , y ) U : y x } I 2 = { ( x , y ) U : x + β y 2 α } I 2 + = { ( x , y ) U : x + β y 2 α } .
We also define
R 1 = I 1 I 2 R 2 = I 1 I 2 + R 3 = I 1 + I 2 + R 4 = I 1 + I 2 .
Figure 4 illustrates the relative positions between the four regions R i , the three isoclines I j , and the critical curve L C . This figure is crucial for understanding the proof of Theorem 4 and provides a clear insight into the dynamics of the orbits in the neighborhood of the interior fixed point ( x , x ) , as will be shown subsequently.
Theorem 4. 
Let α ( 0 , 1 ) , β > 0 , and σ ( 0 , 1 ) . Then, the following relations hold:
1. 
F R 1 I 1 ,
2. 
F R 2 I 2 + ,
3. 
F R 3 I 1 + ,
4. 
F R 4 I 2 .
Proof. 
We prove (1) and (2) as the other proofs are similar with reversed inequalities.
1.
Let ( x , y ) R 1 = I 1 I 2 , so that 0 y x and x + β y 2 α . Then,
y ( 1 σ 2 ) x ( 1 σ 2 ) x ( e α x β y 2 σ 2 ) = f ( x , y ) σ 2 x ,
which after rearrangement yields g ( x , y ) f ( x , y ) . That is, F U ( x , y ) I 1 .
2.
First, note that if 0 x 1 x 2 1 , then f ( x 1 , y ) f ( x 2 , y ) and g ( x 1 , y ) g ( x 2 , y ) . Since the graph of I 2 is decreasing, it suffices to prove the claim on the left boundary of I 1 I 2 + . That is, it suffices to prove the result for ( x , y ) where x x and either x = y or x + β y 2 = α .
In the first case,
f ( x , y ) + β g ( x , y ) 2 = f ( x , x ) + β g ( x , x ) 2 = x e α x β x 2 + β x 2
Since the derivative of x e α x β x 2 + β x 2 is ( 1 x ) e α x β x 2 + 2 β x ( 1 x e α x β x 2 ) , it is a strictly increasing function of x for x x , and we have that
x e α x β x 2 + β x 2 x e α x β ( x ) 2 + β ( x ) 2 = x + β ( x ) 2 = α .
In the second case, f ( x , y ) = x and g ( x , y ) y , so
f ( x , y ) + β g ( x , y ) 2 x + β y 2 α .
Thus, not only do points in U proceed counterclockwise around ( x , x ) under the mapping F U , but in each loop around ( x , x ) they cannot skip any of the R i in their itinerary. Consequently, we have the following corollary.
Corollary 1. 
Let α ( 0 , 1 ) , β > 0 , and σ ( 0 , 1 ) . Then, F U has no points with minimal period 2 or 3.
Now, we turn to the stability of the two fixed points and their bifurcations.
Theorem 5. 
The origin is always a saddle point. The interior fixed point ( x , x ) is locally asymptotically stable if and only if
σ 2 ( 2 α 1 ) 1 + σ 2 < x .
Proof. 
We use a Jacobian matrix. We calculate
J ( 0 , 0 ) = e α 0 σ 2 1 σ 2
which has eigenvalues e α > 1 and 1 σ 2 ( 0 , 1 ) , from which it follows that the origin is a saddle point. For the interior fixed point, we have
J ( x , x ) = 1 x 2 β ( x ) 2 σ 2 1 σ 2 .
Using x α = β ( x ) 2 , we may also express this as
J ( x , x ) = 1 x 2 x 2 α σ 2 1 σ 2 .
From the Jury conditions [42], a necessary and sufficient criteria for local asymptotic stability of a fixed point is given by
| Tr ( J ( x , x ) | < det ( J ( x , x ) ) + 1 < 2 .
We claim that the first of these inequalities is always satisfied. Indeed, for the parameter values given, the trace is always positive since x < α < 1 , so to verify the first inequality we simply observe that the stated inequality is equivalent to
2 x σ 2 < 2 x σ 2 + x σ 2 + 2 β ( x ) 2 σ 2
which is equivalent to the true statement
0 < x σ 2 + 2 β ( x ) 2 σ 2 .
Using the second expression for J ( x , x ) , the second inequality requires that
1 x σ 2 + x σ 2 2 x σ 2 + 2 α σ 2 < 1 .
Rearranging, this is equivalent to
σ 2 + 2 α σ 2 < x + x σ 2 .
Isolating the x term yields
( 2 α 1 ) σ 2 1 + σ 2 < x
as claimed. □
Example 2. 
Let α = 3 4 and σ = 1 2 . Then, x = 3 2 + 2 1 + 3 β . From Theorem 5, x is stable if and only if 1 10 < x . Solving 1 10 < 3 2 + 2 1 + 3 β yields β > 65 , so that the fixed point loses stability at β = 65 . Indeed, at this value, we have x = 1 10 and
J 1 10 , 1 10 = 9 / 10 13 / 10 1 / 4 3 / 4 .
The characteristic polynomial is given by det J 1 / 10 , 1 / 10 λ = 0 which yields
λ 2 33 20 λ + 1 = 0 .
The solutions to this quadratic yields the eigenvalues λ = 33 40 ± i 511 40 . A calculation shows that | λ | = 1 for both eigenvalues so that the fixed point loses stability as β passes through 65.
Corollary 2. 
For σ ( 0 , 1 ) and 0 < α 1 2 , the fixed point ( x , x ) is locally asymptotically stable for all β > 0 .
Proof. 
The bound
σ 2 ( 2 α 1 ) 1 + σ 2 < x
is trivially satisfied since x > 0 and the left-hand side cannot be positive. □
Corollary 3. 
Let α ( 1 2 , 1 ) and σ ( 0 , 1 ) . Then, the fixed point ( x , x ) is locally asymptotically stable if and only if
β < 1 4 α 2 α + σ 2 σ 2 ( 2 α 1 ) 2 1 4 α .
Proof. 
The fixed point is locally asymptotically stable if and only if
σ 2 ( 2 α 1 ) 1 + σ 2 < 2 α 1 + 1 + 4 α β .
Since α > 1 2 , the given inequality is equivalent to
1 + 4 α β < 2 α + 2 α σ 2 σ 2 ( 2 α 1 ) 1 = 2 α + σ 2 σ 2 ( 2 α 1 ) .
Squaring both sides and isolating β gives the desired inequality. □
For fixed α ( 1 2 , 1 ) and σ ( 0 , 1 ) , the mapping undergoes a bifurcation at the critical value β = 1 4 α 2 α + σ 2 σ 2 ( 2 α 1 ) 2 1 4 α . We recall that for smooth maps varying with a parameter, there are three types of bifurcations at a fixed point, depending on how the fixed point loses stability as the parameter changes. A fold bifurcation occurs when one of the eigenvalues of the Jacobian matrix becomes 1, while a period doubling bifurcation occurs when one of the eigenvalues of the Jacobian matrix becomes 1 . A third type, Neimark–Sacker bifurcation, occurs when a pair of complex eigenvalues both have magnitude 1, such as in Example 2. For a comprehensive review of bifurcations, see [43]. We will prove momentarily that in Model (3), the fixed point ( x , x ) has a Neimark–Sacker bifurcation at β = 1 4 α 2 α + σ 2 σ 2 ( 2 α 1 ) 2 1 4 α , if α ( 1 2 , 1 ) and σ ( 0 , 1 ) , but does not undergo a fold (saddle node) bifurcation or flip (period doubling) bifurcation.
Proposition 2. 
For Model (3), let α ( 0 , 1 ) and σ ( 0 , 1 ) . Then, the fixed point ( x , x ) does not undergo a flip or fold bifurcation at any β > 0 .
Proof. 
Fix α ( 0 , 1 ) and σ ( 0 , 1 ) . Let x ( β ) = 2 α 1 + 1 + 4 α β . For brevity, let T ( β ) = Tr ( J ( x , x ) ) and D ( β ) = det ( J ( x , x ) ) .
The characteristic equation of J ( x , x ) is given by λ 2 T ( β ) λ + D ( β ) = 0 . Thus, the eigenvalues λ ± ( β ) of J ( x , x ) are given by
λ ± ( β ) = T ( β ) 2 ± T ( β ) 2 4 D ( β ) .
To have a fold bifurcation, we would need T ( β ) 2 4 D ( β ) to be real and for
T ( β ) 2 4 D ( β ) = T ( β ) 2 1 2
or equivalently for
T ( β ) = D ( β ) + 1 .
Likewise, for a flip bifurcation we would need T ( β ) 2 4 D ( β ) to be real and for
T ( β ) 2 4 D ( β ) = T ( β ) 2 + 1 2
or equivalently for
T ( β ) = D ( β ) + 1 .
However, in the proof of Theorem 5, we proved that | T ( β ) | < D ( β ) + 1 , so we cannot have a fold or flip bifurcation. □
Lemma 3. 
For Model (3), let α ( 1 2 , 1 ) and σ ( 0 , 1 ) . Let β = 1 4 α 2 α + σ 2 σ 2 ( 2 α 1 ) 2 1 4 α and let x ( β ) , x ( β ) denote the interior fixed point. Then, the eigenvalues λ ± of J ( x , x ) satisfy | λ ± | = 1 .
Proof. 
Fix α ( 1 2 , 1 ) and σ ( 0 , 1 ) . Let β = 1 4 α 2 α + σ 2 σ 2 ( 2 α 1 ) 2 1 4 α and let x ( β ) = 2 α 1 + 1 + 4 α β . As in the proof of Proposition 2, let T ( β ) = Tr ( J ( x , x ) and D ( β ) = det ( J ( x , x ) . The eigenvalues λ ± ( β ) of J ( x , x ) are given by
λ ± ( β ) = T ( β ) 2 ± T ( β ) 2 4 D ( β ) .
Since x < α < 1 , we have 0 < T ( β ) < 2 , and consequently we have T ( β ) 2 4 < 1 . We note D ( β ) = 1 , so in a neighborhood of β the eigenvalues λ ± ( β ) are complex and are given by
λ ± ( β ) = T ( β ) 2 ± i D ( β ) T ( β ) 2 4
with | λ ± ( β ) | = T ( β ) 2 4 + D ( β ) T ( β ) 2 4 = D ( β ) . Thus, both eigenvalues have magnitude 1 at β = β . □
As β passes through the bifurcation value β , the fixed point loses stability, and Lemma 3 shows the possibility of a Neimark–Sacker bifurcation. One of the main features of Neimark–Sacker bifurcations is the introduction of closed invariant curves. In a supercritical Neimark–Sacker bifurcation, the fixed point loses stability and stable invariant curves appear, while in a subcritical Neimark–Sacker bifurcation, there are unstable closed invariant curves that vanish at the same time the fixed point loses stability. Whether such a bifurcation is supercritical or subcritical can be determined by converting the map to the Neimark–Sacker normal form. However, this normal form changes depending on whether two nondegeneracy conditions are satisfied. The main theorem for the generic case is given below.
Theorem 6 
([43]). Suppose a two-dimensional discrete-time system
x F ( x , β ) , x R 2 , β R
with smooth F has a fixed point x ( β ) , where D F ( x ) has eigenvalues
λ ± ( β ) = r ( β ) e ± i θ ( β ) ,
where r ( β ) = 1 and θ ( β ) = θ 0 . Let the following conditions be satisfied:
  • (C.1) r ( β ) 0
  • (C.2) e i k θ 0 1 for k = 1 , 2 , 3 , 4
Then, by the introduction of a complex variable z and a parameter transformation γ = r ( β ) 1 , the mapping x F ( x , β ) can be transformed into the Poincare normal form
z e i θ ( γ ) ( 1 + γ + d ( γ ) + i b ( γ ) | z | 2 ) z + O ( | z | 4 )
with parameter γ. If d ( 0 ) < 0 , then the bifurcation is supercritical, and if d ( 0 ) > 0 , then the bifurcation is subcritical.
Lemma 4. 
For Model (3), conditions (C.1) and (C.2) in Theorem 6 are satisfied for α ( 1 2 , 1 ) and σ ( 0 , 1 ) .
Proof. 
Fix α ( 1 2 , 1 ) and σ ( 0 , 1 ) , and let x ( β ) = 2 α 1 + 1 + 4 α β denote the fixed point and β = 1 4 α 2 α + σ 2 σ 2 ( 2 α 1 ) 2 1 4 α denote the bifurcation value. Note that ( x ) ( β ) < 0 and that
r ( β ) = | λ ± ( β ) | = D ( β ) = ( 1 x ) ( 1 σ 2 ) 2 x σ 2 + 2 α σ 2 .
It then follows that
r ( β ) = ( x ) ( β ) ( σ 2 3 ) > 0 .
In particular, r ( β ) > 0 so that condition (C.1) is satisfied.
Since
e i k θ 0 = λ + k ( β ) = T ( β ) 2 ± i 1 T ( β ) 2 4 k ,
we only have e i k θ 0 = 1 for k { 1 , 2 , 3 , 4 } if T ( β ) { 2 , 1 , 0 , 2 } , which is not the case as T ( β ) ( 0 , 2 ) . Thus, condition (C.2) is satisfied. □
Lemmas 3 and 4 are sufficient to establish the following proposition.
Proposition 3. 
Let α ( 1 2 , 1 ) and σ ( 0 , 1 ) . Then, the fixed point ( x , x ) undergoes a Neimark–Sacker bifurcation at β = 1 4 α 2 α + σ 2 σ 2 ( 2 α 1 ) 2 1 4 α .
A prototype of the stability region in the parameter space is shown in Figure 5. For α , σ , and β within the shaded region, the interior fixed point ( x , x ) is locally asymptotically stable, corresponding to the region referenced in the conditions of Theorem 5. The conditions in Corollary 2 define the region on the left side of the concave portion such that α 1 2 , while the conditions in Corollary 3 define the region below the concave portion such that α > 1 2 . The Neimark–Sacker bifurcation mentioned in Proposition 3 occurs precisely on the surface of the concave part.
While we have established the existence of a Neimark–Sacker bifurcation, the computations thus far have not yet shown if it is supercritical or subcritical. In principle, one can calculate the coefficient d ( 0 ) from Theorem 6, but the calculation for Model (3) is lengthy and impractical, as it does not allow one to determine the sign of d ( 0 ) for general parameter values (see Appendix A), due to the size of the coefficients involved. Instead, we will present a geometric argument using stable and unstable manifolds in the next section.

5. Stable and Unstable Manifolds

For a fixed point p of an invertible mapping F : X X , we denote the stable manifold of p by
W s ( p ) = { x X : lim n F n ( x ) = p }
and likewise its unstable manifold by
W s ( p ) = { x X : lim n F n ( x ) = p } .
For the mapping F U , the stable manifold of the origin is the portion of the y-axis that lies in U, that is, W s ( 0 , 0 ) = { ( x , y ) U : x = 0 } .
Some numerically computed images of the unstable manifold at the origin are shown in Figure 6. We will now proceed to show these plots are typical.
In J ( 0 , 0 ) , the eigenvector σ 2 + ( e α 1 ) σ 2 corresponds to the unstable direction in the linearized map with eigenvalue e α . Since W u ( 0 , 0 ) lies tangent to this eigenvector, W u ( 0 , 0 ) is initially in R 1 . We also recall the cycling nature of the R i , the direction of the mapping in each R i , and that the manifold cannot self-intersect. These facts ensure the argument of W u ( 0 , 0 ) is strictly increasing when viewed as centered at ( x , x ) .
If the curve remains in R 1 , then it must monotonically approach ( x , x ) . Otherwise, it must eventually intersect I 2 so that it enters R 2 . Again, if it remains in R 2 , then it must monotonically approach ( x , x ) , otherwise it must intersect I 1 . Eventually, if it does not stay in some R i , it will intersect the isoclines an infinite number of times as it spirals counterclockwise. Viewing ( x , x ) as the center and looking at polar coordinates ( r , θ ) , we can look at the return map that gives the radius r from ( x , x ) generated by intersections with any fixed θ value. Since the manifold does not self intersect, these return maps are strictly decreasing. As they are also bounded, we conclude the following theorem.
Theorem 7. 
Let α ( 0 , 1 ) , β > 0 , and σ ( 0 , 1 ) . Then, the unstable manifold W u ( 0 , 0 ) is asymptotic to either the fixed point ( x , x ) or a closed curve that encloses ( x , x ) .
In Figure 6, we illustrate the concept underlying Theorem 7. As is clearly seen in the presented four examples, all the orbits behave in a counterclockwise manner, either converging to the fixed point when the absolute values of both eigenvalues are less than one or orbiting around the fixed point when both eigenvalues exceed absolute value 1. This is precisely the main idea behind the proof of Theorem 7.
Let us now fix α and σ and see how the unstable manifold evolves as β increases. First, we prove the following lemma.
Lemma 5. 
Let α , σ ( 0 , 1 ) . Let p ( β ) denote the discriminant of the characteristic polynomial of J x , x . Then p ( β ) < 0 .
Proof. 
Recall that we have p ( β ) = T ( β ) 2 4 D ( β ) , where
D ( β ) = det ( J ( x , x ) = 1 x 2 x 2 α σ 2 1 σ 2 = 1 x σ 2 + σ 2 x 2 x σ 2 + 2 α σ 2 = ( 1 σ 2 + 2 σ 2 α ) + x ( σ 2 1 ) ,
T ( β ) = Trace ( J ( x , x ) ) = 2 x σ 2
and
x = 2 α 1 + 1 + 4 α β .
Note also that x β < 0 .
Then,
p ( β ) = 1 2 T ( β ) x β x β ( σ 2 1 ) = x β 1 + 1 2 x + 1 2 σ 2 + σ 2 + 1 = x β 1 2 x + 3 2 σ 2 < 0 .
Initially, when β = 0 , it is straightforward to prove the set R 1 is forward invariant, so that the unstable manifold W u ( 0 , 0 ) forms a heteroclinic connection between the origin and ( x , x ) . In this case, x = α , and the eigenvalues of J ( x , x ) are 1 α and 1 σ 2 , i.e., p ( 0 ) > 0 . The heteroclinic connection for β = 0 persists under perturbations of β and is tangent to one of the eigenvectors at ( x , x ) . This idea is illustrated in Figure 7. In this example, the tangency is clearly observed. Let us examine some properties of these eigenvectors in more detail.
Lemma 6. 
Let α , σ ( 0 , 1 ) . Let β be such that the eigenvalues of J ( x , x ) are real. Then, any eigenvalue λ of J ( x , x ) satisfies 0 < λ < 1 σ 2 .
Proof. 
Since λ is an eigenvalue of J ( x , x ) , then
1 x λ 2 x 2 α σ 2 1 σ 2 λ
has vanishing determinant. Since σ 2 > 0 and 2 x 2 α < 0 , the product
( 1 x λ ) ( 1 σ 2 λ )
must be negative. Since x < α < 1 , if λ 0 , then both terms are positive. If λ 1 σ 2 , then both terms are negative or zero. This implies 0 < λ < 1 σ 2 . □
If λ is a real eigenvalue of J ( x , x ) , then the line through ( x , x ) that is parallel to the corresponding eigenvector has a slope of σ 2 1 σ 2 λ . From the previous lemma, this slope is always negative. The larger of the two eigenvalues, given by λ + = T ( β ) 2 + T ( β ) 2 4 D ( β ) , corresponds to the line with a steeper slope. It is this line that W u ( 0 , 0 ) is asymptotic to.
A calculation shows
lim β p ( β ) = ( 2 σ 2 ) 2 4 1 + σ 2 2 σ 2 α = σ 2 4 ( σ 2 8 α ) .
If α σ 2 8 , then by continuity the eigenvalues of J ( x , x ) remain real for all β . For α > σ 2 8 , eigenvalues transition at a critical value β ¯ from real to complex. At such values, the curve W u ( 0 , 0 ) transitions from directly looping into ( x , x ) and instead starts an infinite spiral, intersecting the isoclines an infinite number of times.
If W u ( 0 , 0 ) is asymptotic to a closed curve, then that curve is necessarily stable. Consequently, the curve cannot be an unstable curve that emerged from a subcritical Neimark–Sacker bifurcation. For β ¯ < β < β , this spiral approaches ( x , x ) and the heteroclinic orbit is maintained. At β , the heteroclinic orbit breaks as a supercritical Neimark–Sacker bifurcation occurs, and a stable closed curve bifurcates from the fixed point for β > β .
Example 3. 
Let α = 3 4 and σ = 1 2 , so that x = 3 2 + 2 1 + 3 β . We have
T ( β ) = 7 4 x
and
D ( β ) = 9 8 5 4 x .
The eigenvalues of J ( x , x ) are real provided
T ( β ) 2 4 D ( β ) > 0
which yields
1 4 ( x ) 2 + 3 8 x 23 64 > 0
which for positive x implies that the eigenvalues are real if x > 2 3 4 . Solving
3 2 + 2 1 + 3 β > 2 3 4
for positive β yields
β < 216 529 80 529 2 .
Thus, β ¯ = 216 529 80 529 2 . Since β = 65 (see Example 2), the dynamics for α = 3 4 and σ = 1 2 may be described as follows. For 0 < β < 216 529 80 529 2 , the unstable manifold W u ( 0 , 0 ) is tangent to the real eigenvector and convergence to the fixed point is guaranteed. For 216 529 80 529 2 < β < 65 , the eigenvalues are now complex, and the unstable manifold infinitely spirals around the fixed point, but still converges to it. For β > 65 , where the Neimark–Sacker bifurcation occurs, a new forward invariant closed curve appears and the unstable manifold spirals towards this curve instead. See also Figure 6.
In all cases, since W u ( 0 , 0 ) forms a barrier which points cannot cross, this unstable manifold dictates global behavior on F U . In particular, in the absence of a closed invariant curve, the fixed point is globally asymptotically stable.

6. Conclusions

In this paper, we studied the local properties, the global properties, and the bifurcation of a nonlinear planar system of coupled difference equations consisting of the evolution of a species of Ricker type. By studying the model when the trait equation is coupled from the population dynamic equation, we were able to conclude that, under certain restrictions on the parameters, the model is forward invariant and globally attracting in a subset of the phase space, bounded by the coordinate axes, the horizontal line y = 1 , and a critical curve L C , which is the image under the mapping that rules the model where the Jacobian vanishes. We show that in this subset of the phase space, there are four regions where the orbits behave in a cyclic nature (counterclockwise direction) and converge towards the asymmetric unstable manifold of the origin, which will attract all the orbits.
Moreover, in this set, the model possesses two symmetric fixed points: the origin, which is always a saddle fixed point, and an interior fixed point that is locally asymptotically stable for certain values of the parameters. This non-trivial fixed point undergoes a Neimark–Sacker bifurcation for a critical parameter value, a phenomenon that does not occur in the original one-dimensional Ricker model. By using a geometric argument with the unstable manifold at the origin, we show that this bifurcation is supercritical, and thus a stable asymmetric closed curve bifurcates from the fixed point once we cross the critical value. Finally, in the absence of the referred closed invariant curve, the symmetric non-trivial fixed point is globally asymptotically stable.
We have utilized tools from both analytic and geometric approaches to achieve our goals. The analytic tools are similar to those existing in the literature for these types of models, including the relative position of critical curves using singularity theory, the dynamics of sets around the isoclines, and local stability using Jury’s conditions.
Typically, the study and classification of the Neimark–Sacker bifurcation are conducted using the normal forms of the equations. Nondegeneracy conditions are usually computed by taking certain derivatives, as explained in Section 4. However, in Model (3), we were unable to use d ( 0 ) to classify the bifurcation, as is mentioned in the end of Section 4. Therefore, we employed a geometric argument, which, to the best of our knowledge, is the first time it has been used in this context. This geometric argument allowed us to demonstrate the global stability of the positive fixed point, meaning that every orbit of a point in the interior of the first quadrant will eventually converge to the fixed point. This is a central conclusion in real-world applications, particularly in ecology.
We should mention that the tools employed in this model are extendable to other discrete evolutionary models as well, such as the evolutionary Beverton–Holt model and the logistic model.

7. Future Work

In future work, we plan to study the dynamics of the model in seasonal environments, specifically focusing on how the discrete model behaves under periodic changes. This involves analyzing how the system’s behavior is influenced by parameters that fluctuate regularly, reflecting seasonal variations in the environment. In particular, we aim to understand how these periodic parameter changes affect the long-term dynamics of the model, i.e., the study of the model
x n + 1 = x n e α n β y n 2 x n y n + 1 = ( 1 σ 2 ) y n + σ 2 x n ,
when the sequence of parameters α n is p -periodic, i.e., α n + p = α n for all n and some integer p > 1 .
We believe that the orbits in the periodic model (4) will exhibit similar behavior, with counterclockwise dynamics around the periodic cycle in phase space. Furthermore, since it will be challenging to explicitly write the periodic solution, the geometric method developed for classifying the Neimark–Sacker bifurcation will be a useful tool. As observed in simulations, if the period of the cycle equals the period of the system, the same number of cycling orbits circulate counterclockwise around the periodic points. Moreover, when the eigenvalues of the Jacobian of the composition map, evaluated at the periodic orbit, have absolute values less than one, convergence can be observed. Thus, the geometric approach will be a valuable method for classifying the nature of the invariant manifolds and the dynamics around the isoclines.
Understanding these dynamics can provide insights into the long-term behavior and global stability of the population under study. In particular, it will be interesting to determine whether the periodic evolutionary model exhibits attenuation, resonance, or neither attenuation nor resonance, a phenomenon observed by R. Sacker [44] in the periodic one-dimensional Ricker model.
We also plan to use similar techniques to study a generalization of Model (3) for n competing species, resulting in a model with 2 n equations. Our real challenge will be generalizing the geometric approach to higher dimensions, as it may become impossible to explicitly compute the fixed points after a certain number of species. Hence, a geometric approach will be crucial.
An interesting study would also investigate whether our tools are applicable to a model similar to Model (3) but with several traits in the canonical equation for evolution.

Author Contributions

Conceptualization, R.L. and B.R.; methodology, B.R.; software, R.L. and B.R.; validation, R.L. and B.R.; formal analysis, R.L. and B.R.; investigation, R.L. and B.R.; Writing—original draft preparation, review, and editing: R.L. and B.R. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Appendix A

The following theorem, which summarizes several theorems found in [43], gives a way to calculate the coefficient d ( 0 ) in Theorem 6.
Theorem A1 
([43]). Let
x F ( x , β ) , x R 2 , β R
be as in Theorem 6. Let q and p be eigenvectors of D F ( x ( β ) ) and D F T ( x ( β ) ) at λ + ( β ) and λ ( β ) , respectively, normalized so that p , q = 1 . Let
H ( z , z ¯ ) = p , F ( x + z q + z ¯ q ¯ , β ) x .
Then, the Taylor expansion at ( z , z ¯ ) = ( 0 , 0 ) is of the form
e i θ 0 z + 2 j + k 3 1 j ! k ! g j k z j z ¯ k + O ( | z | 4 )
and
d ( 0 ) = R e g 20 g 11 ( 1 2 e i θ 0 ) 2 ( e 3 i θ 0 e 2 i θ 0 ) + | g 11 | 2 e i θ 0 1 + | g 02 | 2 2 ( e 3 i θ 0 1 ) + g 21 2 e i θ 0 .
We now show the calculation for Model (3). The complex eigenvectors of J ( x , x ) and J T ( x , x ) at λ ± ( β ) = e ± i θ 0 are
q = q 1 q 2 = c 1 λ + ( β ) + σ 2 1 σ 2
and
p = p 1 p 2 = c 2 σ 2 λ ( β ) + x 1 = c 2 σ 2 ( λ + ( β ) + σ 2 1 )
respectively, for some c 1 , c 2 C . The last equality in the p vectors follows from the trace yielding λ + ( β ) + λ ( β ) = 2 x σ 2 .
Using the standard complex inner product and the fact that λ ( β ) ¯ = λ + ( β ) , we obtain
p , q = p 1 q 1 ¯ + p 2 q 2 ¯ = c 1 ¯ c 2 σ 2 λ ( β ) λ + ( β ) = c 1 ¯ c 2 σ 2 i 4 ( 2 x σ 2 ) 2 .
Thus, we may take any c 1 , c 2 such that c 1 ¯ c 2 = i σ 2 4 ( 2 x σ 2 ) 2 .
Let
A ( z , z ¯ ) = p 1 ( x + z q 1 + z ¯ q 1 ¯ ) ,
B ( z , z ¯ ) = β q ¯ 2 z ¯ + q 2 z + x 2 q ¯ 1 z ¯ + α q 1 z x
and
C ( z , z ¯ ) = p 2 σ 2 q ¯ 1 z ¯ q ¯ 2 z ¯ + q 1 z q 2 z + q ¯ 2 z ¯ + q 2 z + x .
We calculate
H ( z , z ¯ ) = p 1 p 2 , F x + z q 1 + z ¯ q 1 ¯ x + z q 2 + z ¯ q 2 ¯ x x = p 1 p 2 , ( x + z q 1 + z ¯ q 1 ¯ ) e α x z q 1 z ¯ q 1 ¯ β ( x + z q 2 + z ¯ q 2 ¯ ) 2 x σ 2 q ¯ 1 z ¯ + q 1 z + x + 1 σ 2 q ¯ 2 z ¯ + q 2 z + x = A ( z , z ¯ ) e B ( z , z ¯ ) + C ( z , z ¯ ) .
Observe that
A ( 0 , 0 ) = p 1 x A z ( 0 , 0 ) = p 1 q 1 , A z ¯ ( 0 , 0 ) = p 1 q ¯ 1 A z z ( 0 , 0 ) = A z z ¯ ( 0 , 0 ) = A z ¯ z ¯ ( 0 , 0 ) = 0 B ( 0 , 0 ) = 0 B z ( 0 , 0 ) = q 1 2 β x q 2 , B z ¯ ( 0 , 0 ) = q ¯ 1 2 β x q ¯ 2 B z z ( 0 , 0 ) = 2 β q 2 2 , B z z ¯ ( 0 , 0 ) = 2 β q 2 q ¯ 2 , B z ¯ z ¯ ( 0 , 0 ) = 2 β q ¯ 2 2 C ( 0 , 0 ) = p 2 x C z ( 0 , 0 ) = p 2 q 2 + p 2 q 1 σ 2 p 2 q 2 σ 2 , C z ¯ ( 0 , 0 ) = p 2 q ¯ 2 + p 2 q ¯ 1 σ 2 p 2 q ¯ 2 σ 2 C z z ( 0 , 0 ) = C z z ¯ ( 0 , 0 ) = C z ¯ z ¯ ( 0 , 0 ) = 0 .
Then,
g 20 = H z z ( 0 , 0 ) = ( A z z + 2 A z B z + A B z z + A B z 2 ) e B + C z z | ( z , z ¯ ) = ( 0 , 0 ) = p 1 2 β q 2 2 x + x 2 β q 2 x + q 1 2 2 q 1 2 β q 2 x + q 1 ,
g 11 = H z z ¯ ( 0 , 0 ) = ( A z z ¯ + A z B z ¯ + A z ¯ B z + A B z z ¯ + A B z B z ¯ ) e B + C z z ¯ | ( z , z ¯ ) = ( 0 , 0 ) = p 1 q 1 2 β x 1 x q ¯ 2 + x 2 q ¯ 1 + 2 β q 2 x q ¯ 2 2 β x 2 1 + x 1 q ¯ 1 ,
g 02 = H z ¯ z ¯ ( 0 , 0 ) = ( A z ¯ z ¯ + 2 A z ¯ B z ¯ + A B z ¯ z ¯ + A B z ¯ 2 ) e B + C z ¯ z ¯ | ( z , z ¯ ) = ( 0 , 0 ) = p 1 4 β x 1 x q ¯ 2 q ¯ 1 + 2 β x q ¯ 2 2 2 β x 2 1 + x 2 q ¯ 1 2
and
g 21 = H z z z ¯ ( 0 , 0 ) = z ¯ ( A z z + 2 A z B z + A B z z + A B z 2 ) e B + C z z | ( z , z ¯ ) = ( 0 , 0 ) .
Since A z z z ¯ = B z z z ¯ = C z z z ¯ = 0 it follows
g 21 = 2 A z B z z ¯ + A z ¯ B z z + A z ( B z z ) 2 + 2 A B z z ¯ B z e B + B z g 20 = p 1 2 β q 1 q 2 2 x 1 q ¯ 2 + x q 1 4 3 x + q 2 8 β 3 q 2 3 x 4 q 1 3 x 2 + p 1 2 β q 2 2 q ¯ 1 + 4 β 2 q 2 2 x 2 2 q ¯ 2 + q 2 + q 1 q 2 2 + x 2 2 3 x .
For specific α , σ , these values allow us to calculate the coefficient d ( 0 ) . As an example, let α = 3 4 and σ = 1 2 . Then, β = 65 and x ( 65 ) = 1 10 . Then,
J ( x , x ) = 9 / 10 13 / 10 1 / 4 3 / 4
and λ ± = 33 40 ± i 511 40 . We take as eigenvectors of J ( x , x ) and J T ( x , x )
q = c 1 3 40 + i 511 40 1 4
and
p = c 2 1 4 3 40 i 511 40
respectively. For the normalization, we must choose c 1 , c 2 so that c 1 ¯ c 2 = 80 i 511 . Taking c 1 = 1 and c 2 = 80 i 511 , we have
g 20 = 1227 400 + 3209 i 400 511
g 11 = 65 i 4 511
g 02 = 1227 400 + 3209 i 400 511
g 21 = 636107 16000 + 1674581 i 16000 511
from which we calculate that
d ( 0 ) = 1382119 89600 .
As this value is negative, the bifurcation is supercritical.

References

  1. Ricker, E. Stock and recruitment. J. Fish. Res. Board Can. 1954, 11, 559–623. [Google Scholar] [CrossRef]
  2. May, R.M.; Oster, G.F. Bifurcations and Dynamic Complexity in Simple Ecological Models. Am. Nat. 1976, 110, 573–599. [Google Scholar] [CrossRef]
  3. Costantino, R.F.; Desharnais, R.A.; Cushing, J.M.; Dennis, B. Chaotic Dynamics in an Insect Population. Science 1997, 275, 389–391. [Google Scholar] [CrossRef]
  4. Zimmer, C. Life After Chaos. Science 1999, 284, 83–86. [Google Scholar] [CrossRef]
  5. Ackleh, A.S.; Cushing, J.M.; Salceanu, P.L. On the dynamics of evolutionary competition models. Nat. Resour. Model. 2015, 28, 380–397. [Google Scholar] [CrossRef]
  6. Cushing, J.M. An evolutionary Beverton–Holt model. In Theory and Applications of Difference Equations and Discrete Dynamical Systems; Springer: Berlin/Heidelberg, Germany, 2014; Volume 102, pp. 127–141. [Google Scholar] [CrossRef]
  7. Cushing, J.M. Difference equations as models of evolutionary population dynamics. J. Biol. Dyn. 2019, 13, 103–127. [Google Scholar] [CrossRef]
  8. Cushing, J.M. A Darwinian Ricker equation. In Progress on Difference Equations and Discrete Dynamical Systems; Springer: Cham, Switzerland, 2020; Volume 341, pp. 231–243. [Google Scholar] [CrossRef]
  9. Cushing, J.M.; Stefanko, K. A Darwinian dynamic model for the evolution of post-reproduction survival. J. Biol. Syst. 2021, 29, 433–450. [Google Scholar] [CrossRef]
  10. Cushing, J.M. The evolutionary dynamics of a population model with a strong Allee effect. Math. Biosci. Eng. 2015, 12, 643–660. [Google Scholar] [CrossRef]
  11. Rael, R.C.; Vincent, T.L.; Cushing, J.M. Competitive outcomes changed by evolution. J. Biol. Dyn. 2011, 5, 227–252. [Google Scholar] [CrossRef]
  12. Darwin, C. The Origin of Species; Avenel Books: London, UK, 1859. [Google Scholar]
  13. Weibull, J.W. Evolutionary Game Theory; MIT Press: Cambridge, MA, USA, 1995. [Google Scholar]
  14. Vincent, T.L.; Brown, J.S. Evolutionary Game Theory, Natural Selection, and Darwinian Dynamics; Cambridge University Press: Cambridge, UK, 2005. [Google Scholar]
  15. Ch-Chaoui, M.; Mokni, K. A discrete evolutionary Beverton–Holt population model. Int. J. Dyn. Control 2023, 11, 1060–1075. [Google Scholar] [CrossRef]
  16. Mokni, K.; Ch-Chaoui, M. Asymptotic stability, bifurcation analysis and chaos control in a discrete evolutionary Ricker population model with immigration. In Advances in Discrete Dynamical Systems, Difference Equations and Applications; Springer: Cham, Switzerland, 2023; Volume 416, pp. 363–403. [Google Scholar] [CrossRef]
  17. Elaydi, S.; Kang, Y.; Luís, R. The effects of evolution on the stability of competing species. J. Biol. Dyn. 2022, 16, 816–839. [Google Scholar] [CrossRef] [PubMed]
  18. Mokni, K.; Elaydi, S.; Ch-Chaoui, M.; Eladdadi, A. Discrete evolutionary population models: A new approach. J. Biol. Dyn. 2020, 14, 454–478. [Google Scholar] [CrossRef]
  19. Elaydi, S.; Kang, Y.; Luís, R. Global asymptotic stability of evolutionary periodic Ricker competition models. J. Differ. Equ. Appl. 2024, 30, 1091–1116. [Google Scholar] [CrossRef]
  20. Rosenkranz, G. On global stability of discrete population models. Math. Biosci. 1983, 64, 227–231. [Google Scholar] [CrossRef]
  21. Cull, P. Stability in one-dimensional models. Sci. Math. Jpn. 2003, 58, 349–357. [Google Scholar]
  22. Singer, D. Stable orbits and bifurcation of maps of the interval. SIAM J. Appl. Math. 1978, 35, 260–267. [Google Scholar] [CrossRef]
  23. Smith, H.L. Planar competitive and cooperative difference equations. J. Differ. Equ. Appl. 1998, 3, 335–357. [Google Scholar] [CrossRef]
  24. Hirsch, M.W.; Smith, H. Monotone dynamical systems. In Handbook of Differential Equations: Ordinary Differential Equations; Elsevier: Amsterdam, The Netherlands, 2005; Volume II, pp. 239–357. [Google Scholar]
  25. Kulenović, M.R.S.; Moranjkić, S.; Nurkanović, M.; Nurkanović, Z. Global asymptotic stability and Naimark-Sacker bifurcation of certain mix monotone difference equation. Discret. Dyn. Nat. Soc. 2018, 2018, 7052935. [Google Scholar] [CrossRef]
  26. Balreira, E.C.; Elaydi, S.; Luís, R. Global stability of higher dimensional monotone maps. J. Differ. Equ. Appl. 2017, 23, 2037–2071. [Google Scholar] [CrossRef]
  27. Smith, H.L. Global stability for mixed monotone systems. J. Differ. Equ. Appl. 2008, 14, 1159–1164. [Google Scholar] [CrossRef]
  28. Whitney, H. On singularities of mappings of euclidean spaces. I. Mappings of the plane into the plane. Ann. Math. 1955, 62, 374–410. [Google Scholar] [CrossRef]
  29. Mira, C.; Gardini, L.; Barugola, A.; Cathala, J. Chaotic Dynamics in Two-Dimensional Noninvertible Maps; Nonlinear Sciences Series A; World Scientific: Singapore, 1996. [Google Scholar]
  30. Mira, C. Chaotic Dynamics; World Scientific: Singapore, 1987. [Google Scholar]
  31. Gardini, L. Some global bifurcations of two-dimensional endomorphisms by use of critical lines. Nonlinear Anal. 1992, 18, 361–399. [Google Scholar] [CrossRef]
  32. Cabral Balreira, E.; Elaydi, S.; Luís, R. Local stability implies global stability for the planar Ricker competition model. Discret. Contin. Dyn. Syst. Ser. B 2014, 19, 323–351. [Google Scholar] [CrossRef]
  33. Ryals, B.; Sacker, R.J. Global stability in the 2D Ricker equation. J. Differ. Equ. Appl. 2015, 21, 1068–1081. [Google Scholar] [CrossRef]
  34. Ryals, B.; Sacker, R.J. Global stability in the 2D Ricker equation revisited. Discret. Contin. Dyn. Syst. Ser. B 2017, 22, 585–604. [Google Scholar] [CrossRef]
  35. Ryals, B. A note on a parameter bound for global stability in the 2D coupled Ricker equation. J. Differ. Equ. Appl. 2018, 24, 240–244. [Google Scholar] [CrossRef]
  36. Elaydi, S. Global dynamics of discrete dynamical systems and difference equations. In Difference Equations, Discrete Dynamical Systems And Applications; Springer: Cham, Switzerland, 2019; Volume 287, pp. 51–81. [Google Scholar] [CrossRef]
  37. Luís, R. Open Problems and Conjectures in the Evolutionary Periodic Ricker Competition Model. Axioms 2024, 13, 246. [Google Scholar] [CrossRef]
  38. Dercole, F.; Rinaldi, S. Analysis of Evolutionary Processes: The Adaptive Dynamics Approach and Its Applications; Princeton Series in Theoretical and Computational Biology; Princeton University Press: Princeton, NJ, USA, 2008. [Google Scholar]
  39. Lande, R. Natural selection and random genetic drift in phenotypic evolution. Evolution 1976, 30, 314–334. [Google Scholar] [CrossRef]
  40. Lande, R. A Quantitative Genetic Theory of Life History Evolution. Ecology 1982, 63, 607–615. [Google Scholar] [CrossRef]
  41. Abrams. Modelling the adaptive dynamics of traits involved in inter- and intraspecific interactions: An assessment of three methods. Ecol. Lett. 2001, 4, 166–175. [Google Scholar] [CrossRef]
  42. Jury, E. On the roots of a real polynomial inside the unit circle and a stability criterion for linear discrete systems. IFAC Proc. Vol. 1963, 1, 142–153. [Google Scholar] [CrossRef]
  43. Kuznetsov, Y.A. Elements of applied bifurcation theory. In Applied Mathematical Sciences, 4th ed.; Springer: Cham, Switzerland, 2023; Volume 112, pp. xxvi+703. [Google Scholar] [CrossRef]
  44. Sacker, R.J. A note on periodic Ricker maps. J. Differ. Equ. Appl. 2007, 13, 89–92. [Google Scholar] [CrossRef]
Figure 1. The graph of L C separates the plane into two distinct regions. All points in the interior of the shaded region have two distinct preimages. Points in the unshaded region have no preimage.
Figure 1. The graph of L C separates the plane into two distinct regions. All points in the interior of the shaded region have two distinct preimages. Points in the unshaded region have no preimage.
Symmetry 16 01139 g001
Figure 2. Graphs of L C , L C 1 = F 1 ( L C ) , and a line L with slope σ 2 1 σ 2 and its image under F are shown. The image of the line is symmetric about the point on L C 1 , where each point to the lower right of L C 1 may be identified with a point to the upper left, as they have the same dynamics. To illustrate this, four points A , B , C , D are shown on F ( L ) , with their preimages under F marked on L. For instance, the two points marked F 1 ( A ) both map to A and share the same future thereafter. The values used in this graph are α = 1 2 , β = 2 , σ = 3 4 . Lemmas 1, 2, and Theorem 2 have shown that this depiction is typical for all lines with this slope for any α ( 0 , 1 ) , β > 0 , and σ 2 ( 0 , 1 ) .
Figure 2. Graphs of L C , L C 1 = F 1 ( L C ) , and a line L with slope σ 2 1 σ 2 and its image under F are shown. The image of the line is symmetric about the point on L C 1 , where each point to the lower right of L C 1 may be identified with a point to the upper left, as they have the same dynamics. To illustrate this, four points A , B , C , D are shown on F ( L ) , with their preimages under F marked on L. For instance, the two points marked F 1 ( A ) both map to A and share the same future thereafter. The values used in this graph are α = 1 2 , β = 2 , σ = 3 4 . Lemmas 1, 2, and Theorem 2 have shown that this depiction is typical for all lines with this slope for any α ( 0 , 1 ) , β > 0 , and σ 2 ( 0 , 1 ) .
Symmetry 16 01139 g002
Figure 3. The region U bounded by x = 0 , y = 0 , y = 1 , and the curve L C , where the parameters values are α = 1 2 , β = 2 , and σ = 3 4 . This region is forward invariant, globally attracting, and the restriction of F to U is injective, as proved in Theorem 3. Lemma 2 shows that the general shape of region U is typical.
Figure 3. The region U bounded by x = 0 , y = 0 , y = 1 , and the curve L C , where the parameters values are α = 1 2 , β = 2 , and σ = 3 4 . This region is forward invariant, globally attracting, and the restriction of F to U is injective, as proved in Theorem 3. Lemma 2 shows that the general shape of region U is typical.
Symmetry 16 01139 g003
Figure 4. The figure depicts a generic representation of the regions R i as determined by the isoclines and the curve L C . In the plot shown, the parameter values are α = 3 4 , β = 20 , and σ = 1 2 , though Lemma 2, and the equations for the isoclines show that this plot is typical for all α ( 0 , 1 ) , β > 0 , and σ ( 0 , 1 ) .
Figure 4. The figure depicts a generic representation of the regions R i as determined by the isoclines and the curve L C . In the plot shown, the parameter values are α = 3 4 , β = 20 , and σ = 1 2 , though Lemma 2, and the equations for the isoclines show that this plot is typical for all α ( 0 , 1 ) , β > 0 , and σ ( 0 , 1 ) .
Symmetry 16 01139 g004
Figure 5. A prototype of the stability region in the parameter space of the interior fixed point ( x , x ) of Model (3).
Figure 5. A prototype of the stability region in the parameter space of the interior fixed point ( x , x ) of Model (3).
Symmetry 16 01139 g005
Figure 6. Numerically computed the unstable manifolds of the origin for α = 3 4 and σ = 1 2 . The β values in the four plots are 0.1 , 1, 20, and 100, respectively. The fixed point ( x , x ) for α = 3 4 and σ = 1 2 is locally asymptotically stable for β < 65 . In the first three images, ( x , x ) is locally asymptotically stable and the unstable manifold W u ( 0 , 0 ) is asymptotic to the fixed point ( x , x ) . For small β , it converges directly to it, while for larger β it spirals into it. In the last image, these spirals converge to an invariant curve instead of to the unstable fixed point.
Figure 6. Numerically computed the unstable manifolds of the origin for α = 3 4 and σ = 1 2 . The β values in the four plots are 0.1 , 1, 20, and 100, respectively. The fixed point ( x , x ) for α = 3 4 and σ = 1 2 is locally asymptotically stable for β < 65 . In the first three images, ( x , x ) is locally asymptotically stable and the unstable manifold W u ( 0 , 0 ) is asymptotic to the fixed point ( x , x ) . For small β , it converges directly to it, while for larger β it spirals into it. In the last image, these spirals converge to an invariant curve instead of to the unstable fixed point.
Symmetry 16 01139 g006
Figure 7. Numerically computed unstable manifold W u ( 0 , 0 ) for α = 3 4 , β = 1 10 , and σ = 1 2 , along with lines parallel to the eigenvectors at ( x , x ) (shown as dotted lines). The isoclines are also shown.
Figure 7. Numerically computed unstable manifold W u ( 0 , 0 ) for α = 3 4 , β = 1 10 , and σ = 1 2 , along with lines parallel to the eigenvectors at ( x , x ) (shown as dotted lines). The isoclines are also shown.
Symmetry 16 01139 g007
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Luís, R.; Ryals, B. Invariant Sets, Global Dynamics, and the Neimark–Sacker Bifurcation in the Evolutionary Ricker Model. Symmetry 2024, 16, 1139. https://doi.org/10.3390/sym16091139

AMA Style

Luís R, Ryals B. Invariant Sets, Global Dynamics, and the Neimark–Sacker Bifurcation in the Evolutionary Ricker Model. Symmetry. 2024; 16(9):1139. https://doi.org/10.3390/sym16091139

Chicago/Turabian Style

Luís, Rafael, and Brian Ryals. 2024. "Invariant Sets, Global Dynamics, and the Neimark–Sacker Bifurcation in the Evolutionary Ricker Model" Symmetry 16, no. 9: 1139. https://doi.org/10.3390/sym16091139

APA Style

Luís, R., & Ryals, B. (2024). Invariant Sets, Global Dynamics, and the Neimark–Sacker Bifurcation in the Evolutionary Ricker Model. Symmetry, 16(9), 1139. https://doi.org/10.3390/sym16091139

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop