Selection by differential mortality rates

The Theory of Evolution through Natural Selection [3, 12, 15] is central to the understanding of patterns observed in biological systems. The body of theory had gone through many revisions [17], which are natural, if we recall that even the mechanisms of inheritance were still unknown when it was first conceived.

Given the immense scope of Evolution Theory, there are variations in its formulations, depending on the biological phenomenon approached. Nevertheless, it is possible to identify a central invariant core that is present in almost any evolutionary explanation of biological observations. Authors in [20,21] indicate three fundamental components:

There is variation in morphological, physiological, and behavioural traits among members of a species (the principle of variation).

Thus, the concept of fitness (directly related to the number of offspring) plays a fundamental role in Evolution. This concept has evolved and has been used with different meanings [11] depending on the theoretical context. For example, in the framework of Population Genetics, fitness usually means the rate of increase/decrease in frequency of certain genotype in each generation [11, 19]. In studies of altruistic behaviour, the concept of individual fitness is expanded to include the fitness of close relatives, leading to inclusive fitness [18]. Some authors [10, 11] argue further that fitness should ultimately be expressed in terms of gene fitness.

In all of these formulations, the fitness value is influenced by many factors: degree of relatedness, form of reproduction, reproduction rate, mortality rates and others. In this article, we will focus on the mortality factors, that is, how selection operates when individuals have distinct mortality rates.

There are many examples of traits that influence mortality rates of organisms. Mutations affecting senescence mechanisms [5], resistance to toxin produced by prey [7, 27] and the degree of dominance of a gene that provides insecticide resistance [6] are examples of such traits. In fact, any adaptation for defence against predators or parasites (speed, thick shells, and quick reflexes) can be related to the increase of fitness through the reduction of mortality rates.

Many methods can be used to model evolution: Game theory [28], population genetics [16, 19] and dynamical systems related to Game Theory [25] are examples. To model selection through differential mortality rates, we are going to use an aspect-space approach [2]. Such approach has been used, for example, to model the evolution of host-specificity behaviour in host-parasite systems [1] and shows theoretical consistency with other approaches [2]. It consists, basically, of describing a population as distributed in an aspect space (or trait space) and writing the dynamics that contain the three fundamental conditions for evolution to occur [20, 21]. For the details on the construction of models using this approach we refer to [1, 2].

Here we consider scenarios in which all phenotypes in the aspect space have the same reproduction rates but possess distinct mortality rates. We deduce some conditions, depending on the distribution of mortality rates, for the permanence of the population. We also deduce some properties of equilibrium solutions, investigating the particular cases of directional, stabilising and disruptive selection. Finally, we show some differences and similarities of the process of selection through differential mortality rates when compared with the process of differential reproduction rates [2].

We begin by stating the basic assumptions of the model:

Individuals are distributed in a phenotype space Ω = [−M,M] ⊂ ℝ.

Hypothesis 1 simply means that we are describing the evolution of a quantitative trait. Given possible discontinuities in the morphogenesis [23], it is possible that similar quantitative traits lead to dissimilar phenotypes, leading to discontinuities in the mortality function. Hypothesis 2 and 3 are related to the discussion in consideration here, i.e., the focus of selection is on mortality rates and not on differential reproduction rates.

Hypothesis 4 simply states that there is imperfect heredity of traits. Variation on traits might come from recombination of genes, mutations or even changes in environmental factors. The change may be considered to be “small” in the timescale of reproduction and has no bias to either increasing or decreasing the trait described. Additive genetic models [29] can display this kind of behaviour as long as a quantitative trait is influenced by a sufficient large number of alleles. In [13] the authors estimate that a particular quantitative trait in Drosophila melanogaster (number of sensory bristles) are influenced at least by 38 distinct loci, a large enough number to generate unbiased variation through mutation.

Finally, hypothesis 5 includes intra-specific competition. Although models that include unlimited population growth can display a certain type of selection represented by distinct fractions of types [30], a more realistic selection may be obtained with the inclusion of competition and limited population growth.

