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 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, b) 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 r j and the dynamics of selection of strategies is described by the replicator dynamics (Maynard Smith 1982;Cressman 1992;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ṅ j = n j r j to relative frequencies q j = n j /n (where n = w j=1 n j and w is the number of strategies), which leads tȯ wherer = w j=1 q j r j is the average payoff in the population. However, instead of the Malthusian parameter describing the payoff we can explicitly consider the individual fertility f j and mortality d j 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 Table 1 List of important symbols n -population size n j -number of individuals carrying the jth strategy τ -interaction rate f j ( f i j )-fertility payoff of the jth strategy (of the jth strategy at age i ) s j (s i j = 1 − τ d i j )-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 a i j -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 p i j -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) type of behaviour. These can be described by average demographic outcomes per interaction f j and d j and these focal interactions will occur at the rate τ f . Other interactions, not related to the analyzed trait, can be described by average fertility f b and mortality d b , occurring at rate τ b . Products of interaction rates and demographic payoffs will constitute the respective vital rates: game fertility rate τ f f j and mortality rate τ f d j , background fertility rate τ b f b and mortality rate τ b d b . Later the focal game interaction rate τ f can be set to 1 by timescale adjustment and the background fertility and mortality rates become In addition we can add density dependent juvenile recruitment survival (Argasinski and Kozłowski 2008;Broom 2012, 2018a, b). To do this we should multiply fertilities by the logistic suppression coefficient (1 − n/K ) (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: wheref = w j=1 q j f j andd = w j=1 q j d j , the details of which appear in Argasinski and Broom (2012), Argasinski and Broom (2018a, b).
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 q i j = n i j /n j , where n j = i n i j , for the j-th subpopulation) and an additional system describing the dynamics of proportions between those subpopulations p j = n j / z n z . 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 f i j and d i j 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, are the mean fertility and mortality, respectively, in the subpopulation of the j-th strategy carriers andf andd 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 q i j , p j 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 m + 1 age classes, n i is the size of the ith age class and f i is fertility and s i is survival, respectively, in this class. Thus n 0 (t + 1) = i n i (t) f i and the transition between subsequent age classes is n i (t + 1) = s i−1 n i−1 (t). 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, agestructure 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 i−1 z=0 s z 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 dl n(t + dt, l + dl) = n(t, l) + ∂n ∂t dt + ∂n ∂l dl = s(l)n(t, l) where τ d(l) 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 dl = dt we obtain the McKendrick von Foerster equation which should be completed by boundary conditions n(t, 0) = ∞ 0 n(t, l)τ f (l)dl 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.

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 p j but for each strategy subpopulation we have a respective age structure described by parameters a i j = n i j /n j (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 p = [p 1 , . . . , p w ] as in the clas-sical replicator models but also on the age of the opponents, thus the set of vectors of the age structures for all strategies a = [a 1 , . . . , a w ] (where a j = [a 0 j , . . . , a m j ] 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 a + 1. 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 a i j = n i j / m z=0 n z j and a scaling parameter n.
is the respective averaged value. If the growth rates τr j (t) are nearly constant, then for the chosen timescale described by interaction rate τ changes of the strategy frequencies during a single time unit are , thus they are sublinear. Here τ should be as big as possible to minimize the number of equations, but small enough that payoff function arguments p j (and similarly others) should change their values only slightly (i.e. are small enough, but not necessarily infinitesimal), so that the resulting changes of τ f i j and τ d i j are negligible. Then the discretization is acceptable and we obtain: where a 0 j = 1 − m i=1 a i j and the Malthusian parameter describing the growth of the jth strategy is a)) describes aggregated outcomes of the game rounds occurring during a time unit. Therefore it is distinct from the survival probability of a single round 1 − d i j ( p, a) which should be used in trade-off functions when only survivors of the game round can reproduce Broom 2012, 2018a, b), leading to fertility (1 − d i j ( p, a)) f i j . In addition, due to nearly linear behaviour within a single time unit the system (11,12) is equivalent to the large Leslie matrix (7) a)) and then parameter τ describes the fraction of individuals that played the single game round. Parameter τ always acts as the multiplier of game payoffs f i j and d i j (thus the resulting survival rate is 1 − τ d i j ). 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 f i j and s i j which will be interpreted as the vital rates.
Assume the absence of density dependence. Since the r.h.s. of our system (11) is the negative function of a i j , the following attracting nullcline manifold exists (for constant mortalities s j this is an attracting steady state): Note thatâ 0 j will satisfy the general form forâ i j 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 m +1 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, f i j and s i j describe, respectively, the fertility and survival of the j-strategist in age class i. Two perspectives are possible (see Fig. 2): (a) Firstly we consider the impact of the age structure in sub-populations strategically homogenous on selection of the strategies, denoted as system S a . This can be described by coordinates: . 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 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 0 j , . . . , a m j for j = 1, . . . , w the age structure of the j-strategists p 1 , . . . , p w the strategy frequencies in the whole population, where a i j = n i j / z n z j and p j = z n z j /n . (b) Secondly we consider how selection within each age class affects the overall age structure, denoted as system S b . It can be described by coordinates: . . , m strategy frequencies in age class i a 0 , . . . , a m the age structure of the population, where p i j = n i j / z n i z and a i = z n i z /n . Thus in both cases we will have a core system describing the whole population (strategic composition in S a and age structure in S b ) completed by the respective subsystems describing the age structure of the subpopulation of strategy carriers (for S a ) or the strategic age class composition (for S b ). 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: S a to S b : and S b to S a : 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 S a we have the following system of differential equations (see "Appendix C"): For system S b we have (see "Appendix D" for a detailed derivation): The expanded form of the above system will bė Note that Eq. (33) is equivalent to (25) and in both cases Recall 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 τ f i j and In contrast to the basic replicator Eq.
(2), parameter τ cannot easily be removed from systems S a and S b 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 S a and S b . However, for simplicity, in this paper we do not deal explicitly with the background payoffs.

