Modelling the evolution of traits in a two-sex population, with an application to grandmothering

TL;DR

This paper introduces a mathematical model for the evolution of heritable traits affecting life-history timing in two-sex populations, demonstrating how grandmothering influences sex ratios and longevity through fitness landscape analysis.

Contribution

It provides a novel formula for the fitness landscape incorporating two sexes, age structure, and multiple traits, extending the Price Equation framework.

Findings

Grandmothering can increase sex ratios due to longer female post-fertile lifespan.

Longer lifespans lead to higher fitness for both sexes.

Grandmothering alone can drive evolution towards increased longevity.

Abstract

We present a mathematical simplification for the evolutionary dynamics of a heritable trait within a two-sex population. This trait is assumed to control the timing of sex-specific life-history events, such as the age of sexual maturity and end of female fertility, and each sex has a distinct fitness tradeoff associated with the trait. We provide a formula for the fitness landscape of the population and show a natural extension of the result to an arbitrary number of heritable traits. Our method can be viewed as a dynamical systems generalisation of the Price Equation to include two sexes{, age structure} and multiple traits. We use this formula to examine the effect of grandmothering, whereby post-fertile females subsidize their daughter's fertility by provisioning grandchildren. Grandmothering can drive a shift towards higher sex ratios due to lengthening female post-fertile longevity, leading to changes in fitness for both sexes. For males, increased longevity is accompanied by a substantially longer fertile lifespan resulting in higher sex ratios in the fertile ages. Our fitness landscapes show a net increase in fitness for both males and females at longer lifespans, and as a result, we find that grandmothering alone provides an evolutionary trajectory to higher longevities.

Tables

Table 1. Table 1: Parameter values

Symbol	Definition	Value
$m$	Benefit received by fertile female from grandmothering	$3$
${\tau }_T\left(L\right)$	Age of frailty	$\min (2L,75)$
${\tau }_f\left(L\right)$	Age of female sexual maturity	$L/2.5+2$
${\tau }_m$	Age of male sexual maturity	15
${\hat{\tau }}_f$	Age female fertility ends	45
${\hat{\tau }}_m\left(L\right)$	Age male fertility ends	$\min (2L,75)$ (equal to ${\tau }_T$)
$\varphi \left(L\right)$	Male fertility-longevity tradeoff	Function of $L$, see Eq. (9)
$\mu \left(L\right)$	Mortality rate	$1/L$
$\overline{b}\left(L\right)$	Birth rate	$4.522/L-0.023$

Equations

$$\frac{\partial u}{\partial t}+\frac{\partial u}{\partial a}=-\mu (a,x)u.$$

$$u(0,x,t)=\frac{1}{to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}{\text{leavevmode}\text{set@color}S}_f(t)to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}{\text{leavevmode}\text{set@color}S}_m(t)}{\int }_{-\infty }^{\infty }{S}_f(y,t){\int }_{-\infty }^{\infty }{S}_m(to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}\text{leavevmode}\text{set@color}z,t)N\left(\frac{to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}\text{leavevmode}\text{set@color}z+y}{2},x,{\epsilon }^2\right)dto0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}\text{leavevmode}\text{set@color}zdy.$$

$$S_f(x,t)$$

$$S_m(x,t)$$

$$u(0,x)$$

$$=\frac{2}{to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}{\text{leavevmode}\text{set@color}S}_{f}to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}{\text{leavevmode}\text{set@color}S}_m}{S}_f(x)\ast S_m(x)[2x]\ast N(x,0,{\epsilon }^2)[x],$$

$$S_f(x)$$

$$S_m(x)$$

$$F(x)={\int }_{{\tau }_f}^{{\hat{\tau }}_f}b(a,x)\exp \left(-{\int }_0^a\mu (s,x)ds\right)da,$$

$$M(x)={\int }_{{\tau }_m}^{{\hat{\tau }}_m}\varphi (a,x)\exp \left(-{\int }_0^a\mu (s,x)ds\right)da.$$

$$u_{n+1}(x)=2\left(\frac{F(x)u_n(x)}{⟨F(x),u_n(x)⟩}\right)\ast \left(\frac{M(x)u_n(x)}{⟨M(x),u_n(x)⟩}\right)[2x]\ast N(x,0,{\epsilon }^2)[x],$$

