# Fitness versus Longevity in Age-Structured Population Dynamics

## I Introduction

The goal of this paper is to understand the role of mutations on the evolution of fitness and age characteristics of individuals in a simple age-structured population dynamics model. While there are many classical models to describe single-species population dynamics, consideration of age-dependent characteristics is a more recent development. Typically, age characteristics of a population are determined by studying rate equations which include age-dependent birth and death rates. Here we will study an extension of age-structured population dynamics in which the characteristics of an offspring are mutated with respect to its parent. In particular, an offspring may be more “fit” or less fit than its parent, and this may be reflected in attributes such as its birth rate and/or its life expectancy.

In our model, we characterize the fitness of an individual by a single heritable trait – the life expectancy $`n`$ – which is defined as the average life span of an individual in the absence of competition. This provides a useful fitness measure, as a longer-lived individual has a higher chance of producing more offspring throughout its life span. We allow for either deleterious or advantageous mutations, where the offspring fitness is less than or greater than that of the parent, respectively (Fig. 1). This leads to three different behaviors which depend on the ratio between these two mutation rates. When advantageous mutation is favored, the fitness distribution of the population approaches a Gaussian, with the average fitness growing linearly with time $`t`$ and the width of the distribution increasing as $`t^{1/2}`$. Conversely, when deleterious mutation is more likely, a steady-state fitness distribution is approached, with the rate of approach varying as $`t^{2/3}`$. When there is no mutational bias, the fitness distribution again approaches a Gaussian, but with average fitness growing as $`t^{2/3}`$ and the width of the distribution again growing as $`t^{1/2}`$.

In all three cases, the average population age reaches a steady state which, surprisingly, is a decreasing function of the average fitness. Thus within our model, a more fit population does not lead to an increased individual lifetime. Qualitatively, as individuals become more fit, competition plays a more prominent role and is the primary mechanism that leads to premature mortality.

In the following two sections, we formally define the model and outline qualitative features of the population dynamics. In Secs. IV-VI, we analyze the three cases of deleterious, advantageous, and neutral mutational biases in detail. We conclude in Sec. VII. Various calculational details are provided in the Appendices.

## II The Model

Our model is a simple extension of logistic dynamics in which a population with overall density $`N(t)`$ evolves both by birth at rate $`b`$ and death at rate $`\gamma N`$. Such a system is described by the rate equation

$$\dot{N}=bN\gamma N^2,$$
(1)   
with steady-state solution $`N_{\mathrm{}}=b/\gamma `$. Our age-structured mutation model incorporates the following additional features:

1. Each individual is endowed with given life expectancy $`n`$. This means that an individual has a rate of dying which equals $`1/n`$.
2. Death by aging occurs at a rate inversely proportional to the life expectancy.
3. Individuals give birth at a constant rate during their lifetimes.
4. In each birth event, the life expectancy of the newborn may be equal to that of its parent, or the life expectancy may be increased by or decreased by 1. The relative rates of these events are $`b`$, $`b_+`$, and $`b_{}`$, respectively.

Each of these features represent idealizations. Most prominently, it would be desirable to incorporate a realistic mortality rate which is an increasing function of age. However, we argue in Sec. VII that our basic conclusions continue to be valid for systems with realistic mortality rates.

To describe this dynamics mathematically, we study $`C_n(a,t)`$, the density of individuals with life expectancy $`n1`$ and age $`a`$ at time $`t`$. We also introduce $`P_n(t)=_0^{\mathrm{}}C_n(a,t)𝑑a`$, which is the density of individuals with life expectancy $`n`$ and any age at time $`t`$. Finally, the total population density is the integral of the population density over all ages and life expectancies,

$$N(t)=\underset{n=1}{\overset{\mathrm{}}{}}_0^{\mathrm{}}C_n(a,t)𝑑a=\underset{n=1}{\overset{\mathrm{}}{}}P_n(t).$$
(2)   
According to our model, the rate equation for $`C_n(a,t)`$ is

$$\left(\frac{}{t}+\frac{}{a}\right)C_n(a,t)=\left(\gamma N(t)+\frac{1}{n}\right)C_n(a,t).$$
(3)   
The derivative with respect to $`a`$ on the left-hand side accounts for the continuous aging of the population. On the right-hand side, $`\gamma NC_n`$ is the death rate due to competition, which is assumed to be independent of an individual’s age and fitness. As discussed above, the mortality rate is taken as age independent, and the form $`C_n/n`$ guarantees that the life expectancy in the absence of competition equals $`n`$. Because birth creates individuals of age $`a=0`$, the population of newborns provides the following boundary condition for $`C_n(0,t)`$,

$$C_n(0,t)=bP_n(t)+b_+P_{n1}(t)+b_{}P_{n+1}(t).$$
(4)   
Finally, the condition $`P_0=0`$ follows from the requirement that offspring with zero life expectancy cannot be born.

## III Basic Population Characteristics

Let us first study the fitness characteristics of the population and disregard the age structure. The rate equation for $`P_n(t)`$ is found by integrating Eq. (3) over all ages and then using the boundary condition Eq. (4) to give

$$\frac{dP_n}{dt}=\left(b\gamma N\frac{1}{n}\right)P_n+b_+P_{n1}+b_{}P_{n+1}.$$
(5)   
This describes a random-walk-like process with state-dependent hopping rates in the one-dimensional fitness space $`n`$. Notice the hidden non-linearity embodied by the term $`\gamma NP_n`$, since the total population density $`N(t)=_nP_n(t)`$. From Eq. (5), we find that $`N(t)`$ obeys a generalized logistic equation

$$\frac{dN}{dt}=(b+b_++b_{}\gamma N)N\underset{n=1}{\overset{\mathrm{}}{}}\frac{P_n}{n}b_{}P_1.$$
(6)   
The three different dynamical regimes outlined in the introduction are characterized by the relative magnitudes of the mutation rates $`b_+`$ and $`b_{}`$. The main features of these regimes are:

* Subcritical case. Here $`b_+<b_{}`$, or deleterious mutations prevail. The fitness of the population eventually reaches a steady state.
* Critical case. Here $`b_+=b_{}`$, or no mutational bias. The average fitness of the population grows as $`t^{2/3}`$ and the width of the fitness distribution grows as $`t^{1/2}`$.
* Supercritical case. Here $`b_+>b_{}`$, or advantageous mutations are favorable. The average fitness grows linearly in time and the width of the fitness distribution still grows as $`t^{1/2}`$.

In all three cases, the total population density $`N`$ and the average age $`A`$, defined by

$$A=\frac{1}{N}\underset{n=1}{\overset{\mathrm{}}{}}𝑑aaC_n(a)$$
(7)   
saturate to finite values. The steady state values of $`N`$ and $`A`$ may be determined by balance between the total birth rate $`Bb+b_++b_{}`$ and the death rate $`\gamma N`$ due to overcrowding. For example, in the critical and supercritical cases, there are essentially no individuals with small fitness, so that the last two terms in Eq. (6) may be neglected. Then the steady-state solution to this equation is simply

$$N=\frac{B}{\gamma }.$$
(8)   
This statement also expresses the obvious fact that in the steady state the total birth rate $`B`$ must balance the total death rate $`\gamma N`$. (For fit populations, the death rate due to aging is negligible.) Similarly, the average age may be inferred from the condition it must equal the average time between death events. Thus

$$A=\frac{1}{\gamma N}=\frac{1}{B}.$$
(9)   
The behavior of the average age in the subcritical case is more subtle and we treat this case in detail in the section following.

## IV the subcritical case

When deleterious mutations are favored ($`b_{}>b_+`$), the random-walk-like master equation for $`P_n`$ contains both the mutational bias towards the absorbing boundary at the origin, as well as an effective positive bias due to the $`1/n`$ term on the right-hand side of Eq. (5). The balance between these two opposite biases leads to a stationary state whose solution is found by setting $`\dot{P}_n=0`$ in Eq. (5). To obtain this steady state solution, it is convenient to introduce the generating function

$$F(x)=\underset{n=1}{\overset{\mathrm{}}{}}P_nx^{n1}.$$
(10)   
Multiplying Eq. (5) by $`x^{n1}`$ and summing over $`n`$ gives

$$0=(b\gamma N)F\underset{n=1}{\overset{\mathrm{}}{}}\frac{P_n}{n}x^{n1}+b_+xF+b_{}\left(\frac{F}{x}P_1\right).$$
(11)   
The term involving $`P_n/n`$ is simplified by using

$$\underset{n=1}{\overset{\mathrm{}}{}}\frac{P_n}{n}x^{n1}=\frac{1}{x}_0^xF(y)𝑑y.$$
(12)   
Multiplying Eq. (11) by $`x`$ and differentiating with respect to $`x`$ gives

$$\frac{F^{}(x)}{F(x)}=\frac{\gamma Nb+12b_+x}{b_+x^2(\gamma Nb)x+b_{}},$$
(13)   
where the prime denotes differentiation.

As in the case of the master equation for $`P_n`$, this differential equation for $`F`$ has a hidden indeterminacy, as the total population density $`N=F(x=1)`$ appears on the right-hand side. Thus an integration of Eq. (13), subject to the boundary condition $`F(1)=N`$, actually gives a family of solutions which are parameterized by the value of $`N`$. While the family of solutions can be obtained straightforwardly by a direct integration of Eq. (13), only one member of this family is the correct one. To determine this true solution, we must invoke additional arguments about the physically realizable value of $`N`$ for a given initial condition.

An upper bound for $`N`$ may be found from the steady-state version of Eq. (6),

$$(B\gamma N)N=\underset{n=1}{\overset{\mathrm{}}{}}\frac{P_n}{n}+b_{}P_1.$$
(14)   
Since the right-hand side must be non-negative, this provides the bound $`\gamma N<B`$. On the other hand, we may obtain a lower bound for $`N`$ by considering the master equation for $`P_n`$ in the steady state. For $`n\mathrm{}`$, we may neglect the $`P_n/n`$ term in Eq. (5) and then solve the resulting equation to give $`P_n=A_+\lambda _+^n+A_{}\lambda _{}^n`$, where

$$\lambda _\pm =\left[\gamma Nb\pm \sqrt{(\gamma Nb)^24b_+b_{}}\right]/2b_{}.$$
(15)   
For $`P_n`$ to remain positive, $`\lambda _\pm `$ should be real. Hence we require $`\gamma N>b+2\sqrt{b_+b_{}}`$. We therefore conclude that $`N`$ must lie in the range

