Discrete Evolutionary Population Models: A New Approach

In this paper, we apply a new approach to a special class of discrete time evolution models and establish a solid mathematical foundation to analyse them. We propose new single and multi-species evolutionary competition models using the evolutionary game theory that require a more advanced mathematical theory to handle effectively. A key feature of this new approach is to consider the discrete models as non-autonomous difference equations. Using the powerful tools and results developed in our recent work [E. D’Aniello and S. Elaydi, The structure of ω-limit sets of asymptotically non-autonomous discrete dynamical systems, Discr. Contin. Dyn. Series B. 2019 (to appear).], we embed the non-autonomous difference equations in an autonomous discrete dynamical systems in a higher dimension space, which is the product space of the phase space and the space of the functions defining the non-autonomous system. Our current approach applies to two scenarios. In the first scenario, we assume that the trait equations are decoupled from the equations of thepopulations. This requires specialized biological and ecological assumptions which we clearly state. In the second scenario, we do not assume decoupling, but rather we assume that the dynamics of the trait is known, such as approaching a positive stable equilibrium point which may apply to a much broader evolutionary dynamics. ARTICLE HISTORY Received 23 August 2019 Accepted 5 March 2020


Introduction
Population dynamics can be traced back to Malthus' 1798 essay [1], where he communicates that an unrestrained population will grow exponentially unless the amounts of food that is available becomes limited. The classical theories of population dynamics have a long history of being very useful in studying biological systems and ecological interactions and their outcomes. Determining the dynamics of a single or multi-species interactions is, however, not a straightforward undertaking when the evolution plays a crucial role in shaping the behaviour of these interactions.
The fundamental principle of evolution was first pioneered by the naturalist Charles Darwin in his 1859 book On the origin of species [13] where he communicated his reasonings verbally. Darwinian evolution theory is founded on three axioms known as the axioms of natural selection which include variation, inheritance, and competition. Variation is where individuals within a population have different traits or phenotype, whereas inheritance is where offspring inherits a mixture of both parents' traits. Competition is where more offsprings are produced than can survive, so offspring with traits better matched to the environment will survive and reproduce more effectively than others [13,30]. Based on these three axioms, evolution theory asserts that a population will accumulate the traits that allow for more successful competition, survivability, and better reproduction. In other words, the theory of evolution underlines the key role of selection of life traits or behavioural strategies against certain criteria [20,27,31]. Since its inception, numerous efforts have been made to quantitatively formalize Darwin's theory of evolution by natural selection as a mathematical game using concepts from the well established field of game theory.
Game theory is a mathematical tool that was first developed by von Neumann & Morgernstern (1944) and mainly concerned with the rational choice between decision-makers. It was later made prominent by [25] through his concept of Nash Equilibrium [25], where the optimal outcome of a game is one where no player has an incentive to deviate from their chosen strategy after considering an opponent's choice. Game theory has a long tradition in the economic and social sciences and recently has been successfully applied to the study of evolutionary dynamics in nature [26]. In classical game theory, a game has players, strategies, payoffs and sets of rules. In biological systems, however, strategies are not the outcome of rational choices but rather are due to the behaviour shaped through evolution by natural selection. Evolution by natural selection can be considered an evolutionary game in the way that it has players, strategies, strategy sets, and payoffs. The players are the individual organisms. Strategies are their heritable phenotypic traits and evolve by means of natural selection. Payoffs are their individual fitness, and the rules are provided by the environment [31]. The mean strategy of a population changes over time, moving in the direction that will increase its fitness.
Evolutionary game theory (EGT) originated as an application of game theory to biological problems that are mainly concerned with the strategies (i.e. heritable phenotypes) of individual species. EGT provides a mathematical framework for understanding and modeling Darwinian evolution by natural selection. It was first introduced by Maynard Smith and Price [22,23] through the solution concept of an evolutionarily stable strategy (ESS) by bringing together the evolution and ecological processes. An ESS is a strategy (or set of strategies) such that, if all the members of a population adopt it, then no mutant strategy could invade the population under the influence of natural selection. In other words, at an ESS, the strategies are ecologically stable in the sense that they can co-exist together at positive population sizes, and they are evolutionarily stable when they cannot be invaded by rare alternative strategies. This has been fundamental to the development of EGT which has contributed to the understanding of ecological dynamics. EGT has advanced greatly since then as it is well evidenced through a myriad of publications on both the theoretical and application aspects (see for examples [7,24] and reference therein).
In this paper, we consider discrete population dynamical models which are governed by maps or difference equations. These equations describe typically autonomous, discrete time dynamics with population vital rates (coefficients) as the only temporal change. These coefficients can change in time through density effects or because of evolutionary processes according to Darwinian principles resulting in non-autonomous difference Equation [10]. The theory of non-autonomous discrete dynamical systems has been developed by many authors (see for example [12] and references within). The work in [12], which is the foundation of the new approach presented in this paper, focuses on the topological and the dynamical properties of the omega limit sets in non-autonomous discrete dynamical systems that are asymptotic to autonomous systems. The theoretical results were then applied to non-autonomous triangular maps and to some classical population models such as Ricker's model. This paper is organized as follows. In Section 2, we give a brief review of EGT methodology and set up single-and multi-species discrete evolutionary population models building upon prior works in [9][10][11]17,20,31]. In Section 3, we give a brief exposition of the theory of non-autonomous difference equations and the construction of the associated autonomous skew product discrete dynamical systems based on our recently published work [12]. We then use it to develop the mathematical foundation of the results in this paper as described in the new Theorem 3.2. In Section 4, we apply these results to establish the global stability of the equilibrium points of populations of single and multi-species. Here, we only consider a special type of evolutionary dynamics in which the mean trait equation is decoupled from the equation(s) of the species. The theoretical results are applied to Beverton-Holt, Ricker, and Leslie-Gower evolutionary models. In Section 5, we consider evolutionary hierarchical competition models of populations of multi-species. In hierarchical models, species are classified based on their dependence or independence from the other competing species, [5,6,16]. In the dynamical systems literature, hierarchical systems are represented by maps and are called triangular maps. These maps have a lower triangular Jacobian matrix [4]. Analogous results are obtained for evolutionary predator-prey models [3] as presented in Section 6 that would improve the results therein. In Section 7, we present general evolutionary models and discuss their dynamics using the new approach. Finally, in Section 8, we investigate the cases when the equilibrium points of the trait equation are unstable and either a saddle-node bifurcation or a period-doubling bifurcation occur. In the case of the saddle-node bifurcation and exchange of stability occurs and a new asymptotically stable equilibrium point is born. On the other hand, in the case of periodic-doubling bifurcation, the fixed point loses its stability and a stable new periodic cycle of period 2 is born. This period doubling bifurcation will lead to chaos [15]. We summarize and conclude with Section 9.