$$G(x)=\{\begin{array}{l}\begin{array}{rl}\frac{(x-23{)}^a}{b+(x-23{)}^c}& \text{if }x\ge 23,\\ 0& \text{if }x<23,\end{array}\end{array}$$

$$\varphi (L)=\exp \left({\int }_{L_0}^L\delta (s)ds\right),$$

$$k(L)=\left(\int_{\tau_{m}}^{\hat{\tau}_{m}}\exp\left(-\int_{0}^{a}\mu(s,L)\,ds\right)\,da\right)\Bigg{/}\left(\int_{\tau_{f}}^{\hat{\tau}_{f}}\exp\left(-\int_{0}^{a}\mu(s,L)\,ds\right)\,da\right),$$

$$u_{n+1}(x)=2\left(\frac{F(x)u_n(x)}{⟨F(x),u_n(x)⟩}\right)\ast \left(\frac{M(x)u_n(x)}{⟨M(x),u_n(x)⟩}\right)[2x].$$

$$\Lambda (u_{n+1}(x))=\frac{1}{2}\left(\Lambda \left(\frac{F(x)u_n(x)}{⟨F(x),u_n(x)⟩}\right)+\Lambda \left(\frac{M(x)u_n(x)}{⟨M(x),u_n(x)⟩}\right)\right).$$

$$u_n(x;{\mu }_n,\alpha )=\{\begin{array}{l}\frac{1}{2\alpha }\text{if}x\in [{\mu }_n-\alpha ,{\mu }_n+\alpha ],\\ 0\text{otherwise}.\end{array}$$

$${\mu }_{n+1}=\frac{1}{2}\left(\frac{{\int }_{{\mu }_n-\alpha }^{{\mu }_n+\alpha }xF(x)dx}{{\int }_{{\mu }_n-\alpha }^{{\mu }_n+\alpha }F(x)dx}+\frac{{\int }_{{\mu }_n-\alpha }^{{\mu }_n+\alpha }xM(x)dx}{{\int }_{{\mu }_n-\alpha }^{{\mu }_n+\alpha }M(x)dx}\right).$$

$$\mu_{n+1}=\Phi(\mu_{n}):=\mu_{n}+\frac{\alpha^{2}}{6}\left(\frac{\frac{d}{dx}(F(x)M(x))}{F(x)M(x)}\right)\Bigg{|}_{x=\mu_{n}}.$$

$$\frac{d}{dx}(F(x)M(x))\Big{|}_{x=\mu^{*}}=0.$$

$$\frac{d}{dx}\Phi(\mu^{*})=1-\frac{\alpha^{2}}{6}\left(\frac{\frac{d}{dx}(F(x)M(x))}{(F(x)M(x))^{2}}\right)\Bigg{|}_{x=\mu^{*}}+\frac{\alpha^{2}}{6}\left(\frac{\frac{d^{2}}{dx^{2}}(F(x)M(x))}{F(x)M(x)}\right)\Bigg{|}_{x=\mu^{*}}.$$

$$\frac{-12F(x)M(x)}{\alpha^{2}}<\frac{d^{2}}{dx^{2}}(F(x)M(x))\Big{|}_{x=\mu^{*}}<0.$$

$$u(0,\mathbf{x})=\frac{1}{to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}{\text{leavevmode}\text{set@color}S}_{f}to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}{\text{leavevmode}\text{set@color}S}_m}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\cdots {\int }_{-\infty }^{\infty }{S}_f(\mathbf{y})S_m(to0.0pt\text{leavevmode}\text{set@color}\mathbf{z}\text{hss}\text{leavevmode}\text{set@color}\overline{})N\left(\frac{\mathbf{y}+to0.0pt\text{leavevmode}\text{set@color}\mathbf{z}\text{hss}\text{leavevmode}\text{set@color}\overline{}}{2},\mathbf{x},{\epsilon }^2\right)d\mathbf{y}dto0.0pt\text{leavevmode}\text{set@color}\mathbf{z}\text{hss}\text{leavevmode}\text{set@color}\overline{},$$

$$N(\mathbf{x},\mathbf{y},{\epsilon }^2)=(2\pi {\epsilon }^2{)}^{-n/2}\exp \left(\frac{-(\mathbf{x}-\mathbf{y}{)}^T(\mathbf{x}-\mathbf{y})}{2{\epsilon }^2}\right).$$

$$u(0,\mathbf{x})$$

$$=\frac{2^n}{to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}{\text{leavevmode}\text{set@color}S}_{f}to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}{\text{leavevmode}\text{set@color}S}_m}{S}_f(\mathbf{x})\otimes S_m(\mathbf{x})[2\mathbf{x}]\otimes N(\mathbf{x},0,{\epsilon }^2)[\mathbf{x}],$$