Mixed PDE-ODE versions of systems S a and S b
We can derive mixed PDE-ODE versions of systems S a and S b , 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 d j (t, l) and f j (t, l) 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 S a we have with a j (t, with a(t,

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(Argasinski , 2013(Argasinski , 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 P j ∈ [0, 1] is carried by x j females and y j 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 P = y/(x + y), primary sex ratio (average strategy of females)P pr = w j=1 x j x P j and vectors G and M where: the gene frequencies, M j = y j x j + y j the sex ratios in the carrier subpopulations.
Then P = w j=1 G j M j andP pr = w j=1 x j x P j = w j=1 G j (1 − M j ) 1 − P P j . 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 equationṡ The 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 M j and the threshold between growth and decline is the adult sex ratio P = w j=1 G j M j . The parameters M j are determined by the Tug of War dynamics (52) describing the impact of female carriers producing newborns according to the carried strategy P j and randomly drawn female partners of male carriers producing newborns according to the average strategy of femalesP pr .

The extension of the sex ratio model to the age structured case
We will extend this system in the following way (see Fig. 3).
System S a will be applied to extend the gene pool dynamics to the system with explicit age structure for each subpopulation of carriers (described by a i j for the jth 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 S b 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 S b , 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 M j Eq. (52) in the original model. Assume that survival, described by s i f for females and by s i m 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 are Analogous parameters for the whole population arē We also have P = w j=1 G j i a i j M i j , 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 is The 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 1 − P op /P op describes the number of partners and 1 − M op j /M op j the number of female carriers ("sisters") of the average male carrier of the focal strategy gene. Therefore, the male operational fertility payoff f op m describes the fertility of their partners with the average strategy and "sisters" carrying the same gene. The gene payoff f g 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":Ġ  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: 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 S a . Then each age structure subspace is completed by subspaces describing carriers' sex ratios, according to system S ḃ where average survival probabilities arē are the fractions of sexually active females and males among the P j 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 S m j and S f j 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 bŷ Sex ratios among the jth strategy carriers of particular ages converge tô Note that when we assume that there are no differences in survival probabilities between sexes (s i f = s i m ) then the system (65-69) reduces to the simplified version. Equation (66) will be independent of parameters M i j and the Eq. (68) will converge to a constant value over the whole life cycle (M i j = M 0 j for all i). Therefore all strategies will have the same age structure. In effect the bracketed term s j −s , 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 f g (P j , a, G, M) −f (a, G, M) 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 S f j = 1 − S m j , we have that operational sex ratios (56) are M op j = S m j and P op =S m (S f = 1 − P op ). Therefore for the operational sex ratio P op = 0.5 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 S f j =S f and S m j =S m 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 P op ) will have a greater value of (73). Since all individuals of the same sex suffer the same mortality, the values of parameters S f j and S m j 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 P 1 = 0.05, P 2 = 0.55, P 3 = 0.95 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 c = 8 until age d = 15 while males are active from age a = 8 until age b = 20. The initial population size was n(0) = 40 with a carrying capacity K = 10 000. Initial conditions are G 1 (0) = 0.9, G 2 (0) = G 3 (0) = 0.05, M 0 1 = 0.7 and M 0 2 = M 0 3 = 0.1. 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 M i 1 = 0.9 and M i 2 = M i 3 = 0.8. 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. 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.

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 Trajectories show that P op is the threshold between growth and decline of the gene frequency depending on the value of M op j . 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 M op j passes the population's operational sex ratio P op , which is the threshold between growth and decline. The second is when the average operational sex ratio P op passes the value of 0.5 which inverts the strategic situation, since the opposite sex is in the minority when this happens the derivation of the fertility-survival trade-off functions used in situations when only survivors of the interaction can reproduce 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(Argasinski , 2013(Argasinski , 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).

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 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 andSpreij 2009), state based games McNamara 1991, 1999;Argasinski and Rudnicki 2020) or "large games" with a distinction between strategy sets and population states (Wieczorek and Wiszniewska 1998;Wieczorek 2004Wieczorek , 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 Fig. 6 Trajectories of age specific sex ratios. The pattern caused by the assumed changes in survival probabilities is clearly visible 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 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. s j −s and the gene frequency growth rates f g −f (1 − n/K ) + s j −s 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 those events. However we can imagine a timescale where interactions are sufficiently rare, so that only a small fraction of individuals play a single round of the game. Then the dynamics is technically linear and survival is described by the first term of the Taylor series of the exponential function. Thus assume that time interval dt = dl = 1 is small enough that a small fraction τ d(t, l)dt of individuals will die due to aggregated outcomes of independent game rounds. The remaining 1− τ d(t, l)dt survivors will be moved to the next age compartment. Then from Eq. (9) we have that n(t + 1, l + 1) = (1 − τ d(t, l)) n(t, l) = s(t, l)n(t, l) which describes the move from point t, l to point t + 1, l + 1. Assume that during a unit interval all surviving individuals from age l will be moved to age l + 1 while all individuals from l + 1 will be moved from this age to the next age step or die. Therefore during a single time unit we have linear movement occurring at incoming rate s(t, l − 1)n(t, l − 1) and removal rate -n(t, l), since all individuals will be removed during a single time unit. Therefore this linear process can be well approximated by the first order Taylor expansion for t = 1 where the bracketed term can be interpreted as the first derivative; therefore the bracketed term constitutes derivative dn/dt, leading to Note that (74) can be presented in the form n(t + 1, l) = s(t, l − 1)n(t, l − 1) which leads to the Leslie matrix (7). When we change notation to the numbered age classes describing age increments and assume that aggregated survival rate s i (t) and fertility rate τ f i (t) can change in time, we obtain the systeṁ It is reasonable to assume that s 0 = 1, (other values are equivalent to s 0 = 1 with rescaled fertilities f i ) and s m = 0. However, to be compatible with the replicator dynamics and game theoretic machinery, the dynamics should be expressed in terms of phenotype frequencies. Let us change the coordinates to the frequencies a i = n i /n (with n = m i=0 n i ) describing the age structure. The system (75,76) can be presented in the form of the Malthusian equations: Therefore, this system can be presented as a system of frequency dependent replicator equationsȧ i = a i (r i −r ) and a single equation on the scaling parameterṅ = nr . Since m i=0 a i = 1 and s m = 0 we have the average Malthusian growth rate as Then we can formulate a system of frequency dependent replicator equations by transforming Eq. (78) To add density dependence we should multiply the fertility rate by the logistic suppression coefficient 1 − n K .
(2) Frequency and density dependence and the choice of time unit determining the discretization step When the above system describes the dynamics of the age structure of the subpopulation of carriers of some strategy competing with other strategies (indexed by subscripts) then the parameters s i j (t) = 1 − τ d i j ( p(t), a(t)) and τ f i j ( p(t), a(t)) (thus r j = τr j = τ m i=0 a i j f i j (1 − n/k) − d i j ) are game payoffs depending on strategy frequencies p j = n j / k n k (and their age distributions, but now we limit our reasoning to strategy frequencies p). We want to choose the longest possible time step to reduce the number of equations. However the frequencies will change in time, so the discretization step cannot be too big, since the aggregated payoffs during unit interval will depend on the changes of τ d i j ( p(t)) and τ f i j ( p(t)) during that time interval. The time unit should be short enough that these vital rates will not change significantly and the number of individuals will change nearly linearly within each age class. Thus we should analyze how much the strategy frequencies p j can change during unit time and how this affects the vital rates. For small t = 1 we have a change of n j = n j (t)τr j (t) t = i n i j (t)τ f i j ( p(t)) 1 − n K − d i j ( p(t)) t (positive or negative) for each j, leading to the change n = j n j = n(t)τr (t) t for the population size. Then p j (t + 1) = n j (t) n(t) + n + n j n(t) + n = n j (t) n(t) n(t) n(t) + n + n j n(t) + n = n j (t) n(t) 1 − n n(t) + n + n j n(t) + n = p j (t) + n j − p j (t) n n(t) + n .
Thus the pace of increment is nearly linear or slower since τ/(1 + τr (t)) < τ . We can set the timescale by adjusting parameter τ in the formulae f i j = f i j ( p(t) + p(τ )) − f i j ( p(t)) and d i j = d i j ( p(t) + p(τ )) − d i j ( p(t)) to make their values small enough that τ f i j and τ d i j are negligible. For the density factor we have that (1 − n(t + 1)/K ) = (1 − (n(t) + n)/K ) and the change of fertility rate is −τ f i j ( p(t)) n/K , thus it depends on τ/K and is negligible for even big time steps.

Appendix B: The stationary age distribution and the Euler-Lotka equation in the continuous case
From (80) (recall that τ is hidden in the fertilities f i ), for the non-infant age classes (i > 0) the stationary points for the age structure of this system are: (note the similarity to the Euler-Lotka equation). The stable age structure is a unique vector of frequencies among age classes, conditional on the average Malthusian growth rate of the population. Now let us prove the equivalence with the Euler-Lotka equation. After substitution of the stable age frequencies from Eq. (85) into Eq. (77) we obtain: Frequency equilibrium implies that per capita growth rates in all age classes are equal to the average growth rater (â). Thus equality will also be satisfied for the growth rate of the 0 age class, leading to which is the Euler-Lotka equation.

Appendix C: Derivation of system S a
We start from the Malthusian system describing the dynamics of age classes in the subpopulation of carriers of the j-th strategy: According to (11) the above system can be transformed into the frequency replicator dynamics of age classes: