Difference equations as models of evolutionary population dynamics

ABSTRACT We describe the evolutionary game theoretic methodology for extending a difference equation population dynamic model in a way so as to account for the Darwinian evolution of model coefficients. We give a general theorem that describes the familiar transcritical bifurcation that occurs in non-evolutionary models when theextinction equilibrium destabilizes. This bifurcation results in survival (positive) equilibria whose stability depends on the direction of bifurcation. We give several applications based on evolutionary versions of some classic equations, such as the discrete logistic (Beverton–Holt) and Ricker equations. In addition to illustrating our theorems, these examples also illustrate other biological phenomena, such as strong Allee effects, time-dependent adaptive landscapes, and evolutionary stable strategies.


Introduction
Difference equations have a long history of use as discrete time models of population dynamics. These equations describe typically autonomous, discrete time dynamics and assume that the only temporal change in vital rates are due to a dependence on population density (giving rise to nonlinear difference equations). There are, however, numerous other mechanisms and circumstance that cause vital rates to change over time. For example, the population's physical and biological environments can change or fluctuate, either regularly (e.g. due to daily, monthly, or seasonal cycles) or randomly (due to stochastic fluctuations of weather, resource availability, etc.). An individual's vital behaviour and activities (reproduction, survival, resource acquisition, etc.) can similarly change and fluctuate. Such explicit dependencies on time can be modelled, for example, by making equation coefficients periodic in time, in the case of regular fluctuations, or random in time, in the case of stochastic fluctuations. The result in the first case is a nonautonomous, periodically forced difference equation, while in the second case is a stochastic difference equation. These types of difference equations have received increasing mathematical attention in recent years.
There is another reason why model coefficients might not remain constant in time. Birth rates, survival probabilities, and behavioural characteristics that effect reproduction and survival are subject to change due to Darwinian evolution. A modelling methodology that included the principles of Darwinian evolution could help to provide explanations for why populations possess the traits and characteristics they do, analyse the ability of a population to adapt and survive in the face of environmental changes, understand the development of pathogen immunities, and so on. In this paper, we describe a methodology (called evolutionary game theory) that a modeller can use to account for evolutionary changes to the coefficients in any difference equation population model of interest. This method is developed in the book [22] (also see [16]). It is based on the three axioms of Darwinian evolution: variability, heritability, and differential fitness.
In order to set the groundwork, we begin in Section 2 with an account of population modelling using difference equations and give general theorems that focus mainly on the basic question of survival verses extinction, or mathematically on the bifurcation that occurs when the extinction equilibrium destabilizes. This primary and basic bifurcation results in the creation of positive (survival) equilibria whose stability, at least near the bifurcation point, is dependent on the direction of bifurcation. While we give some general results about positive equilibria outside the neighbourhood of the bifurcation point (Theorems 2.2 and 3.1), we do not in this paper consider the well-known propensity of difference equations to display further bifurcations and even cascades to chaos. These dynamic phenomena are model specific and our focus here is on the primary bifurcation that is quite general and concerned with the fundamental biological question of extinction versus survival.
Given a difference equation model for the dynamics of a population, evolutionary game theory assumes that model coefficients are dependent on a suite of phenotype traits (strategies). From a definition of fitness based on the population equation, the method sets up a system of difference equations that includes both the population dynamics and the dynamics of the trait means, which together define the Darwinian dynamics. These equations are based on trait-driven variability [22] within a (sexually or asexually reproducing) population and account for the coupled dynamics of both the population and its mean traits. In Section 3, we generalize the fundamental bifurcation Theorem 2.2 to Darwinian equations. The results in Section 3 generalize those in [7] for the case of a single trait (m = 1).
In Section 4 we give several examples to illustrate the modelling methodology and the results in Section 3. These examples are based on classic difference equations for population dynamics, namely, the discrete logistic ( Beverton-Holt) and Ricker equations. These examples also illustrate several other biological phenomena. Another important concept in evolutionary theory is that of an evolutionary stable strategy (ESS). 1 This is based on invasion-driven variability, in that the question is whether a population can be invaded by a mutant population which is reproductively isolated from the original population. The reproductive isolation assumption is appropriate for asexually reproducing populations, but is more problematic for sexually reproducing populations unless some additional assumptions are made, such as strong assortative mating or a long term period of geographic isolation. If a stable equilibrium of the Darwinian equations lies on a global maximum of the adaptive landscape (see Remark 3.2 in Section 3), then it is an ESS, otherwise it is not. This is the ESS maximum principle in [22]. We illustrate fitness landscapes and ESS traits using the examples in Section 4 and see how they are not fixed in time while the mean trait attempts to move uphill to a peak (which can allow it to move from one peak to another) and how traits from a stable equilibria are not always ESS traits. The examples in Section 4 also illustrate the relevance of unstable bifurcations, which can (and typically do in population models) lead to strong Allee effects, i.e. a multi-attractor scenarios in which both an extinction and a positive (survival) equilibrium are stable.

Difference equation models of population dynamics
As a population dynamic model, a difference equation predicts the population density x t+1 at time t+1 from a knowledge of the density x t at time t. When additions and subtractions from the population are due only to births and deaths (i.e. there is an absence of immigration/emigration, seeding/harvesting, etc.), individuals present at time t+1 were either born during the time interval or they were present at time t and survived the time unit. Thus, to model g(x t ) we add together models for newborns and survivors. If b denotes the per capita number of newborns born during the time interval that survive to the end of the time unit, then the total number of newborns at time t+1 is bx t . If s denotes the per capita probability of surviving one time unit, the sx t is the number of survivors at time t + 1. This leads to the difference equation where b > 0, 0 ≤ s < 1.
In carefully derived models, attention is paid to the unit of time. Often the unit of time is taken to be such that an individual can reproduce no more than once during one time unit (for example, a maturation period). In this case, s > 0 implies there exist survivors who live to reproduce more than once. Such a population is called iteroparous and generations overlap. This is in contrast to a semelparous population (or monocarpic, for plants) when s = 0 and an individual has (at most) one reproductive event before it dies. In this case, generations do not overlap. If b and s remain constant in time, then the population goes extinct if and grows exponentially if r > 1. There are, however, numerous reasons why the vital rates b and/or s would not remain constant in time. These include their dependence on population density, evolutionary adaptation, and time fluctuating environments (such as daily, monthly, or annual variations in the environment). Suppose instead that b and/or s are at any time t dependent on population density x t at that time. Then b and s are changed by factors β(x t ) and/or σ (x t ), respectively, and the difference equation becomes nonlinear To be biologically meaningful, it is required that bβ(x) be nonnegative and sσ (x) remain between 0 and 1. A1: Assume b > 0 and 0 ≤ s ≤ 1. Assume β and σ are twice continuously differentiable on an open interval containing the half line x ≥ 0 and satisfy β(x) ≥ 0 and 0 ≤ sσ (x) ≤ 1 for all x ∈ .
We will consider the case when the birth rate is positive in the absence of any density effects, i.e. β(0) > 0. We can assume without loss in generality that β(0) = 1 and, in this way, b retains its interpretation of the density-free or inherent birth rate. Similarly, so that s is the inherent survival rate, we assume σ (0) = 1.
A2: Assume in A1 that β(0) = σ (0) = 1. A negative derivative ∂ x β(x) < 0 represents a negative effect of population density on the birth rate in the sense that an increase in population density from x results in a lower birth rate. A positive derivative ∂ x β(x) > 0 represents a positive density effect on the birth rate in the sense that an increase in population density from x results in a higher birth rate. Negative and positive density effects can also be present in the survival factor σ (x).
It is typically assumed in population models that negative density effects occur, if not at all density levels, at least at all high densities x > > 0. Such negative effects are generally considered to occur because of crowding effects, competition for resources, susceptibility to diseases, etc. Positive density effects at low population densities x 0 are called component Allee effects, which have recently received considerable attention in both the theoretical and biological literature [2,9,19]. Many mechanisms for such positive density effects have been found in a wide variety of species, including protection against predation by increased population densities, an increased probability in successfully mating in sexually reproducing species, group foraging, etc.
Equilibria (fixed points) of (2) are roots of the algebraic equation Population extinction is represented by the root x e = 0, which is an equilibrium for all values of b and s. Other biologically feasible equilibria are positive roots x > 0, whose existence depend (in general) on b and s. Suppose we fix s and consider positive equilibria x > 0 as a function of b. An equilibrium pair (b, x) satisfies of the equation The graph G of this equation, i.e. the graph of  which positive equilibria exist, how many positive equilibria exist at a specific value of b, and bifurcation points where the number of equilibria changes. The (local) stability properties indicated in Figure 1 follow from Theorem 2.1 below. From this graph, one can also locate intervals where multiple stable equilibria exist and hysteresis loops and their accompanying tipping points (fold or tangent or blue-sky bifurcations) occur. This includes a strong Allee effect, which is defined as a multiple attractor scenario in which the extinction equilibrium is stable and there also exists (at least one) non-extinction attractor, such as a stable positive equilibrium. By the linearization principle, a positive equilibrium x > 0 is (locally asymptotically) stable or unstable if is less or greater than 1, respectively. Suppose the point (b, x) lies on the graph G. Then from (3)  If b > b 0 (or b < b 0 ) for the bifurcating positive equilibrium points (b, x) near the bifurcation point (b 0 , 0), then the bifurcation is said to be forward (or backward). That is to say, the bifurcation is forward if the graph G is increasing, and backward if it is decreasing, at (b 0 , 0). By Theorem 2.1 near the bifurcation point the bifurcating positive equilibria are stable if the bifurcation is forward and are unstable if it is backward. Thus, the direction of bifurcation determines the stability properties of the bifurcating positive equilibria near the bifurcation point. For conciseness, we simply say a forward bifurcation is stable and a backward bifurcation is unstable.
A forward (respectively backward) bifurcation occurs if db/dx| x=0 > 0 (respectively < 0) and the calculation shows that the direction of bifurcation is, if κ = 0, determined by a weighted average κ of the survival and birth rate sensitivities ∂ x σ (0) and ∂ x β(0) at low population levels. This, together with Theorem 2.1, gives us the following result concerning the bifurcation at the extinction equilibrium.