Let v(τ,y) denote the number of individuals at time τ who belong to genotype y, and let V denote the diffusion coefficient, which in our case is a mutation parameter. In the equation below, we introduce the term rV containing the reproduction within the diffusive process, because at the genetic level, mutation occurs before reproduction. So new-borns of parents in genotype class x_i may belong to the nearby genotype classes, x_i±1. This makes a substantial difference also with standard reaction-diffusion equations in population theory.

The above assumptions lead us to the following partial differential equation (1) ∂v∂τ=rV∂2v∂y2+rv−m(y)v−αv∫Ωv dy {{\partial v} \over {\partial \tau}} = rV{{{\partial^2}v} \over {\partial {y^2}}} + rv - m(y)v - \alpha v\int_\Omega v\;dy defined for τ > 0 and y on an interval of real line.

In this section, we analyse what happens with the solution as time evolves when describing the equilibrium solution.

We start by rescaling the variables as follows τ=tr,y=xV,u(t,x)=αVrv(τ,y). \tau= {t \over r},\quad y = x\sqrt V,\quad u(t,x) = {{\alpha \sqrt V} \over r}v(\tau,y).

In this way, we end up with the following dimensionless equation (2) ∂u∂t=∂2u∂x2+[1−μ(x)]u−u∫Du dx {{\partial u} \over {\partial t}} = {{{\partial^2}u} \over {\partial {x^2}}} + \left[ {1 - \mu (x)} \right]u - u\int_Du\;dx in which μ(x)=m(xV)/r \mu (x) = m(x\sqrt V)/r , for all x ∈ D = [−L,L]. For this equation we are considering both homogeneous Dirichlet and Neumann boundary conditions.

We could easily check that ū ≡ 0 is an equilibrium solution of the previous equation. Assuming there is also a non-zero equilibrium solution, say ū(x), then it turns out that (3) u¯″(x)=[μ(x)+N−1]u¯(x),N=∫Du¯(x)dx. {\bar u^{''}}(x) = \left[ {\mu (x) + N - 1} \right]\bar u(x),\quad N = \int_D\bar u(x)dx.

By integrating the previous equation on the domain D, we obtain (4) N(1−N)=∫Dμ(x)u¯(x)dx−∫Du¯″(x)dx. N(1 - N) = \int_D\mu (x)\bar u(x)dx - \int_D{\bar u^{''}}(x)dx.

The last equality implies the following:

For non-zero solutions we have the following result:

We can also state the following result about the shape of a non-zero equilibrium solution ū(x):

Previous results show us that ū(x) and μ(x) have opposite monotonicity. Thus, by Chebyshev’s inequality for integrals we can conclude that ∫Dμ(x)u¯(x)dx≤12L∫Dμ(x)dx∫Du¯(x)dx=N2L∫Dμ(x)dx. \int_D\mu (x)\bar u(x)dx \le {1 \over {2L}}\int_D\mu (x)dx\int_D\bar u(x)dx = {N \over {2L}}\int_D\mu (x)dx.

Applying last inequality to Eq. (4) it turns out that (6) N≥1−12L∫Dμ(x)dx. N \ge 1 - {1 \over {2L}}\int_D\mu (x)dx.

It is not difficult to check that the same inequality is true when μ(x) is an even function whose restriction on [0,L] is a monotone function.

Eq. (6) implies that the existence of a non-zero equilibrium solution relies on how much the value of the integral changes according to changes in the values of L. When the order of change of the integral L^p with respect to L is less than one, 0 < p < 1, the inequality becomes more restrictive as L increases. On the other hand, when the integral of the mortality rate increases as L^1+s, s > 0, then the inequality becomes less restrictive as L increases.

From (14) and (15), the stability of a non-zero solution relies on the value of the unknown quantity N. In order to overcome this problem, we can consider the instability conditions of the zero equilibrium solution. In this case, Eq. (15) becomes (16) 12L∫Dμ(x) dx<1−θ12. {1 \over {2L}}\int_D\mu (x)\;dx < 1 - \theta_1^2.

Looking at Eq (14) we see that the condition for stability of non-zero equilibrium solutions is less restrictive than the one of the zero solution, namely 1−N−θ12<12L∫Dμ(x) dx<1−θ12. 1 - N - \theta_1^2 < {1 \over {2L}}\int_D\mu (x)\;dx < 1 - \theta_1^2.

