A bifurcation theorem for Darwinian matrix models and an application to the evolution of reproductive life-history strategies

We prove bifurcation theorems for evolutionary game theoretic (Darwinian dynamic) versions of nonlinear matrix equations for structured population dynamics. These theorems generalize existing theorems concerning the bifurcation and stability of equilibrium solutions when an extinction equilibrium destabilizes by allowing for the general appearance of the bifurcation parameter. We apply the theorems to a Darwinian model designed to investigate the evolutionary selection of reproductive strategies that involve either low or high post-reproductive survival (semelparity or iteroparity). The model incorporates the phenotypic trait dependence of two features: population density effects on fertility and a trade-off between inherent fertility and post-reproductive survival. Our analysis of the model determines conditions under which evolution selects for low or for high reproductive survival. In some cases (notably an Allee component effect on newborn survival), the model predicts multiple attractor scenarios in which low or high reproductive survival is initial condition dependent.


Introduction
Individual life-history traits or strategies affecting survival and reproduction, and how they affect population level dynamics and basic issues such as population survival or extinction, play a prominent role in life-history theory [21,29,33]. The question of why organisms adopt one strategy over another is a very old question that traces its roots back to Aristotle and Linnaeus [23]. In an influential paper from last century, L. Cole [6] focussed on the strategies of semelparity and iteroparity. Individuals in a semelparous population do not survive reproduction whereas those in an iteroparous populations have multiple reproductive episodes. Cole concluded that evolution should favour a semelparous strategy, a conclusion that became known as 'Cole's Paradox' since iteroparous populations are ubiquitous in nature. Numerous mechanisms were explored in response to this paradox and have been offered to account for the evolution of these reproductive strategies, including juvenile survival [8,31]; density dependence [1,5,18,22,34,35]; bet-hedging against environmental stochasticity [3,4,18,33,36,37,41], and many others [19,29,30,33].
Central to life-history theory are trade-offs [21,29]. Trade-offs are often viewed as an allocation issue, specifically the allocation of food or energy resources to metabolism growth, or reproduction, etc. For example, resources allocated to fertility are not available to apply towards post-reproductive survival and vice versa, in which case fertility and post-reproductive survival are negatively correlated.
Further complications are that the binary descriptions of semelparity and iteroparity are often not suitable for many populations and circumstances. Hughes [23] argues in favour of a 'continuum of modes of parity', so that all individuals in a population may not have the same life-history traits or have plasticity with regard to traits that are affected by environmental circumstances. With this point of view in mind, our goal here is to derive and analyse some evolutionary game theoretic models in which fertility and survival depend on continuous phenotypic traits that are subject to the Darwinian evolution and, using example models, to investigate the evolutionary stability or instability of different life-history strategies. We will focus on the trait dependence of two key features: population density effects on fertility and a trade-off between inherent fertility and post-reproductive survival. In addition to the standard negative effects of population density on fertility and/or survival, our model allows for positive effects on newborn survival (Allee component effects [7]).
In Section 2, we derive a discrete time Darwinian model for dynamics of a population and a mean phenotypic trait on which fertility, survival, and density effects depend. In Section 4, we analyse the equilibrium and stability properties of the model using some theorems proved in Section 3 concerning equilibrium existence and stability for a general Darwinian model. The results in that section extend the local bifurcation and stability results for non-evolutionary model equations for the dynamics of structured populations found in [16] to an evolutionary setting. They also generalize the results for the Darwinian, structured population dynamic models studied in [10,17] and [28] by allowing for a more general bifurcation parameter μ. In these earlier results, either a linearly appearing parameter or the inherent population growth rate (which does not appear explicitly as a model parameter) is used as the bifurcation parameter, both of which constrain the applicability of these otherwise general theorems. Our results in Section 3 allow the use of any model parameter as a bifurcation parameter. In Section 4, we apply these theorems to the model derived in Section 2 and, in addition, supplement the analysis with some numerical simulation studies. The biological punch lines obtained from this analysis are discussed in Section 5.

