Towards a replicator dynamics model of age structured populations.

We present a new modelling framework combining replicator dynamics, the standard model of frequency dependent selection, with an age-structured population model. The new framework allows for the modelling of populations consisting of competing strategies carried by individuals who change across their life cycle. Firstly the discretization of the McKendrick von Foerster model is derived. We show that the Euler-Lotka equation is satisfied when the new model reaches a steady state (i.e. stable frequencies between the age classes). This discretization consists of unit age classes where the timescale is chosen so that only a fraction of individuals play a single game round. This implies a linear dynamics and individuals not killed during the round are moved to the next age class; linearity means that the system is equivalent to a large Bernadelli-Lewis-Leslie matrix. Then we use the methodology of multipopulation games to derive two, mutually equivalent systems of equations. The first contains equations describing the evolution of the strategy frequencies in the whole population, completed by subsystems of equations describing the evolution of the age structure for each strategy. The second contains equations describing the changes of the general population's age structure, completed with subsystems of equations describing the selection of the strategies within each age class. We then present the obtained system of replicator dynamics in the form of the mixed ODE-PDE system which is independent of the chosen timescale, and much simpler. The obtained results are illustrated by the example of the sex ratio model which shows that when different mortalities of the sexes are assumed, the sex ratio of 0.5 is obtained but that Fisher's mechanism, driven by the reproductive value of the different sexes, is not in equilibrium.

Introduction

Among the most important approaches to the modelling of evolutionary processes are life history optimization and evolutionary games. Classical life history theory (Stearns 1992; Roff 1992) relies on optimization models, where there are no interactions among individuals and no density dependence:

“Life history evolution usually ignores density and frequency dependence. The justification is convenience, not logic, or realism” (Stearns 1992).

On the other hand, in classical game theoretic models there is no age or stage structure. Payoffs describe the averaged lifetime activity of an individual, which can be found for example in Cressman (1992):

“...an individual’s strategy is fixed over its lifetime or, alternatively, the life history of an individual is its strategy.”

Thus the synthesis of these perspectives can be very fruitful for theoretical insight (McNamara 2013). Methods used in life history optimization are closely related to classical demographic methods such as Bernadelli–Lewis–Leslie matrices (Caswell 2001). However, how to construct a general description of the relationships between demographic structure and population dynamics is still an unsolved problem (Caswell 2011). More precise than matrix models are continuous approaches arising from Lotka’s renewal equation (Lotka 1911, Diekmann et al. 2020a, b) and the McKendrick von Foerster model (McKendrick 1926). The combination of demography with a game theoretic perspective focused on frequency dependent selection, advocated by McNamara (2013), can be very useful since demographers are interested in the patterns produced by heterogeneity in the populations (Vaupel et al. 1979; Vaupel and Yashin 1983; Hougaard 1984; Vaupel and Yashin 1985). The game theoretic structure can explain the mechanisms shaping those patterns. The first papers combining both approaches are Garay et al. (2016) devoted to the particular biological problem of sib cannibalism, Li et al. (2015) and Lessard and Soares (2018) containing the approach incorporating age structure into a matrix game. These results show that after introduction of the age structure, matrix notation becomes very complicated and makes analysis difficult even in the case of two competing strategies and few age classes. In addition, previous works do not study the interplay between game dynamics and demographic structure in detail, assuming a fixed demographic structure. However, the game interactions described by demographic payoffs should affect the demographic structures of subpopulations of carriers of different strategies. In addition Li et al. (2015) assumes that payoffs are described by a standard payoff matrix, thus the same actions performed in different ages/stages will generate the same payoffs. However, we can expect that outcomes of individual actions may vary for different ages due to different experiences and physical condition of the playing individuals.