Modeling evolution as a game
To set the stage for the mathematical modeling framework and discussion presented in this paper, we begin with a brief review of the main concepts of the EGT methodology used in this paper. We follow Vincent and Brown's EGT methodology [31] and Cushing's mathematical framework for deriving discrete time evolutionary population models [10]. The development of the evolutionary and ecological dynamics is informed by Darwin's axioms. The central tenet in EGT modeling process lie in defining evolution as change in strategy (heritable phenotypes) frequency and ecological processes as the per capita growth rates of these strategies. We start by defining an evolutionary discrete model for the population dynamics of the species of interest. Let x = (x 1 , x 2 , . . . , x n ) ∈ R n + denote the vector of population numbers or densities and u = (u 1 , u 2 , . . . , u n ) denote the vector of mean strategies u i (heritable phenotypic) in U used by the species x i , where U is a one-dimensional set of biologically feasible strategies (phenotypic traits). The difference equation describing the population dynamics for species i from one time period to the next is written as a function of the strategies and the population density, is the individual fitness function (as defined by per-capita growth rate).
Model parameters in a difference equation population model may depend on a phenotypic trait v subject to the axioms of Darwinian evolution. The modeling methodology in [31] assumes the trait v is normally distributed throughout the population at all times with a constant variance and, consequently, the trait distribution is determined by the population mean trait u. Vincent and Brown [31] developed the concept of the fitness generating function or G-function (for short) as a single mathematical expression to describe the fitness function for individuals using strategy u i when it is substituted for the focal individual's strategy v in U. In other words, G(x, v, u) gives the expected per capita growth rate of a focal individual using strategy v in U when the population is in state (x, u). As a result, when v is replaced by a mean strategy value u i , the G-function generates the fitness H i for population i, that is, G(x, v, u)| v=u i = H i (x, u). G-function can also be interpreted as an adaptive landscape which is a plot of the fitness of a species with mean strategy v, given the population vector x and strategy vector u. For a discrete time model, fitness is given as the logarithm of the per capita population growth rate, i.e. ln G(x, v, u) and the fitness gradient is with respect to the trait v, i.e. ∂ ln G(x,v,u) ∂v | v=u i . The dynamical equation for the change in trait equation is described by where the constant of proportionality σ 2 is the variance of the distribution of strategies (phenotypic traits) in species x i about the mean phenotypic trait u i . σ 2 also known as the evolutionary speed. The phenotypic trait equation states that the change in the mean trait is proportional to the fitness gradient (with respect to an individuals trait). This means that larger variation in strategies produces more rapid evolution. The fitness gradient moves the strategy uphill, in the direction of the positive gradient on the adaptive landscape, which in accordance with Fisher's fundamental theorem of natural selection [18]. Evolutionary or Darwinian dynamics is then modeled by coupling the ecological population dynamics in terms of a G-function and the evolutionary dynamics equation together to give the following general system of difference equations The phenotypic trait Equation (1b) is often called the canonical equation of evolution, Lande's Equation [20], Fisher's equation for additive genetic variance [18] or the Breeder's Equation [10,31].