This means that Eq. (2) could have a stable zero solution and, if it exists, a stable non-zero equilibrium solution as well.

Now, let us revert to the original variables in order to infer the influence of the original parameters in the instability conditions. In this case, Eq. (16) is rewritten as (17) 12M∫Ωm(x) dx<r−π2VrbM2 {1 \over {2M}}\int_\Omega m(x)\;dx < r - {{{\pi^2}Vr} \over {b{M^2}}} in which b = 1 in case of Dirichlet boundary conditions and b = 4 for Neumann conditions. From Eq. (17), increasing the mutation rate V leads to a more restrictive condition for the establishing of the population. For this to occur, the greater the value of V either the smaller the mortality rate m(x) or the bigger the value of r must be. Of course, this is true due to the boundary conditions imposed on the equation. When the mutation rate is small compared to the size M of the domain then what really counts is whether the mean value of the mortality is smaller than the natality rate r.

Let us assume that the natality rate is greater than the mortality rate, m(x) ≤ (1 − β)r for all x ∈ Ω. In this case, we could naively imagine the existence of an establishing population. However, performing the calculations under this condition, the criterion of instability of the zero solution becomes (18) π2VbM2<β. {{{\pi^2}V} \over {b{M^2}}} < \beta.

Thus, as we can clearly see, the mutation rate V plays an important role for the existence of a non-zero equilibrium solution.

In the next subsections, we study the circumstances under which the selection by differential mortality rate could lead to the three most basic types of evolution: directional, disruptive and stabilising selection. We remark that condition (16) is always satisfied when we consider Neumann boundary conditions.

4.1

Directional selection

Directional selection occurs when an extreme phenotype is selected over the others on the phenotype space. In our context, it means that the mortality rate function is monotone over the domain D = [−L,L]. The simplest function that describes this type of selection is given by μ(x)=h2L(L−x) \mu (x) = {h \over {2L}}(L - x) and for this function it turns out that ∫Dμ(x)dx=hL. \int_D\mu (x)dx = hL.

Thus, the criterion (16) for establishing a population becomes h<2(1−θ12) h < 2(1 - \theta_1^2) which in turn, for sufficient large values of L, gives just h < 2.

In other words, when the adaptive value is increasing linearly with respect to the measurement of the trait then a population is experiencing directional selection if the maximum of the mortality is less than the double of the natality rate (recall that for the dimensionless equation r = 1).

Now, let us consider the shape of the non-zero equilibrium solution. First of all, let us observe that, in the context of directional selection, equilibrium solutions on bounded domains are a kind of transient solution. As we have previously discussed, according to Proposition 5 if μ(x) is a decreasing function the equilibrium solution must be an increasing function in case of Neumann boundary conditions. Moreover, the extreme point −L is a local minimum, L is a local maximum of the solution and there is no other critical point in the domain.

In case of Dirichlet boundary conditions, there is a critical point in the interior of the domain. If x̄ is such that μ(x̄) = 1 − N, using Eq. (5), we can show that the critical point must be greater than x̄.

In Figure 2 we have the non-zero equilibrium solution of Eq. (2) with mortality rate given by μ(x) = 0.15(5 − x). The graph on the left is a numerical simulation of the equilibrium solution for Neumman boundary condition whereas, on the right, we see a numerical simulation for Dirichlet boundary condition. Clearly, the shapes of these solutions show that an extreme phenotype is favoured over the others.

Fig. 1

Fig. 2

Fig. 3

4.2

Disruptive selection

Now, let us consider the case in which the mortality rate is defined as μ(x)=max{h(a2−x2)a2,0},a≤L. \mu (x) = \max \left\{{{{h({a^2} - {x^2})} \over {{a^2}}},0} \right\},\quad a \le L.

This function describes a situation in which the extreme traits have greater adaptive value over the mean trait. That is, the mortality rate decreases with respect to distance to the mean trait. The zero phenotype has the greatest mortality rate μ(0) = h and phenotypes satisfying a ≤ |x| ≤ L have mortality rate zero.