$$f(\mathbf{x})\otimes g(\mathbf{x})[\mathbf{x}]={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\cdots {\int }_{-\infty }^{\infty }f(y_1,y_2,\dots ,y_n)g(x_1-y_1,x_2-y_2,\dots ,x_n-y_n)d\mathbf{y}.$$

$$S_f(\mathbf{x},t)$$

$$S_m(\mathbf{x},t)$$

$$u_{t+1}(\mathbf{x})=2^n\left(\frac{F(\mathbf{x})u_t(\mathbf{x})}{⟨F,u_t⟩}\otimes \frac{M(\mathbf{x})u_t(\mathbf{x})}{⟨M,u_t⟩}\right)[2\mathbf{x}],$$

$$\Lambda (u_{t+1}(\mathbf{x}))=\frac{1}{2}\left(\Lambda \left(\frac{F(\mathbf{x})u_t(\mathbf{x})}{⟨F,u_t⟩}\right)+\Lambda \left(\frac{M(\mathbf{x})u_t(\mathbf{x})}{⟨M,u_t⟩}\right)\right),$$

$$\Lambda (f(\mathbf{x}))=\left(\begin{array}{c}{\int }_{{\mathbb{R}}^n}{x}_{1}f(\mathbf{x})d\mathbf{x}\\ {\int }_{{\mathbb{R}}^n}{x}_{2}f(\mathbf{x})d\mathbf{x}\\ ⋮\\ {\int }_{{\mathbb{R}}^n}{x}_{n}f(\mathbf{x})d\mathbf{x}\end{array}\right).$$

$$u_t(\mathbf{x},{\mathbf{\mu }}_t,\alpha )=\{\begin{array}{l}\begin{array}{rl}& \frac{1}{(2\alpha {)}^n}\text{if }x\in A,\\ & 0\text{otherwise,}\end{array}\end{array}$$

Full text

11footnotetext: School of Mathematics and Statistics, University of Sydney, NSW 2006, Australia. M.Chan@maths.usyd.edu.au22footnotetext: Department of Anthropology, University of Utah, Salt Lake City, UT 84112, USA. Hawkes@anthro.utah.edu33footnotetext: School of Mathematics and Statistics, University of Sydney, NSW 2006, Australia. PKim@maths.usyd.edu.au

Key-words: Population dynamics, Evolutionary dynamics, Sexual conflict, Grandmother Hypothesis, Human evolution.

1 Introduction

The modelling of a two-sex population with heritable traits has often proved to be a formidable problem due to the interplay of many seemingly crucial evolutionary mechanisms. Our study is motivated by a specific example of this problem: the evolution of longevity in ancestral populations, where sex-specific life histories introduce biases in the mating sex ratio and each sex has a fitness tradeoff associated with longevity. The tradeoffs and life histories vary with longevity, which is expected from Charnov’s model of life history variation in female mammals [4], and when models include both males and females, result in a sexual conflict, wherein each sex is optimised at different longevities, but both optima cannot be obtained simultaneously [13, 3]. This leads to the selected longevity of the population being a compromise between the two optima. Although an observation shared by previous studies is that the sexual conflict plays a significant role in the evolutionary dynamics, it remains a major obstacle in gaining a deeper understanding of the problem as it has been difficult to explicitly quantify [8, 9, 11, 10, 2]. This further adds to the confusion on exactly which components of the problem generate the sexual conflict and how strongly such components contribute to the sexual conflict.

Instead of reasoning about this problem in a discrete sense, for example with an agent-based approach, we show that it is possible, and more convenient, to consider continuous distributions along trait-space. This approach allows for a simplified computation of the evolutionary dynamics and the explicit calculation of the fitness landscape of the population. Moreover, the male and female fitnesses across trait-space can be separately calculated, which provides a way to visualise the sexual conflict. Since the calculation of the fitness landscape is explicit and inexpensive, comparisons of fitness landscapes generated from different sex-specific life-history strategies and tradeoffs can be made. In Section 5.2, we generalise this result for an arbitrary number of heritable traits that affect the sex-specific tradeoffs and life-history parameters. A common approach to such problems is discrete formulations aiming to generalise Price’s equation (for examples see Batty et al. [1] and Grafen, A. [7]), which is a mathematical statement of how trait frequencies in a population change over time, given the initial population frequencies over trait space and the fitnesses of the traits (see Frank [6] for a detailed discussion and also Price’s original papers [15, 16]). We take an alternative dynamical-systems approach; thus, our result can be thought of as a dynamical-systems generalisation of Price’s equation to include two sexes, age structure and multiple traits.

In Section 3, we apply the result to determine whether grandmothering, whereby post-fertile females increase the fertility rate of their daughters by provisioning their grandchildren, can drive the evolution of increased longevity by allowing for a transition between equilibria in the trait-space. This question stems from the Grandmother Hypothesis, which proposes increased longevity was selected for when ancestral populations began relying on foods weaned juveniles could not acquire effectively for themselves. Under such circumstances, mothers would have to feed their offspring longer, but with subsidies from grandmothers, they could have next offspring sooner; and longer-lived grandmothers, without infants themselves, could support more grandchildren.

Whether grandmothering alone can propel a great ape-like population to higher longevities has been the subject of several modelling studies recently, many of which involve computationally expensive simulations [9, 11, 10, 14]. Similar to Kim et al. [11], Kim et al. [10] and Chan et al. [2], who considered the problem with comparable assumptions, we find here that grandmothering alone can drive the evolution of increased longevity. The key difference from previous models is that our approach produces an explicit fitness landscape which is inexpensive to compute. This allows for a straightforward comparison between fitness landscapes of both sexes and the population compromise that evolves with and without grandmothering.

2 Model

2.1 Problem setup

We consider a two-sex population in which each individual possesses a heritable longevity trait value $x$. Females are fertile from $\tau_{f}(x)$ to $\hat{\tau}_{f}(x)$, while males are fertile from $\tau_{m}(x)$ to $\hat{\tau}_{m}(x)$, where $\hat{\tau}_{f}(x)$ and $\hat{\tau}_{m}(x)$ denote the end of fertility for females and males respectively for a particular longevity trait value $x$. We make no assumptions on the form of $\tau_{f}(x)$, $\tau_{m}(x)$, $\hat{\tau}_{f}(x)$ or $\hat{\tau}_{m}(x)$, only that they are piecewise continuous functions. Both males and females have equal age-specific mortality rates $\mu(x)$, leading to equivalent age-profiles for both sexes.

In each time interval, fertile males compete over every fertile female for a chance at paternity. The competitive success rate for a male with a specific trait value $x$ is governed by the male longevity-fertility tradeoff function $\phi(x)$, representing the relative probability that a male will out-compete others for a chance at paternity. Likewise, the female birth rate at each time $t$ is governed by the female longevity-fertility tradeoff function $b(x)$, which accounts for varying interbirth intervals as longevity increases. Offspring inherit the mean value of their parents’ longevity trait values, with a probability to mutate according to a normal distribution with mean 0 and variance $\varepsilon^{2}$. We assume that there is equal probability for an offspring to be male or female, leading to an equal sex ratio. Finally, to ensure that the population converges to a finite equilibrium, we assume that the number of offspring with a specific trait value $x$ entering the population is regulated by a competition factor equal to the total number of births at time $t$.

2.2 Model formulation

We let $u(a,x,t)$ denote the density of individuals with age $a$, longevity $x$ at time $t$, and model the age and mortality dynamics via the McKendrick von-Foerster model

[TABLE]

The mating and mutation dynamics are addressed in the boundary condition, given by

[TABLE]

The boundary condition gives the total density of offspring with longevity $x$ entering the system at any time $t$. We give an explanation of its form below. The functions $S_{f}(x)$ and $S_{m}(x)$ respectively denote the density of births by fertile females with longevity $x$ and paternities by males with longevity $x$, that is,

[TABLE]

The functions \hbox to0.0pt{\hskip 0.85414pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorS}}_{f}(t):=\int_{-\infty}^{\infty}S_{f}\,dx and \hbox to0.0pt{\hskip 0.85414pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorS}}_{m}(t):=\int_{-\infty}^{\infty}S_{m}\,dx represent the total density of births and a measure of total male competitiveness at time $t$ respectively. The last factor in the integrand $N(x,\mu,\varepsilon^{2})$ denotes the normal density function with mean $\mu$ and variance $\varepsilon^{2}$; this factor governs the mating and mutation dynamics, where the parameter $\varepsilon$ in Eq. (2) denotes the mutation rate, which we assume to be small. The expression {S_{m}(\hbox to0.0pt{\hskip 0.04514pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorz}},t)}/{\hbox to0.0pt{\hskip 0.85414pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorS}}_{m}(t)} represents the probability of a female mating with a male with longevity \overline{\hbox{}}$$z, or equivalently, the probability of a male with longevity \overline{\hbox{}}$$z securing a mate. The total paternities by males with longevity \overline{\hbox{}}$$z is then given by S_{m}(\hbox to0.0pt{\hskip 0.04514pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorz}},t)S_{f}(y,t)/\hbox to0.0pt{\hskip 0.85414pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorS}}_{m}(t). Including the mutation factor N((\hbox to0.0pt{\hskip 0.04514pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorz}}+y)/2,x,\varepsilon^{2}) and integrating over $y$ gives the total density of offspring with longevity $x$ entering the system at any time. Finally, we divide by \hbox to0.0pt{\hskip 0.85414pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorS}}_{f}(t) to regulate the number of offspring entering the system, which serves to keep the total population density finite as described in Section 2.1. This process results in the boundary condition Eq. (2).