$$b+2\sqrt{b_+b_{}}\gamma N<B.$$
(16)   
While the foregoing suggests that $`N`$ lies within a finite range, we find numerically that the minimal solution, which satisfies the lower bound of Eq. (16), is the one that is generally realized (Fig. 3). This selection phenomenon is reminiscent of the corresponding behavior in the Fisher-Kolmogorov equation and related reaction-diffusion systems, where only the extremal solution is selected from a continuous range of a priori solutions.

To understand the nature of this extremal solution in the present context, notice that with the bounds on $`N`$ given in Eq. (16), $`\lambda _+`$ lies within the range $`[\mu ,1)`$, where

$$\mu \sqrt{\frac{b_+}{b_{}}}$$
(17)   
is the fundamental parameter which characterizes the mutational bias. Consequently the steady-state fitness distribution decays exponentially with $`n`$, namely $`P_n\lambda _+^n`$. When the total population density attains the minimal value $`N_{\mathrm{min}}=(b+2\sqrt{b_+b_{}})/\gamma `$, $`\lambda _+`$ also achieves its minimum possible value $`\lambda _+^{\mathrm{min}}=\mu `$, so that the fitness distribution has the most rapid decay in $`n`$ for the minimal solution. Based on the analogy with the Fisher-Kolmogorov equation, we infer that there are two distinct steady-state behaviors for $`P_n`$ as a function of the initial condition $`P_n(0)`$. For any $`P_n(0)`$ with either a finite support in $`n`$ or decaying at least as fast as $`\mu ^n`$, the extremal solution $`P_n\mu ^n`$ is approached as $`t\mathrm{}`$. Conversely, for initial conditions in which $`P_n(0)`$ decays more slowly than $`\mu ^n`$, for example as $`\alpha ^n`$, with $`\alpha `$ in the range $`(\mu ,1)`$, the asymptotic solution also decays as $`\alpha ^n`$. Correspondingly, Eq. (5) in the steady state predicts a larger than minimal population density $`N=(b+b_{}\alpha +b_+\alpha ^1)/\gamma `$.

We also find that the extremal and the non-extremal solutions exhibit different relaxations to the steady state. For those initial conditions which evolve to the extremal solution, the deviation of $`N`$ and indeed each of the $`P_n`$ from their steady state values decay as $`t^{2/3}`$, while for all other initial conditions, the relaxation to the steady state appears to follow a $`t^1`$ power law decay. The power-law approach to the steady state is surprising, since the overall density obeys a logistic-like dynamics, $`\dot{N}=bN\gamma N^2`$, for which the approach to the steady state is exponential. These results are illustrated in Fig. 4 which shows the asymptotic time dependence of $`N(t)`$ based on a numerical integration of Eq. (5) with the fourth-order Runge-Kutta algorithm. The demonstration of the $`t^{2/3}`$ relaxation to the extremal solution relies on a correspondence to the transient behavior in the critical case. This is presented in Appendix B.

A disconcerting feature of the numerical calculation for $`N(t)`$ is the small disagreement between the numerically observed values of the steady-state population density and the expected theoretical prediction (Fig. 5). This discrepancy arises from the finite computer precision which causes very small values of $`P_n`$ to be set to zero. To confirm this, we changed the computer precision from $`10^{100}`$ to the full machine precision of $`10^{308}`$ (Fig. 5). As the precision is increased, $`N`$ saturates to progressively higher values and approaches the theoretical prediction. A similar precisions-dependent phenomenon has been observed in the context of traveling Fisher-Kolmogorov wave propagation.

For the relevant situation where the density $`N`$ takes the minimal value, we may rewrite Eq. (13) as

$$\frac{F^{}}{F}=\frac{2\mu }{1\mu x}+\frac{1}{b_{}(1\mu x)^2}$$
(18)   
Integrating from $`x=1`$ to $`x`$ and using $`F(1)=N`$ gives

$$F(x)=N\left(\frac{1\mu }{1\mu x}\right)^2\mathrm{exp}\left\{\frac{1}{b_{}\mu }\left(\frac{1}{1\mu x}\frac{1}{1\mu }\right)\right\}.$$
(19)   
One can now formally determine $`P_n`$ by expanding $`F(x)`$ in a Taylor series. For example,

$`P_1`$ $`=`$ $`N(1\mu )^2\mathrm{exp}\left\{{\displaystyle \frac{1}{b_{}(1\mu )}}\right\},`$   
$`P_2`$ $`=`$ $`N(1\mu )^2\left(2\mu +{\displaystyle \frac{1}{b_{}}}\right)\mathrm{exp}\left\{{\displaystyle \frac{1}{b_{}(1\mu )}}\right\}.`$   
For many applications, however, there is no need to deal with these unwieldy expressions. As we now discuss, the overall fitness or age characteristics of the population can be obtained directly from the generating function without using the explicit formulae for the $`P_n`$.

### A Fitness characteristics

Consider the average fitness of the population

$$n=\frac{1}{N}\underset{n=1}{\overset{\mathrm{}}{}}nP_n,$$
(20)   
which can be expressed in terms of the generating function as

$$n=\frac{1}{N}\frac{dF}{dx}|_{x=1}+1.$$
(21)   
From Eq. (19) we thereby obtain the average fitness