Another problem is that game theoretic models operate in abstract terms of costs and benefits measured in units of fitness mostly without deeper insight into their meaning or interpretation. This problem was analyzed in Argasinski and Broom (2012) where relationships between classical demography and evolutionary games are described in detail. This approach was later clarified in Argasinski and Broom (2018a, 2018b) by definition of the vital rates (birth and death rates) as the product of the interaction rates, describing the distribution of interactions (game rounds) in time and demographic game payoffs describing the number of offspring and the probability of death during a single interaction. The main conclusion there is that instead of excess above average fitness, models should be described explicitly by mortality and fertility, which are basic opposite forces shaping population dynamics (Doebeli et al. 2017). These results are significant progress in ecological realism, emphasizing the role of background mortality and fertility or the turnover of individuals (Argasinski and Kozłowski 2008). However, that approach is still very primitive. Mortality is described as an exponential decay of the population, which implies that the length of an individual’s lifetime is potentially unbounded, and there is no aging and no age specific payoffs. The goal of this paper is to fill this gap and develop a mathematical structure combining selection of individual strategies with an age structured population which will allow us to overcome the problems arising from increasing complexity of the models shown in Li et al. (2015) and simplifications ignoring the age dependence of payoffs resulting from certain actions and feedbacks driving the interplay between game dynamics and demography (the fixed age structure assumption). For practical reasons we will develop a high dimensional ODE system consisting of relatively simple equations, which can be generated by a simple loop and solved in every popular numerical platform.

The classical approach to evolutionary games and replicator dynamics

In the following subsections, we describe the state of the art in relation to our problem. A list of existing (and indeed new, see later) parameters are described in Table 1. Traditionally, in evolutionary game theory the payoff obtained by the jth strategy is proportional to its Malthusian growth rate and the dynamics of selection of strategies is described by the replicator dynamics (Maynard Smith 1982; Cressman 1992; Hofbauer and Sigmund 1988, 1998; Weibull 1995; Nowak 2006; Broom and Rychtár 2013; McNamara and Leimar 2020). We can derive this by rescaling the Malthusian equations for competing strategies to relative frequencies (where and w is the number of strategies), which leads towhere is the average payoff in the population. However, instead of the Malthusian parameter describing the payoff we can explicitly consider the individual fertility and mortality of a j-strategist. The explicit distinction between fertility and mortality was proposed also by Doebeli et al. (2017) as the cornerstone of a mechanistic model of natural selection. Note that in real life organisms are involved in different types of interactions with others or elements of the environment. Game theoretic models are focused on the outcomes of the particular interactions (such as fights as in the Hawk Dove game) responsible for selection of the analyzed trait or type of behaviour. These can be described by average demographic outcomes per interaction and and these focal interactions will occur at the rate . Other interactions, not related to the analyzed trait, can be described by average fertility and mortality , occurring at rate .

Table 1

List of important symbols

n -population size

n_j-number of individuals carrying the jth strategy

τ-interaction rate

f_j (fji)-fertility payoff of the jth strategy (of the jth strategy at age i )

s_j (sji=1-τdji)-survival payoff of the jth strategy (of the jth strategy at age i )

f¯_j, f¯^i, f¯-average fertility for the jth strategy, ith age class, whole population

s¯_j, s¯^i, s¯-average survival for jth strategy, ith age class, whole population

r¯=τr~-Malthusian parameter, product of the interaction rate and the game payoffs

m+1-number of age classes

w-number of strategies

K-carrying capacity, maximal population load

aji-frequency of individuals at age i among j-strategists

p_j-frequency of j-strategists in the population

a^i-proportion of individuals in the ith age class

pji-frequency of j-strategists in the ith age class

P_j-sex ratio strategy (fraction of males in the brood of the female)

x (x_j)-number of females (carrying the jth strategy)

y (y_j)-number of males (carrying the jth strategy)

G_j=x_j+y_j/∑l=1wx_l+y_l-frequency of the jth strategy gene

P=y/(x+y) -secondary sex ratio (proportion of males)

P¯_pr-primary sex ratio (average strategy of females)

Sjf=∑l=cdajl(1-Mjl)-proportion of active females among the jth strategy carriers

Sjm=∑l=abajlMjl-proportion of active males among the jth strategy carriers

M_j (Mji)-sex ratio of the population of the jth strategy carriers (of jth strategy carriers at age i)

k-number of offspring in the brood of a female

Mjop=Sjm/Sjm+Sjf-operational sex ratio of jth strategists

P_op=S¯^m/S¯^m+S¯^f-operational sex ratio in the population