To simplify Eq. (2), we let \hbox to0.0pt{\hskip 0.04514pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorz}}=2z-y to obtain

[TABLE]

where $\ast$ denotes the convolution operator, i.e. $f(x)\ast g(x)[x]=\int_{-\infty}^{\infty}f(y)g(x-y)\,dy$.

Since the mortality and ageing dynamics in Eq. (1) and the evolutionary dynamics in Eq. (4) are on different timescales, we simplify the system by assuming that the population is always at a stable age distribution. Thus, the density of births and paternities for specific longevity values, given above by $S_{f}$ and $S_{m}$, can be simplified to

[TABLE]

For convenience, we define

[TABLE]

The evolutionary dynamics of the system, given in Eq. (4), is then simplified to the following integrodifference equation

[TABLE]

where the angled brackets denote the $L^{2}$ inner product $\langle f(x),g(x)\rangle:=\int_{-\infty}^{\infty}f(x)g(x)\,dx$. Since $\frac{F(x)u_{n}(x)}{\langle F(x),u_{n}(x)\rangle}$ and $\frac{M(x)u_{n}(x)}{\langle M(x),u_{n}(x)\rangle}$ both integrate to 1, is it instructive to view them as probability density functions. The term $2\left(\frac{F(x)u_{n}(x)}{\langle F(x),u_{n}(x)\rangle}\right)\ast\left(\frac{M(x)u_{n}(x)}{\langle M(x),u_{n}(x)\rangle}\right)[2x]$ can then be interpreted as another probability density function which lies “in-between” $\frac{F(x)u_{n}(x)}{\langle F(x),u_{n}(x)\rangle}$ and $\frac{M(x)u_{n}(x)}{\langle M(x),u_{n}(x)\rangle}$ assuming $u_{n}(x)$ is bell-shaped, see Figure 1a. Convolving with $N(x,0,\varepsilon^{2})[x]$, the term responsible for mutations, “smooths” this function. The convolution of probability density functions results in another probability density function, thus $u_{n+1}(x)$ integrates to 1 and can be viewed as a probability density function for all $n$. We provide an illustration of one hypothetical iteration of Eq. (7) in Figures 1a-1b, where one iteration encapsulates the mating, mutation and birth dynamics for one generation. The process illustrated in Figures 1a-1b is a continuous perspective on the discrete problem described in Section 2.1. Its main advantage is that the computation of convolutions of functions is much faster and more elegant than the agent-based approach of simulating discrete individuals. We note that this approach does take into account the biases in the mating sex ratio, even though there is no age-structure in Eq. (7).

The function $F(x)$ gives the expected number of births by a female with longevity $x$ in her lifetime, while $M(x)$ gives a value proportional to the fraction of available mating opportunities a male with longevity $x$ will have in his lifetime. We let $F(x)$ and $M(x)$ be measures of fitness for females and males respectively, although we note that $M(x)$ is technically not a measure of reproductive success, instead it is a measure of the expected proportion of matings a male with longevity $x$ will have in his lifetime. We note that a special case of $F(x)$ and $M(x)$ is when the tradeoff functions are independent of age $a$, in which case $F(x)=b(x)k_{1}(x)$ and $M(x)=\phi(x)k_{2}(x)$, where $k_{1}(x)$ and $k_{2}(x)$ are the expected number of years survived in the fertile ages for females and males respectively with longevity trait value $x$.

A natural question is, to which longevity value $x$ does the population converge for given parameter values and tradeoff functions of the model? Through the stability analysis in Section 5.1, we show that the answer is simply the longevity value $x$ that maximises the product $F(x)M(x)$. In fact, $F(x)M(x)$ represents the fitness landscape of the scenario outlined in Section 2.1, and so one can examine the shape of $F(x)M(x)$ and infer the optimal longevity value $x$ of the population without performing a simulation of the PDE. The magnitude of the sexual conflict in the system can be inferred from a comparison of $F(x)$, $M(x)$ and $F(x)M(x)$, that is, a comparison of the female, male and combined fitness landscapes. Furthermore, we show that the speed of convergence to the optimum of the fitness landscape is proportional to the slope of the fitness landscape $F(x)M(x)$. We generalise these results for an arbitrary number of traits in Section 5.2.

3 Including grandmothering

In this section, we consider the same problem as in Section 2.1, except that we allow for grandmothering, i.e. provisioning from post-fertile females. We let $L$, an individual’s life expectancy, be the only heritable trait in the population. Similar to Chan et al. [2], we assume that post-fertile females take care of all their matrilineal grandchildren and that grandmothers boost their daughters’ birth rates. The birth rate $b(a,L)$ in Eq. (6) is then of the form $B(G(L))\bar{b}(a,L)$, where $\bar{b}(a,L)$ is the base birth rate without the help of grandmothers, $G(L)$ is the proportion of fertile females covered by a post-fertile mother for a fixed longevity $L$, assuming a stable age distribution, and $B(x)=mx+(1-x)$ represents the boost in fertility. Following Chan et al. [2], we let the base birth rate $\bar{b}(L)$ be independent of age and require it to be equal to 0.3/year at $L=14$ and 0.11/year at $L=34$. Also guided by Chan et al. [2], we set $m=3$, so that grandmothers boost their daughters’ birth rates by a factor of 3, thus increasing the birth rate to 0.33/year for an individual with a life expectancy of 34 years, who is also supported by a post-fertile mother.

We use the parameter settings of Kim et al. [10] for $\tau_{f}(L)$, $\tau_{m}$, $\hat{\tau}_{f}$, $\hat{\tau}_{m}(L)$ and $\mu(L)$. These parameters are listed in Table 1. We include the assumption from Kim et al. [10] that individuals are frail, and thus exit the population, at the age $\tau_{T}=\min(2L,75)$ to prevent individuals living to unrealistic ages. Moreover, males are fertile from sexual maturity $\tau_{m}$ until frailty $\tau_{T}$. We note that the monotonically increasing nature of the age of female sexual maturity $\tau_{f}(L)$ is based on Charnov’s model of mammalian life history, which predicts an increase in maturation age with greater longevities [4]. The male age of sexual maturity is kept constant for simplicity.

Since we could not find a straightforward analytic way to calculate $G(L)$ given arbitrary parameter values, we approximate $G(L)$ by using an agent-based model (ABM). The ABM runs as follows: First a longevity value $L$ is chosen and is fixed for the entire simulation. All individuals have the same longevity $L$ and there is zero probability for mutations. Females are fertile from $\tau_{f}(L)$ to $\hat{\tau}_{f}$. Due to longevity being fixed and the assumption in Section 2.1 that every female is guaranteed to mate in every time step, there is no need for the ABM to track males; tracking males is only essential if $L$ is allowed to evolve. At each time step, every fertile female gives birth to one offspring with probability $1-\exp(-\bar{b}(x))$ if she does not have a surviving post-fertile mother or with probability $1-\exp(-3\bar{b}(x))$ if she does. If there are over 200 newborns, then newborns are randomly removed from the population until only 200 remain. Individuals die with probability $1-\exp(-\mu(x))$ and any individuals past the age of $\min(2L,75)$ are removed from the population. Figure 2a shows 100 samples of $G(x)$ for each integer $L$ past 22, since there are no grandmothers unless $L\geq 23$ due to $\hat{\tau}_{f}=\min(2L,45)$. We fit the function

[TABLE]

through the means of each data set collected in Figure 2a and use the method of least squares to obtain $a=1.164$, $b=43.83$ and $c=1.36$.

Following Kim et al. [11], we let the male fertility-longevity tradeoff $\phi(L)$, the relative probability for a male with longevity $L$ to secure a mate, to decrease exponentially with longevity by requiring it to satisfy $\phi^{\prime}(L)=\delta(L)\phi(L)$. The function $\phi(L)$ is then uniquely defined, up to a scaling factor, by

[TABLE]

where we choose $L_{0}=20$ and $\delta(L)=-0.4\exp(-0.087L)$.

To examine the life expectancy that is selected in the population with and without grandmothering, we compute $F(L)M(L)$. As shown in Figure 2b, the fitness landscape is the same before life expectancies of $L=23$, since $L=23$ is the earliest life expectancy at which the population has any post-fertile females. However, for life expectancies greater than 23, grandmothering alters the fitness landscape and enables a transition to a unique optimum corresponding to a higher longevity.

We confirm that the population converges to the maxima of the fitness landscape in Figure 3a, where we show the expected trait value of the population over time, calculated via Eq. (7). Moreover, to elucidate the role of the sexual conflict in the system, Figure 3b provides a way to visualise the magnitude of the sexual conflict by comparing the female fitness function $F(L)$, male fitness function $M(L)$ and the two-sex fitness landscape $F(L)M(L)$. From Figure 3b, male fitness increases concurrently with longevity $L$; however, this is counterbalanced by a decrease in female fitness, which creates the unique maximum point in the combined two-sex fitness landscape. One can observe the optima preferred by the males and females through their respective fitness landscapes, and also the compromise adopted by the population overall that lies between these two optima, as shown by the two-sex fitness landscape. We note that the male fitness landscape is the same regardless of whether there is grandmothering in the population, since we have made the assumption that the only function of grandmothers in the population is to boost their daughter’s birth rate.

4 Discussion

We have presented a computationally inexpensive method to compute the fitness landscape of the problem presented in Section 2.1 and have applied it to determine if grandmothering in a population can enable an evolutionary trajectory to higher longevities. Previous studies of this problem have included both short- and long-time scale components, which result in computationally expensive simulations and hinders inference from the model. Our method essentially compacts the short-time scale components of the problem, such as mortality and ageing, by assuming a stable age-distribution, while focusing on the evolution of traits in the population, which we assume to be of a much longer-time scale. The resultant equation describing the evolutionary dynamics, given by Eq. (7), can be viewed as a dynamical systems generalisation of Price’s Equation that includes two-sexes and multiple traits.

A difficulty encountered by previous studies in the literature is assessing the role of the sexual conflict in the problem. This is because the sexual conflict is a manifestation of the many different components; specifically, the mating sex ratio (which is itself determined by the male and female fertile ages, and the mortality rate) and the male and female tradeoffs with longevity. Although it is tempting to analyse the sexual conflict through the mating sex ratio, it is not the correct approach; we stress that since a given mating sex ratio value can be achieved through a variety of different fertile ages and mortality rates, there is not a unique correspondence between mating sex ratio values and fitness landscapes. This is evident from the form of the female fitness landscape $F(L)$ and male fitness landscape $M(L)$ (see Eq. (6)). More specifically, the function $k(L)$, denoting the mating sex ratio at longevity $L$ and given by

[TABLE]

is not guaranteed to be invertible. This implies that the same mating sex ratio value may be obtainable at different longevity $L$ values. Consequently, although the sexual conflict is clearly present in the system, it is not an aspect that can be accurately represented through a single indicator and this hinders its analysis. In this study, we provide the next best approach: a method to compute the female and male fitness landscapes separately, and also the combined fitness landscape. Through this, one can observe how each sex contributes to the resulting evolutionary dynamics of the overall problem and how the male and female fitness landscapes vary with the mating sex ratio. In particular, when the female and male tradeoffs are quantitatively known, our method allows one to assess how much contribution a change in mating sex ratio contributes to the fitness landscape of a population.