$$n=\frac{2\mu }{1\mu }+\frac{1}{b_{}(1\mu )^2}+1.$$
(22)   
As one might anticipate, the average fitness diverges as $`\mu 1`$ from below, corresponding to the population becoming mutationally neutral. To determine the dispersion of the fitness distribution we make use of the relation

$$n(n1)=\frac{1}{N}\underset{n=1}{\overset{\mathrm{}}{}}n(n1)P_n=\frac{1}{N}\frac{d^2(xF)}{dx^2}|_{x=1}.$$
(23)   
Substituting Eqs. (19) and also Eq. (22) then gives

$`n^2=1+{\displaystyle \frac{6\mu }{(1\mu )^2}}+{\displaystyle \frac{3(1+\mu )}{b_{}(1\mu )^3}}+{\displaystyle \frac{1}{b_{}^2(1\mu )^4}}.`$   
Thus the dispersion $`\sigma ^2=n^2n^2`$ in the fitness distribution is

$$\sigma ^2=\frac{2\mu }{(1\mu )^2}+\frac{\mu +1}{b_{}(1\mu )^3}.$$
(24)   
As the mutational bias vanishes, $`\mu 1`$, the average fitness and the dispersion diverge as $`nb_{}^1(1\mu )^2`$ and $`\sigma \sqrt{2/b_{}}(1\mu )^{3/2}`$. Thus these two moments are related by $`\sigma n^{3/4}`$. As we shall see in Sec. VI, this basic relation continues to hold in the critical case.

### B Age characteristics

In the steady state, we solve Eq. (3) to give the concentration of individuals with age $`a`$ and fitness $`n`$

$$C_n(a)=P_n\left(\gamma N+\frac{1}{n}\right)\mathrm{exp}\left[\left(\gamma N+\frac{1}{n}\right)a\right].$$
(25)   
The average age of the population is