Products of interaction rates and demographic payoffs will constitute the respective vital rates: game fertility rate and mortality rate , background fertility rate and mortality rate . Later the focal game interaction rate can be set to 1 by timescale adjustment and the background fertility and mortality rates become and . In addition we can add density dependent juvenile recruitment survival (Argasinski and Kozłowski 2008; Argasinski and Broom 2012, 2018a, b). To do this we should multiply fertilities by the logistic suppression coefficient (where the carrying capacity K is interpreted as the maximal environmental load, Hui 2006). Since fertilities but not mortalities are so scaled, the turnover of generations will not be suppressed at the equilibrium as it is in the classical logistic model (which leads to an immortal and childless population at equilibrium K). This gives the following variant of the replicator equations:where and , the details of which appear in Argasinski and Broom (2012), Argasinski and Broom (2018a, 2018b).

It was shown (Argasinski 2006) that every single population system described by the replicator Eq. (1) can be divided into the product of subsystems describing the dynamics in arbitrary chosen disjoint subpopulations (described by a frequencies , where , for the j-th subpopulation) and an additional system describing the dynamics of proportions between those subpopulations . This is useful when indiviuals differ not only by strategies but also by another second trait such as sex, age or developmental stage. Then, for example, we can decompose the population into subpopulations of carriers of the same strategy and describe the dynamics of the second trait among them. Then the dynamics in each subpopulation will have the form (1) and will depend on the excess of the strategy payoff from the average payoff in this subpopulation. Therefore, the same operation can be carried out for Eq. (2), and we obtain the system:where and are the fertility and mortality, respectively, of the i- th type (such as age or sex) in the subpopulation of the j-th strategy carriers, and are the mean fertility and mortality, respectively, in the subpopulation of the j-th strategy carriers and and are the respective values in the global population.

Note that we can decompose the initial population with respect to the second trait and describe the dynamics of strategic composition among individuals in the same age or sex class. Then the equations will describe variables , and n.

The classical approach to the modelling of age structured populations

Now we focus on age structured models (age classes will be indexed by superscripts). The classical approach to the modelling of age structured populations is related to Bernadelli–Lewis–Leslie matrices (Bernadelli 1941; Lewis 1942; Leslie 1945; Charlesworth 1994; Caswell 2001), following the matrix equation:where there are age classes, is the size of the ith age class and is fertility and is survival, respectively, in this class. Thus and the transition between subsequent age classes is . When the time unit equals the time step between age classes the above system is a good model of age structure. This age-structured growth model suggests a steady-state, or stable, age-structure and growth rate. The growth rate can be calculated from the characteristic polynomial of the Bernadelli–Lewis–Leslie Matrix called the Euler–Lotka equation (Caswell 2001):where r is the intrinsic growth rate of the population and describes survival to age i. We note here that in reality r will not be an independent parameter, and moreover will change in time as the distributions of the sizes of age classes change. An equilibrium distribution over the age classes in turn will allow us to define r in terms of the other model parameters. A simple ODE generalization of this system with continuous time but discrete age structure can be obtained by application of the delayed differential equations (Caswell 2001) where survival rates may describe aggregated exponential survival between respective age classes (Diekmann et al. 2017). However this approach may not work if the mortality function depends on the actual population state (as in game theory). Here the mortality rate may be unknown since it will depend on the trajectory of the dynamics during the age class. Then we can consider the continuous time limit of an infinite number of infinitely small age classes where population structure becomes a function n(t, l) of time t and continuous age l describing a moment in the lifetime of an individual. Then we can imagine the Taylor expansion analogous to the transition equation describing a small time step dt leading to ageing dlwhere is the continuous time mortality rate (at age l) similarly to the game models but without the distinction between the focal game and the background interactions. Since we obtain the McKendrick von Foerster equationwhich should be completed by boundary conditions and initial age distribution n(0, l).

The paper structure

In this paper we derive the discretization of the McKendrick von Foerster model allowing for the derivation of frequency dependent models. This is motivated by the fact that the discretized approach can be easily numerically solved by basic ODE solvers from popular numerical platforms. Thus the developed methodology does not need advanced knowledge in numerical analysis. Another advantage is that it will be compatible with standard game theoretic notation based on matrix games. Using the derived discretization we build two approaches to modelling selection among competing strategies with life cycles in an asexual population. One is focused on the impact of age structures of strategies on selection, while the second shows the impact of selection dynamics on the age structure of the whole population. The models obtained are generalized to mixed PDE-ODE models with continuous non-discretized age structures to outline the direction of future development. This framework is illustrated by a sex ratio example combining the two approaches, allowing us to model the sexually reproducing population. Our intention is to build a simple ready to use modelling methodology which can be extended in the future. However, we believe that even after clarification of the PDE based approach and development of simple solvers of coupled integro-differential PDE-ODE systems, our approach will still be useful for practical reasons arising from the simplicity of the methods based on matrix payoffs, which are much simpler to derive than continuous payoff functions. Thus it can be, for example, easily used for building initial toy models.

Results

Presenting the McKendrick von Foerster model as a system of ODE’s

In this section we will build the submodel describing the age structure dynamics of a subpopulation of carriers of some strategy competing with other strategies. Demographic vital rates will be outcomes of interactions between carriers of different strategies, interpreted as rounds of evolutionary games as in Argasinski and Broom (2018a). Thus as in replicator dynamics models we have the state of the population described by strategy frequencies but for each strategy subpopulation we have a respective age structure described by parameters (frequencies of individuals of age i among j-th strategy carriers). Demographic payoffs determining the vital rates will depend not only on the strategy frequencies as in the classical replicator models but also on the age of the opponents, thus the set of vectors of the age structures for all strategies (where describes the age structure of the j-th strategy carriers subpopulation) should be another argument of the payoff functions.

A major technical difference between the McKendrick von Foerster model and replicator dynamics is that the first is a PDE (or system of PDE’s as for example in Rundnicki and Mackey 1994) and the second is a system of ODE’s. The simple combination of both approaches will lead to a mathematically elegant but technically intractable system due to the lack of a general theory for mixed PDE-ODE systems and software for numerical computation. This methodology should be developed in the future, however before that, we need a useful approach based on existing solutions. To solve this problem we can approximate the continuous system by a large number of ODE’s describing unit interval age classes consisting of all individuals of age from a to . The discrete structure will allow us to use standard matrix payoff functions. The chosen time unit should be as long as possible to reduce the number of equations. Since we want to model frequency dependent selection, the mortality and fertility payoffs will depend on the trajectory of the population state. Therefore we cannot use simplified delayed differential equations since we do not know the trajectories during the time delay interval. Instead we can assume that the unit of a timescale described by interaction rate is short enough that the changes of the population state are small enough with respect to the population size (e.g. 50 births in a population of 30000), that the resulting changes of frequency dependent birth and death rates will be negligible.

Following “Appendix A” we see that Eq. (10) can be discretized and approximated by the replicator dynamics (see Fig. 1 for the intuitive presentation of the discretization scheme for frequency dependent vital rates). In particular for the jth strategy we describe the system in frequencies and a scaling parameter n. Assume that is the game payoff component of the growth rate (then ) and is the respective averaged value. If the growth rates are nearly constant, then for the chosen timescale described by interaction rate changes of the strategy frequencies during a single time unit are (where ), thus they are sublinear. Here should be as big as possible to minimize the number of equations, but small enough that payoff function arguments (and similarly others) should change their values only slightly (i.e. and are small enough, but not necessarily infinitesimal), so that the resulting changes of and are negligible. Then the discretization is acceptable and we obtain:where and the Malthusian parameter describing the growth of the jth strategy isIt is important that age class survival describes aggregated outcomes of the game rounds occurring during a time unit. Therefore it is distinct from the survival probability of a single round which should be used in trade-off functions when only survivors of the game round can reproduce (Argasinski and Broom 2012, 2018a, b), leading to fertility . In addition, due to nearly linear behaviour within a single time unit the system (11,12) is equivalent to the large Leslie matrix (7) with survival and then parameter describes the fraction of individuals that played the single game round. Parameter always acts as the multiplier of game payoffs and (thus the resulting survival rate is ). Since in the next sections we will focus on the derivation of the dynamics, where the structure of the vital rates is not so important, for simplicity we can incorporate the interaction rate into the birth and death vital rates and skip it in the notation. Therefore, below, will be hidden inside functions and which will be interpreted as the vital rates.

Fig. 1

Schematic presentation of the discretization of the continuous age dynamics. The assumed unit time step between age classes is associated with a change of the population state, which may induce change of the frequency dependent payoffs. However, while the resulting changes of the vital rates are negligible, values of payoffs can be approximated by their initial values at the beginning of the transition between age classes

Assume the absence of density dependence. Since the r.h.s. of our system (11) is the negative function of , the following attracting nullcline manifold exists (for constant mortalities this is an attracting steady state):Note that will satisfy the general form for in Eq. (14). In addition the Euler–Lotka equation is satisfied (for a derivation and proof, see “Appendix B”). In the density dependent case, the age structure attractor (14) will change with the growth of the population. Now we can use the derived submodel for derivation of the full model.

The extension to multipopulation replicator dynamics

Now we can incorporate the above model into a multipopulation evolutionary game (Argasinski 2006). Recall that we have w strategies and age classes indexed from 0 to m. Assume that p describes the strategy (phenotype) fraction and a describes the frequency of the age class. As before, and describe, respectively, the fertility and survival of the j-strategist in age class i. Two perspectives are possible (see Fig. 2):

Fig. 2

The difference between two alternative formulations of the problem: system a describes the evolution of the gene pool according to age structures of carrier subpopulations, system b describes the evolution of the global age structure driven by strategy selection in age classes

(a) Firstly we consider the impact of the age structure in sub-populations strategically homogenous on selection of the strategies, denoted as system . This can be described by coordinates:

for the age structure of the j-strategists

the strategy frequencies in the whole population,

where and .

(b) Secondly we consider how selection within each age class affects the overall age structure, denoted as system . It can be described by coordinates:

for strategy frequencies in age class

the age structure of the population,

where and .

Thus in both cases we will have a core system describing the whole population (strategic composition in and age structure in ) completed by the respective subsystems describing the age structure of the subpopulation of strategy carriers (for ) or the strategic age class composition (for ).

Now we describe the transition of coordinates between the formulations. First we define the auxiliary canonical coordinates without subclasses:Now following Argasinski (2006) we define transitions between systems:

to :and to :Now let us derive systems of equations operating in both coordinate systems. In the following we use the within group averaging terms:We also use two global averages, which can each be written in two ways:For system we have the following system of differential equations (see “Appendix C”):givingFor system we have (see “Appendix D” for a detailed derivation):The expanded form of the above system will beNote that Eq. (33) is equivalent to (25) and in both casesRecall that for simplicity we assumed that functions act as the vital rates with interaction rate hidden inside. When we insert it back it would appear as and . In contrast to the basic replicator Eq. (2), parameter cannot easily be removed from systems and by simple timescale adjustment. A similar situation occurs with the background payoff components and , which simply cancel out in (2) but are still present in the population size equation. This will not be the case for systems and . However, for simplicity, in this paper we do not deal explicitly with the background payoffs.

Mixed PDE-ODE versions of systems and

We can derive mixed PDE-ODE versions of systems and , where the age profile is a continuous function, which are simpler and more mathematically elegant. The advantage is that they are independent of the timescale since the interaction rate will simply cancel out (see “Appendix E” for derivations). Thus the previous simplifying assumption about skipping it is obsolete in this case. Payoffs and are now continuous functions of the lifetime l and the strategic composition at time t. In addition the distinction between aggregated age class survival and game round survival discussed below Eq. (13) is not necessary since PDE versions of both systems will be driven by game payoffs only. Therefore for system we havewith , , , and .

For system we havewith , , , and , where and .

A sex ratio example

Now we will show how the methods presented in the previous sections can be used to extend the simpler age independent model to the age dependent case and how they can be used to model a sexually reproducing population. We will show this methodology by example of the synthetic sex ratio model (Argasinski 2012, 2013, 2017) combining simple explicit genetics (similar to the more advanced approaches as in Karlin and Lessard 1986) with rigorous strategic analysis. We will use the formulation of the model focused on selection of genes encoding sex ratio strategies (Argasinski 2013). Below we outline the basic details of this model. The introduction of the life cycle perspective to theoretical studies on the sex ratio is important, since data show the huge impact age specific mortalities can have on the dynamics of age specific sex ratios (for example see Orzack et al. 2015 for data showing the changes of the human sex ratio from conception to death).