A limiting assumption of our model is that fertile males compete over every fertile female for a chance at paternity regardless of the mating sex ratio. However, a shifting mating sex ratio will influence the availability of partners due to the asymmetry in the male and female fertile ages. Hence, mating strategies are likely to change along with the mating sex ratio (see Schacht & Bell [17], Coxworth et al. [5] and Loo et al. [12] for examples). Thus, an avenue for future work is including such changes in mating strategies in response to a shifting mating sex ratio via the evolution of longevity in a population. This could in turn have important consequences on the evolutionary trajectory the population follows towards higher longevities.

Acknowledgements:

MHC and PSK were supported by the Australian Research Council, Discovery Project (DP160101597).

5 Appendix

5.1 Stability

We are interested in the equilibria of the system and how they are affected by the shape of $F(x)$ and $M(x)$. Specifically, we wish to find the expected value of $u_{n}(x)$ as $n\to\infty$. To examine this, we consider the limiting case $\varepsilon\to 0$, which represents a scenario where the probability for mutations in $x$ is zero. This yields the system

[TABLE]

Since Eq. (7) is a regular perturbation problem with respect to $\varepsilon$, we expect that Eq. (11) is an accurate approximation of Eq. (7). One can check that this has an infinite number of fixed points in the form $u^{*}(x)=\delta(x-\mu)$, for any constant $\mu$. However, assuming that the initial condition $u_{0}(x)$ is strictly positive, bounded and integrable, the sequence of functions $u_{n}(x)$ only converges to $u^{*}(x)$ for particular values of $\mu$. To determine such values of $\mu$, we compute the first moment of both sides of Eq. (11), which we denote by $\Lambda$. It is well known that if $X$ and $Y$ have the density functions $f$ and $g$ respectively, then the density function for $\frac{1}{2}(X+Y)$ is $2(f(x)\ast g(x))[2x]$. We have that $\Lambda(2(f(x)\ast g(x))[2x])=\text{\bf{E}}(\frac{1}{2}(X+Y))=\frac{1}{2}(\Lambda f+\Lambda g)$. Thus,

[TABLE]

Due to the fixed points being in the form of $u^{*}(x)=\delta(x-\mu)$, for any constant $\mu$, to assess stability we fix $u_{n}(x)$ in the following box function form

[TABLE]

where $\alpha>0$ is a small parameter and $\mu_{n}$ is the expected value of $u_{n}(x;\mu,\alpha)$ by construction. This acts as a perturbation of $\delta(x-\mu)$ for $\alpha$ small. Substituting this into Eq. (12), we obtain

[TABLE]

We require that $\mu_{n}$ such that it is a fixed point of Eq. (12). Since $\alpha$ is arbitrarily small, we take the first order approximation of $F(x)$ and $M(x)$ about $\mu_{n}$, that is, we let $F(x)=F(\mu_{n})+F^{\prime}(\mu_{n})(x-\mu_{n})$ and $M(x)=M(\mu_{n})+M^{\prime}(\mu_{n})(x-\mu_{n})$. After some algebraic manipulation, Eq. (14) simplifies to the recurrence relation

[TABLE]

Thus, fixed points $\mu^{*}$ of Eq. (15) satisfy the condition

[TABLE]

To determine the stability of $\mu^{*}$ satisfying the above condition, we examine $\Phi(\mu^{*})$. We require $|\Phi^{\prime}(\mu^{*})|<1$ for local stability (see Theorem 1.5 in [18]). We have that

[TABLE]

Thus, fixed points satisfying Eq. (16) are locally stable when

[TABLE]

We note that $\alpha$ can be chosen such that $-12/\alpha^{2}$ is arbitrarily negative. The conditions given in Eq. (16) and Eq. (18) indicate that $F(x)M(x)$ represents the total fitness landscape. As $n$ increases, the first moment of $u_{n}(x)$ moves in the direction of increasing gradient in $F(x)M(x)$. We note that it is possible for the solution to be trapped in a local optima, i.e. the population converges to the nearest peak on the fitness landscape.

5.2 Generalisation to $n$ traits

In this section we generalise the results in Section 5.1 to $n$ traits, where we emulate the derivation, except in $n$-dimensions. Let $\mathbf{x}=(x_{1},x_{2},\dots,x_{n})^{\prime}$ be the vector of $n$ traits. The PDE in Eq. (1) remains unchanged, except that the functions $u$ and $\mu$ map from $\mathbb{R}^{n}$ to $\mathbb{R}$. The boundary condition given in Eq. (2) becomes

[TABLE]