Remark 2.2:
Backward bifurcations typically result in a strong Allee effect [5], as illustrated in Figure 1 when κ > 0. This is because a typical assumption is that the birth rate drops to 0 as density increases without bound, i.e. lim x→+∞ β(x) = 0. This implies, together with equation (3) that defines the graph G, that b increases without bound as x → +∞ which in turn implies b has a minimum at a critical point less than b 0 (see in Figure 1). By Theorem 2.1 there are also stable positive equilibria in a neighbourhood of this critical point and, consequently, there occurs a strong Allee effect on an interval of b values less than b 0 .

Remark 2.3: Theorem 2.1 remains valid if instead of b
we use the inherent population growth rate r = b+s or the inherent reproduction number R 0 = b/(1 − s) as bifurcation parameters. In this case, the bifurcation at x = 0 occurs at r = R 0 = 1. See [4].
As is well known, positive equilibria of difference equations are not always stable and other bifurcations, such as period doubling cascades to chaos, can occur along the graph G. The occurrence of secondary bifurcations is highly dependent on the specific properties of the nonlinearities in β(x)x and σ (x)x.

Difference equation models of Darwinian dynamics
The coefficients b and s in the basic equation (1) can change in time not only through density effects, as modelled by equation (2), but can change in time because of evolutionary processes according to Darwinian principles. In this section, we describe one way in which evolutionary adaptation can be incorporated into population model Equations (2), namely by the method of evolutionary game theory or Darwinian dynamics [22]. This method assumes the existence of a vector v = col (v i ) of m ≥ 1 phenotypic (real valued) traits v i that are subject to Darwinian evolution, i.e. traits that determine an individual's fitness, have variation throughout the population, and have a heritable component. The traits v i are assumed (normally) distributed throughout the population with means u i , variances σ 2 i , and covariances σ ij = σ ji , i = j, as given by the variance-covariance matrix Fitness of an individual with trait vector v is defined to be ln r where Here an individual's fitness can depend on the traits of other individuals as represented by the mean u i (which is why the methodology is sometimes referred to as evolutionary game theory). The reason we assume the inherent birth and survival rates b and s do not depend on u is that, by definition, they are rates of an individual in the absence of density effects, i.e.
To retain the interpretation as inherent rates for b(v) and s(v) we assume The dynamics of the population and its mean phenotypic traits are governed by the equations (see (5.16), page 138, in [22]). (The variance-covariance matrix V is assumed constant in time.) Here ∇ v is the gradient (column) vector of r with respect to v. The trait equation (5b) is often called the canonical equation for evolution [10], Lande's equation [12,13], Fisher's equation for additive genetic variance [1], or the breeder's equation [15]. If V is the zero matrix, then there are no trait dynamics (i.e. evolution does not take place) and the Darwinian equations (5) reduce to just the population dynamics equation (5a) with fixed u (which is an equation of type (2) in Section 2 ). We will assume H3. The variance-covariance matrix V is invertible. An equilibrium (x, u) of the Darwinian equations (5) satisfies the algebraic equations We are interested in two types of equilibria: extinction equilibria in which x = 0 and positive equilibria in which x > 0. A pair (x, u) = (0, u e ) is an extinction equilibrium if and only if u = u e is a critical trait, i.e. satisfies the equation Our goal is to study the existence and (local asymptotic) stability of both types of equilibria and to obtain a generalization of Theorem 2.2 for the Darwinian equations (5) when r is given by (4). To do this, we will place some additional assumptions on the inherent birth and death rates. We begin by assuming This requirement is motivated by the biological assumption that there exists a trait at which the inherent birth rate b has an extreme value (either a minimum or maximum) and that there is a trade-off between reproduction and survival, so that s has the opposite extreme (maximum or minimum) at the same trait. Trade-offs such as this are central and important in the study of life history strategies [18]. Of course, assumption (6) allows other biological scenarios as well (for example, that b and s have the same type of extremum properties at v c and instead are positively correlated). Since we place no scale or units on the trait vector v we can assume, without loss in generality, that v c = 0 ∈ U. In keeping with the pointof-view taken for the non-evolutionary model equations in Section 1, we will use b(0) as a bifurcation parameter. Toward this end we write b(v) = b 0b (v),b(0) = 1 and use b 0 as the bifurcation parameter. This motivates the following assumption. H4. In H1 we assume Under H4, b 0 and s 0 are the inherent birth and survival rates for individuals with trait v = 0. From (4) and H2 we find that We include the bifurcation parameter b 0 in the argument list of and in the Darwinian equations (5) where we have made use of the notation Under assumption H4, the equilibrium equations (7) is where x, u, u) with respect to u, I is the m × m identity matrix, and τ denotes matrix transpose. Evaluated at an extinction equilibrium (x, u) = (0, 0) the Jacobian becomes the block diagonal matrix where by H2  In case (b) the population is not threatened with extinction from the extinction equilibrium (0, 0), and so we focus on case (a) when the extinction equilibrium loses stability as b 0 increases through the critical value b * 0 . It is often assumed that the covariances σ ij among the traits are equal to 0 or in any case are small [22]. It is proved in [21] (also see [20]) that if H(b 0 , 0, 0) is negative definite and V is diagonally dominant (i.e. The destabilization of the extinction equilibrium at b 0 = b * 0 that occurs in Theorem 3.1(a) and in Corollary 1 suggests that a bifurcation will occur that creates nonextinction equilibria and, possibly, positive equilibria. This is the subject of Theorem 2.2 below.
Let a superscript 0 denote evaluation at the extinction equilibrium (x, u) = (0, 0) and 0, 0, 0) , etc. and Under the assumptions that H5. κ = 0 and H 0 is nonsingular, we prove the following theorem in the appendix.
(a) For b 0 near b * 0 there exist equilibria of the form