Single-species evolutionary models
We first start by considering the general difference equation of a single-species population given by where x(t + 1) and x(t) are the populations size or densities in successive generations and r(x(t)) is the density-dependent per capita population growth rate of the population from one time period to the next. Difference equation models of population dynamics are generally represented by x(t + 1) = f (x(t)), where f is called a map. They have been extensively investigated analytically and their global dynamics are studied using the following fundamental Theorem.
To give the model (in 2) an evolutionary dimension, we follow the EGT mathematical framework described in [10] and let the growth rate r be a function of the phenotypic trait of an individual, v, and the mean trait of the population, u, subject to the axioms of Darwinian evolution. The single-species model is then given by the coupled population and trait equations. This modeling methodology assumes the trait v to be normally distributed throughout the population at all times with a constant variance and, consequently, the trait distribution is determined by the population mean trait u. Letting r = r(x, v, u), we get Throughout this paper, we will assume that the fitness gradient is independent of the population density x, in other words: for all x, v, u in their domains. This condition leads to the decoupling of the trait equation from the population equation. Hence Equation (3) become ∂v | v=u . Next, we illustrate this modeling approach using Beverton-Holt and the Ricker evolutionary models.

Example 2.2:
The Beverton-Holt evolutionary model was first studied by Cushing [9] and is given by: where coefficients b (inherent growth rate) and c (intraspecific competition) are assumed to be functions of a phenotypic trait that are subject to Darwinian dynamics. Using the same set of assumptions from [9] which we state again here: b depends only on the individual's trait v (and not on the traits of others in the population), and c is a function of the difference between traits v and u which is maximized (or minimized) when v = u.
Using assumption (4), the Beverton-Holt evolutionary model can now be decoupled as follows:

Example 2.3:
The second example is the Ricker Evolutionary model [14,15] given by: where coefficients α (inherent growth rate) and c (intraspecific competition) are assumed to be functions of a phenotypic trait that are subject to Darwinian dynamics. Under the same assumptions as in Beverton-Holt example, the Ricker evolutionary model can now be decoupled as follows:

Multi-species evolutionary models
The above model can be generalized to multiple species. Let x 1 , x 2 , . . . , x n be n interacting species. Then an evolutionary competition model may be given by where i = 1, 2, . . . , n, x = (x 1 , x 2 , . . . , x n ) ∈ R n + . As in Section 2.2, we will assume throughout this paper that the fitness gradient of each species is independent of x i . In other words, we will assume that for all x i , v i , u i in their domains. This main assumption leads to the decoupling of the trait equations from the population equations. Hence we have the following system of difference equations give two examples to illustrate this modeling procedure: the Leslie-Gower and Ricker evolutionary competition model of two species.