where $b(\mathbf{x})$, $\phi({\mathbf{x}})$, $\tau_{f}({\mathbf{x}})$, $\tau_{m}({\mathbf{x}})$, $\hat{\tau}_{f}({\mathbf{x}})$, $\hat{\tau}_{m}({\mathbf{x}})$, $S_{f}(\mathbf{x})$ and $S_{m}({\mathbf{x}})$ are now functions which map from $\mathbb{R}^{n}$ to $\mathbb{R}$. The vectors $\mathbf{y}:=(y_{1},y_{2},\dots,y_{n})^{\prime}$ and \hskip 0.27779pt\hbox to0.0pt{\leavevmode\hbox{\set@color\mathbf{z}}\hss}{\hskip-0.27779pt\leavevmode\hbox{\set@color\overline{\hbox{}}}}:=(\hbox to0.0pt{\hskip 0.04514pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorz}}_{1},\hbox to0.0pt{\hskip 0.04514pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorz}}_{2},\dots,\hbox to0.0pt{\hskip 0.04514pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorz}}_{n})^{\prime} are integration dummy variables. Furthermore, $N\left(\mathbf{x},\mathbf{y},\varepsilon^{2}\right)$ is the $n$ dimensional normal density function given by

[TABLE]

Lastly, \hbox to0.0pt{\hskip 0.85414pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorS}}_{f} denotes $\int_{\mathbb{R}^{n}}S_{f}(\mathbf{x})\,d\mathbf{x}$ and similarly for \hbox to0.0pt{\hskip 0.85414pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorS}}_{m}.

Letting \hskip 0.27779pt\hbox to0.0pt{\leavevmode\hbox{\set@color\mathbf{z}}\hss}{\hskip-0.27779pt\leavevmode\hbox{\set@color\overline{\hbox{}}}}=2\mathbf{z}-\mathbf{y}, we obtain

[TABLE]

where $\otimes$ denotes the $n$-dimensional convolution operator, that is,

[TABLE]

Letting $\varepsilon\to 0$ and assuming that the population is at stable age distribution for all time, that is,

[TABLE]

we obtain

[TABLE]

where $\langle F,u_{t}\rangle$ denotes $\int_{\mathbb{R}^{n}}F(\bm{x})u_{t}(\mathbf{x})\,d\mathbf{x}$ and similarly for $\langle M,u_{t}\rangle$.

Computing the first moment of both sides leads to

[TABLE]

where we have defined

[TABLE]

Since the fixed points of Eq. (24) are in the form $u^{*}=\delta(x_{1}-\mu_{1})\delta(x_{2}-\mu_{2})\dots\delta(x_{3}-\mu_{3})$, we examine stability by perturbing this fixed point by setting $u_{t}(\mathbf{x})$ to the following box function form

[TABLE]

where $A$ is the $n$-dimensional box centered at coordinates given by $\bm{\mu}_{t}$ with length $2\alpha$ along each dimension. Thus, the volume of $A$ is equal to 1. Substituting this into Eq. (25), we obtain

[TABLE]

where

[TABLE]

and similarly for the other integrals in Eq. (28). The superscript on $\bm{\mu}_{t}$ denotes the index, for example, $\mu^{n}_{t}$ denotes the $n$th element of $\bm{\mu}_{t}$. Since $\alpha$ is arbitrarily small, we approximate $F(\mathbf{x})$ and $M(\mathbf{x})$ to first order, that is, we let $F(\mathbf{x})=F(\bm{\mu})+\frac{\partial}{\partial x_{1}}F(\mu_{1})(x_{1}-\mu_{1})+\frac{\partial}{\partial x_{2}}F(\mu_{2})(x_{2}-\mu_{2})+\dots+\frac{\partial}{\partial x_{n}}F(\mu_{n})(x_{n}-\mu_{n})$ and $M(\mathbf{x})=M(\bm{\mu})+\frac{\partial}{\partial x_{1}}M(\mu_{1})(x_{1}-\mu_{1})+\frac{\partial}{\partial x_{2}}M(\mu_{2})(x_{2}-\mu_{2})+\dots+\frac{\partial}{\partial x_{n}}M(\mu_{n})(x_{n}-\mu_{n})$. Substituting this into Eq. (28) and simplifying, we obtain

[TABLE]

Thus, fixed points of Eq. (30) are also critical points of $F(\mathbf{x})M(\mathbf{x})$. In order to assess local stability, we examine the Jacobian matrix $\bm{J}$ of Eq. (30) evaluated at the critical points $\bm{x^{*}}$ of $F(\mathbf{x})M(\mathbf{x})$. It is given by

[TABLE]

where $\bm{H}(F(\mathbf{x})M(\mathbf{x}))$ denotes the Hessian matrix of $(F(\mathbf{x})M(\mathbf{x})$ and $\bm{I}$ denotes the identity matrix. For critical points to be locally stable, we require that all eigenvalues $\lambda$ of the matrix given in Eq. (31) to satisfy $|\lambda|<1$. Since we have that

[TABLE]

if all eigenvalues $\lambda$ of Eq. (31) satisfy $|\lambda|<1$, then all eigenvalues of \bm{H}(F(\mathbf{x})M(\mathbf{x}))\Big{|}_{\mathbf{x}=\mathbf{\bm{x^{*}}}} are negative. Thus, locally stable fixed points of Eq. (30) correspond to local maxima of $F(\mathbf{x})M(\mathbf{x})$. From Eq. (30), we can see that the population converges to locally stable fixed points along the gradient of $F(\mathbf{x})M(\mathbf{x})$, i.e. the direction of steepest ascent on the fitness landscape.