For this function, it turns out that 12L∫Dμ(x) dx=2ah3L {1 \over {2L}}\int_D\mu (x)\;dx = {{2ah} \over {3L}} and thus the condition for an establishing population becomes ah<3L2(1−θ12), ah < {{3L} \over 2}\left({1 - \theta_1^2} \right), which in terms of the dimensional variables turns into ah<32(MV−π2VbM). ah < {3 \over 2}\left({{M \over {\sqrt V}} - {{{\pi^2}\sqrt V} \over {bM}}} \right).

In the last expression, the value of b depends on the type of boundary condition.

According to Eq. (4), we have 0 < N < 1 and, hence, 1 − N > 0. When we consider Neumann boundary conditions for Eq (2), by Proposition 2, both −L and L are critical points of the equilibrium solution. Since μ(−L) = μ(L) = 0 < 1 − N, it follows that these points are local minima of ū(x). On the other hand, since μ(x) is an even function the equilibrium solution ū must be an even function as well. Indeed, it cannot be odd, because it must be non-negative. Thus, it follows that x = 0 is a critical point of ū. According to Proposition 2 we claim that this critical point is a local minimum.

The same statement is true when we consider Dirichlet boundary conditions. That is, for a non-constant equilibrium solution, x = 0 is a local minimum for the solution ū(x). From Eq. (5) it follows that ū(−L) > 0 and ū(x) < 0 for x near zero. Thus, observing that μ(x) is increasing on (−L,0) we conclude that ū(x) has a local maximum on (−L,0). Hence, by symmetry, ū(x) has a local maximum on (0,L) as well.

In Figure 4, we see numerical solutions illustrating this behaviour.

Fig. 4

Fig. 5

4.3

Stabilising selection

Let us consider now the situation in which the mean phenotype has the greatest adaptive value over all the others. This situation can be described by the following mortality rate function: μ(x)=min{hx2a2,h}. \mu (x) = \min \left\{{{{h{x^2}} \over {{a^2}}},h} \right\}.

For this function, it is not difficult to check that ∫Dμ(x) dx=2Lh−2ah3. \int_D\mu (x)\;dx = 2Lh - {{2ah} \over 3}.

Now, assuming that a is a small fraction of the domain, say a = β L, then the criterion for the instability of the zero solution (16) becomes: h(1−β3)<1−θ12. h\left({1 - {\beta\over 3}} \right) < 1 - \theta_1^2.

Since μ(x) is even, according to Proposition 3, the non-constant equilibrium solution has a local maximum at x = 0. Further, according to Proposition 2, for Neumann boundary conditions, both −L and L are local minima of ū(x). In Figure 6, we see numerical solutions illustrating the typical behaviour of an equilibrium solution under evolution by stabilising selection.

Fig. 6

4.4

Survival of the flattest

We finally address another important case of selection: the survival of the flattest. That is, we consider the case in which the fitness landscape has two peaks: one low and flat and another one higher but narrower.

In terms of mortality rates, we are interested in discussing the shape of possible equilibrium solutions when the mortality function μ(x) has two local minima: one high and flat, and the other lower but narrower (see the left frame of Figure 7).

Fig. 7

As shown in previous sections, according to Proposition 2, the existence of local maximum points for an equilibrium solution relies on whether μ(x) is smaller than 1 − N or not. For this case, we have three possible scenarios: both local minima are less than 1 − N; both local minima are greater than 1 − N; just the lower local minimum is less than 1 − N. Here, we consider mainly the first scenario.

If we look at Eq. (5) we can see that the difference between μ(x) and 1−N controls the sign of the derivative of ū(x). It means that what is important here are both the size of the intervals in which μ(x) is greater or less than 1 − N as well as the difference between the mortality rate and 1 − N in each one of these intervals. Thus, at least hypothetically, we could have a non-zero equilibrium solution with either one local maximum, peaked near one of the local minima of μ(x), or two local maximum points peaked near both the local minimum points of μ(x).