We have a population consisting of x females and y males. All of them are carriers of a single gene encoding one from a finite number w of competing sex ratio strategies which are expressed by females (strategy is carried by females and males and describes the fraction of male newborns in the brood of a female). Then the population state can be expressed by the population’s sex ratio , primary sex ratio (average strategy of females) and vectors G and M where:

the gene frequencies,

the sex ratios in the carrier subpopulations.

Then and . Strategy genes are inherited from mother or father with probability 0.5. The sex specific payoff functions describe the impact of direct reproductive success (offspring of the focal female or offspring of partners of the focal male) and the per capita normlized contribution of the same strategy carriers of the opposite sex. Therefore, the payoffs of male and female carriers and the average gene carrier are:where k is the number of offspring per female. The average payoffs are:We can obtain the system describing the dynamics of gene frequencies and the sex ratios in the carrier subpopulations:leading to the following system of equationsThe above system can be regarded as an example of multi-level selection since the fate of a gene is determined by the actual composition of the carrier subpopulation described by the carriers’ sex ratio and the threshold between growth and decline is the adult sex ratio . The parameters are determined by the Tug of War dynamics (52) describing the impact of female carriers producing newborns according to the carried strategy and randomly drawn female partners of male carriers producing newborns according to the average strategy of females .

The extension of the sex ratio model to the age structured case

We will extend this system in the following way (see Fig. 3).

Fig. 3

The extension of the phase space of the sex ratio model to the age structured case. The gene pool phase space is completed by respective subspaces describing the age structures among carriers of the particular genes, as in system . Then each age structure subspace is completed by subspaces describing carriers’ sex ratios, according to system

System will be applied to extend the gene pool dynamics to the system with explicit age structure for each subpopulation of carriers (described by for the j-th gene) of the particular gene. This means that each Eq. (51) will be transformed to the form (21) and completed by the respective subsystem (20) describing the age structure of the subpopulation of carriers of the particular gene. In addition, for the age structure of each strategy we will apply system to describe the dynamics of the sex ratios within each age class. Thus for each strategy, the respective subsystem (20) will be the core subsystem (28) of system , and it will be completed by the respective subsystems (26,27), describing the dynamics of strategy carriers’ sex ratios in particular age classes. This structure will be the generalization of the Eq. (52) in the original model. Assume that survival, described by for females and by for males, depends only on sex and age. Males are active in the age classes from a to b and females from c to d, and fractions of sexually active female and male individuals carrying the j -th strategy areAnalogous parameters for the whole population areWe also have , and the primary sex ratio is:Thus this is the average strategy of active females describing the proportion of males among all newborns or zygotes. The operational sex ratio among active carriers of strategy j and the equivalent average value for the population isThe equations on G should be updated according to the additional assumptions on age limits of sexual activity (age classes from a to b for males and c to d for females). We should also derive the respective forms of per capita fertility payoffs described in the new coordinates. For derivation of the dynamics we need the following operational male fertility payoff of active males, average per capita gene fertility payoff and the average fertility in the whole population (the detailed derivation is in “Appendix F”):Note that describes the number of partners and the number of female carriers (“sisters”) of the average male carrier of the focal strategy gene. Therefore, the male operational fertility payoff describes the fertility of their partners with the average strategy and “sisters” carrying the same gene. The gene payoff describes the aggregated fertility of all female carriers and all partners of male carriers. Thus we will obtain the following general system derived in “Appendix G”:where describes the average survival of the carrier of the jth strategy determined by the actual carriers sex ratio. Then and . Thus the general Eqs. (20–22) have become Eqs. (60, 61, 64) through the sequences: (21) (49) (60) , (20) (28) (61), (26) (62), (27) (63). Figure 3 shows how the phase space of the original model was extended to the age structured case. After substitution of the payoff functions (see “Appendix G”) we obtain the system:

where average survival probabilities areare the fractions of sexually active females and males among the gene carriers, and the respective averages. Thus the selection mechanism is seriously altered by the age structure. The above system shows that differences in mortality between sexes and different ages of sexual activity can significantly affect the selection of individual strategies. Equation (67) contain the terms and describing the fractions of sexually active individuals and are the equivalent of the Tug of War dynamics (52). The dynamics of the age structure of each strategy is attracted bySex ratios among the jth strategy carriers of particular ages converge toNote that when we assume that there are no differences in survival probabilities between sexes () then the system (65–69) reduces to the simplified version. Equation (66) will be independent of parameters and the Eq. (68) will converge to a constant value over the whole life cycle ( for all i). Therefore all strategies will have the same age structure. In effect the bracketed term , describing the excess of the mortality payoff from average mortality, will vanish in Eq. (65) and selection of the genes will be driven by the excess fertility payoff bracket describing the Fisherian mechanism driven by the difference in reproductive value between the sexes;which is equivalent to the Shaw–Mohler formula (Shaw and Mohler 1953). If we assume that both sexes are mature in the same age classes then , we have that operational sex ratios (56) are and (). Therefore for the operational sex ratio the above formula equals zero for all strategies. When we additionally assume that the sex specific survivals for different ages are the same, the system replicates the results of the original model . In the general case if and then obviously the operational sex ratios (56) are equal. In other cases, the strategies with the greater fraction of the sex which is in the minority among active individuals (according to the operational sex ratio ) will have a greater value of (73). Since all individuals of the same sex suffer the same mortality, the values of parameters and are determined by the allocation of sexes at birth, determined by their encoded strategy. Due to the constant brood size k, an increase of female newborns leads to a decrease of male newborns and vice versa. Therefore, this allocation will determine operational sex ratios and selection should act accordingly to differences in operational sex ratios, similarly to (51).

We can see this in Fig. 4 depicting a numerical simulation for the case of three competing strategies , , with 25 age classes plus infant age class 0. For simplicity we assumed that age class survivals will be the same with only one change at some arbitrary age, different for males and females. For females we have survival 0.95 until age 10 and 0.80 in subsequent ages. For males we have 0.88 until age 15 and then 0.72 subsequently. By definition survival in the last age classes is zero. Females are fertile from age until age while males are active from age until age . The initial population size was with a carrying capacity K = 10 000. Initial conditions are , , and . We start from a very young population where adult age classes have frequencies 0.001 leading to a 0.025 proportion of non-infant individuals and sex ratios are and . These exaggerated conditions show the initial dynamics of the growing cohort leading to the interesting patterns depicted in Fig. 5 depicting the age structure and Fig 6 showing the dynamics of age specific sex ratios.

Fig. 4

Panel a dynamics of gene frequencies, panel b operational sex ratios, for strategies , , and primary and operational sex ratios of the population and , panel c population size. Trajectories show that is the threshold between growth and decline of the gene frequency depending on the value of . This is shown by the example of strategy 0.05, where bumps in the marked areas are caused by two types of events. The first is when the strategy’s operational sex ratio passes the population’s operational sex ratio , which is the threshold between growth and decline. The second is when the average operational sex ratio passes the value of 0.5 which inverts the strategic situation, since the opposite sex is in the minority when this happens

Fig. 5

Trajectories of age classes. The initial behaviour is caused by huge differences in the initial sex ratios. The assumed changes in age specific survivals slightly affect the trajectories

Fig. 6

Trajectories of age specific sex ratios. The pattern caused by the assumed changes in survival probabilities is clearly visible

Figure 7 shows the delayed convergence to the respective Euler–Lotka manifolds. At the global equilibrium excess fertility payoffs (73) (the difference between the payoff and the mean) are not equal to zero because they must balance the nonzero values of the excess survival payoffs for growth rates to be equal (Fig. 8). Therefore, the classical Fisherian equilibrium focused only on fertility payoffs is not reached here. A question arises about the interplay between fertility and survival and how it leads to the primary sex ratio of 0.5. In addition the operational sex ratio is far from 0.5. Therefore in this case the Fisherian mechanism is not enough to explain the origins of the primary sex ratio of 0.5. Figure 4 shows that the mechanism driven by the operational sex ratios still works but all values are rescaled, and we also have different mortalities for different strategies. The interplay between the Fisherian mechanism, driven by fertility and differences in reproductive value between the sexes, and age structure, driven by survival differences between the sexes, needs an explanation which will be the subject of future work.