Remark 3.2:
The inequality κ > 0 (respectively κ < 0) means, as it does for the nonevolutionary equations in Section 2, that negative (respectively positive) density effects are dominant at those low-level population density points near the bifurcation point (b 0 , x, v, u) = (b * 0 , 0, 0, 0). According to the formula for κ, these density effects are measured by the weighted average (1 − s 0 )∂ 0 x β + s 0 ∂ 0 x σ of the derivatives (sensitivities) of the birth and survival rates with respect to population density at the bifurcation point. With regard to the existence of positive equilibria outside a neighbourhood of the bifurcation point, it is shown in [17] that the bifurcating branch of b 0 values and positive equilibria (x e , u e ) connects to either another bifurcation point located at an extinction equilibrium or is unbounded. With regard to stability, it is well known for non-evolutionary model equations (2) that positive equilibria can destabilize (or stabilize), causing further (secondary) bifurcations, and even routes to chaos, so this is certainly a possibility for the Darwinian equations (5). These issues are model specific and we do not pursue them here.

Remark 3.4:
As pointed out in Section 2, backward bifurcations in population models usually give rise to strong Allee effects. We expect that to be true for EGT model equations as well (when κ < 0 in Theorem 3.2), although we do not study this here, since it concerns positive equilibria outside of a neighbourhood of the bifurcation point (see Remark 3.3). An example in which this occurs is given in Section 4.2 (cf. Figure 2).  Another concept of a stable trait involves its resistance to invasion by mutant trait which is reproductively isolated from the original population. By a mutant trait is meant one possessed by a similar population that is reproductively isolated (and hence a separate species) from the original population. If the trait u e of a stable positive equilibrium (x e , u e ) also possesses this invasibility property, then the trait u e is called an evolutionarily stable trait or an ESS [22].
One need not construct multiple species models to find an ESS, however. It is shown in [22] that a trait u e from a (locally asymptotically) stable positive equilibrium of the Darwinian equations (7) is an ESS if the global maximum of the fitness function ln r(x e , v, u e ), as a function of v, occurs at v = u e (this is called the ESS maximum principle). Thus, an analysis of a Darwinian model (7) is often accompanied by a plot of ln r(x e , v, u e ) (the adaptive landscape) to see if a global maximum occurs at v = u e . The equilibrium equations satisfied by (x e , u e ) imply the trait component u e is a critical point on the adaptive landscape, but a maximum does not, of course, necessarily occur at a critical point. In fact, examples are given in [22], and reference cited therein, in which u e is located at a minimum of ln r(x e , v, u e ), an occurrence which plays a role in theories of speciation (also see [10]). This might seem surprising since the trait equation (5b) implies that u t climbs uphill on the adaptive landscape. The explanation lies in the fact that the adaptive landscape ln r(x t , v, u t ) is not fixed in time. Because of this, the path of u t can sometimes be complicated and non-intuitive. An example appears in Section 4.3.