$`A`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle _0^{\mathrm{}}}𝑑aaC_n(a)`$ (26)   
$`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{P_n}{\gamma N+n^1}},`$ (27)   
where the second line is obtained by using Eq. (25). This expression can be rewritten directly in terms of the generating function by first noticing that

$`{\displaystyle _0^1}x^\nu F(x)𝑑x`$ $`=`$ $`{\displaystyle _0^1}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}P_nx^{n+\nu 1}dx`$ (28)   
$`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{P_n}{n+\nu }}.`$ (29)   
Thus we re-express Eq. (26) in a form which allows us to exploit Eq. (28). After several elementary steps, we obtain

$`A`$ $`=`$ $`{\displaystyle \frac{1}{\gamma N}}{\displaystyle \frac{1}{N}}{\displaystyle \frac{1}{(\gamma N)^2}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{P_n}{n+(\gamma N)^1}}`$ (30)   
$`=`$ $`{\displaystyle \frac{1}{\gamma N}}{\displaystyle \frac{1}{N}}{\displaystyle \frac{1}{(\gamma N)^2}}{\displaystyle _0^1}𝑑xx^{\frac{1}{\gamma N}}F(x).`$ (31)   
This expression should be compared with the result for the critical and supercritical cases, namely $`A=(\gamma N)^1=B^1`$ (see Eq. (9)). In the subcritical case, $`\gamma N<B`$ and the above two expressions $`A_{\mathrm{min}}=B^1`$ and $`A_{\mathrm{max}}=(\gamma N)^1`$ provide lower and upper bounds for the average age. This is proved in Appendix A. Fig. 6 shows the surprising feature of Eq. (30) that the average age decreases as the population gets fitter! We also see that the average age of the least fit population $`(\mu 0)`$ is twice that of the increasingly fit populations in the critical and supercritical cases. We now demonstrate this fact. To provide a fair comparison we take the total birth rate rate to be equal to unity in both cases and also choose $`b=0`$ for simplicity. For fit populations (critical and supercritical cases), the average age is simply $`A=B^1=1`$. For the least fit population $`\mu 0`$, and correspondingly $`N0`$. In this limit, we may write Eq. (26) as,

$$A=\frac{1}{N}\underset{n=1}{\overset{\mathrm{}}{}}\frac{P_n}{\gamma N+n^1}\frac{1}{N}\underset{n=1}{\overset{\mathrm{}}{}}nP_n=n.$$
(32)   
On the other hand, from Eq. (22) we have

$$n1+\frac{1}{b_{}}2.$$
(33)   
The relation $`A=n`$ is natural for the least fit population, as the total density is small and competition among individuals plays an insignificant role. Thus the average age may be found by merely averaging the intrinsic life expectancy of the population. Intriguingly, in this limit the average individual in the least fit population lives twice as long as individuals in relatively fit populations.

It is also worth noting that in the limit of a minimally fit population ($`\mu 0`$) we can expand the generating function in Eq. (19) systematically. We thereby find that the density $`P_n`$ exhibits a super-exponential decay, $`P_n=Ne^1/(n1)!`$.

## V the supercritical case

When advantageous mutations are favored, the master equation for $`P_n`$, Eq. (5), can be viewed as a random walk with a bias towards increasing $`n`$. Because there is no mechanism to counterbalance this bias, the average fitness grows without bound and no steady state exists. As in the case of a uniformly biased random walk on a semi-infinite domain, the distribution of fitness becomes relatively localized in fitness space, with the peak drifting towards increasing $`n`$ with a velocity $`V=b_+b_{}`$. Since the fitness profile is non-zero only for large $`n`$ in the long time limit, it is appropriate to adopt continuum approach. We therefore treat $`n`$ as continuous, and derive the continuum limits of Eqs. (5) and (6). For the time evolution of the fitness distribution $`P(n,t)`$, we obtain the equation of motion

$$\left(\frac{}{t}+V\frac{}{n}\right)P=\left(B\gamma N\frac{1}{n}\right)P+D\frac{^2P}{n^2}.$$
(34)   
This is just a convection-diffusion equation, supplemented by a birth/death term. Here the difference between advantageous and deleterious mutations provides the drift velocity $`V=b_+b_{}`$, and the average mutation rate $`D=(b_++b_{})/2`$ plays the role of diffusion constant. For the total population density we obtain

$$\frac{dN}{dt}=(B\gamma N)N_0^{\mathrm{}}𝑑n\frac{P(n,t)}{n}.$$
(35)   
To determine the asymptotic behavior of these equations, we use the fact that the fitness distribution becomes strongly localized about a value of $`n`$ which increases as $`Vt`$. Thus we replace the integral in Eq. (35) by its value at the peak of the distribution, $`N/Vt`$. With this crude approximation, Eq. (35) becomes

$$\frac{dN}{dt}=\left(B\gamma N\frac{1}{Vt}\right)N.$$
(36)   
Thus we conclude that the density approaches its steady state value as

$$\gamma NB\frac{1}{Vt}.$$
(37)   
This provides both a proof Eq. (8), as well as the rate of convergence to the steady state.

We now substitute this asymptotic behavior for the total population density into Eq. (34) to obtain

$$\left(\frac{}{t}+V\frac{}{n}\right)P=\left(\frac{1}{Vt}\frac{1}{n}\right)P+D\frac{^2P}{n^2}.$$
(38)   
Notice that the birth/death term on the right hand side is negative (positive) for subpopulations which are less (more) fit than average fitness $`Vt`$. This birth/death term must also be zero, on average, since the total population density saturates to a constant value. Moreover, this term must be small near the peak of the fitness distributions where $`nVt`$. Thus as a simple approximation, we merely neglect this birth/death term and check the validity of this assumption a posteriori. This transforms Eq. (38) into the classical convection-diffusion equation whose solution is

$$P(n,t)=\frac{N}{\sqrt{4\pi Dt}}\mathrm{exp}\left[\frac{(nVt)^2}{4Dt}\right].$$
(39)   
This basic results implies that the fitness distribution indeed is a localized peak, with average fitness growing linearly in time, $`n=Vt`$, and width growing diffusively, $`\sigma =\sqrt{2Dt}`$. We now check the validity of dropping the birth/death term in Eq. (38). Near the peak, $`|nVt|\sqrt{Dt}`$, so that the birth/death term is of order $`t^{3/2}\times P`$. On the other hand, the other terms in Eq. (38) are of order $`t^1\times P`$, thus justifying the neglect of birth/death term.

We now turn to the age characteristics. Asymptotically, the density of individuals with given age and fitness changes slowly with time because the overall density reaches a steady state. Consequently, the time variable $`t`$ is slow while the age variable $`a`$ is fast. Physically this contrast reflects the fact that during the lifetime of an individual the change in the overall age characteristics of the population is small. Thus in the first approximation, we retain only the age derivative in Eq. (3). We also ignore the term $`C_n/n`$, which is small near the peak of the asymptotic fitness distribution. Solving the resulting master equation and using the boundary condition of Eq. (4) we obtain

$`C_n(a,t)`$ $``$ $`P_n(t)\gamma Ne^{\gamma Na}`$ (40)   
$`=`$ $`{\displaystyle \frac{\gamma N^2}{\sqrt{4\pi Dt}}}\mathrm{exp}\left[\gamma Na{\displaystyle \frac{(nVt)^2}{4Dt}}\right].`$ (41)   
Integrating over the fitness variable, we find that the age distribution $`C(a,t)=𝑑nC_n(a,t)`$ has the stationary Poisson form

$$C(a)=\gamma N^2e^{\gamma Na}.$$
(42)   
From this, the average age is $`A=(\gamma N)^1=B^1`$ in agreement with Eq. (9). As discussed in Sec. III B, the surprising conclusion is that the average age in the supercritical case is always smaller than that in the subcritical case.

## VI the critical case

We now consider the critical case where the rates of advantageous and deleterious mutations are equal. Without the $`1/n`$ term and with $`\gamma Nb=b_++b_{}`$, Eq. (5) becomes the master equation for an unbiased random walk on the semi-infinite range $`n0`$. Due to the $`1/n`$ term, the system has a bias towards increasing $`n`$ which vanishes as $`n\mathrm{}`$ (see Fig. 2). Thus we anticipate that the average fitness will grow faster than $`t^{1/2}`$ and slower than $`t`$. Hence we can again employ the continuum approach to account for the evolution of the $`P_n`$. In this limit, the corresponding master equation for $`P(n,t)`$ becomes

$$\frac{P}{t}=\left(B\gamma N\frac{1}{n}\right)P+D\frac{^2P}{n^2}.$$
(43)   
Numerically, we find $`nt^{2/3}`$, while the dispersion still grows as $`t^{1/2}`$, that is, as $`\sigma \sqrt{t}`$. Thus these two quantities are related by $`\sigma n^{3/4}`$, as derived analytically for the subcritical case.

To provide a more quantitative derivation of the above scaling laws for $`n`$ and $`\sigma `$, as well as to determine the fitness distribution itself, we examine the equation for $`P(n,t)`$. First note that the total population density still obeys Eq. (35), as in the supercritical case. Under the assumption that the fitness distribution is relatively narrow compared to its mean position, a result which we have verified numerically, we again estimate the integral on the right-hand side of Eq. (35) to be of the order of $`N/n`$. This leads to

$$\gamma NB\frac{1}{n}.$$
(44)   
Substituting this into Eq. (43) yields

$$\frac{P}{t}=\left(\frac{1}{n}\frac{1}{n}\right)P+D\frac{^2P}{n^2}.$$
(45)   
Given that the peak of the distribution is located near $`nn`$, it proves useful to change variables from $`(n,t)`$ to the comoving co-ordinates $`(y=nn,t)`$ to determine how the peak of the distribution spreads. We therefore write the derivatives in the comoving coordinates

$`{\displaystyle \frac{}{t}}={\displaystyle \frac{}{t}}|_{\mathrm{comov}.}{\displaystyle \frac{dn}{dt}}{\displaystyle \frac{}{y}},{\displaystyle \frac{}{n}}={\displaystyle \frac{}{y}},`$   
and expand the birth/death term in powers of the deviation $`y=nn`$

$`{\displaystyle \frac{1}{n}}{\displaystyle \frac{1}{n}}={\displaystyle \frac{y}{n^2}}{\displaystyle \frac{y^2}{n^3}}+\mathrm{}`$   
Now Eq. (45) becomes

$$\frac{P}{t}\frac{dn}{dt}\frac{P}{y}=\frac{y}{n^2}P\frac{y^2}{n^3}P+D\frac{^2P}{y^2}.$$
(46)   
Let us first assume that the average fitness grows faster than diffusively, that is, $`n\sqrt{t}`$. With this assumption, the dominant terms in Eq. (46) are

$$\frac{dn}{dt}\frac{P}{y}=\frac{y}{n^2}P.$$
(47)   
These terms balance when $`n(ty)^1yn^2`$. Using this scaling and balancing the remaining sub-dominant terms in Eq. (46) gives $`y\sqrt{t}`$. The combination of these results yields $`nt^{2/3}`$. This justifies our initial assumption that $`n\sqrt{t}`$. Now we write $`n=(ut)^{2/3}`$, with $`u`$ of order unity, to simplify Eq. (47) to

$$\frac{P}{y}=\frac{3y}{2u^2t}P.$$
(48)   
In terms of $`n=y+n`$ the solution is

$$P(n,t)=N\sqrt{\frac{3}{4\pi u^2t}}\mathrm{exp}\left\{\frac{3\left[n(ut)^{2/3}\right]^2}{4u^2t}\right\}.$$
(49)   
Thus the fitness distribution is again Gaussian, as in the supercritical case, but with the average fitness growing as $`n=(ut)^{2/3}`$. Finally, we determine $`u=\sqrt{3D}`$ by substituting $`n=(ut)^{2/3}`$ in Eq. (46) and balancing the sub-dominant terms.

The age distribution in the critical case can be obtained in similar manner as in the supercritical case. The approximations that were invoked to determine the age distribution in the supercritical case still apply. Consequently, the asymptotic form for $`C_n(a)`$ is still given by Eq. (40), and this gives the same expression for $`C(a)`$ after integrating over $`n`$, as in Eq. (42). Hence the average age is again $`B^1`$, as in Eq. (9).

## VII summary and discussion

We have introduced an age-structured logistic-like population dynamics model, which is augmented by fitness mutation of offspring with respect to their parents. Here fitness is quantified by the life expectancy $`n`$ of an individual. We found unusual age and fitness evolution in which the overall mutational bias leads to three distinct regimes of behavior. Specifically, when deleterious mutations are more likely, the fitness distribution of the population approaches a steady state which is an exponentially decaying function of fitness. When advantageous mutations are favored or when there is no mutational bias, a Gaussian fitness distribution arises, in which the average fitness grows as $`n=Vt`$ and as $`n=(3D)^{1/3}t^{2/3}`$, respectively.

Paradoxically, the average age of the population is maximal for a completely unfit population. Conversely, individuals are less long-lived for either positive or no mutational bias, even though the average fitness increases indefinitely with time. That is, a continuous “rat-race” towards increased fitness leads to a decrease in the average life span. As individuals become fit, increased competition results in their demise well before their intrinsic lifetimes are reached. Thus within our model, an increase in the average fitness is not a feature which promotes longevity.

Our basic conclusions should continue to hold for the more realistic situation where the mortality rate increases with age . The crucial point is that old age is unattainable within our model, even if individuals are infinitely fit. When the mutational bias is non-negative, old age is unattainable due to keen competition among fit individuals, while if deleterious mutations are favored, age is limited by death due to natural mortality. In either case, there are stringent limits on the life expectancy of any individual. To include an age-dependent mortality into our model, we may write the mortality term $`f_n(a)C_n(a,t)`$ instead of $`n^1C_n(a,t)`$ in Eq. (3), where $`f_n(a)`$ is the mortality rate for individuals of age $`a`$. Similarly, the term $`n^1P_n`$ in Eq. (5) should be replaced by $`𝑑af_n(a)C_n(a,t)`$. However, these generalized terms play no role for large $`n`$, since $`f_n(a)`$ is a decreasing function of $`n`$ and old age is unattainable.

We gratefully acknowledge partial support from NSF grant DMR9632059 and ARO grant DAAH04-96-1-0114.

## A Bounds for the average age

The upper bound, $`A<(\gamma N)^1`$, follows immediately from Eq. (30), so we just prove $`A>A_{\mathrm{min}}`$. We have

$`A_{\mathrm{min}}`$ $`=`$ $`B^1={\displaystyle \frac{1}{b+b_{}(1+\mu ^2)}}`$   
$`A`$ $`=`$ $`{\displaystyle \frac{1}{b+2b_{}\mu }}{\displaystyle \frac{1}{(b+2b_{}\mu )^2}}{\displaystyle _0^1}𝑑xx^{\frac{1}{b+2b_{}\mu }}\left({\displaystyle \frac{1\mu }{1\mu x}}\right)^2\mathrm{exp}\left\{{\displaystyle \frac{1}{b_{}\mu }}\left({\displaystyle \frac{1}{1\mu x}}{\displaystyle \frac{1}{1\mu }}\right)\right\}.`$   
Using these expressions and performing elementary transformations we reduce the inequality $`A>A_{\mathrm{min}}`$ to

$$_0^1\frac{dx}{b_{}(1\mu x)^2}x^{\frac{1}{b+2b_{}\mu }}\mathrm{exp}\left\{\frac{1}{b_{}\mu }\left(\frac{1}{1\mu x}\frac{1}{1\mu }\right)\right\}<\frac{b+2b_{}\mu }{b+b_{}(1+\mu ^2)}.$$
(A1)   
Let us now introduce the variable

$$v=\frac{1}{b_{}\mu }\left(\frac{1}{1\mu x}\frac{1}{1\mu }\right),$$
(A2)   
so that $`dv=dx/b_{}(1\mu x)^2`$, which varies in the range $`[0,V]`$, with $`V=\frac{1}{b_{}(1\mu )}`$. This simplifies Eq. (A1) to

$$_0^V𝑑ve^v\left[\frac{1V^1v}{1V^1\mu v}\right]^{\frac{1}{b+2b_{}\mu }}<\frac{b+2b_{}\mu }{b+b_{}(1+\mu ^2)}.$$
(A3)   
We now use the inequality

$$\left[\frac{1p}{1q}\right]^\nu <e^{(qp)\nu }$$
(A4)   
which holds for $`0<q<p<1`$ and $`\nu >0`$. This inequality is readily proven by taking the logarithm on both sides and using the expansion $`\mathrm{ln}(1u)=_{k1}u^k/k`$. Now we apply Eq. (A4) to the integrand in (A3) and then replace the upper limit $`V`$ in the integral by $`\mathrm{}`$ to give

$`{\displaystyle _0^V}𝑑ve^v\left[{\displaystyle \frac{1V^1v}{1V^1\mu v}}\right]^{\frac{1}{b+2b_{}\mu }}`$ $`<`$ $`{\displaystyle _0^V}𝑑v\mathrm{exp}\left\{vv{\displaystyle \frac{b_{}(1\mu )^2}{b+2b_{}\mu }}\right\}`$   
$`<`$ $`{\displaystyle _0^{\mathrm{}}}𝑑v\mathrm{exp}\left\{v{\displaystyle \frac{b+b_{}(1+\mu ^2)}{b+2b_{}\mu }}\right\}`$   
$`=`$ $`{\displaystyle \frac{b+2b_{}\mu }{b+b_{}(1+\mu ^2)}}.`$   
This completes the proof.

The lower bound $`A_{\mathrm{min}}`$ turns out to be very accurate in the case when mutations are slightly deleterious. To see this let us write $`b_+=1,b_{}=(1+ϵ)^2`$, where $`ϵ1`$. Replacing $`x`$ by the transformed variable $`v=ϵ^1(1+ϵx)^1`$ recasts the integral Eq. (30) as

$`ϵ^2{\displaystyle _0^{\frac{1}{ϵ(1+ϵ)}}}𝑑ve^v\left(1{\displaystyle \frac{ϵ^2v}{1ϵv}}\right)^{\frac{1}{b+2+2ϵ}}.`$ (A5)   
We now expand the integrand,

$`\left(1{\displaystyle \frac{ϵ^2v}{1ϵv}}\right)^{\frac{1}{b+2+2ϵ}}=1{\displaystyle \frac{ϵ^2v}{b+2+2ϵ}}{\displaystyle \frac{ϵ^3v^2}{b+2+2ϵ}}+𝒪(ϵ^4),`$   
replace the upper limit in the integral Eq. (A5) by $`\mathrm{}`$, and compute the resulting simple integrals explicitly to obtain a series expansion in $`ϵ`$ for the average age. Together with analogous expansions for $`A_{\mathrm{max}}`$ and $`A`$ we have

$`A_{\mathrm{max}}`$ $`=`$ $`{\displaystyle \frac{1}{b+2+2ϵ}}`$ (A6)   
$`A_{\mathrm{min}}`$ $`=`$ $`A_{\mathrm{max}}{\displaystyle \frac{ϵ^2}{(b+2+2ϵ)^2}}+{\displaystyle \frac{ϵ^4}{(b+2+2ϵ)^3}}+𝒪(ϵ^6)`$ (A7)   
$`A`$ $`=`$ $`A_{\mathrm{max}}{\displaystyle \frac{ϵ^2}{(b+2+2ϵ)^2}}+{\displaystyle \frac{ϵ^4}{(b+2+2ϵ)^3}}+{\displaystyle \frac{2ϵ^5}{(b+2+2ϵ)^3}}+𝒪(ϵ^6).`$ (A8)   
Thus the difference between the exact value and $`A_{\mathrm{min}}`$ is of order $`ϵ^5`$.

## B Transient behavior of the total density

Numerically, we found that in the subcritical case the total population density approaches the steady state value $`N_{\mathrm{}}=(b+2\sqrt{b_+b_{}})/\gamma `$ from below with a deviation that vanishes as $`t^{2/3}`$. We now explain this behavior by constructing a mapping between this transient behavior in the subcritical case and the transient behavior in the critical case. We start with the basic rate equation, Eq. (5). We may remove the term $`\gamma NP_n`$ through the transformation

$$Q_n(t)=P_n(t)\mathrm{exp}\left\{\gamma _0^t𝑑t^{}N(t^{})\right\},$$
(B1)   
which simplifies Eq. (5) to

$$\frac{dQ_n}{dt}=\left(b\frac{1}{n}\right)Q_n+b_+Q_{n1}+b_{}Q_{n+1}.$$
(B2)   
Next, the steady state behavior $`P_n\mu ^n`$ suggests replacing the $`Q_n`$ by $`R_n(t)=\mu ^nQ_n(t)`$. This also removes the asymmetry in the birth terms and gives

$$\frac{dR_n}{dt}=\left(b\frac{1}{n}\right)R_n+b_{}(R_{n1}+R_{n+1}),$$
(B3)   
where we use the shorthand notation $`b_{}=\sqrt{b_+b_{}}`$.

One cannot use the continuum approximation to determine the steady-state solutions for $`P_n`$ or $`Q_n`$. However, the continuum approximation is appropriate for the $`R_n`$. Then Eq. (B3) reduces to

$$\frac{R}{t}=\left(b+2b_{}\frac{1}{n}\right)R+b_{}\frac{^2R}{n^2},$$
(B4)   
which is very similar to Eq. (45). Hence we expect that the distribution of $`R_n`$ is peaked around $`n(3b_{})^{1/3}t^{2/3}`$. It proves convenient to make this scaling manifest. To this end we change variables once more,

$$S_n(t)=R_n(t)\mathrm{exp}\left\{(b+2b_{})t+\left(\frac{9t}{b_{}}\right)^{1/3}\right\},$$
(B5)   
to get

$$\frac{S}{t}=\left(\frac{1}{n}\frac{1}{n}\right)S+b_{}\frac{^2S}{n^2}.$$
(B6)   
Repeating the procedure of Sec. V we determine the asymptotic solution to Eq. (B6) as

$$S_n(t)\frac{1}{\sqrt{4\pi b_{}t}}\mathrm{exp}\left\{\frac{(nn)^2}{4b_{}t}\right\}.$$
(B7)   
To find the asymptotics of the total population density let us compute $`Q_n(t)`$. First, (B1) can be expressed as

$$\underset{n=1}{\overset{\mathrm{}}{}}Q_n(t)=N(t)\mathrm{exp}\left\{\gamma _0^t𝑑t^{}N(t^{})\right\}.$$
(B8)   
On the other hand,

$$\underset{n=1}{\overset{\mathrm{}}{}}Q_n(t)=\underset{n=1}{\overset{\mathrm{}}{}}\mu ^nR_n(t)=\mathrm{exp}\left\{(b+2b_{})t\left(\frac{9t}{b_{}}\right)^{1/3}\right\}\underset{n=1}{\overset{\mathrm{}}{}}\mu ^nS_n(t).$$
(B9)   
In the last sum, the factor $`\mu ^n`$ suggests that only terms with small $`n`$ contribute significantly. Although the asymptotic expression (B7) is formally justified only in the scaling region, where $`|nn|\sqrt{b_{}t}`$, the continuum approach typically provides a qualitatively correct description even outside this region. Therefore we take Eq. (B7) to estimate $`S_n`$ for small $`n`$. We find

$$\underset{n=1}{\overset{\mathrm{}}{}}\mu ^nS_n(t)\mathrm{exp}\left\{C\left(\frac{9t}{b_{}}\right)^{1/3}\right\},$$
(B10)   
where we use $`nt^{2/3}`$, as in the critical case, and $`C`$ is a constant. By substituting Eq. (B10) into Eq. (B9) we obtain

$$\underset{n=1}{\overset{\mathrm{}}{}}Q_n(t)\mathrm{exp}\left\{(b+2b_{})t(1+C)\left(\frac{9t}{b_{}}\right)^{1/3}\right\}.$$
(B11)   
Combining Eqs. (B8) and (B11) we arrive at the asymptotic expansion

$$\gamma _0^t𝑑t^{}N(t^{})=(b+2b_{})t(1+C)\left(\frac{9t}{b_{}}\right)^{1/3}+\mathrm{},$$
(B12)   
which implies

$$\gamma N(t)=b+2b_{}\mathrm{const}\times t^{2/3}+\mathrm{}$$
(B13)