Fig. 7

A plot of the convergence to the respective Euler–Lotka manifolds (dashed lines) for arbitrarily chosen age classes for strategy 0.05. The convergence is delayed by some inertia caused by the age dynamics

Fig. 8

Plots of the excess fertility payoffs , excess mortality payoffs and the gene frequency growth rates from the gene pool dynamics (65). Fertility payoffs are not equal as in the classical theory, and the same situation is true for mortality payoffs, but the right hand sides of the equations are zero. This shows that the explanation for the primary sex ratio being 0.5 needs an explicit consideration of the interplay between fertility and mortality

Discussion

In this work we presented a new modelling framework combining evolutionary dynamics with demographic structure. This approach can be a useful tool in the development of the synthesis between evolutionary game theory and life history theory. We started with the derivation of the ODE discretized approximation of the McKendrick von Foerster model of age structured populations and its critical manifold equivalent to the Euler–Lotka equation. This was extended to the explicit case of multiple competing strategies and transformed into two types of age structured replicator dynamics. The first focused on the selection of strategies when each strategy is described by a subsystem describing the dynamics of the age structure. The second focused on the age structure of the whole population and a subsystem of the strategies within each age class. These led to huge ODE systems which are equivalent to systems of Bernadelli–Lewis–Leslie matrices. Another complication is that for the discretized age structure we need age class survival functions which will describe the aggregated outcomes of all interactions (game rounds) that have happened during a single time unit. This survival function is distinct from game round survival which can be used for the derivation of the fertility-survival trade-off functions used in situations when only survivors of the interaction can reproduce (Argasinski and Broom 2012, 2018a, b). In addition we have outlined the PDE versions of the obtained systems to indicate a future direction of research.

Both approaches are combined in the illustrative example of a sex ratio model. This is an extension of the dynamic sex ratio model (Argasinski 2012, 2013, 2017). It shows that when we assume different mortalities for both sexes, the classical Fisherian explanation based on the differences of reproductive values of offspring is not enough to explain convergence to the primary sex ratio of 0.5. The excess fertility payoff does not converge to 0 which would be equivalent to an equal reproductive values for both sexes, but its non-zero value is equal to the value of the excess survival payoff. The question of how this mechanism works in detail should be explained in future research. The new model provides a theoretical framework that can be used to explain the mechanisms shaping the patterns observable in data collected on age specific sex ratios from conception to death, as in Orzack et al. (2015) and Orzack (2016).

The obtained results clearly show that a life cycle perspective plays a crucial role in evolutionary processes. In the classical approaches to evolutionary game theory individuals cannot change their properties during their lifetime. Thus their life history is a memoryless process, and survival of a single interaction does not change the state of the individual. This is caused by the fact that the classical approaches to evolutionary games are focused on the strategies interpreted as patterns of behaviour, not on the individual itself. The exception to this rule is the state based approach (Houston and McNamara 1999; Argasinski and Rudnicki 2020). The explicit description of the life cycle and the different payoffs at different ages leads to a more complicated game theoretic structure. In particular a mixed PDE-ODE approach will lead to more complex payoff functions based on continuous distributions of ages for different strategies. This will need more sophisticated methods, such as models with function valued traits (Oechssler and Riedel 2001; Dieckmann et al. 2006, van Veelen and Spreij 2009), state based games (Houston and McNamara 1991, 1999; Argasinski and Rudnicki 2020) or “large games” with a distinction between strategy sets and population states (Wieczorek and Wiszniewska 1998; Wieczorek 2004, 2005), as opposed to basic two person matrix games. The modelling framework proposed in this paper can also be a useful tool in the research on animal personalities The combination of game theoretic analysis with an explicit age structure will allow us to analyze the relationships between behavioural strategies (such as aggression or cowardice) and life history traits (such as allocation of energy into growth or reproduction). This is important because life history trade-offs are shaped by external mortality which is the outcome of interactions with the environment. On the other hand, the demographic outcomes of interactions such as mortality are affected by phenotypic traits such as growth shaped by life history strategies (Wolf and Weissing 2010; Wolf and McNamara 2012). This constitutes life-history-behavioural feedback.