However, more interesting is the role played by the mutation rate V in switching the maximum value of the equilibrium solution from one valley of μ(x) to the other. To illustrate this behaviour, consider the mortality rate function on the left in Figure 7. This function is given by μ(x)=min{1.1,1.1+0.06(x+2)(x+8)}+min{0,(x−4)(x−6)}. \mu (x) = \min \left\{{1.1,1.1 + 0.06(x + 2)(x + 8)} \right\} + \min \left\{{0,(x - 4)(x - 6)} \right\}.

We also set M = 10, r = 1 and α = 1. For this mortality function, the flattest valley is on the interval [−10,0].

To measure the influence of V on the equilibrium solution ū(x), let us define p by (19) p=∫−M0u¯(x)dx∫−MMu¯(x)dx. p = {{\int_{- M}^0\bar u(x)dx} \over {\int_{- M}^M\bar u(x)dx}}.

In words, p defined as Eq. (19) measures the proportion between the number of individuals on the left over the total number of individuals. Thus, letting the mutation rate V vary we capture its influence over the equilibrium solution. In Figure 7 right, we can see that as V increases the equilibrium solutions switch from the right, i.e. a low p means more population on the rightmost half of the domain, to the left, i.e. a large p means more population on the leftmost part of the domain of the aspect space. Thus for small values of V the equilibrium solution is peaked near the minimum point of the narrower valley. Further, increasing V we clearly see the transition from the narrower valley to the flatter one, shown in Figures 8, 9.

Fig. 8

Fig. 9

In Figure 8 we have the graphs of two equilibrium solutions when V = 0.4 and V = 0.52. For these values of V, both equilibrium solutions have a global maximum near x = 5. However, the solution on the right has also a local maximum at x = − 5.

In Figure 9 we have the graphs of two equilibrium solutions when V = 0.54 and V = 0.7. For these values of V, both equilibrium solutions have a global maximum near x = −5. However, the solution on the left has also a local maximum at x = 5.

In this paper we presented some results relative to evolutionary selection through differential mortality rates. The question of population permanence has been discussed, given fixed mortality distributions and boundary conditions. Two main effects may push the population to extinction: high mutation rates leading to population dissipation or high mortality rates resulting in insufficient reproduction.

The first effect is closely related to the classical problem of critical patch size for population permanence under diffusion [4,22]. In this case, the coefficient of phenotype change, V, plays the role of a diffusion/dispersal constant. Patch size is related with how narrow the adaptation peak occupied by the biological population is. Under very narrow adaptation peaks, only species with low mutation rates can survive. This effect is present in all permanence conditions calculated in this paper. The second effect is the more obvious restriction on average mortality when compared with the reproduction rate, which can be explicitly observed in inequality (17). There is evidence that pathogens can be controlled by causing an ‘error threshold catastrophe’, by pushing mutation rates using mutagens [9, 26].

The model has shown coherence in the results of all three particular cases analysed: directional, stabilising and disruptive selection, with results mimicking qualitatively the expected and observed biological results, see Chapter 12 of [15]. Selection of quasispecies was another characteristic displayed by the model, which is also observed in other theoretical approaches [24], experiments [8] and in an aspect space model with differential reproduction rates [2].

However, an important difference in the quasispecies selection based on differential mortality rates was found. The transition from the selection of one quasispecies to another one is smooth. While replicator and differential reproduction models lead to abrupt transitions from one type to another, selection through differential reproduction rates can display coexistence of types and smooth transitions from one type to another (as shown in Figure 7). Also, in this type of selection, increased mutation rates favoured the ‘selection of the flattest’, being coherent with the results obtained in other models and experiments [24, 25].

There was also an important difference observed between the dynamics of selection through mortality rates and selection based on differential reproduction rates. In the scenario with just one adaptive peak, regardless of mutation rates, selection of the fittest is observed in the model with differential mortality rates. The same cannot be said to occur in classical approaches [14] or even in the aspect space model with focus on reproduction rates [2]. In such models, the mutation rate must be below a certain threshold for the fittest to be selected, a result that is coherent with biological observations.

Since mortality is just one of the many components of fitness in real systems, the results obtained by the model should not be considered as a contradiction, since a mutation rate below a certain threshold is still needed for the selection of the fittest when its other components are included. Finally, we point out that the analytical results concerning the aspect space model are original and represent the important development in understanding the dynamics of such models.