Example 2.4:
Consider the Leslie-Gower competition model of two-species, .
where a and b are the intrinsic population growth rates and c ij intraspecific (for i = j) or interspecific (i = j) competition coefficients. Using the EGT methodology of Vincent and Brown [31], and assuming that a, b and c ij to be functions of a phenotypic trait that are subject to Darwinian dynamics, we define the Leslie-Gower evolutionary competition model by We will assume that a and b are functions of the corresponding individual traits v 1 and v 2 , respectively, the intraspecific competition parameters c 11 and c 22 are functions of the difference in traits v i − u i , i = 1, 2, respectively, and the interspecific competition parameters c 12 and c 21 are constants. We further assume that c ij (0, 0) = 0, i, j = 1, 2. Applying assumption (9), the Leslie-Gower evolutionary competition model can be uncoupled as follows:

Example 2.5:
Consider the Ricker competition model of two species [15], where α and β are the intrinsic growth rate for species x and y, respectively, and c ij intraspecific (for i = j) or interspecific (i = j) competition coefficients. Using the EGT methodology of Vincent and Brown [31], and assuming that α, β and c ij to be functions of a phenotypic trait that are subject to Darwinian dynamics, we define the evolutionary Ricker competition model by Making the same assumptions as in the Leslie-Gower competition model and applying assumption (9), we get the uncoupled system for the Ricker evolutionary competition model

Nonautonomous difference equations
In this section, we develop the mathematical foundation of the results in this paper and give a brief exposition of the theory of non-autonomous difference equations and the construction of the associated autonomous skew product discrete dynamical systems based on our recently published work [12]. Let {f t : t ∈ Z + } be a sequence of functions f t : R n + → R n + . Then this sequence of functions would generate the nonautonomous difference equation We will assume that the map {f t } ∞ t=0 are continuous and converges uniformly to a continuous function f : R n + → R n + . The map f generates the autonomous difference equation: Equation (14)  Then we may define a skew-product semi-dynamical system as follows. Figures 1 and 2). Then it may be shown [12] that the system (R n + × F, Z + , π) ≡ (R n + × F, π) is a discrete semidynamical system. Moreover, one may extend this dynamical system to the closure of F, In the sequel, we will restrict our study on nonautonomous population models in order to avoid certain pathological examples in which all orbits of the nonautonomous system converge to an unstable fixed point of the limiting equation. The following example from Cushing [8] illustrates this situation.  To avoid the above scenarios, we put conditions on the nonautonomous system as well as on the limiting autonomous system.
Let R n + denote the cone of nonnegative vectors in R n and let int( R n + ) and ∂(R n + ) denote the interior and the boundary of R n + , respectively. Assume A 1 : f and f t : R n + −→ R n + are continuos for all t ∈ Z + , f t converges uniformly to f as t → ∞. Then x(0) ∈ R n + implies solutions of the nonautonomous difference equation satisfies x(t) ∈ R n + , for all t ∈ Z + . (That is to say R n + is forward invariant). The same is true for solutions of the limiting equation A key assumption is . Then it is always true that x(0) ∈ int(R n + ) implies solutions of the nonautonomous difference Equation (14)

Applications
In this section, we apply Theorem 3.2 to one and two-species competition models. Let us start with the single-species evolutionary models.