A discrete-time Darwinian model with a fertility-survival trade-off
Let x(t) denote the density of a population, at times t = 0, 1, 2, , 3, . . ., subject only to birth and death processes and in the absence of immigration or emigration. We take the unit of time to be an individual's maturation period and build a population dynamic model by accounting for the population at time t + 1 by only two means: newborns (individuals not present at time t) and those reproducing individuals at time t that survive to time t + 1. We assume that an individual's fertility is dependent upon available resources which are allocated between the production of newborns and processes that promote post-reproductive survival. Let f denote the fraction allocated to reproduction (0 ≤ f ≤ 1), b > 0 denote the per capita production of newborns per unit resource (per unit time), and let s n (0 < s n ≤ 1) denote the probability that a newborn born at time t survives to the next census t + 1. Then the density of newborns at time t + 1 equals bfs n x(t). In this model, we will assume that the fertility/survival trade-off is described by a probability of post-reproductive survival given by the formula s a (1 − f ) where s a (0 ≤ s a ≤ 1) is the maximal possible survival probability. Then is the population growth rate.
In addition to the fertility/survival trade-off, we want to study the effects that population density has on fertility and survival. We do this by introducing density factors into the growth rate where β, α and γ describe the effects that population density has on fertility, newborn survival, and post-reproductive survival, respectively. We assume H1: β, α and γ ∈ C 2 ( 1 ,R 1 HereR 1 + is the set of nonnegative real numbers and 1 is an open interval of real numbers that containsR 1 + . The reason for the normalizations at x = 0 is so that b, s n and s a retain their biological interpretations as inherent rates (i.e. rates in the absence of density effects).
A familiar factor used throughout the literature to describe negative density effects (i.e. a decreasing function of x > 0) is We refer to the coefficient c as the intraspecific competition coefficient since it accounts for the effect on an individual's fertility due to the presence of other individuals in the population (due, say, to competition for resources). We are also interested here in possible positive density effects (at least a low densities), which are called component Allee effects. Many mathematical expressions for such a factor can be found throughout the recent literature. In Section 4.2, we use for this purpose. (This expression seems to have been first introduced into population models by Jacobs [24].) The scalar difference equation together with (1), represents a general model for the dynamics of a population with the fertility versus survival trade-off and density-dependent features described above. Note that if all resource allocation goes towards reproduction, i.e. f = 1, then post-reproductive survival equals 0. In this case, the population is strictly semelparous in that no individual survives reproduction. We added the adjective 'strictly' since a population whose individuals have a very small post-reproductive survival rate could, for all practical purposes, still be considered semelparous. Only if f is significantly larger that 0 would the population be clearly designated as iteroparous, the cut-off point being subjectively chosen by the modeller. See [23] for a treatment of the difficulties with the binary classification of semelparity and iteroparity. Whatever language is used, we consider the fraction f in the model (4)-(1) to be a basic life-history strategy, measuring the fraction of resources allocated to fertility versus post-reproductive survival.
To introduce evolutionary changes in the components of the model (4)-(1) by means of evolutionary game theory, we consider the components in the per capita growth rate r to be associated with an individual in the population, at least one of which depends on a quantified phenotypic trait v of the individual. In this methodology, it is assumed the trait v is, at all times, normally distributed throughout the population with a constant variance. It follows that the trait distribution is determined by the population mean trait u. We allow model coefficients in r to depend on both the trait v of the individual and, in order to capture the influence of other individuals, on the population mean trait u. Thus, we write r = r (x, u, v). Evolutionary game theory describes the coupled dynamics of x and u by the Darwinian equations [27,40] Here ln r is taken to be the measure of population fitness. The coefficient σ 2 ≥ 0 is the speed of evolution (it is proportional to the assumed constant variance of the trait throughout the population). If σ 2 = 0 (there is no trait variability), then u(t) = u 0 remains constant (i.e. there is no evolution) and the Darwinian equations reduce to the population Equation (4). A study of the asymptotic dynamics of (5) and (6) begins with a consideration of equilibria (x e , u e ) with x e ≥ 0. These are found by solving the equations We distinguish two types of equilibria: extinction equilibria (x e , u e ) = (0, u 0 ) where u 0 solves and positive equilibria with x e > 0. These represent extinction and survival equilibrium states, respectively. Note that the first equilibrium equation for positive equilibria reduces to 1 = r(x, u, v)| v=u . Also of interest, of course, are the stability properties of equilibria. A stable extinction equilibrium represents a threat to survival whereas a stable positive equilibrium (x e , u e ) represents the survival of a population with density x e > 0 and mean trait u e > 0.
In the next section, we prove some equilibrium existence and stability results for a general class of Darwinian models (which allow for structured populations and multiple evolving traits). Existence and stability of equilibria depend, of course, on the model coefficients and the approach we take is that of bifurcation theory, using any chosen parameter appearing in the population equation as a bifurcation parameter. The stability of extinction equilibria and the creation of positive equilibria upon destabilization of extinction equilibria are studied in Section 3. In Section 4, we apply these results to a model derived from (1)-(5)-(6) that focuses on the question of the evolutionary selection of post-reproductive survival.

Preliminaries
Let R p denote the space of real p-dimensional (column) vectors, R p + denote the cone of positive vectors, andR p + denote the closure of R p + (i.e. the cone of nonnegative vectors). Letx ∈R m + denote the demographic vector of a population whose individuals have been classified into m different classes such as age, size, life cycle stages, disease state, etc. Let v = col(v i ) ∈ R n denote a n-dimensional column vector of phenotypic traits which are subject to Darwinian evolution (whose axioms are variability, heritability, and differential fitness). The phenotypic trait v i is assumed to have a heritable component and be normally distributed throughout the population, with a constant variance, throughout time. The trait distribution is therefore determined by its population mean trait u i . Letû = col(u i ) ∈ R n denote the vector of population mean traits.
A discrete time, matrix model for the dynamics of the structured population is defined by a population projection matrix P = [p ij ] whose entries describe class specific, per capita (individual) rates of fertility, survival, and transitions of individuals between classes. The projection matrix advances the demographic population vector from one census time t to the next t + 1 by means of matrix multiplication. In density-dependent models the entries p ij are functions ofx and in Darwinian models they are also functions ofv and/orû, so that we write P(μ, where μ denotes a model parameter that we have designated as a bifurcation parameter.
We make the following assumption about P and its entries p ij . Let μ be an open interval in R 1 for the bifurcation parameter μ and n an open set in R n for the trait vectorsû andv. Let the domain D of p ij be an open sent in μ × R m × R n × R n that contains μ ×R m + × n × n . Let0 m and0 n denote the zero vectors in R m and R n , respectively.
The matrix P(μ,0 m ,û,v) is called the inherent or intrinsic projection matrix, since entries describe fertility, survival and transition rates in the absence of density effects. The assumption A1 is designed so that its entries p ij (μ,0 m ,û,v) have interpretations as inherent (i.e. density free) fertility, survival, and transition rates. For this reason, they are assumed independent of the population mean traitû and to depend only on μ andv, as the vanishing of the gradient indicates.
The equations for both the population and mean trait dynamics provided by evolutionary game theory are [27,40] where M = (σ ij ) is a (symmetric) n × n variance-covariance matrix for trait evolution and is the gradient of fitness ln r(μ,x,û,v) with respect to the n-vectorv. Here we use the notation for partial differentiation. The diagonal entry σ 2 i := σ ii of M is the variance of the trait v i from the mean u i throughout the population at all times and σ ij = σ ji for i = j is the covariance of the i th and j th traits. Introducing the simplifying notations we re-write Equations (9)-(10) aŝ These equations form an m × n system of difference equations whose state variable is the pair (x,û) ∈ m × n . If (x e ,û e ) ∈ m × n is a constant solution (equilibrium) for whichx e ∈ R p + , then we refer to (x e ,û e ) as a positive equilibrium. An equilibrium (x e ,û e ) = (0 m ,û e ) is called an extinction equilibrium. The existence and stability properties of equilibria depend on the designated parameter μ. We refer to [μ, (x e ,û e )] as an equilibrium pair. If (x e ,û e ) is a locally asymptotically stable (or unstable) equilibrium of Equations (11) and (12) then we say the equilibrium pair [μ, (x e ,û e )] is stable (or unstable).
In addition to the assumption A1 on the projection matrix P, we will make use of the following assumptions on the variance-covariance matrix M.
Under assumption A2, the pair [μ, (x e ,û e )] is an equilibrium pair of Equations (11) and (12) if and only if μ, x =x e , and u =û e satisfy the algebraic equationŝ Our goal is to solve these equations for (x,û) = (x e (μ),û e (μ)) as a function of μ. We first consider extinction equilibria.

Existence and stability of extinction equilibria
Clearlyx = 0 solves the first equilibrium Equation (13) for all μ andû. In order for (x,û) = (0 m ,û e ) to be an extinction equilibrium for a value of μ, the trait componentû =û e must satisfy the equation From this assumption, we get the following lemma.

Lemma 3.1:
Under assumptions A1-A3, the Darwinian Equations (11) and (12) have an In order to investigate the stability of an equilibrium (x e ,û e ) by means of the linearization principle we consider the eigenvalues of the Jacobian associated with Equations (11) and (12) evaluated at (x e ,û e ). We write the (m + n) × (m + n) Jacobian in block form where J is the m × m Jacobian of P(μ,x,û,û)x with respect tox ; is the m × n matrix whose n columns are ϒ is the n × m matrix whose i th row is the transpose of the gradient is the n × n matrix I n + MH(μ,x,û,û) where I n is the n × n identity matrix and H is the n × n matrix

Remark 3.1:
A mathematical consequence of assumption A1 is that a partial derivative with respect to a component u k ofû of any entry p ij (μ,0 m ,û,v) of the inherent projection matrix P(μ,0 m ,û,û), and hence of r(μ,0 m ,û,v) (or their partial derivatives with respect to v i ), is identically equal to 0 for all μ,û, andv. For example When evaluated at an extinction equilibrium the Jacobian is block triangular (0 n×m is the n × m matrix of zeros) and thus its eigenvalues are the eigenvalues of the diagonal blocks P(μ,0 m ,û 0 ,û 0 ) and I n + MH(μ,0 m ,û 0 ,û 0 ). Note that by A1 which we see is the Hessian of fitness ln r(μ,x,û,û) with respect tov evaluated at (μ,x,û,û) = (μ,0 m ,û 0 ,û 0 ).
Putting Lemmas 3.2 and 3.3 together, we obtain the following theorem.

Existence and stability of positive equilibria
Our goal in this section is to solve the Equations (13) and (14) for positive equilibria and study their stability properties as they depend on μ. We will do this using local bifurcation theory in a neighbourhood of the extinction equilibrium (x e ,û e ) = (0 m ,û 0 ) and for μ near μ 0 given by assumption A4. Our approach is to solve the trait Equation (14) for u as a function of x and μ, substitute the result into the population Equation (13), and solve the resulting equation for x as a function of μ by means of results in [16]. The proof of the following theorem appears in the Appendix.
and η(ε),ẑ(ε) andξ(ε) are twice continuously differentiable on the interval |ε| < ε 0 witĥ When ε = 0, the equilibrium described in Theorem 3.2 equals the extinction equilibrium (0 m ,û 0 ) and for that reason the branch of equilibrium pairs in Theorem 3.2 is said to (transcritically) bifurcate from the extinction equilibrium at μ = μ 0 . Sincew R ∈ R m + we see that the bifurcating equilibrium pairs are positive for ε > 0. Whether they exist for μ greater than or less than (but near) μ 0 depends on the sign of κ/d. If κ/d > 0 and the bifurcating positive equilibrium pairs exist for μ > μ 0 , we say the bifurcation is forward. If κ/d < 0 and the bifurcating positive equilibrium pairs exist for μ < μ 0 , we say the bifurcation is backward. If the positive equilibrium pairs are stable or unstable, then we say the bifurcation is stable or unstable, respectively. The next theorem, whose proof is given in the Appendix, describes the stability properties of the bifurcating positive equilibrium guaranteed by Theorem 3.2 and relates them to the direction of bifurcation.  Remark 3.2: Theorem 3.3(a) is false without the assumption that P(μ 0 ,0 m ,û 0 ,û 0 ) is primitive. It is well known that if P(μ 0 ,0 m ,û 0 ,û 0 ) is imprimitive then the direction of bifurcation does not in general determine the stability properties of the bifurcating positive equilibrium [17,38,39]. In the imprimitive case, the bifurcation in general involves the bifurcation of positive periodic cycles and other invariant sets in addition to positive equilibria; the bifurcation is not well understood except for lower dimensional cases and some other special cases.

Evolutionarily stable strategies (ESS)
Another concept in evolutionary population dynamics is that of an evolutionary stable strategy or trait (an ESS) by which is meant a trait that cannot be replaced by a population with a mutant trait. For populations modelled by the Darwinian Equations (5) and (6), the ESS maximum principle [40] states that the trait component u e of a positive equilibrium (x e , u e ) is an ESS if and only if (x e , u e ) is a (locally asymptotically) stable equilibrium of (5) and (6) and, in addition, a global maximum of the adaptive landscape L(x e , u e , v) := ln r(x e , u e , v) occurs at v = u e . Thus, that an equilibrium (x e , u e ) of the Darwinian equations is (locally asymptotically) stable does not necessarily imply that u e is an ESS trait. The trait Equation (6) implies that u t moves in the uphill direction on the landscape L(x t , u t , v), and the trait equilibrium Equation (7) implies that v = u e is a critical point on the landscape L(x e , u e , v). It is possible, however, that v = u e is located at a local, but not global maximum on L(x e , u e , v) (see [13] for an example). It is also possible that v = u e is located at a minimum of L(x e , u e , v), as the example in Section 4.1 illustrates. That this can happen is explained by the fact that the adaptive landscape L(x t , u t , v) changes with time; see Figure 2 in Section 4.1.

Two Darwinian models with reproductive-survival trade-offs
A central role in studies of life-history theory is played by reproductive strategies and central to those are trade-offs between reproductive effort and post-reproductive survival. The allocation of more resources towards higher reproductive effort often comes at the cost of fewer resources placed into processes affecting post-reproductive survival. An extreme case is semelparity, the life-history strategy in which individuals reproduce only once in their life time and die afterward. In an influential paper, Cole [6] argued that evolution should favour this strategy which, because iteroparity (possessing more than one reproductive episode per lifetime) is pervasive throughout nature, became known as Cole's paradox. Many critiques of Cole's argument and resolutions of the paradox (i.e. explanations for how iteroparity can arise by evolutionary or environmental processes) have been offered. These are based on a variety of mechanisms and circumstances such as fertility-survival tradeoffs, nonlinear density effects, fluctuating environments and bet-hedging, spacial effects, and combinations of these [29,33]. Here we investigate the first two of these mechanisms by using model equations of the form (1)-(5)- (6).
The goal is to determine what circumstances lead to equilibrium states with low postreproductive survival and which lead to high post-reproductive survival. As we will see, there is no general answer to this question, even for the relatively simply low dimensional model (1)-(5)- (6). Which strategy is selected by evolution is model dependent, specifically, is significantly dependent on the type of density factors used and on the manner in which they depend on the evolving trait.
We illustrate the diversity of possible outcomes by means of two examples selected to illustrate two different model predictions. Both models are based on an m = 1 dimensional population equation and a single n = 1 phenotypic trait. The first example is an evolutionary version of the discrete logistic (Beverton-Holt) equation with a type of density trait dependence called symmetric. The second example includes an Allee component with a hierarchical (or asymmetric) trait dependence. In both examples, we assume that fertility has a Gauss-type distribution with respect to the trait v Thus, individuals with trait v = 0 attain maximal fertility b, but as a result have a postreproductive survival rate of 0 (i.e. are semelparous).
We also assume that the effect that population density has on an individual's fertility and survival depends on both the individual's trait v and the population mean trait u. This assumed dependence on u reflects the fact that density effects involve intraspecific interactions. More specifically we assume (as is often done in Darwinian models [40]) that the density factors depend on the difference v−u, that is to say, the effects that intraspecific interactions have on an individual depends on how different its trait v is from the typical individual in the population. Finally, in these examples, we use the maximal fertility rate as a bifurcation parameter μ = b.

A logistic model
In (1)-(5)-(6), we take density factors β(x) ≡ 1 and This describes a situation in which the only density effects in play are negative effects on newborn and adult survival probabilities, both in the same way. An expression commonly used for the trait dependence of the competition coefficient is [40] c This implies that the maximal value of the intraspecific competition coefficient, namely c 0 , occurs when v = u, i.e. when the individual's trait is equal to that of a typical individual it is likely to encounter, as represented by the population mean. Note that this modelling assumption implies c(v − u) symmetric in the difference v−u. Under these modelling assumptions, we have and A1 and A2 are easily seen to hold on the domain D = R 1 + × x × R 1 × R 1 where x is the half line x > −1/c 0 . The Darwinian Equations (1)-(5)-(6) become A3 holds forû 0 = 0 and (x, u) = (0, 0) is an extinction equilibrium for all values of b > 0. It is not difficult to see, for all initial conditions x 0 ≥ 0, that bs n − s a ≤ 0 implies x(t + 1) < s a x(t) and hence that x(t) → 0 as t → ∞. Therefore, we assume from now on that bs n > s a . and hence A4 holds. Finally, from r(b, 0, 0, 0) = bs n we find that A5 holds with μ 0 = b 0 = 1/s n (since d = s n > 0). Theorems 3.1 implies that the extinction equilibrium (0, 0) loses stability as μ = b increases through b 0 := 1/s n . Theorem 3.3 in turn implies, provided σ 2 is sufficiently small, the forward bifurcation of stable positive equilibria for b > b 0 .
We can obtain more results for (20) and (21) by noting the scalar difference Equation (21) is uncoupled from Equation (20). Straightforward manipulations of inequalities show that for all b > s a /s n provided σ 2 < 1/w 0 . It follows, under this condition on σ 2 , that |u(t + 1)| < |u(t)| and consequently all solutions (21) tend to 0 as t → +∞. Thus, for any initial condition u(0), the Equation (20) is asymptotically autonomous with limiting equation For this (discrete logistic) equation It is well known that x = 0 is globally asymptotically stable if bs n < 1. If bs n > 1 the positive equilibrium x = (bs n − 1)/c 0 is globally asymptotically stable for x(0) > 0. Theorems dealing with asymptotically autonomous difference equations (for example see Theorem 4 in [9]) imply the following results. (20) and (21). A take away message from this theorem is that, under the modelling assumption made in this example, the only possible positive equilibria have mean trait component u = 0, i.e. are semelparous equilibria, which is in accord with Cole's assertion [6]. However, is u = 0 and ESS trait? That is to say, does the adaptive landscape

Theorem 4.1: Consider the Darwinian model Equations
From these, we see that w 1 w 0 > bs n − s a bs n − 1 (22)   Figure 1(b)).
implies the landscape has a local minimum at v = 0 and in this case the equilibrium trait is not an ESS. Near the bifurcation point bs n 1 the semelparous equilibrium is an ESS, but if bs n is large enough that (22) holds then it is not an ESS. See Figure 1 for numerical simulations that illustrate this loss of ESS status. In Figure 1, it is perhaps counter-intuitive that the trait component of a positive equilibrium is located at a minimum on the adaptive landscape, since the trait equation in a Darwinian model implies that u(t) moves uphill on the landscape. The explanation lies in the fact that the landscape changes over time, which can (and does in Figure 1(b)) allow for uphill movement that ends in a valley. See Figure 2.

A model with an allee component
Our second example is a Darwinian model (1)-(5)-(6) that is again designed to investigate trait dependent fertility-survival trade-offs and density effects, but under different assumptions from those made in Section 4.1. This model will, unlike the model Section 4.1, allow for ESS traits that yield high post-reproductive survival (iteroparity). We consider a model the incorporates density effects in adult fertility β(x) and newborn survival α(x) (and ignores density effects on adult survival γ (x) ≡ 1). Specifically, we assume population density affects adult fertility in a negative way, again given by the factor (2) β (x) = 1 1 + cx but affects newborn survival in a positive way by using a factor α(x) of the form (3), namely α (x) = 1 + ρax q 1 + ax q , a ≥ 0 and 1 < ρ < 1/s n where we take q ≥ 1 to be an integer. This factor is the fractional increase of the inherent newborn survival rate s n as a function of increasing population density. Such a so-called component Allee effect can be the result of numerous mechanisms [7], an example of which is increased protection of newborns from predators by adults. Under this assumption with a > 0, newborn survival probability increases from s n to ρs n as population density x increases without bound. We will refer to a as the Allee coefficient since it measures how quickly newborn survival increases with x, as can be seen by the fact that newborn survival reaches the midpoint s n (1 + ρ)/2 at x = 1/a 1/q ; this midpoint is reached more quickly for large values of a > 0. In the Darwinian model (5)-(6) equations, we assume in this example that both the competition and Allee coefficients are trait dependent. In contrast to the example in Section 4.1, we assume the competition coefficient has an asymmetric trait dependence [25]. Specifically, we take This assumption can be viewed as describing a hierarchical dependence of intra-specific competition on the trait in that individuals with larger traits are subject to less intraspecific competition. Since c(0) = c 0 , the coefficient c 0 measures the effect of competition experienced by the typical individual with trait v = u.
Finally, we assume the Allee coefficient is an increasing function of the population mean trait u so that newborns of all traits are (equally) benefitted by a population with larger mean trait. Specifically, in this example, we take a = a 0 u 2 . In summary, these assumption imply roughly that as an individual's trait v > 0 increases it obtains higher post-reproductive survival probably in trade-off with decreased inherent fertility, which is to some degree compensated for by both a lower competition coefficient and a higher Allee coefficient. These would suggest the possibility of positive equilibrium states that have trait component different from 0 (in contrast to the example in Section 4.1) and hence the possibility of significant post-reproductive survival.
With these modelling choices, we have in the Darwinian Equations (5) and (6). We again choose μ = b as the bifurcation parameter to be used in Theorems 3.2 and 3.3, for which we need to verify assumptions A1-A5.
First, we note that if bs n ≤ s a then r(b, x, u, v)x ≤ s a x for all x ≥ 0 and for all v, u. Since s a < 1 it follows in this case that the population goes asymptotically extinct regardless of the trait dynamics, i.e. x(t) → 0 as t → +∞ for all initial conditions x(0) and u(0). Therefore, we assume from now on that bs n > s a .
Under this added assumption, we turn our attention to assumptions A1-A5. With the parameter interval μ taken to be the half line μ = b > s a /s n and the trait interval n = 1 to be R 1 , it is straightforward to verify that assumptions A1, A2, and A3 hold if σ 2 > 0 (evolution occurs) and that (x e , u e ) = (0, 0) is an extinction equilibrium for all b ∈ μ . The 1 × 1 Hessian is nonsingular for b ∈ μ which implies assumption A4 holds. With assumptions A1-A5 verified, we obtain the following theorem from Theorems 3.1(a), 3.2, and 3.3(a). (5) and (6) with (23), under the assumption that σ 2 is sufficiently small, the extinction equilibrium (x e , u e ) = (0, 0) destabilizes as b increases through b 0 = 1/s n and a forward bifurcation of (locally asymptotically) stable positive equilibria occurs.

Remark 4.1:
We can obtain an upper bound for how small σ 2 must be in Theorem 4.2 from a consideration of the proof of Theorem 3.1 in the Appendix. The assumption of small σ 2 is only needed so as to insure that the eigenvalue values of lower right block I n + MH(μ,0 m ,û 0 ,û 0 ) in the block diagonal Jacobian J(μ,0 m ,û 0 ,û 0 ) are all less than 1 in absolute value. Then it is the eigenvalues of the upper left block P(μ,0 m ,û 0 ,û 0 ) that determine stability (by linearization). For the Equations (5) and (6) Figure 3 illustrates the bifurcation described in Theorem 4.2, for selected parameter values, when b is greater that 1/s n (hence the extinction equilibrium is unstable) but is close to 1/s n . In the x, u-phase plane, the orbit with initial condition (x(0), u(0)) = (5, 0) is shown approaching the equilibrium (x 1 , u 1 ) = (0.096, 0.093), which is near the extinction equilibrium (0, 0) as predicted by Theorem 4.2. However, Figure 3 also illustrates another phenomenon, namely that in this sample case there exists a second positive equilibrium (x 2 , u 2 ) = (6.789, 1.086), which is approached by the orbit with initial condition (x(0), u(0)) = (6, 0). This example case in Figure 3 suggests a multi-attractor scenario in which the asymptotic, equilibrium states of orbits are initial condition dependent. This is not precluded by Theorem 4.2 since that theorem only concerns a local bifurcation of positive equilibrium near the extinction equilibrium.
Also shown in Figure 3 are the adaptive landscapes at each of these two equilibria and how the trait component of each equilibrium lies at a global maximum of its   Figure 3 except that b has been changed to b = 3.2. Since b < b 0 = 1/s n = 10/3 the extinction equilibrium (x 1 , u 1 ) = (0, 0) is stable, as is illustrated by the orbit with initial condition (x(0), u(0)) = (8, 0). However, the orbit with initial condition (x(0), u(0)) = (9, 0) approaches a positive equilibrium (x 2 , u 2 ) = (6.036, 1.065), whose trait component is an ESS as the accompanying graph of the adaptive landscape L(v) shows.
respective adaptive landscape. Thus, both equilibria have ESS trait components. Note that the post-reproductive survival probability s a (1 − f (u 1 )) = 0.008 at the equilibrium (x 1 , u 1 ) is quite low, while it is significantly higher at the equilibrium (x 2 , u 2 ), namely s a (1 − f (u 2 )) = 0.623. Therefore, in this sample case, whether evolution selects for a low or high post-reproductive survival rate is initial condition dependent. Figure 4 shows another interesting phenomenon that can occur in the model (5)-(6) with (23). It shows two-phase plane orbits in a case when b < 1/s n and the extinction equilibrium is stable. One orbit does approach the extinction equilibrium, but the other approaches a positive equilibrium. This suggests that a strong Allee effect [7] occurs in this example, which means survival is initial condition dependent. The trait component of the positive equilibrium is an ESS (see Figure 4) and yields a high post-reproductive survival of s a (1 − f (u 2 )) = 0.611. Strong Allee effects arise in population models most commonly due to a backward bifurcation [11]; however, in this case, the bifurcation is forward.
A rigorous mathematical treatment of positive equilibria outside a neighbourhood of the bifurcation point remains an open problem. In the next Section 5, we offer a plausible, but unproved, explanation of the multi-attractor cases suggested by Figures 3 and 4.

Concluding remarks
In Section 4, we applied the general bifurcation theorems in Section 3 for Darwinian matrix models to a low dimensional model (5) and (6) designed to study the evolution post-reproductive survival and to address Cole's paradox that evolution should favour semelparity. The models focus on two trait dependent mechanisms: inherent fertilityversus-survival trade-offs and density effects. In both models, fertility is assumed to have a Gaussian distribution with respect to the phenotypic trait v.
The model studied in Section 4.1 assumes logistic density effects on newborns and adults and a trait dependent intra-specific competition coefficient that is a symmetric function of the different v−u. This model supports Cole's assertion in that it predicts that the only survival equilibria have an ESS semelparous trait. In contrast, the model studied in Section 4.2 includes an Allee component in the density effect on juvenile survival and a trait dependent intra-specific competition coefficient that is a asymmetric (monotone) function of the different v−u. The fundamental bifurcation that occurs upon destabilization of the extinction equilibrium creates positive equilibria with low post-reproductive survival (Theorem 4.2). These equilibria again support Cole's contention that evolution favours semelparity. However, simulations of the model show that it is possible for the trait component of positive equilibria lying outside a neighbourhood of the bifurcation point to correspond to a high post-reproductive survival rate, which runs counter to Cole's contention (Figures 3 and 4). These simulations also suggest that these iteroparous equilibria can occur simultaneously with either semelparous equilibria or with the extinction equilibrium, in multiple attractor scenarios. While it remains an open mathematical question to establish these phenomena rigorously, a hint towards a possible mathematical explanation can be found in the bifurcation analysis of the non-evolutionary version of the model.
When the mean trait u is held fixed and not allowed to evolve (σ 2 = 0), the two Equations (5) and (6) reduce to just the scalar difference Equation (5) whose (positive) equilibrium equation is For any fixed mean trait value u, a plot of x as a function b as defined by this equation, yields a equilibrium bifurcation diagram. An example is shown in Figure 5, using the same parameter values as in Figures 3 and 4, but with u fixed at u = 1 (which is near the equilibrium trait component for the positive equilibrium obtained from the Darwinian model used in Figures 3 and 4). The bifurcation diagram in Figure 5 shows the destabilization of the extinction equilibrium x = 0 as b increases through b 0 = 3 + e/3 = 3.91 and multiple positive equilibria on an interval of b values less than b 0 where x = 0 is stable (as in Figure 3) and an interval of b values greater than b 0 on which there exists a multiple positive equilibria where x = 0 is unstable (as in Figure 4). It is known from general theorems that the positive equilibria on the decreasing segment in the bifurcation diagram are unstable, but that those on increasing segments, at least near the two tangent bifurcations, are (locally asymptotically) stable [12,14]. The hysteresis loop in this bifurcation diagram is due to the interplay of the negative and positive density effects described in the fertility and Allee factors β and α, respectively. While this is not the bifurcation diagram for the Darwinian model (5)-(6)- (23) in which u dynamically evolves, it does suggest that perhaps these same alternates do occur in the Darwinian model, as is suggested by Figures 3  and 4, and is due to the presence of a hysteresis in the bifurcation diagram. The results in Section 4, and those in other studies of other low dimensional Darwinian models [13,15,20], show that the answer to the question of whether evolution favours semelparity or iteroparity is much more complicated than Cole's simple assertion that it favours semelparity. In these models, whether evolution selects in favour of low or high post-reproductive survival depends crucially on the details of how the evolving trait affects the birth versus survival trade-off and on the specifics of the density effects, with some circumstances favouring one reproductive strategy and others favouring the other strategy. Moreover, this selection can be initial condition dependent. Other complications are introduced by a nonlinear trait dependent inherent fertility versus survival trade-off [20] and multiple peaks in the fertility distribution f (v) [13,20], neither of which we considered here, although there is evidence for them in nature [26,32]. evaluated at (x e ,û e ) = (0 m ,û 0 ) the first term vanishes by assumption A1. This implies that the Jacobian of ∇ˆvr(μ,x,û,û) with respect toû is nonsingular, when evaluated at μ = μ 0 and (x,û) = (0 m ,û 0 ), equals the Hessian H(μ,0 m ,û e ,û e ) and hence, by A4, is nonsingular. The implicit function theorem then implies that there exists a solutionû =υ(μ,x) of Equation (14), which is twice continuously differentiable in a neighbourhood of μ = μ 0 andx =0 m , that satisfiesυ(μ 0 ,0 m ) = u 0 . Substituting this solution into Equation (13), we are then tasked to solve the equation x = P μ,x,υ μ,x ,υ μ,x x forx ∈ R m + as a function of μ. To do this we apply Theorem 3.1 in [16] to the equation whereP(μ,x) = [p ij (μ,x)] is defined byp ij (μ,x) := p ij (μ,x,υ(μ,x),υ(μ,x)). To apply this theorem it is required that: (a)P(μ,x) is nonnegative and irreducible; (b) the entries inP(μ,x) are twice continuously differentiable in a neighbourhood of (μ,x) = (μ 0 ,0 m ); (c)P(μ 0 ,0 m ) = P(μ 0 ,0 m ,û 0 ,û 0 ) has dominant eigenvalue value r(μ 0 ,0 m ,û 0 ,û 0 ) = 1; (d)w T L ∂ μP (μ 0 ,0 m )w R = 0.
Following the same calculations found in the proof of Theorem 3 and Lemma 1 in [17], we obtain (19).

Proof of Theorem 3.3.:
To obtain the stability properties of the positive equilibrium pairs in Theorem 3.2 by the linearization principle, we consider the Jacobian J(μ,x,û,û) evaluated at the