Examples
To devise an evolutionary population model (7) one needs to specify the density factors β(x) and σ (x) and how they, and the inherent fertility and survivorship b and s, depend on the adaptive trait v and the population mean trait u. In this way, r(b 0 , x, v, u)

An EGT discrete logistic model
If β (x) = 1 1 + cx and s 0 = 0 (i.e. the population is semelparous and generations do not overlap) then the fitness function ln(r(b 0 , x, v, u)) in the Darwinian Equations (7) is If no evolution occurs (V is the zero matrix), then the population equation reduces to . This is the famous discrete logistic (or Beverton-Holt) equation 2 whose dynamics are well known: if b < 1 the extinction equilibrium is (globally asymptotically) stable and if b > 1 the extinction equilibrium is unstable and there exists a (globally asymptotically) stable survival equilibrium. If evolution does occur (V is not the zero matrix), then equations (7) become a Darwinian version of the discrete logistic equation. The coefficient c(v, u) measures the negative effect that population density has on the birth and survival rates of an individual with rate v when the population has mean trait u. This negative effect can be attributed to intraspecific competition among individuals for resources necessary for reproduction and survival. An increased (decreased) value of c(v, u) means that an individual with trait v suffers a greater (lesser) negative effect (at a fixed population size x).
A common assumption is that the competition coefficient depends on the difference v − u, that is to say, how different an individual with trait v is from the average individual with trait u. The resulting Darwinian discrete logistic model, with competition coefficient c(v − u) was studied in [6] with a one dimensional m = 1 trait vector (i.e. a single adaptive trait) where it was shown that if evolution is not too fast, then b 0 < 1 implies the extinction equilibrium (x, u) = (0, 0) is (globally asymptotically) stable and b 0 > 1 implies the extinction equilibrium is unstable and there exists a (globally asymptotically) stable survival equilibrium (x, u), x > 0. These results are also implied by Theorem 3.2 with m = 1 but only in a neighbourhood of the bifurcation point (0, 0) when b 0 = 1. Theorem 3.2 also applies when m > 1, but the global results obtained in [6] remain an open question in this case.
Here is an example with two traits (m = 2): In this example, the inherent birth rate is binormally distributed as a function of two traits v 1 and v 2 and the competition coefficient c(v − u) is based on an assumption of (asymmetric) hierarchical competition. Individuals with larger traits v i experience less intraspecific competition; c 0 is the competition intensity when an individual's trait v equals the mean trait u. We assume no trait covariances, so that It is clear that H1-H3 hold with trait interval U = R. From we get the Darwinian equations we see that the origin v 1 = v 2 = 0 is the only critical trait and H4 holds. By Theorem 3.1 it follows that the extinction equilibrium (x, u 1 , u 2 ) = (0, 0, 0) destabilizes as b 0 increases through 1.
Calculations show which is invertible and negative definite, and that It follows that H5 holds and Theorem 3.2(a) implies the forward bifurcation of positive equilibria at b 0 = 1. From we see that Theorem 3.2(b) applies provided both σ 2 i < 2w 2 i . That is to say, if the variance of the (constant) traits in the population is less than twice the variance in the birth rate as functions of the trait vector components. In this case, it follows that a forward and stable bifurcation of survival equilibria (x e , u e ) occurs at b 0 = 1. We next ask whether the trait vectors u e associated with these equilibria are ESS traits or not.
Near the bifurcation point, the bifurcating positive equilibria are approximated by the expansions in Theorem 3.2(a) with b * 0 = 1 and i.e. by The positive valued function r(b 0 , x e , v, u e ) of v is bounded and tends to 0 as |v| → ∞. Straightforward differentiation analysis shows that there exists a unique critical point (i.e. vector v at which ∇ v r(b 0 , x e , v, u e ) = 0). It follows that r(b 0 , x e , v, u e ), and therefore the adaptive landscape ln r(b 0 , x e , v, u e ), has a global maximum at this critical point. Since the trait equation implies the equilibrium v = u e is a critical point of the adaptive landscape, it follows that a global maximum occurs at u e , which is therefore an ESS.
For this particular Darwinian version of the discrete logistic, we draw the following conclusions. If evolution is not too fast (0 < σ 2 i < 2w 2 i ), a forward stable bifurcation of positive equilibria (x, u e ) occurs in a neighbourhood of the extinction equilibrium at b 0 = 1. (We conjecture that in fact there exists a globally asymptotically stable equilibrium for each b 0 > 1, as for the non-evolutionary version when σ 2 i = 0). Therefore, the population is threatened with extinction for b 0 1 and survives and equilibrates for b 0 1. In the latter case, the mean trait u e is an ESS. Moreover, the components of u e are positive, which means evolution does not optimize the (inherent) birth rate of the average individual (i.e. individuals with trait v = u e ). Instead, evolution decreases the inherent birth rate of the average individual in favour of increasing its competitive advantage against other individuals (in the sense that the competitive density coefficient c(v − u), the competitive effect on individuals with trait v = u e , is greater when u = u e than when u = 0).

An EGT Ricker/Allee model
In the example in Section 4.1 only negative density effects are included. This is why the bifurcation at b 0 = 1 is forward. As pointed out in Remark 3.2 in Section 3 positive density effects are required at low population densities for a backward bifurcation, which will occur if they outweigh any negative density effects at low densities (as measured by κ < 0). In this section, we consider a model equation with positive density effects included. We will based this model on the famous Ricker equation in which a negative density effect on the birth rate is described by an exponential factor exp(−cx).
As with the discrete logistic model in Section 4.1, we take in the Darwinian equation (7). In this example, however, we allow iteroparity (s 0 > 0) and assume that both the birth and survival rate are optimal at trait v = 0 as given by , 0 ≤ s 0 < 1.
In this example, inherent birth and survival rates both attain their maxima at trait v = 0.
In the term β describing the effect of density on the birth rate we include the factor exp(−cx) with the same trait dependent competitive coefficient c given by (10) as used in Section 4.1. However, we also include a positive density effect on fertility by including an additional factor a(x) that is increasing as a function of x, so that The factor a(x) is called a component Allee effect [2]. In this example, we use a (x) = exp p αx 1 + αx in which the coefficient α measures the benefit an individual attains from a population of density x and p > 0 measures the sensitivity of the response to low density population increases, i.e. ∂ 0 x a = αp. We assume that α is a function of v − u which attains its maximum at v = u. The reason for this is that the component Allee effect is assumed to be due to an intraspecific cooperative mechanism and hence the maximum benefit to an individual is attained when its trait is the same as most other individuals, as measured by the mean trait u. Specifically, we take Finally, we again assume no trait covariances, so that V is given by (11). All these model prescriptions result in the fitness function ln (r(b 0 , x, v, u)) with with c(v − u) and α(v − u) given by (10) and (12) respectively. A straightforward calculation shows that v = 0 is a critical trait for all values of b 0 (and is the only critical trait). Therefore, assumptions H1-H4 hold for the Darwinian equations defined by this r(b 0 , x, v, u). Furthermore, straightforward calculations show that which is invertible and negative definite, and that It follows that assumption H5 holds if p = c 0 . Therefore if p = c 0 then Theorem 3.2(a) implies the bifurcation of positive equilibria at b 0 = 1 − s 0 . From we see that Theorem 3.2(b) applies provided  of the bifurcation point, does not address this possibility. However, we can corroborate that it occurs in this model, for selected parameter values, by means of numerical simulations. See Figure 2 for an example that shows initial condition dependent survival.
In addition, for this example the plot of the adaptive landscape in Figure 3 for the case in Figure 2(A) shows that the positive equilibrium in this example is an ESS since u e is located at a (rather sharply defined) global maximum. Note that there is also a lower peak located near u = 0 where the inherent vital rates b(v) and s(v) are maximized.
As pointed out in Remark 3.5, adaptive landscapes in Darwinian models are not fixed in time. This is illustrated in Figure 4 which shows several landscapes at some select times, including the initial landscape. The chosen times show the trait vector u t making an interesting sojourn, starting at a location where it would appear poised to climb a peak located near u = 0. As time goes on, however, another peak appears and at the times shown u t crosses a saddle between the two peaks (all the while moving upward on the landscape, as the model equation for the trait dynamics prescribes) so as to ultimately climb the higher peak shown in Figure 4. It is not intended here to imply that this particular sequence of landscapes and the mean trait trajectory always (or even frequently) occurs in this example. Indeed, each solution gives rise to its own sequence of landscapes and trajectory paths. This example was chosen simply to illustrate the point of temporally changing landscapes and to illustrate one interesting possibility, namely, that of adaptive traits crossing a valley on the landscape while climbing uphill according to the Darwinian trait equation.

A model with non-ESS equilibrium traits
In This model adds survival (making the population iteroparous and have overlapping generations) with a binormally distributed survival rate s(v) with 0 < s 0 < 1. The distribution for b(v) is designed so that the inherent birth rate is maximized at several trait vectors, namely v = col(0, 0), col(1, 0), col(0, 1) and col (1,1). This leads to the fitness function ln(r(b 0 , x, v, u)) with ).  The purpose of this example is to show that a stable positive equilibrium (x e , u e ) does not necessarily correspond to an ESS trait u e . This is illustrated in Figure 5 where an equilibrium lies at a peak on the adaptive landscape that is not a global maximum. This stable equilibrium trait u e is located very near u = 0 where both the inherent birth rate and the survival rate are maximal, but it is not an ESS trait.
The example in Figure 5 turns out to have multiple positive equilibria (Theorem 3.2 deals only with positive equilibria in a neighbourhood of the bifurcation point (x, u) = (0, 0) at b * 0 = 1 − s 0 . Other positive equilibria may well exist outside this neighbourhood). Figure 6 shows the adaptive landscape corresponding to a different positive equilibrium, approached by an orbit with an initial condition different from that used in Figure 5, and that the equilibrium trait u e ≈ col(1.10, 0.95) does, in this case, correspond to a global maximum and therefore is an ESS trait. This trait is near u ≈ col(1, 1) which is where the inherent birth rate is maximized, but the inherent survival rate is not.