Single-species evolutionary models
Example 4.1: Beverton-Holt evolutionary model [9] x(t where b(u) > 0 is twice differentiable on its domain.
The survival (positive) equilibrium is given by ( 1+c 0 x , which has the equilibrium point (x * , u * ). It is a well known fact for the autonomous Beverton-Holt model represented by the map, the equilibrium point x * is globally asymptotically stable. Hence by Theorem 3.2, the equilibrium point (x * , u * ) is globally asymptotically stable in R + × W. On the other hand, if b(u * ) ≤ 1, then (0, u * ) is the only equilibrium point of Equations (16), and (17), which is globally asymptotically stable in R + × W.

Example 4.2: Ricker evolutionary model
where α(u) > 0 is twice differentiable on its domain.
The model has two equilibria, the extinction equilibrium (0, u * ) and the survival equilibrium (x * , u * ), where x * = α(u * ) c 0 . Note that u * is any value such that α (u * ) = 0. The following statements hold true.
To prove these two statements, consider the autonomous Ricker model It is well known [14,15], that f (x) > x when x < x * and f (x) < x when x > x * . The Schwarzian derivate of f can be written as follows where P(x) = −c 2 0 x 2 + 4 c 0 x − 6 a second degree polynomial, with discriminant < 0. Thus sign(Sf (x)) = sign(−c 2 0 ), and thus Sf (x) < 0 for all x ∈ R + . Hence by Theorem 2.1, the equilibrium point x * of the autonomous Ricker model is globally asymptotically stable. Now, the maps {f t : t ∈ Z + } converges uniformly to the map f (x) = e α(u * )−c 0 x x, which has the equilibrium point (x * , u * ), and by Theorem 3.2, the equilibrium point (x * , u * ) is globally asymptotically stable on the interior of R + × W.

Example 4.6: A variation of the Leslie-Gower evolutionary model ([2]).
A variation of the above model was given in [2]. Here, the authors assumed that a and b are constants and are not dependent on the individual or the mean trait. However, the interspecific parameters are dependent on the traits of both species. The authors considered only the case where only one of the species evolutionary adapts while the other species does not adapt. Hence we have u 1 (t) = u * 1 for all t, and σ 2 1 = 0. The Leslie-Gower evolutionary model now becomes.
We assume further that species y selected a bad adaptation that led to its extinction. That is, lim t→∞ y(t) = 0 This yields We assume a > 1. One may easily show that lim t→∞ . Equation (24) converges to the limiting equation where If |f (u * 2 )| < 1, the equilibrium point u * 2 of the trait equation is asymptotically stable. Hence there is an open neighborhood U of u * 2 such that lim t→∞ u 2 (t) = u * 2 . By Theorem 8.1,

Example 4.7:
Consider now the Ricker competition evolutionary model of two-species.
Now if |1 + α (u * 1 )| < 1 and |1 + β (u * 2 )| < 1, then there exists neighborhoods G 1 of u * 1 and G 2 of u * 2 such that lim t→∞ u 1 (t) = u * 1 and lim t→∞ u 2 (t) = u * 2 if u 1 (0) ∈ G 1 and u 2 (0) ∈ G 2 . Thus the nonautonomous system (25) and (26) is asymptotic to the limiting system There are four equilibrium points The following result on the global stability of E * 4 may be found in [5] which was improved by [28]. But before stating this result, we need to introduce a few definitions.

Definition 4.8:
The set of singular points S is defined as the set of all points (x, y) ∈ R 2 + , for which the determinant of the Jacobian matrix is equal to zero.

Hierarchical competition models
By hierarchical models of n species x 1 , x 2 , . . . , x n , we mean that species x 1 is 'a silverback' species that gets first choice of resources and whose growth is limited by its own intraspecific competition, while the last species x n is an "inferior" species that gets what resources are left after all the other species, and the growth of species x i depends on all the species . , x 1 . These models have been investigated by [6,16,19], to mention only a few. Such models may be represented by the difference system where F : R n f n (x 1 , x 2 , . . . , x n )). In dynamical systems, these maps are called triangular maps [5] since their Jacobian matrix JF(X) is a lower triangular given by The main result on hierarchical models is the following result that describes their global dynamics. But before stating the theorem, we need to make the following assumptions:

Corollary 5.2: Under the above assumptions, if the map F has a locally asymptotically stable equilibrium point, then it is globally asymptotically stable.
Recall that by a Kolmogorov map, we mean a map of the form F(x 1 , x n g(x 1 , x 2 , . . . , x n )). This is to ensure that the origin is an equilibrium point and, more importantly, these are types of models that are considered here. Let us illustrate these results by the following example.

Example 5.3: Consider the 3-species Ricker competition model
This system has seven equilibrium points, the origin, three equilibria on the three axes, three planar equilibria, and one coexistence (positive) equilibrium point. We will assume here that the coexistence equilibrium x * = (x * , y * , z * ) is locally asymptotically stable.

Theorem 5.4 ([16])
: If x * is locally asymptotically stable, where x * is the coexistence equilibrium point of (30), then it is globally asymptotically stable in the interior of R 3 + . The corresponding evolutionary hierarchical model is given by Let u * 1 , u * 2 , u * 3 be the equilibrium points of the equations in (32), respectively, when α (u * Then there exist open neighborhoods U 1 of u * 1 , U 2 of u * 2 , and U 3 of u * 3 such that lim t→∞ u i (t) = u * i , i = 1, 2, 3, whenever u i (0) ∈ U i , i = 1, 2, 3. Now under Assumption (33), and using Theorem 5.4, we conclude that lim t→∞ (x(t),

A predator-prey model
The following example is motivated by the interaction of sperm whales and their food resources (e.g. squids), Ackleh et al. [3] posed the following question: In a predator-prey systems, does evolution in the prey to resist toxicants have the potential to positively affect the predator dynamics and density levels? To address this question, we will recall the dynamics of a non-evolutionary Predator-Prey Model: where φ(n) is the growth rate of prey at density level n, f (p)p is the fraction of prey consumed in the presence of p predators, b(n)n is the conversion factor when n preys are present, and s is the inherent predator survival probability when no prey is present. The following theorem summarizes the dynamics of this equation. Next, we consider the evolutionary model associated with (34), (35). We will make the following assumptions: • Only prey is assumed to evolve and the prey life span is significantly shorter than the predator • The predator is impacted by prey evolution, indirectly, through changes in prey density and directly through assumed tradeoff between toxicant resistance and the prey's ability to escape predation. • There exists a phenotypic trait u measuring the mean level of toxicant resistance in the prey. The higher the trait u, the better the prey resistance to toxicants.
The model below (36)-(38) is an example that fulfils these assumptions. where where c(u) is the searching capacity of the predator.
The trait equation can be simplified to The trait equation has solutions u = 0 and u = u * The following theorem addresses the stability of the positive equilibrium point (x * , u * ) of the evolutionary model. n(t + 1) = n(t)g(n(t), p(t), u * ) (40)

General evolutionary models
Let x 1 , x 2 , . . . , x n be n species. Then a general evolutionary model may be given by where i = 1, 2, . . . , n, x = (x 1 , x 2 , . . . , x n ) ∈ R n + . Here, it is assumed that each species x i is influenced by mean traits of its own and of all other species. Moreover, the mean trait of each species is affected by the mean traits of all the species. As in Section 2.2, we will assume throughout this paper that the fitness gradient of each species is independent of x i . In other words, we will assume that for all x i , v i , u i in their domains. This main assumption leads to the decoupling of the trait equations from the population equations. This yields the following system of difference equations where i = 1, 2, . . . , n, u = (u 1 , u 2 , . . . , u n ) is the vector mean trait of species x 1 , x 2 , . . . , x n , and The main result that we need here is the following theorem ( [12], Theorem 4.1.)

Theorem 7.1 ([12]):
Assume (44) and x * ∈ R n + is an equilibrium point of the limiting equation Then Let us revisit the evolutionary predator-prey model [3]. We modify the assumptions made in [3] as follows: • Both prey and predator are assumed to evolve and the prey life span is significantly shorter than the predator. • The prey and the predator are impacted by the evolution of the other, indirectly, through changes in density and directly through assumed tradeoff between toxicant resistance and the prey's ability to escape predation. • There exists phenotypic traits u 1 and u 2 measuring the mean levels of toxicant resistance in the prey and the predator, respectively. The higher the traits, the better resistance to toxicants.
Hence we have the following system The dynamics of this model is determined by applying Theorem (52).

Asymptotically periodic non-autonomous difference equations
To this end, we have investigated the case when the trait equation has stable equilibrium points. In this section, we will investigate the cases when the equilibrium points of the trait equation are unstable and either a saddle-node bifurcation or a period-doubling bifurcation occur. In the case of the saddle-node bifurcation and exchange of stability occurs and a new asymptotically stable equilibrium point is born. On the other hand, in the case of periodic-doubling bifurcation, the fixed point loses its stability and a stable new periodic cycle of period 2 is born. This period doubling bifurcation will lead to chaos [15].

Theoretical development
Let us assume that {f t : t ∈ Z + } be a sequence of functions f t : R n + → R n + , that converges uniformly to the p-periodic system G = {g t : t = 0, 1, 2, . . . , p − 1}, where g t : R n + → R n + . Then we have two equations, a non-autonomous difference equations and a non-autonomous periodic difference equation: and y(t + 1) = g t (y(t)), t = 0, 1, 2, . . . , (p − 1)) We extend the skew-product semi-dynamical system (R n + × F, Z + , π) ≡ (R n + × F, π) to the closure of F, F = F {G} by letting π((x, g i ), t) = ( t,i (x), g (i+t) mod p ) (x), f ). (see Figure 1, 6). In this setting, we are going to make two assumptions similar to those we had in Section 4. We assume A 1 : g i and f t : R n + −→ R n + are continuos for all t ∈ Z + , f t converges uniformly to G as t → ∞. Then x(0) ∈ R n + implies solutions of the nonautonomous difference equation x(t + 1) = f t (x(t)), x = (x 1 , x 2 , . . . , x n ) ∈ R n + .
satisfies x(t) ∈ R n + , for all t ∈ Z + . That is to say R n + is forward invariant. The same is true for solutions of the limiting system y(t + 1) = g t (y(t)), t = 0, 1, 2, . . . , (p − 1)) ( and A 2 f t : int(R n + ) −→ int(R n + ) Then it is always true that x(0) ∈ int(R n + ) implies solutions of the nonautonomous difference Equation (54) satisfies x(t) ∈ int(R n + ), for all t ∈ Z + . The main result that here is the following theorem ( [12], Theorem 4.1.)

Theorem 8.1 ([12]):
Assume A 1 and A 2 and that the periodic system G has a globally asymptotically stable cycle c p of period p or a divisor of p. Then (i) if c p ∈ int(R n + ), and if it is globally asymptotically stable on int(R n + ) as a periodic cycle of the limiting Equation (15), then all solutions of the nonautonomous difference Equation (14) with x(0) ∈ int(R n + ) tend to c p . (ii) if c p ∈ ∂(R n + ), and if it is globally asymptotically stable on int(R n + ), then all solutions of the nonautonomous difference Equation (14) with x(0) ∈ int(R n + ) tend to c p .
Proof: Let = g p−1 • g p−2 • · · · • g 1 • g 0 . Then an equilibrium point of the composition map is a periodic point of period p or of a divisor of p. The dynamics of the periodic system G = {g t : t = 0, 1, 2, . . . , p − 1} is determined by the single map . Hence, applying Theorem 8.1, the conclusion of the theorem follows.
Let us illustrate the effectiveness of this theorem by the following example of a twopeaked adaptive landscape model.
• The equilibrium u * 2 is unstable. • The equilibrium u * 3 is asymptotically stable if σ 2 < 2 3 . A period-doubling bifurcation occurs if σ 2 > 2 3 . This forces the non-autonomous system to have a period-doubling bifurcation as well.

Conclusion
In this paper, we applied a new approach to a special class of discrete time evolution singleand multi-species models. The main tool of our mathematical analysis was the embedding the non-autonomous evolutionary models into autonomous difference systems via the construction of a skew-product discrete dynamical systems as developed in Elaydi et al. [12]. This enabled us to utilize the well-developed theory of autonomous dynamical systems. While the theory may be extended to continuous evolutionary systems, we restricted our study to discrete-time evolutionary models only.
Our analysis, however, applies to special cases of evolutionary dynamics and does not apply to a more general setting. It is still an open question of how to extend our techniques and methodology to the cases when the trait equation and the species equations are not decoupled. This clearly requires the development of new tools to establish global stability of equilibrium points, including using the method of Liapunov functions which was effectively used in [9], but only for the simple Beverton-Holt single-species evolutionary model. Finding the right Liapunov function has, however, eluded researchers in discrete-time modeling of populations. Future work will broaden the study presented in this paper and will try to implement the ideas alluded to in the last sentence. Another direction, not explored here, is to apply our methodology to investigate models of infectious disease that are of paramount importance for all the living species. Such study may prove to be a breakthrough in this important field.