Concluding remarks
In Section 2, we present a basic method for deriving difference equations as models of population dynamics. This general method is simply a bookkeeping method that accounts for all sources of individuals present at discrete census times by predicting the population density at time t+1 from a knowledge of the population density at time t. For closed populations in which the only source of new individuals is by means of a per capita birth rate and the only loss is from the probability of dying (or otherwise being removed), the individuals present at time t+1 are newborns and survivors. The general model equation (2) results from the assumption that per capita birth and survival rates are population density dependent. Simple methods of analytic geometry and the stability linearization principle lead to the general bifurcation diagrams in Figure 1, which account for the existence of positive equilibria (as a function of the inherent birth rate), the destabilization of the extinction equilibrium at a critical value of the inherent birth rate, and the stability properties of positive equilibria near bifurcation points.
Evolutionary game theory extends the nonlinear population model Equation (2) to an evolutionary setting by assuming equation coefficients depend on a vector of phenotypic traits subject to Darwinian evolution and providing an expanded system of difference equations (7) for the population dynamics plus the dynamics of the mean phenotypic traits. Under assumptions H1-H5, Theorems 3.1 and 3.2 extend the fundamental bifurcation that occurs when extinction equilibria are destabilized to the Darwinian model equations (7).
The bifurcation Theorem 3.1 has also been extended to Darwinian models for matrix population equations, i.e. systems of difference equations that model the dynamics of structured populations [3,8]. These extensions do not include Theorem 3.2 as a special case, however, since they utilize a different bifurcation parameter.
Some open questions remain. The correspondence between the direction of bifurcation and the stability of the bifurcating positive equilibria provided by Theorem 3.2(b) assumes the speed of evolution is not too fast (ρ[I + VH 0 ] < 1). What is the nature of the bifurcation when the speed of evolution is fast? In this case, the extinction equilibrium loses stability because the trait equation destabilizes (see the Jacobian (9)), which will result in new extinction equilibria (possibility with a non-constant trait component if the destabilization is because an eigenvalue leaves the complex unit circle somewhere other than 1). A study of these extinction states, their stability in the Darwinian model (7), and the nature of bifurcations when they destabilize is an open problem. Another open problem is the existence of positive equilibria outside of a neighbourhood of the bifurcation point and in particular of the bifurcating branch itself. The global extent of the bifurcating branch is described in [17], using general global bifurcation theorems from nonlinear functional analysis, but only for the case of a single phenotypic trait (m = 1). This global existence theorem is important in establishing strong Allee effects from backward bifurcations [5]. This also involves establishing stability and instability results of positive equilibria outside a neighbourhood of the bifurcation point (as inspired by Theorem 2.1 and Figure 1 for the non-evolutionary case).

Notes
1. Strategy and trait are synonymous in this theory. The word strategy derives from game theory and from viewing this method of modelling evolution as a 'game' in the sense that an individual's fitness is the 'pay-off' from its strategy (trait) and that of others, represented by the mean. However, evolution differs from classic game theory in that an individual does not rationally choose a strategy, but instead biologically inherits a trait. 2. The equation seems to have been first extensively utilized as a population model by Leslie [14] in his study of logistic growth for age structured populations.
Making use again of ∇ 0 v r = 0 and ∇ 0 u r = 0, we find from the first equation that From the second equation we have From H4 and H2 follows that ∂ b 0 ∇ v r = 0, that ∇ u ∇ v r is the zero matrix and that Thus, (A2) evaluated at these points becomes so that ∂ b 0 υ(b 0 ) = u 1 .
(b) We consider the Jacobian of the Darwinian equations (7) evaluated at the positive equilibrium (x, u) = (ξ(b 0 ), υ(b 0 )), which we denote by J (b 0 ) := J (b 0 , ξ (b 0 ) , υ (b 0 )) , and calculate its eigenvalues to lowest order in b 0 − b * 0 . From (9) we have x ∇ v r I + VH 0 whose eigenvalues are 1 and those of I + VH 0 , the latter of which all lie in the complex unit circle by assumption. The eigenvalues of the Jacobian evaluated at the bifurcating equilibria approach these eigenvalues as b 0 approaches b * 0 and therefore m of them lie inside the complex unit circle for b 0 close to b * 0 . Let λ(b 0 ) denote the eigenvalue that approaches 1, i.e. λ(b * 0 ) = 1. Our goal is to calculate ∂ b 0 λ(b * 0 ) in order to determine whether λ(b 0 ) is greater or less than 1 for b near b * 0 . Let w L (b 0 ) and w R (b 0 ) denote the left and right eigenvectors of The Fredholm alternative implies the right side is orthogonal to w L (b * 0 ) so that Left and right eigenvectors of J(b * 0 ) associated with eigenvalue 1 have the forms Since ∂ 0 b 0 r = 1 and, by the product and chain rules, we find from (8) that It follows that the numerator and denominator in (A3) equal −1 and 1 respectively. Thus ∂ 0 b 0 λ(b * 0 ) = −1 and For b 0 − b * 0 ≈ 0 it follows that λ < 1 for b 0 > b * 0 and the forward bifurcating positive equilibria are therefore stable while λ > 1 for b < b * 0 and the backward bifurcating positive equilibria are unstable.