Michel Hulot

Loreau

Céline Hauzy

Mathias Gauduchon
email: mathias.gauduchon@univmed.fr

Florence D Hulot
email: florence.hulot@u-psud.fr

Michel Loreau
email: michel.loreau@mcgill.ca

Density-dependent dispersal and relative dispersal affect the stability of predator-prey metacommunities

Keywords: dispersal behaviours, spatial synchrony, top-down control, regional variability, persistence

Density-dependent dispersal and relative dispersal affect the stability of predator-prey metacommunities

INTRODUCTION

Understanding the mechanisms that govern the stability of metacommunities is a key challenge to understand the response of ecosystems to fragmentation. Metacommunities are regional sets of local communities linked by dispersal, i.e., by the movement of individuals from a patch (emigration) to another (immigration) [START_REF] Leibold | The metacommunity concept: a framework for multi-scale community ecology[END_REF]). Both local and regional processes affect the stability of metacommunities. For instance, in predator-prey metacommunities, local processes such as local density regulation mechanisms drive local predator-prey population fluctuations, which may lead to local extinctions when communities are isolated [START_REF] Gause | The Struggle for Existence[END_REF][START_REF] Lotka | Elements of Physical Biology[END_REF][START_REF] Rosenzweig | Graphical representation and stability conditions of predator-prey interactions[END_REF][START_REF] Volterra | Variations and fluctuations in the numbers of individuals in animal species living together[END_REF]. But dispersal, which is a regional process, may allow regional persistence of the prey and the predator [START_REF] Huffaker | Experimental studies on predation: dispersion factors and predatorprey oscillations[END_REF] because it allows recolonization of vacant patches [START_REF] Levins | Extinction. Some mathematical questions in biology[END_REF] or rescue effects [START_REF] Brown | Turnover Rates in Insular Biogeography: Effect of immigration on extinction[END_REF]. This positive effect of dispersal on the regional stability of predator-prey metacommunities requires that dynamics are spatially asynchronous. If local populations fluctuate strongly and synchronously, the whole metapopulation may go extinct because of simultaneous extinction of local populations. But if fluctuations of local densities are asynchronous, recolonizations or rescue effects can occur and the persistence and stability of the metacommunity is improved [START_REF] Comins | Prey-Predator Models in Spatially Heterogeneous Environments[END_REF][START_REF] Hassell | Spatial structure and chaos in insect population dynamics[END_REF]Jansen 2001).

But dispersal and synchrony between communities are not independent processes. On the one hand, dispersal may increase spatial synchrony because exchanges of individuals between communities tend to reduce differences in population sizes, especially when dispersal is high [START_REF] Bjornstad | Spatial population dynamics: analyzing patterns and processes of population synchrony[END_REF]. On the other hand, weak dispersal may also synchronise dynamics in predator-prey metacommunities with nonlinear density regulation [START_REF] Bjornstad | Cycles and synchrony: two historical 'experiments' and one experience[END_REF]Jansen 1999). Thus, both dispersal and local density regulation can synchronise community dynamics. Since there is empirical evidence for asynchronous spatial dynamics in predator-prey metacommunities (Holyoak and Lawler 1996a;Holyoak and Lawler 1996b;[START_REF] Janssen | Metapopulation dynamics of a persisting predator-prey system in the laboratory: Time series analysis[END_REF]Nachman 1991;[START_REF] Van De Klashorst | A demonstration of asynchronous local cycles in an acarine predator-prey system[END_REF], some mechanisms maintain persistent differences between local densities and counteract the synchronising effect of dispersal. These mechanisms have been explored using a number of models, and involve either stochastic or deterministic processes. Asynchrony between patches of the same quality can be due to independent stochastic local factors [START_REF] Taylor | Metapopulations, Dispersal, and Predator-Prey Dynamics: An Overview[END_REF], such as demographic stochasticity in small populations (Nachman 1987a) and stochasticity of dispersal between patches or of environmental variations in population parameters [START_REF] Crowley | Dispersal and the stability of predator-prey interactions[END_REF][START_REF] Reeve | Environmental variability, migration, and persistence in host-parasitoid systems[END_REF]Reeve 1990). These models show that an intermediate dispersal intensity of the prey and the predator insures metacommunity persistence [START_REF] Crowley | Dispersal and the stability of predator-prey interactions[END_REF]Nachman 1987b;Nachman 1991;[START_REF] Reeve | Environmental variability, migration, and persistence in host-parasitoid systems[END_REF]Zeigler 1977 reviewed in Holyoak andLawler 1996b). Dispersal must be sufficiently high to allow recolonization of extinct patches, but sufficiently low to maintain the desynchronising effect of stochastic local processes. Spatial asynchrony can also emerge from the interaction between local dynamics and dispersal (Jansen 1995). In deterministic predator-prey metacommunities with nonlinear density regulation, spatial asynchrony between patches arises when dispersal rates are either low or intermediate. In these models, population dynamics are synchronous when dispersal is high, and temporal variation in regional densities is decreased when dispersal rates are decreased.

In this paper, we investigate the influence of two key aspects of dispersal behaviours on the stability of predator-prey metacommunities: the density dependence in dispersal, and the ratio of predator and prey dispersal rates, which we name "relative dispersal". Densitydependent dispersal, i.e. dispersal changes in response to changes of species densities in the donor patch, has been observed in several experimental systems (see for instance [START_REF] Bernstein | Prey and predator emigration responses in the acarine system Tetranychus urticae-Phytoseiulus persimilis[END_REF][START_REF] French | Density-dependent dispersal in host-parasitoid assemblages[END_REF][START_REF] Hauzy | Intra-and interspecific densitydependent dispersal in an aquatic prey-predator system[END_REF]. It implies a direct effect of local processes (population dynamics) on regional processes (dispersal). A species' population density can affect its own dispersal (intraspecific density-dependent dispersal) or the dispersal of another species (interspecific density-dependent dispersal). Recent models including density-dependent dispersal in predator-prey metacommunities yielded contrasted results. [START_REF] Mccann | Population Outbreaks in a Discrete World[END_REF] suggested that local predator population outbreaks do not depend on interspecific density-dependent dispersal. By contrast, [START_REF] Li | Impact of predator pursuit and prey evasion on synchrony and spatial patterns in metapopulation[END_REF] showed that interspecific density-dependent dispersal can reduce spatial synchrony in both the prey and the predator. Most of these studies, however, have assumed that prey and predator dispersal rates are equal. But empirical data suggest that dispersal abilities of prey and predator can be very different. The large range of predator-prey body-size ratios observed in real food webs [START_REF] Brose | Consumer-resource body-size relationships in natural food webs[END_REF] and the allometric relationship between body size and the scale of species' movements [START_REF] Rooney | A landscape theory for food web architecture[END_REF] or dispersal abilities [START_REF] Jenkins | Does size matter for dispersal distance?[END_REF] suggest that the relative dispersal rate of the predator compared to the prey can vary within a large range. In some ecosystems, scale of movement is strongly positively related to trophic position, suggesting a higher dispersal rate for the predator than for the prey. In other systems, predators can even be less mobile than their prey, suggesting a higher dispersal rate for the prey than for the predator. Differences in the dispersal abilities of the prey and the predator can have important consequences, and yet have received little attention so far.

Here we develop a deterministic metacommunity model of a prey and a predator that includes density-dependent dispersal. We compare the effects of both intra-and interspecific density-dependent dispersal on the dynamics and the stability of the prey and the predator at the regional scale with those of constant dispersal. We study the effects of density-dependent dispersal in two contrasting cases of relative dispersal: either predator dispersal is distinctly higher than prey dispersal or prey dispersal is distinctly higher than predator dispersal. These comparisons allow us to assess the effects of density-dependent dispersal and relative dispersal on the relationship between dispersal and the stability of predator-prey metacommunities.

THE MODEL

We consider a patchy landscape consisting of several localities. Each locality is occupied by a local predator-prey community. Within each community, local dynamics follow a Rosenzweig-MacArthur predator-prey model [START_REF] Rosenzweig | Graphical representation and stability conditions of predator-prey interactions[END_REF]: prey growth is logistic, the predator has a Holling type II functional response and its mortality rate is constant. Local communities are linked by dispersal of both species from one patch to another. A collection of such communities constitutes a predator-prey metacommunity. For the sake of simplicity, and following Jansen (1995Jansen ( , 2001)), we study a 2-patch predator-prey metacommunity in which the two patches are identical. The dynamics of the prey and the predator in the two patches (1 and 2) are given by a system of ordinary differential equations: where N i and P i denote the population densities of the prey and the predator, respectively, in patch i =1, 2. As patches are identical, species parameters are the same in both patches. r and K denote the growth rate and carrying capacity of the prey; b and c are the saturation value and the half-saturation constant of the predator functional response; m is the predator mortality rate. Conversion efficiency from prey to predator has been set to 1. This does not constitute an ecological assumption but arises from a simple mathematical rescaling of N i or

°°°°°°°°°°® 1 • © § 1 • © § 1 1 1 2 2 2 2 2 2 2 2 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 1 1 1 1 1 1 1 1 2 2 2 1 1 1 1 1 1 1 1 1 , , , , 1 ,
P i .
The last two terms of equations describe dispersal from and to each patch. We make the ecological hypothesis that migration is not costly; as a consequence, the emigration term in one patch corresponds exactly to the immigration term in the other. This assumption fits in well with dispersal through corridors that allow fast moves with limited mortality risks.

Dispersal rates are determined locally: the per capita dispersal rates from patch i of the prey D n (N i , P i ) and of the predator D p (N i , P i ) are functions of species densities in the donor patch i only. The dispersal rate of each species may depend on its own local density, which we refer to as "intraspecific density-dependent dispersal", and on the local density of the other species, which we refer to as "interspecific density-dependent dispersal" [START_REF] Hauzy | Intra-and interspecific densitydependent dispersal in an aquatic prey-predator system[END_REF]. Our mathematical analysis deals with the general model, but our numerical simulations investigate the separate effects of intraspecific and interspecific density dependence in dispersal. The dispersal rates then simplify into D n (N i ) or D n (P i ) and D p (P i ) or D p (N i ). In numerical simulations, density dependence in dispersal is modelled through Hill functions, which have been chosen for their convenient plasticity: the per capita dispersal rate from patch i of the prey writes ( ) In this paper, we explore density-dependent dispersal behaviours that mimic those observed in several experiments [START_REF] Bernstein | Prey and predator emigration responses in the acarine system Tetranychus urticae-Phytoseiulus persimilis[END_REF][START_REF] French | Density-dependent dispersal in host-parasitoid assemblages[END_REF][START_REF] Hauzy | Intra-and interspecific densitydependent dispersal in an aquatic prey-predator system[END_REF], i.e., prey dispersal increases with prey and predator densities (positive intra-and interspecific density-dependent dispersal in the prey), and predator dispersal increases with predator density (positive intraspecific density-dependent dispersal in the predator) and decreases with prey density (negative interspecific density-dependent dispersal in the predator). In terms of model parameters, this implies x nn > 0, x np > 0, x pp > 0 and x pn < 0. We conduct both a mathematical and a numerical analysis of the model.

MATHEMATICAL ANALYSIS

Stability of the symmetric equilibrium

We first focus on equilibriums of the model, and more specifically on the symmetric equilibrium, i.e., that for which densities are equal in patches 1 and 2 (N 1 =N 2 and P 1 =P 2 ).

We analyse the stability of this equilibrium using the linearization method developed by [START_REF] Jansen | The dynamics of two diffusively coupled, identical Lotka-Volterra patches Theoretical Aspects of Metapopulation Dynamics[END_REF]. The original method considers a predator-prey metacommunity model with constant dispersal. We extend the method to general forms of density-dependent dispersal and apply it to our model (Appendix A).

The stability of the symmetrical equilibrium depends only on the local parameters r, b, c, m and K. The equilibrium is stable if and only if the prey carrying capacity K is below a threshold value:

m b m b c K ) (
(Appendix A). For higher values of the prey carrying capacity, it is always unstable and a stable limit cycle appears. These results for densitydependent dispersal generalize those obtained for constant dispersal by Jansen (1995) and correspond to the well-known paradox of enrichment in the non-spatial MacArthur-Rosenzweig model. Neither prey nor predator dispersal affects the stability of the symmetric equilibrium and this conclusion is unaffected by density dependence in dispersal. A metacommunity exhibits the same sudden destabilisation of its dynamics when K increases beyond its threshold whether dispersal is constant or density dependent.

Regional densities at the symmetric equilibrium population densities of the prey and the predator at the symmetric equilibrium. From equations (1), we see that emigration from patch i and immigration to patch i cancel each other out for each species, so that the regional densities of the prey and the predator are independent of dispersal. Thus, prey and predator densities at the symmetric equilibrium are unaffected as well by dispersal (Jansen 1995), whether dispersal is constant or density dependent. When the metacommunity has reached equilibrium, local population densities are the same as in isolated communities. Asymptotic population densities at the symmetric equilibrium depend only on local processes. In our two-patch predator-prey metacommunity, when K is below the instability threshold, the asymptotic regional densities are simply equal to the values they would have in an isolated patch.

When

m b m b c K ! ) (
, There is no stable equilibrium any longer, and the prey and the predator show persistent oscillations through time. We used numerical simulations to investigate the effect of density-dependent dispersal on prey and predator asymptotic dynamics at the scale of the metacommunity.

Methods

We studied the various types of dispersal separately, i.e., we studied constant dispersal in both species, intraspecific density-dependent dispersal in either the prey or the predator, and interspecific density-dependent dispersal in either the prey or the predator. For each of these five types of dispersal, we also studied two contrasting relative dispersal rates, one in which the maximal dispersal rate of the predator, d p , is lower than that of the prey, The system of ordinary differential equations (1) was numerically simulated by means of C programs using the GNU Scientific Library. Numerical integrations used the embedded Runge-Kutta Prince-Dormand method. Based on preliminary inspection of numerous simulations, we estimated a time delay before population dynamics entered its asymptotic regime. This "asymptotic time" (t asympt ) was much wider for low dispersal rates;

it was then to 2 x 10 4 time units for high dispersal rates (Min (d n , d p ) > 10 -4 ) and 2 x 10 5 time units for low dispersal rates (Min (d n , d p ) < 10 -4 ). Our results were obtained from population dynamics observed between t asympt and t asympt + 10 4 , except for the first Lyapunov exponent (see below), which was calculated on dynamics between t asympt + 10 4 and t asympt + 2.10 4 .

All the variables and statistics considered in this study are related to the asymptotic dynamical regime and are computed on the attractors reached by population dynamics.

Therefore, we needed a method to detect and follow the system's attractors as some control parameter (here, the maximal dispersal rates) is varied within a study interval. Except for catastrophic changes in the configuration of attractors, smooth variations in control parameters generate smooth "continuation branches", which correspond to the curves (or pieces of curves) shown on figures 2, 3, 4, 5, S3 and S4. Our model is complicated by the possibility of multiple coexisting attractors and chaotic regimes. Instead of using techniques based on numerical continuation of attractors [START_REF] Kuznetsov | Elements of Applied Bifurcation Theory[END_REF]), which proved inefficient and intractable, we used a "sieve" technique of exploration, followed by the ex post facto reconstruction of "continuation branches". For each parameter value, the initial densities of the prey and the predator in the two patches were randomly sampled within ranges [0;20] for the prey and [0;10] for the predator as these ranges were commonly observed in the numerical simulations. This random sampling allowed us to explore the various basins of attraction of the possibly multiple attractors, and hence to detect the main attractors (i.e., those with a wide enough basin of attraction) by means of simulations runs from these initial population densities. The drawback of this method is that only one initial condition is Analyses of dynamics were performed under the Matlab environment (R2008a). We studied two aspects of metacommunity dynamics:

(1) the nature of the attractor, and (2) the degree of spatial synchrony. First, we determined whether population fluctuations were periodic or chaotic by numerically estimating the first Lyapunov exponent. Computation of the Lyapunov exponent in fast oscillating models like ours gives easily rise to an important computational error accumulation. For the sake of numerical accuracy, the first Lyapunov exponent was estimated on a Poincaré section of the orbit, as the exponential increase of deviations from the iterations of the Poincaré map [START_REF] Wolf | Determining Lyapunov Exponents from a Time-Series[END_REF]. In each simulation, we computed the Poincaré section that intersects the trajectory orthogonally at the asymptotic time t asympt . We checked that the 10 4 -time-unit interval we used to compute the Lyapunov exponent provided satisfactory converge of the estimate. A positive value of the first Lyapunov exponent reveals a chaotic attractor, while a negative value indicates a stable equilibrium on the Poincaré map, which corresponds to a stable limit cycle for the population dynamics in 4 . Second, to quantify the effects of dispersal behaviours on spatial synchrony in the asymptotic regime, we used the Pearson moment correlation coefficient between local densities in the two patches. As an alternative synchrony estimate, the correlation of the per capita population growth rates is recommended when population dynamics exhibits long-term transient trends, as often encountered in field or experimental measurements [START_REF] Bjornstad | Spatial population dynamics: analyzing patterns and processes of population synchrony[END_REF]. But as we studied the asymptotic regime here, Pearson correlation between local densities is appropriate. The higher spatial synchrony, the closer to 1 the correlation coefficient. We investigated spatial synchrony for both the prey and the predator.

Because many definitions of stability exist in the literature (see e.g. a review in [START_REF] Loreau | A new look at the relationship between diversity and stability[END_REF], the stability of the metacommunity can be defined in several ways. Two stability properties are regional persistence and the temporal variability of regional densities. Because our model is purely deterministic and does not allow for extinctions, persistence cannot be estimated directly. Persistence is here assessed indirectly by means of the minimal regional densities reached in the asymptotic regime; the lower these minimal densities, the higher the vulnerability of species to stochastic extinction, and the lower the probability of persistence. Results are presented in the supplementary material (supplementary material, fig. S3 andS4). We also measured the temporal variability of regional densities by their coefficient of variation (CV). To disentangle the effects of the temporal mean of regional densities and of their standard deviation, however, we present both the temporal mean and the CV of regional densities. Stability increases when the CV decreases.

In order to investigate the robustness of our results with regard to changes in dispersal parameters, we used a factorial design in which the independent variables were the shape of the density dependence of the dispersal rates (|x αβ |), their half-saturation parameters

Metacommunity dynamics

We first focus on two descriptors of metacommunity dynamics, the nature of the dynamics (chaotic or periodic), and the spatial synchrony between patches. The first

Lyapunov exponent characterises the nature of dynamics of the whole community. The dynamics of the two species cannot be separated; they are either both periodic or both chaotic. Spatial synchrony was calculated for each species separately, but proved to be completely similar between the prey and the predator in all simulations. Therefore we present below results for spatial synchrony of the prey only. For the sake of simplicity, we present results for two contrasted relative dispersal rates but these results are robust to wide variations in relative dispersal rates (supplementary material, fig. S1 andS2). At high dispersal rates, population dynamics are synchronous and periodic for all types of dispersal behaviours (fig. 2). The metacommunity model is then equivalent to a non-spatial predator-prey model. In contrast, at low and intermediate dispersal rates, synchrony depends both on the relative dispersal rate of the prey and the predator and on the type of dispersal.

When dispersal is density-independent (constant) and higher in the prey than in the predator (fig. 2A), spatial synchrony decreases monotonically with dispersal decrease, through both desynchronisation of periodic dynamics and the emergence of chaotic dynamics at intermediate dispersal rate. This decrease is rather low and the minimum value of synchrony is still high. In contrast, when predator dispersal is higher than prey dispersal (fig. 2D), the dynamics can remain perfectly synchronous for all dispersal rates or become very asynchronous, provided that the dispersal rates are not too high. Indeed, for low and intermediate dispersal rates, the metacommunity dynamics shows two coexisting periodic attractors, the first one characterized by a perfect spatial synchrony while the second one shows a strong asynchrony between patches. The synchrony of the later decreases with increasing dispersal rates until the attractor disappears, letting the synchronous attractor alone. This coexistence of synchronous and asynchronous attractors has been observed by Jansen (1995). No chaotic dynamics is observed for this relative dispersal rate.

Density-dependent dispersal in the prey (fig. 2B,2E) affects the dynamics only when prey dispersal is higher than predator dispersal (d n > d p ) (fig. 2B). When predator dispersal is higher than prey dispersal (d p > d n ) (fig. 2E), the dynamics is very similar to the case where dispersal is constant (fig. 2D). When 5A). On the contrary, when prey dispersal is subject to interspecific density dependence, population dynamics are strongly desynchronised whether the dynamics be periodic or chaotic (fig. 2B).

Whatever the relative dispersal rates, intraspecific density-dependent dispersal in the predator does not affect spatial synchrony (fig. 2C,2F). The effect of density-dependence in predator dispersal shows up when it is interspecific and when predator dispersal is higher than prey dispersal (fig. 2F). Indeed, the bistable regime observed for constant dispersal disappears through the loss of the synchronous limit cycle. For low and intermediate dispersal rates, the remaining limit cycle shows strong asynchronous dynamics. The dynamics abruptly become synchronous at intermediate dispersal rate and the transition is accompanied with the emergence of a new chaotic window. Thus, for interspecific densitydependent dispersal in the predator, spatial synchrony dependence upon dispersal rate simplifies into a global increasing relationship.

Regional densities

The temporal mean of regional densities of the prey and the predator are affected by dispersal rates in the same way for all types of dispersal behaviour and of relative dispersal.

For the same range of dispersal rates, prey regional density decreases (fig. 3) while predator regional density increases (fig. 4) in a correlated way. This suggests that dispersal modifies the strength of the impact of predator populations on prey populations at the regional scale.

First, changes in prey and predator regional densities occur for dispersal rate in chaotic windows and for dispersal rate close to chaotic window but at which dynamics are periodic (fig. 3A, 3B, 3C, 3F; fig. 4A,4B,4C,4F). Second, these changes occur for asynchronous limit cycle which exists when predator dispersal is higher than prey's one (fig. 3D,3E,3F;

fig. 4D,4E,4F). While the asynchronous limit cycle deeply differs from the synchronous one in terms of spatial synchrony, it shows identical mean regional densities at low dispersal rates. But when dispersal rates increase, prey regional density slowly decreases while predator regional density increases. These changes in mean regional densities are not accompanied with chaotic dynamics.

Moreover, these mean regional densities can be quantitatively affected by dispersal behaviours and relative dispersal: the mean regional density of the prey decreases and the 4F). As a consequence, density-dependent dispersal and relative dispersal can increase the strength of the impact of predator populations on prey populations at the regional scale, and thus increase the control of the predator on the prey.

Regional variability

The regional variability of prey and predator populations, as measured by the CV of their regional densities, may be strongly affected by density-dependent dispersal through both spatial synchrony and predator control on the prey. We should expect low regional variability when local dynamics are asynchronous and high regional variability when local dynamics are synchronous. Predator control on the prey should also affect the regional variability through changes in mean regional densities: this control should decrease predator regional variability and should increase prey regional variability.

When dispersal is constant, the regional variability of the prey and the predator is modified by their relative dispersal (fig. 3A, 3D; fig. 4A,4D). When prey dispersal is higher than predator dispersal (d n > d p ) (fig. 3A; fig. 4A), the regional variability of the prey and of the predator at low dispersal rates is slightly lower than their regional variability at high dispersal rate. This is caused by desynchronisation of dynamics at low dispersal rates for both species (fig. 2A). But at intermediate dispersal rate, the regional variabilities of the prey and the predator vary in different way. Two mechanisms act in opposite direction for the prey and in the same direction for the predator. The increase in asynchrony tends to lower regional variability of both the prey and the predator, whereas predator control increases regional variability of the prey and decreases the predator's. As a result, the regional variability of the prey increases slightly while the predator's decreases strongly.

When predator dispersal is higher than prey dispersal (d p > d n ) (fig. 3D; fig. 4D), the regional variability of the prey and the predator can decrease when dispersal rates decrease.

Indeed, at intermediate and low dispersal rates, two limit cycles coexist: the synchronous one has high regional variability and the asynchronous one shows low regional variability. Moreover, regional variability of the prey and the predator on this asynchronous limit cycle vary with increasing dispersal rates. The increase in asynchrony tends to lower regional variability of both the prey and the predator, whereas predator control increases regional variability of the prey and decrease the predator's. As a result, the regional variability of the prey increases slightly while the predator's decreases strongly. Thus, the effect of dispersal on a species regional variability depends on its trophic position.

When dispersal rates are high, as in a well-mixed system, the regional variabilities of the prey and the predator are not affected by the type of dispersal behaviour and by the relative dispersal. But when dispersal rates decrease, density-dependent dispersal can affect regional variabilities of both species, depending on their relative dispersal rates (fig. 3B, 3C, 3E, 3F; fig. 4B,4C,4E,4F).

Density-dependent dispersal in the prey can modify regional variability when prey dispersal is higher than predator dispersal (d n > d p ) (fig. 3B; fig. 4B). Intraspecific densitydependent dispersal in the prey can strongly modify regional variability of both species by comparison to constant dispersal. At intermediate dispersal rates, prey and predator regional variabilities vary in opposite way when |x αβ |=4, because of the increase of spatial synchrony 20 (fig. 2B) and of predator control on the prey (fig. 3B; fig. 4B). Prey regional variability increases or decreases compared to constant dispersal, depending on the relative effects of these two mechanisms that act in opposite direction. On the contrary, these two mechanisms both decrease predator regional variability. At low dispersal rates, intraspecific densitydependent dispersal in the prey increases regional variability of both species because of spatial synchrony of periodic dynamics when |x αβ |=4 and when |x αβ |=1 (Appendix B, fig. 5B,C). The effects of interspecific density-dependent dispersal in the prey (fig. 3B; fig. 4B) on the relationship between dispersal rates and regional variability are qualitatively the same as under constant dispersal (fig. 3A; fig. 4A), but they are stronger. Indeed, prey and predator regional variabilities are strongly lower at low dispersal rates than at high dispersal rates because of strong spatial asynchrony at low dispersal rates. The transition between high and low dispersal rates occurs with increase of spatial asynchrony and of predator control on the prey. As a consequence, prey regional variability is increased whereas predator regional variability is decreased compared to constant dispersal.

Only interspecific density-dependent dispersal in the predator modifies regional variability when predator dispersal is higher than prey dispersal (d p > d n ) (fig. 3C, 3F; fig.

4C, 4F). Intraspecific density-dependent dispersal in the predator does not affect prey and predator regional variability (fig. 3C, 3F; fig. 4C,4F), which is consistent with the results on spatial dynamics and regional densities (fig. 2-4C, 2-4F). When predator dispersal depends on prey density, the synchronous limit cycle disappears at low and intermediate dispersal rates, leaving only the asynchronous limit cycle. As a consequence, prey and predator regional variabilities are strongly lower at low dispersal rates than at high dispersal rates, whatever the initial densities are. Moreover, between low and high dispersal rates, prey regional variability results of a balance between spatial synchrony and predator control on the prey, whereas these two mechanisms act in same way to decrease predator regional variability.

DISCUSSION

Our results when dispersal rates of the prey and the predator are constant, such as the potential coexistence of several attractors, are consistent with those obtained by Jansen (1995Jansen ( , 2001)). They further show that density-dependent dispersal in both the prey and the predator and their relative dispersal rates may affect the nature of dynamics (periodic vs chaotic) and the spatial synchrony and stability of species when the symmetric equilibrium of the metacommunity is unstable and populations fluctuate over time.

Dispersal behaviours and spatial synchrony

One important factor that affects metacommunity stability is the spatial synchrony of the population dynamics of the various species. Since a species is more vulnerable to extinction when all its population densities in a metacommunity are low at the same time, spatial synchrony decreases the probability of persistence of that species at the regional scale [START_REF] Briggs | Stabilizing effects in spatial parasitoid-host and predator-prey models: a review[END_REF]. Spatial synchrony in predator-prey metacommunities is a balance between synchronising and desynchronising processes. Processes that desynchronise dynamics include different forms of stochasticity [START_REF] Taylor | Metapopulations, Dispersal, and Predator-Prey Dynamics: An Overview[END_REF]), but also nonlinear trophic interactions [START_REF] Bjornstad | Cycles and synchrony: two historical 'experiments' and one experience[END_REF][START_REF] De Roos | Mobility versus density-limited predator-prey dynamics on different spatial scales[END_REF]Jansen 1999). Local dynamics may be synchronised by environmental forcing (the "Moran" effect), high dispersal rates (Bjornstad et al. 1999, Vasseur and[START_REF] Vasseur | Phase-locking and environmental fluctuations generate synchrony in a predator-prey community[END_REF] and nonlinear trophic interactions when dispersal rates are low, [START_REF] Bjornstad | Cycles and synchrony: two historical 'experiments' and one experience[END_REF][START_REF] De Roos | Mobility versus density-limited predator-prey dynamics on different spatial scales[END_REF]Jansen 1999).

Our results show that density-dependent dispersal behaviours can affect spatial synchrony in metacommunities at low and intermediate dispersal rates. In particular, interspecific density-dependent dispersal can strongly desynchronise dynamics. In a patchy migration followed the gradient of prey density. However, their study was focused on equal dispersal rates in the prey and the predator. Here, we show that this effect depends on the relative dispersal of the two species: interspecific density-dependent dispersal in the prey has a strong impact when prey dispersal is higher than predator dispersal, whereas interspecific density-dependent in the predator has a strong impact when predator dispersal is higher than prey dispersal.

According [START_REF] Briggs | Stabilizing effects in spatial parasitoid-host and predator-prey models: a review[END_REF], spatial asynchrony in predator-prey metacommunities is increased when local population dynamics are uncoupled from immigration, i.e. when a local population receive more individuals when its density is low than when its density is high. Because prey and predator densities do not fluctuate synchronously within patches (predator abundance peaks after prey abundance), the dependence of a species' dispersal on the density of another species (interspecific densitydependent dispersal) could promote their uncoupling and increase spatial asynchrony.

We fund a stabilising effect of interspecific density-dependent dispersal in a patchy

Rosenzweig-MacArthur model with identical patches in agreement with Li et al.'s results (2005) obtained in a similar model. By contrast, in patchy models with consistent differences between patches in parameter values and local dynamics that obey either the Lotka-Volterra model (Murdoch et al. 1992) or a variant with a type-II functional response (Ives 1992), the aggregation of the predator in patches of high prey density can be either synchronising or desynchronising. The desynchronising effect could come either from the absence of density-dependence in the prey growth rate (e.g. logistic growth), which is stabilising, or from the differences between patches.

Our results show also that the synchrony of spatial dynamics is affected by intraspecific density-dependent dispersal in the prey. This dispersal behaviour can fully synchronise population dynamics at low dispersal rates, but can also desynchronise local dynamics at intermediate dispersal rates when the shape of density-dependence in dispersal is sigmoidal. Synchronising and desynchronising effects have been observed also in singlespecies metapopulation models [START_REF] Ylikarjula | Effects of Patch Number and Dispersal Patterns on Population Dynamics and Synchrony[END_REF]. In our metacommunity model, we observed this effect for the prey, but not for the predator. We show that the effects of intraspecific density-dependent dispersal in the prey depend on relative dispersal: these effects are strong when prey dispersal is higher than predator dispersal. Thus, densitydependent dispersal and relative dispersal have interacting effects and play an important part in regulating spatial synchrony in metacommunities.

The interaction between density-dependent dispersal and relative dispersal suggests that density-dependent dispersal behaviours can have different effects on stability in different types of predator-prey systems. Predators often disperse more than their prey [START_REF] Rooney | A landscape theory for food web architecture[END_REF]. In such systems, our results suggest that interspecific densitydependent dispersal in the predator can modify stability whereas density-dependent dispersal in the prey does not. But prey sometimes disperse more than their predators [START_REF] Rooney | A landscape theory for food web architecture[END_REF]. In such situations, density-dependent dispersal in the prey can have a strong effect on the stability of the metacommunity: interspecific density-dependent dispersal can increase spatial synchrony, whereas intraspecific density-dependent dispersal can either decrease or increase it depending on dispersal rates.

Dispersal behaviours and top-down control

The extent to which populations are regulated by predation ("top-down" control) or by the availability of resources ("bottom-up" control) has been intensely debated in ecology (e.g. [START_REF] Hairston | Community structure, population control, and competition[END_REF][START_REF] Murdoch | Community Structure, Population Control, and Competition"-A Critique[END_REF][START_REF] Oksanen | Exploitation Ecosystems in Gradients of Primary Productivity[END_REF]. Most ecologists now agree that both types of control are present in most ecosystems [START_REF] Shurin | A cross-ecosystem comparison of the strength of trophic cascades[END_REF][START_REF] Borer | What determines the strength of a trophic cascade?[END_REF] and that their relative strength depends on local factors such as functional diversity [START_REF] Hulot | Functional diversity governs ecosystem response to nutrient enrichment[END_REF] and coevolution between prey and predators [START_REF] Loeuille | Nutrient enrichment and food chains: can evolution buffer top-down control?[END_REF].

Spatial processes can also modify the strength of these controls. A recent model predicts that nutrient inputs and immigration should generally increase the strength of trophic cascades in open ecosystems [START_REF] Leroux | Subsidy hypothesis and strength of trophic cascades across ecosystems[END_REF].

Moreover, predator-prey metacommunity models reveal that dispersal between local communities can also modify the strength of top-down control. Stochastic predator-prey metacommunity models have shown that the relative abundances of the prey and the predator vary with dispersal rate: the mean of prey abundance decreases and the mean of predator abundance increases at intermediate dispersal rates (e.g. [START_REF] Reeve | Environmental variability, migration, and persistence in host-parasitoid systems[END_REF][START_REF] Zeigler | Persistence and patchiness of predator-prey systems induced by discrete event population exchange mechanisms[END_REF]. We fund this result in our deterministic predator-prey model. It is based on a type-II functional response, which implies that the predator controls the prey at equilibrium. When prey and predator local densities fluctuate over time, we fund that the temporal mean of prey density decreases while the temporal mean of predator density increases at intermediate dispersal rates. These variations are observed both at the regional scale (fig. 34) and at the local scale (results not shown). Thus, results of our deterministic model and of stochastic models suggest that intermediate dispersal rates increase the average strength of the control of the prey by the predator in predator-prey metacommunities. Moreover our results show that intra-and interspecific density-dependent dispersal in the prey.reinforces the increase in the strength of predator control on the prey at intermediate dispersal rates.

The effect of dispersal on the strength of top-down control has some experimental support. In aquatic predator-prey metacommunities, the mean density of the predator (Didinium nasutum) was highest at intermediate dispersal rates whereas the mean density of the prey (Colpidium striatum) was decreased by dispersal (Holyoak and Lawler 1996b). In trophically more complex metacommunities, dispersal between mesocosms increased predator biomass, decreased herbivore biomass and increased primary producer biomass [START_REF] Chase | Habitat isolation moderates the strength of top-down control in experimental pond food webs[END_REF]. This suggests that the trophic cascade was stronger in mesocosms connected by dispersal than in isolated mesocosms. Thus, theory and experiments concur to suggest that dispersal and other spatial flows across ecosystems [START_REF] Leroux | Subsidy hypothesis and strength of trophic cascades across ecosystems[END_REF] generally tend to increase the strength of top-down control food-web metacommunities.

Dispersal behaviours and the stability of metacommunities

Our results show that the synchrony of the prey and the predator are strongly related in metacommunities. Indeed, spatial synchrony is affected in the same way in the prey and in the predator by density-dependent dispersal and relative dispersal. In particular, densitydependent dispersal in a species can modify not only its own spatial synchrony, as in metapopulation models, but also the spatial synchrony of the other species. For instance, interspecific density-dependent dispersal in the prey can increase both its spatial asynchrony and that of the predator.

Metacommunity stability depends not only on spatial asynchrony but also on the strength of prey control by the predator: as discussed above, density-dependent dispersal increases predator control on the prey at intermediate dispersal rates. This generates a decrease in regional variability of the predator, and an increase in that of the prey. Regional persistence shows the same trends: the regional persistence of the predator is increased while that of the prey is decreased, although the latter effect is smaller than for regional variability (See supplementary material, fig. S3 andS4). Thus, the stability of the prey and that of the predator display opposite responses to density-dependent dispersal because of the increase in the average strength of predator control on the prey.

Although the regional stability (as measured by both persistence and variability) of the predator is maximal at intermediate dispersal rates, it is nevertheless lower than that of the prey. Thus, top-down control stabilises the less stable species (the predator) and destabilises the more stable species (the prey). This suggests that the metacommunity as a whole may be stabilised at intermediate dispersal. This conclusion is consistent with the general finding of models including different forms of stochasticity that the persistence of predator-prey metacommunities is enhanced at intermediate dispersal rates [START_REF] Crowley | Dispersal and the stability of predator-prey interactions[END_REF][START_REF] Reeve | Environmental variability, migration, and persistence in host-parasitoid systems[END_REF][START_REF] Zeigler | Persistence and patchiness of predator-prey systems induced by discrete event population exchange mechanisms[END_REF]. The mechanism that is generally thought to generate the humpshaped relationship between the stability of predator-prey metacommunities and dispersal is the balance between the positive effect of low dispersal (allowing recolonization and rescue effects) and the negative effect of high dispersal (synchronisation). Here we suggest that another mechanism could stabilise predator-prey metacommunities at intermediate dispersal rates: the modification of the strength of top-down control, which increases the stability of the most unstable species (the predator).

Our numerical results show a hump-shaped relationship between predator stability and dispersal for both constant and density-dependent dispersal. But they suggest that the hump can become more or less marked when dispersal is density dependent. For instance, the hump is stronger when the prey has intraspecific density-dependent dispersal because the latter destabilises the metacommunity at low dispersal rates (spatial synchronisation) and increases the stability at intermediate dispersal rates (strengthening of predator control on the prey). Although interspecific density-dependent dispersal in the prey stabilises the metacommunity at low dispersal (desynchronisation of local dynamics), the increase in the strength of predator control is sufficient to maintain the hump-shaped relationship between predator stability and dispersal.

Thus, density-dependent dispersal can modify the stability of predator-prey metacommunities by modifying both the spatial synchrony of local populations and the strength of prey control by the predator. We hope that these results will stimulate further work on the effects of dispersal on top-down control and the stability of food-web metacommunities.

Ecosystem fragmentation and conservation management

Ecosystem fragmentation leads to decreased species dispersal between patches.

Therefore studying the relationship between dispersal and metacommunity stability is of particular importance to understand the response of metacommunities to fragmentation. We showed that this response may depend on the existence of density-dependent dispersal in species and on their relative dispersal. If fragmentation decreases dispersal of only one species in the metacommunity because of its specific sensitivity to changes in environmental quality between patches, relative dispersal will be modified. In the case of constant prey dispersal and interspecific density-dependent dispersal for the predator, for instance, our model predicts that decreased predator dispersal only may lead to destabilisation of the metacommunity. Thus our work emphasizes the importance of density-dependent dispersal and relative dispersal of species in understanding the response of ecosystems to fragmentation.

This also suggests that enhancing the dispersal of one species through e.g. corridors may not be the best way to prevent its own extinction. For instance, if prey dispersal is constant and predator dispersal is lower and depends on prey density, the increase in prey dispersal by corridors will not increase its stability. In contrast, an increase in predator dispersal may stabilise both the prey and the predator if the ratio between the dispersal rates of the two species is reversed. Thus, our results emphasize the need to integrate conservation policies at the scale of communities and metacommunities, instead of considering them merely at the scale of populations and metapopulations.

CONCLUSION

Our work shows that density-dependent dispersal, which generates interactions between local and regional dynamics, and relative dispersal of the prey and the predator are key to understand the dynamics and stability of metacommunities. Density-dependent dispersal behaviours and relative dispersal deserve more theoretical and experimental explorations to increase our ability to predict the effects of habitat fragmentation on the stability of ecosystems.

APPENDIX A: STABILITY ANALYSIS OF SYMMETRICAL EQUILIBRIUM

Linearization of two-patch model with density-dependent dispersal

We used the method proposed by [START_REF] Jansen | The dynamics of two diffusively coupled, identical Lotka-Volterra patches Theoretical Aspects of Metapopulation Dynamics[END_REF] and generalised it to densitydependent dispersal. System of equations ( 1), and more generally a two-patch model of predator-prey dynamics, can be written in vectorial form into:

°°°°® ) ( ) ( ) ( ) ( ) ( ) ( X G Y G Y F dt dY Y G X G X F dt dX (2) where 1 • © § 1 1 P N X and 1 • © § 2 2 P N Y
are the vectors of population densities of the prey (N) and of the predator (P) in patches 1 and 2. The local dynamics of the prey and the predator are described by the functions

1 • © § ) ( ) ( ) ( X f X f X F P N and 1 • © § ) ( ) ( ) ( Y f Y f Y F P N
. The dispersal dynamics of both species in patch 1 and 2 are described respectively by the

functions 1 • © § ) ( ) ( ) ( X g X g X G P N and 1 • © § ) ( ) ( ) ( Y g Y g Y G P N .
The symmetric equilibrium of (2), noted (X*, Y*) is those for which X(t)=Y(t) for all t. In order to facilitate the analysis of symmetric equilibrium, we introduce

1 • © § 6 6 6 P N Y X ) ( 2 1
, the vector of average densities of each species over the two patches, and

1 • © § ' ' 6 6 ' P N Y X Y X ) ( 2 1
, the vector of differences between densities in a patch and average densities in the metacommunity. System (2) can then be recast into:

> @ > @ °°°°® ' 6 ' 6 ' 6 ' 6 ' ' 6 ' 6 6 ) ( ) ( ) ( ) ( 2 1 ) ( ) ( 2 1 G G F F dt d F F dt d (3)
We are interested in the equilibrium (Σ*, Δ*) of (3) and more particularly to the symmetric equilibria (Σ*, Δ*=0). Note that for Δ = 0, the dynamics of the mean densities, Σ, is exactly the same as the local dynamics of densities in the absence of migration. In particular, the equilibrium values of Σ are the equilibrium values of the local dynamics determined by F. In order to study the stability near a symmetric equilibrium, the system ( 3) , where F and G correspond to our specific model.

is linearized around Σ = Σ*, Δ = 0: > @ °°°°® ' 6 6 ' u 6 6 ' ' 6 6 6 6 u 6 6 ) 2 . 4 ( ) ( 2 ) ( ) 1 . 4 ( ) (
In our specific model, F and G write: ).Thus, the sign structure of the Jacobian matrix of F and G at Σ* are: CV of regional prey density Mean of regional prey density Mean of regional predator density CV of regional predator density

1 • © § 1 • © § P m N c P N b N c P N b K N N r P N F 1 ) , ( and
1 • © § 6 0 * ) ( * J F and 1 • © § 6 ) ( *

  d n (d p = 0.01 d n ) and one in which the converse is true (d p = 100 d n ). Maximal dispersal rates varied between 10 -4 and 10 2 when d n > d p , and between 10 -6 and 10 0 when d p > d n . Within these ranges, we took 1000 sample values of each parameter (keeping the ratio d n /d p constant), uniformly distributed on a log scale. The parameters of local interactions were the following: r=10, K=25, b=11, c=1, m=10. For the density dependence in dispersal, we focused the numerical analysis on sigmoidal shapes (|x αβ |=4, α=n,p ; β=n,p) and on values for the half-saturation parameters S αβ set close to the mean densities observed when populations are isolated: S nn =S pn =18 and S np =S pp =2.

  simulated for each set of parameter values, so that in the case of multiple coexisting attractors, only one of them is reached. But, provided that parameter values are sampled densely enough in the study interval (we used a fine 1000-values sampling grid for the maximal dispersal rates), the other attractors are very likely to be detected for neighboring values of the parameter. This eventually results in a slightly scattered general portrait of the various "continuation branches". The portrait "gaps" are then filled out by linking pairs of detected attractors, provided that parameter values are not too distant and the attractors are close enough in the state space. Attractor closeness was measured by means of the Hausdorff distance in4 .

(

  S αβ ), and their maximal values for both the prey and the predator (d n and d p ). We compared two different shapes, cyrtoid (|x αβ |=1) and sigmoid (|x αβ |=4). We varied the half-saturation parameters S αβ by 20% above and below the mean densities when population are isolated (S nn =S pn =18 and S np =S pp =2). We also varied the maximal dispersal rates of the prey (d n ) and the predator (d p ) between 10 -6 and 10 2 . Within this range, we took a sample of 100 values for each parameter d n and d p (and up to 200 values for the reference parameter combination |x αβ |=4, S nn =S pn =18 and S np =S pp =2), uniformly distributed on a log scale (see supplementary material, fig. S1 and S2).

  The different shapes (x αβ ) and half-saturation values (S αβ ) of density dependence give qualitatively the same results for each type of dispersal, except for intraspecific density-dependent dispersal in the prey, which is sensitive to |x αβ | (supplementary material, fig. S2). Below, we present the results for |x αβ |=4 and mention when they are different for |x αβ |=1.

  d n > d p , intra-and interspecific densitydependent dispersal in the prey have different effects on spatial synchrony (fig. 2B). Intraspecific density-dependent dispersal synchronises local dynamics when dynamics are periodic (at low and intermediate dispersal rates), whereas spatial synchrony is decreased when dynamics are chaotic. The presence of the synchronising effect of intraspecific density-dependent dispersal does not depend on the shape of density dependence in dispersal (|x αβ |=4: fig. 2B; |x αβ |=1: Appendix B, fig. 5A), but the presence of chaotic dynamics does. When |x αβ |=4 chaotic dynamics arise for a wider range of dispersal rates than when dispersal is constant (fig. 2A, 2B) whereas when |x αβ |=1 chaotic dynamics disappear (Appendix B, fig.

  F and G.The stability of (Σ = Σ*, Δ = 0) is determined by the eigenvalues of matrices ) (Σ = Σ*, Δ = 0) of system (4), or the symmetrical equilibrium (X*, Y*) of system (2), is stable if both matrices ) two-patch Rosenzweig -MacArthur modelIn order to determine the effects of density-dependent dispersal on the stability of the symmetric equilibrium of (1), we analyse the signs of the trace and of the determinant of

Figure 1 .

 1 Figure 1. Per capita dispersal rate D α of species α (α = n or α = p) in response to changes in the density of species β (β = n or β = p). S αβ is the half-saturation parameters. Densityindependent (x αβ =0: solid line) dispersal; positive (black dashed line) and negative (grey dashed line) density-dependent dispersal; cyrtoid (|x αβ |=1: thin dashed lines) and sigmoid (|x αβ |=4: thick dashed lines) shape of density dependence in dispersal.

Figure 2 .

 2 Figure 2. Dynamics and spatial synchrony of the prey in predator-prey metacommunities as functions of its maximal dispersal rate (d n ). Top diagrams illustrate the various cases studied: N and P represent prey and predator densities; dispersal of each species is represented by a wide white arrow; positive or negative effects of densities on dispersal are represented by thin solid arrows with a + or a -sign. Figures in the first column (A and D) present results when dispersal is constant in both species. In the second column (B and E), we present results for (intra-and interspecific) density-dependent dispersal in the prey. In the third column (C and F), we present results for (intra-and interspecific) densitydependent dispersal in the predator. We present results for d n = 100 d p (d n > d p ) and for d p = 100 d n (d p > d n ) in the first (A,B and C), and in the second (D, E and F) row, respectively. Markers on curves (except for A and D) indicate what type of density dependence the curve refers to: intraspecific (dots) or interspecific (triangles). The nature of dynamics obtained using the Lyapunov exponent is represented on curves: chaotic dynamics (with colour thickening) or periodic dynamics (without colour thickening). The specific colours used (gray, cyan and red) as well as the position of the thickening (above or under the curve) only help, in addition to curve markers, to distinguish between the different curves and do not have any further meaning. Notice on figures D, E, F (d p > d n ) the two coexisting attractors (synchronous and asynchronous) for some dispersal rate values. Parameters values are the following. Constant dispersal: x αβ =0; positive density-dependent dispersal: x αβ =4(αβ=nn,np,pp); negative density-dependent dispersal: x pn =-4. Other parameters: S nn =18, S np =2, S pp =2, S pn =18; local interactions: r=10,K=25, b=11, c=1, m=10.

Figure 3 .

 3 Figure 3.Temporal mean of regional density (upper sub-graph) and coefficient of variation (CV, lower sub-graph) of the prey in predator-prey metacommunities as functions of its maximal dispersal rate (d n ). See legend of figure 2.

Figure 4 .

 4 Figure 4. Temporal mean of regional density (upper sub-graph) and coefficient of variation (CV, lower sub-graph) of the predator in predator-prey metacommunities as functions of the maximal dispersal rate of the prey (d n ). See legend of figure 2.

Figure 5 .

 5 Figure 5. The effect of intraspecific density-dependent dispersal in the prey when x nn =1. The spatial synchrony (A), the temporal mean of regional densities (upper sub graph) and the coefficient of variation (CV, lower sub graph) of the prey (B) and of the predator (C) as functions of the maximal dispersal rate of the prey (d n ). We present results for d n = 100 d p (d n > d p ). The curve coding is the same as in figure 2. Notice that, compared to figures 2, 3 and 4, intraspecific density-dependent dispersal does not induce chaotic dynamics any more and that two attractors can now coexist at intermediate dispersal rates. Parameters values are the same as in figure 2, except |x nn |=1.

  αβ are the densities of species β (prey: β stands for "n"; or predators: β stands for "p") for which the per capita dispersal rates D D of species α (α = "n", "p") is half of the maximal per capita dispersal rate d α . Parameters x αβ determine the shape of density

	dependence in dispersal. When x αβ = 0, dispersal is independent of density and occurs at a
	constant per capita rate d α /2; when x αβ > 0, the dispersal rates increase with density;
	whereas when x αβ < 0, they decrease with density. The absolute value of x αβ determines the
	steepness of density dependence: when |x αβ | = 1, dispersal rates are concave functions of
	density up to their saturation level (cyrtoid shape), whereas they are sigmoidal when |x αβ | > 1
	(the inflexion point of the dispersal function then occurs at the threshold density S αβ ); high
	values of |x αβ | model step-like sigmoid curves, whereas low values model smoother
	variations in dispersal with population density.
	n D N	i	d	n	S	nn	nn x N	nn x i N	nn x i	or	( ) i D P d n	n	S	np	np x P i	np x	i P	np x
	and the per capita dispersal rate of the predators writes
	( ) i D P d p	p	S	pp	x	pp P i	x	pp	i P	x	pp	or	( ) i D N p	d	p	S	pn	x	pn N	x i	pn N	x i	pn	.

These formulations imply several assumptions. Species have a maximal per capita dispersal rate d n for the prey, and d p for the predator (fig.

1

). These parameters can reflect the ability of the species to disperse or the distance between patches. The half-saturation parameters S

  mean regional density of the predator increases more strongly when d n > d p and there is intraspecific density-dependent dispersal for the prey (only when |x αβ |=4: fig. 3B; fig. 4B; |x αβ |=1: Appendix B, fig. 5B,C) or interspecific density-dependent dispersal for the prey (fig. 3B; fig. 4B), or when d p > d n and there is interspecific density-dependent dispersal for the predator (fig. 3F; fig. 4F). These changes commonly go hand in hand with the emergence of new chaotic regimes (fig. 3F; fig.

  , where D N and D P are the per capita dispersal rate of the prey and the predator. We allow any form of density-dependant dispersal rates, provided that prey dispersal increases

											G	N	,	P	© §	P N	D D P N	, , N N ( (	) ) P P	1 •
	with prey density (	w D N	w N	!	0	) and predator density (	w D N	w P	!	0	) and that predator
	dispersal increases with predator density (	w D P	w P	!	0	) and decreases with prey density
	(	w D P	w N	0

ACKNOWLEDGEMENTS

We thank Vincent A. A. Jansen, Bernard Cazelles, Silvia de Monte, E. Peter Münger for stimulating discussions and Maurice Sabelis, Louis-Félix Bersier and two anonymous reviewers for providing comments and suggestions on the manuscript. Michel Loreau was supported by a Discovery Grant from the Natural Sciences and Engineering Research Council of Canada and by a team research project from the Quebec Nature and Technologies Research Fund. Céline Hauzy was supported by a grant from the Chancellerie des Universités de Paris.


Aristotelis Kittas
email: akittas@gmail.com

Evolution of the rate of biological aging using a phenotype based computational model

Keywords: evolution, aging, computer simulations, age-structured populations, modelling

Evolution of the rate of biological aging using a phenotype based computational model

Introduction

Aging is probably one of the most familiar, yet least understood aspects in biology [START_REF] Kirkwood | Understanding the odd science of aging[END_REF]. It can be defined, from an evolutionary point of view, as a progressive decline in fitness (the ability to survive and reproduce) with increasing age [START_REF] Partridge | Evolutionary theories of ageing applied to long-lived organisms[END_REF]. It has been used interchangably with the term senescence, which is the lowering of survival rates as time goes by. Cellular senescence is the phenomenon by which normal diploid cells lose the ability to divide, while organismal senescence refers to the aging of organisms. Senescence is a complex process, which may derive from a variety of mechanisms, such as the shortening of telomeres with each cell cycle, oxidative stress (i.e. the imbalance between the production of reactive oxygen and ability to readily detoxify the reactive intermediates) and others, and its role in organismal aging is at present an active area of investigation.

Senescence of the organism can give rise to the Gompertz law of mortality. This law, observed by Gompertz (1825), was proved to be valid for the populations of most industrialised countries, describing the age dynamics of human mortality rather accurately in the age window from about 30 to 80 years of age [START_REF] Thatcher | The long-term pattern of adult mortality and the highest attained age[END_REF]. The Gompertz law of exponential increase in mortality rates with age is observed in many biological species [START_REF] Strehler | Time, Cells and Aging[END_REF][START_REF] Finch | Longevity, Senescence, and the Genome, 1st Edition[END_REF], including humans, rats, mice, fruit flies, flour beetles, and human lice [START_REF] Gavrilov | The Biology of Life Span: A Quantitative Approach[END_REF]. Let S(α) be the probability of surviving from birth to age α, which can be obtained by the cumulative distribution of the stable age distribution. Then, the mortality function μ(α) is:

μ(α) = - d(lnS(α)) dα (1)
Because data normally is given in yearly intervals, equation ( 1) can be approximated by:

μ(α) = ln( S(α) S(α + 1) ) ( 2 )
The Gompertz law states that the mortality rate of certain organisms increases exponentially with age α:

μ(α) ∝ c • e bα (3)
This law, however, is not without exceptions. Some animals, including many reptiles and fish, age extremely slowly and exhibit very long life spans (negligible senescence) [START_REF] Finch | Longevity, Senescence, and the Genome, 1st Edition[END_REF]. Some produce even more offspring as they get older, while their mortality rate actually drops, in disagreement with the Gompertz "law", a phenomenon called negative senescence [START_REF] Vaupel | The case for negative senescence[END_REF]. Another notable exception are Hydras, which do not undergo senescence and are therefore biologically immortal (ageless) [START_REF] Martinez | Mortality patterns suggest lack of senescence in hydra[END_REF].

Aging theories

"There is no shortage of theories of aging" [START_REF] Partridge | A singular view of ageing[END_REF]. Almost every aspect of an organism's phenotype undergoes modification with aging, and this phenomenological complexity has led, over the years, to a bewildering proliferation of ideas about specific cellular and molecular causes [START_REF] Kirkwood | Understanding the odd science of aging[END_REF]). An attempt by Medvedev [START_REF] Mevdevev | An attempt at a rational classification of theories of ageing[END_REF] to rationalise the multiplicity of hypotheses resulted in a listing of more than 300 mechanistic "theories" of aging, that describe how aging occurs. However, there are only a few evolutionary theories about why aging evolved.

Evolutionary biology is concerned with the reasons behind the aging process and the challenge of why aging occurs, in spite of its obvious drawbacks [START_REF] Kirkwood | Understanding the odd science of aging[END_REF]. Aging results in significant loss of Darwinian fitness, giving rise to the question why it has not been eliminated by natural selection. It was proposed [START_REF] Hamilton | The moulding of senescence by natural selection[END_REF][START_REF] Templeton | Population Genetics and Microevolutionary Theory, 1st Edition[END_REF] that if mutations can occur that kill their bearers at a sufficiently advanced age, such mutations are effectively neutral and some will go to fixation, thereby destroying the agelessness of the initial population.

The "mutation accumulation" theory [START_REF] Medawar | An Unsolved Problem of Biology[END_REF] claims that aging results from the non-adaptive process of mutation accumulation. If a deleterious mutation manifests itself at a young age, there will be strong selection pressure to eliminate it because it will affect the fitness of a large majority of the population. However, if that same mutation does not manifest itself until later in life, many of the individuals carrying it will have died before it is expressed. Therefore, mutations killing an individual before it reaches reproductive maturity vanish from the population, while those acting later in life allow it to reproduce and are selected with a much weaker force [START_REF] Charlesworth | Aging: Questions and prospects[END_REF].

Williams in 1957 proposed a similar theory, in which the existence of pleiotropic genes of special sort have opposite effects on fitness at different ages. Their effects are beneficial in early life, when natural selection is strong, but harmful at later ages, when selection is weak. If these genes confer increased reproductive success early in life they would be selected for, despite the fact that they may later cause senescence. This is now known as "antagonistic pleiotropy" theory [START_REF] Williams | Pleiotropy, natural selection, and the evolution of senescence[END_REF], which claims that aging results from an adaptive process of selection favoring certain genes with netpositive effects. The p53 gene is perhaps one of the clearest examples of an antagonistially pleiotropic gene [START_REF] Leroi | What evidence is there for the existence of individual genes with antagonistic pleiotropic effects?[END_REF].

Another major theory is the "disposable soma" theory [START_REF] Kirkwood | Evolution of ageing[END_REF][START_REF] Kirkwood | The evolution of ageing and longevity[END_REF][START_REF] Kirkwood | Evolution of senescence: late survival sacrificed for reproduction[END_REF], which focuses on the idea that cell maintenance (e.g. DNA repair, protein turnover etc) is costly. It suggests that longevity is controlled primarily through genes that regulate the levels of somatic maintenance and repair [START_REF] Kirkwood | Understanding the odd science of aging[END_REF]. The organism should optimally allocate its metabolic resources, chiefly energy, between the maintenance and repair of its soma and its other functions in order to maximise its Darwinian fitness. The necessity for trade-off arises because resources allocated to one function are unavailable for another. Therefore, the allocation of resources to maintenance and repair is determined by evolutionary optimization.

Mutation accumulation and antagonistic pleiotropy provide the backbone for much of the thinking about the evolutionary genetics of aging. However, evidence for mutation accumulation remains limited and controversial [START_REF] Kirkwood | Understanding the odd science of aging[END_REF][START_REF] Shaw | Toward reconciling inferences concerning genetic variation in drosophila melanogaster[END_REF], while evidence for antagonistic pleiotropy, although stronger, is as yet lacking in detail. There have been suggestions [START_REF] Kirkwood | Understanding the odd science of aging[END_REF]) that these theories may all provide a partial explanation for the phenomenon of aging but no single theory has been explicitly proven to explain the aging process.

Computational modelling for biological aging

While there is an abundance of aging theories, the same cannot be said about computer simulations of biological aging [START_REF] Stauffer | The penna model of biological aging[END_REF]. T. J. P. Penna implemented the mutation accumulation theory [START_REF] Penna | A bit string model for biological aging[END_REF] by dividing and individual lifetime into B time intervals, representing the genome with a string of B bits (chronological genome). An analytic solution is also available [START_REF] Coe | Analytical solution of a generalized penna model[END_REF]. The model has been successfull in reproducing the senescence of the pacific salmon [START_REF] Penna | Mutation accumulation and the catastrophic senescence of the pacific salmon[END_REF], and has been modified to account for the antagonistic pleiotropy mechanism [START_REF] Sousa | An unusual antagonistic pleiotropy in the penna model for biological ageing[END_REF][START_REF] Cebrat | Altruism and antagonistic pleiotropy in penna ageing model[END_REF] and spatial distributions on a lattice [START_REF] Makowiec | Penna model of biological aging on a lattice[END_REF]. Speciation in age structured populations has also been of interest, and in specific sympatric speciation [START_REF] Sousa | Sympatric speciation in an age-structured population living on a lattice[END_REF][START_REF] Luz-Burgoa | Computer simulation of sympatric speciation with penna ageing model[END_REF] and parapatric speciation [START_REF] Schwammle | Monte carlo simulations of parapatric speciation[END_REF], in some cases including sexual reproduction, diploid and even triploid organisms [START_REF] De Oliveira | Evolution, ageing and speciation: Monte carlo simulations of biological systems[END_REF][START_REF] Sousa | Evolutionary advantage of diploidal over polyploidal sexual reproduction[END_REF].

A rather obscure model for the evolution of aging is the Heumann-Hotzel model [START_REF] Heumann | Generalization of an aging model[END_REF], which is a generalization of the Dasgupta model [START_REF] Dasgupta | A computer simulation for biological ageing[END_REF]. In this model, each individual carries a "chronological genome" of size x max , with a survival probability per time step G(x) at age x, which can be modified according to some mutations. The model has been modified to account for its weaknesses, namely its incapacity to treat populations with many age intervals and its ability to handle mutations exclusively deleterious (de [START_REF] De Medeiros | The heumann-hotzel model for aging revisited[END_REF], where the catastrophic sensescence and the Gompertz law have been demonstrated. Optimization models of the disposable soma theory have been developed by several groups [START_REF] Cichon | Ageing and typical survivorship curves result from optimal resource allocation[END_REF][START_REF] Mangel | A life-history perspective on short and long-term consequences of compensatory growth[END_REF][START_REF] Abrams | Optimality theory, gompertz' law, and the disposable soma theory of senescence[END_REF][START_REF] Vaupel | The case for negative senescence[END_REF][START_REF] Baudisch | Inevitable aging? Contributions to evolutionarydemographic theory[END_REF][START_REF] Drenos | Modelling the disposable soma theory of ageing[END_REF], which describe how the optimal investment in maintenance is affected by varying the parameters that specify the schedules of reproduction and mortality.

The model

A computational model was constructed to study the evolution of the rate of biological aging. This is a phenotype (in contrast to genotype) centered model and does not focus on a specific genetic procedure that affects the aging process. It does not assume any genetic mechanism of aging at work, just the idea that each individual's ability to survive is age dependent. It focuses on the phenotypic results of the genome i.e. the viability of each individual. The model uses asexual reproduction and assumes a well mixed population (i.e. no spatial effects) and no interaction between the individuals.

The motivation behind this work is to provide a simple model to understand how various parameters, such as the age at maturity or different kinds of mutations affect the evolution of the rate of aging. In this model I assume that the viability of each individual is age dependent and diminishes with age. Given this hypothesis, mutations can produce offspring with increased or decreased viability. Consequently, natural selection acts on these mutations, affecting the overall viability of the population. The model can make predictions about how the mutation probabilities, mutation strength, and age at reproductive maturity shape aging patterns and affect the aging rate and the final state of the system.

I focus in the simplicity of the model rather than how it would be more realistic. As a result, a small number of parameters have been used and a number of simplifying assumptions have been made. It assumes asexual reproduction only, a fixed age at reproductive maturity (which is known to differ from the strategy many organisms use in the face of limited resources), a constant number of offspring and a very specific function which is used to calculate each individual's viability. Consequently the model has its limitations and favors simplicity over realism. The intention behind it is to give a qualitative understanding and some initial predictions on the influence of each individual's viability and how its change affects the aging process and thus and instigate future work on the subject.

Let N 0 be the initial number of individuals. Each individual is characterized by a parameter f that is responsible for the individual's viability F (α), which is age dependent and stays in range from 0 to 1. Parameter f is the measure of aging and the change of its value over time shows the evolution of biological aging. The value of f indicates how an individual's viability changes over its life course; high values of f imply fast aging and low f slow aging (see also Eq. ( 4)). Parameter f can be related to Darwinian fitness under the assumption that individuals with higher survival potential will be able to maintain or increase their numbers in succeeding generations, leading to a better representation of their genome.

The threshold parameter T , which controls the number of hits an individual can withstand is typically set to 2. Simulations have been also been performed with T = 1, T = 4, T = 8, and T = 16. The Gompertz pattern arises typically for small T values (e.g. T ≤ 4), because for higher values stochastic effects are more pronounced and the mortality pattern is more scattered due to noise, deviating significantly from the exponential function. Therefore, a small constant value of T is used to minimize this effect and to reduce the complexity in the study of the model's behavior.

An aging function, which is a logistic function, is used to calculate the individual's viability at age α in the form of:

F (α) = 1 1 + ( 1 B -1) • e fα (4)
where: B is a constant (typically equal to 0.99), f is each individual's characteristic parameter for its viability, which can be modified due to mutations, and F (α) each individual's viability at age α. A time step is completed when all the individuals at the specific time are selected. The simulation begins with N 0 individuals. The algorithm at each time step is as follows:

1. An individual is selected. 2. Its age is increased by 1.

3. F (α) is calculated for the individual's current age α. It then receives a "hit" with probability 1 -F (α). Let H be the number of total hits an individual can withstand. 4. T is the parameter which represents the threshold of hits. If H = T (T typically equal to 2), the individual dies and the next individual is selected. 5. If the individual survives and its age α is α ≥ R (where R is the age at reproductive maturity), it generates 1 offspring. 6. The offspring survives with probability: 1 -

N (t)
Nmax . This is the Verhulst Factor, which accounts for the limited carrying capacity of the environment. 7. If the offspring survives, they receive their parents' genome with mutations that affect their viability (parameter f ). These mutations can have a beneficial effect (f is reduced with probability P b ), dentrimental effect (f is increased with probability P d ), or neutral (f remains constant).

As a result of the aging function (Eq. ( 4)), low values of f correspond to high viability of each individual (and vice versa). The amount that the f is modified in the newborn is parameter M of the simulation (mutation step). To keep the exponent in Eq. ( 4) above zero, the minimum value of f is restricted to be equal to M in each simulation. In this model, the maximum age of each individual does not depend on an external parameter, but emerges naturally from the aging function.

The Verhulst factor is applied here in a different fashion than the Penna or the Heumann-Hotzel model. Instead of acting as a random death probability for all the population, it acts here only on the individuals whose genomes have not been tested by the environment yet, the newborn. From the economics of the population, this is clearly a better choice, since little investment is wasted [START_REF] Lee | Rethinking the evolutionary theory of aging: Transfers, not births, shape senescence in social species[END_REF]. From the biological perspective, although it is not yet the most faithful representation of the real natural processes, it has an advantage since the genome of the newborn is on the average less well-fitted, because of the overwhelming majority of bad over good mutations and random deaths will only occur for a fraction of the population that has more bad mutations than the average [START_REF] Martins | Random deaths in a computational model for age-structured populations[END_REF][START_REF] Niewczas | The influence of the medical care on the human life expectancy in 20th century and the penna ageing model[END_REF].

The quantities I monitor in the simulation are the total population N(t), the mortality function μ(α) when the population reaches a stable age distribution, and the average value f (t) of the population over time, which is a measure of the evolution of aging in the population. Typical values used for the population dynamics in our model are N 0 =100 and N max = 10 4 . All results are an average of 100 independent runs in each case.

Results and Discussion

Firstly, we would like to observe the population dynamics of the model with no mutations present, i.e the parameter f of each lineage remains constant during the simulation. In Fig. 1, N(t) is shown with no mutations for different values of reproduction age R. It is clear that populations with low R become fixed at higher values than the ones where R is high (see also inset of fig. 3). This can be attributed to the fact that as R increases, the average number of individuals that will be given the opportunity to reproduce is decreased. For R = 8 the limiting factor is the carrying capacity of the environment (Verhulst factor) and N(t) quickly approaches N max and settles a little below it, while for R = 30 the reproduction age is principal factor for limiting population growth.

The time it takes for the population to stabilize is also higher for high values of R. Because not everyone has the chance to reproduce, N(t) increases much slower compared to lower values of R. The fluctuations in the population density become more pronounced for high R. For R = 30, when a new generation is born, it will take 30 time intervals for its individuals to reproduce. The effect is obviously less pronounced for lower values of R and disappears when the population density becomes stable. If R is sufficiently high, the population will perish, because very few individuals will reach the age of reproduction and the rate of deaths will surpass the birth rate.

Fig. 2 shows the mortality function μ(α) for different initial values of f 0 , for a stable age distribution. It is clear that the main factor that affects the mortality function is the viability of the individuals, with lower values of f allowing them to survive longer. The model is able to reproduce the Gompertz-law of mortality (eq. 3), and μ(α) is shown to depend exponentially on age a, with the exponent having a linear dependence on the value on f 0 (fig. 2, inset). Simulations were also performed for different values of R, where it is evident that it does not have a dependence on the mortality function (although it can affect the stable state of the system when mutations are present, as is later shown).

Population density N(t) in the presence of mutations is shown in fig. 3. For R = 30, N(t) increases significantly slower than R = 8, just as in the case without mutations. With mutations however, it takes more time to reach a stable population density. There is a significant difference in this case; the population will reach N max regardless of R in contrast to the case with no mutations, where the final population density is greatly affected by R (see inset of fig. 3). While the population has a tendency to die out due to mutations, it is obvious that it stabilizes at a value close to the carrying capacity of the environment. This happens since beneficial mutations increase survival chances, and so even for large R the population evolves to survive to maturity and thus eventually is able to reach its size limit.

Mutations change the viability of the individuals and give the opportunity to the system to reach a different stable state. In fig. 4 I plot f (t) of the population during the course of the simulation. A high value for deleterious mutations is used; 90% of the newborn will not survive as long as their parents (same values of simulation parameters as in fig. 3). A high probability of detrimental mutations tends to reduce significantly the overall viability of the population. This is evident by the initial increase of f (t) and the peak in fig. 4. However, natural selection is at work and the small percentage of individuals with increased viability will reproduce more and carry their healthy genes to their offspring.

Consequently, in the system there are two forces at work: mutations which tend to decrease the viabilty of the population (increase f ) and natural selection, which allows individuals with low f to generate more offspring and so contributes to a decrease of f , namely an increase in the average viability of the population. Therefore, mutations combined with low selection pressure start to move the system towards a low viability with weaker phenotypes and after the time delay selection pressure is tightened and only the small percentage of the healthier phenotypes survive and reproduce, thus increasing the viability of the population. The constant reproductive potential, which the model assumes, favors simplicity over realism. Real organisms produce different number of offspring and this can affect greatly the dynamics of the system.

The effect of R on the final state of the system is evident in fig. 4. High values of R will have as a result a stable state with low f . This happens because for low R individuals with reduced viability will still have a high probability to reproduce, because they reach the reproduction age early in their lives. In this case the selection force on them is very relaxed and f tends to increase and stabilize at higher values. Consequently, the force of natural selection is tightened for higher reproduction ages, e.g. in the case of R=30. Individuals with reduced viability are strongly selected against, because they will not have a chance to reproduce. As a result, the value of f will not have a high tendency to increase compared to lower R values. This explains the lower peak in the case of R = 30 compared to R = 8 in fig. 4 and the lower value of f where the system reaches equilibrium. If, however, the reproduction age is very high the individuals will not have a chance to reproduce, f increases and does not stabilize and the population dies out, as shown below.

In fig. 5-9 I relax the probability of mutations, allowing a large room for neutral mutations, which do not change the viability of individuals. The population density still reaches the environent's carrying capacity regardless of R, as long as the R remains within the limits that allow the individuals to reach reproductive maturity (fig. 5). The mortality function obeys the Gompertz law and does not depend on R (inset of fig. 5). Detrimental mutations still happen much more frequently than benecifial ones. As a direct result f decreases slower and the selection pressure is less evident for different R (compare fig. 4 to fig. 6).

If P d is increased even more (fig. 78) an interesting thing happens, namely the population dies out for low R values. In this case, the relaxed selection pressure at low R can be fatal to the population as f keeps increasing (fig. 8) at a rapid rate at R = 8. Here, detrimental mutations decrease the population's viability because they are not selected against and the population perishes. For high R values, however, f decreases and the system reaches equilibrium, because the strong selection pressure eliminates deleterious mutations and weaker phenotypes. There seems to be a critical value close to R = 16, below which the population dies out and above which it survives and reaches a stable state. Therefore, the system's final state is shown to be very sensitive to the age at reproductive maturity. It should be noted that these dynamics assume a constant age at reproductive maturity in order to keep the model simple, while the age at reproductive maturity changes dynamically for many organisms (e.g. given the strategy that they use in the face of limited resources), and this change could have a great effect on the evolution of the system.

The probability of beneficial or detrimental mutations is therefore shown to affect how the population evolves, as well as R, which plays a critical role in the final state of system. For lower values in the ratio of the probability of detrimental/beneficial mutations (P d /P b ), aging evolves almost as fast regardless of R (fig. 6). Higher values favor significantly species with late age of maturity, in which case there seems to be a critical value in the age at reproductive maturity, below which the population dies out and above which it survives and evolves to a stable state (fig. 8). High absolute values of these probabilities produce peaks in viability of the population f (t), because strong selection pressure takes over after a significant number of weaker phenotypes have appeared due to the high mutation rate (fig. 4). In species with high age at reproductive maturity, natural selection is stronger and quickly eliminates the weaker phenotypes, while in species with low R it is more relaxed and the peaks in f (t) are more pronounced.

The role of mutation step M is examined in fig. 9, which shows f (t) for various M values. The system is again shown to evolve towards a low f of its individuals. However in this case, it will take considerably more time ( f (t) decreases very slow) to reach the final state for M = 0.001, with M controlling the speed with which evolution takes place. The process, however is also controlled by P d and P b . If the mutations are relaxed, the system will evolve to a steady state of individuals with higher potential for survival. If we use P d = 0.9 and P b = 0.1 (same values as fig. 34), the population will die out for M = 0.001. This happens because the small difference in fitness in the newborn is not significant enough to be selected by natural selection and compete with the high tendency of mutations to reduce the population's viability. The carrying capacity of the environment however, is also to be taken into account in conjunction with the value of M. Simulations have been performed with larger populations (N max = 100000), where the population reaches a stable state because a larger number for individuals allow natural selection more "room" to act.

The effect of different mutation probabilities P d and P b is studied in fig. 10. The frequency of beneficial and deleterious mutations and their ratio is again revealed to play a crucial role in the system, controlling not only the evolution rate, but also whether the population survives or perishes. When the probability of beneficial mutations is sufficiently small the population will die out, otherwise the system will evolve to a stable state with better overall viability (see also inset of fig. 10, where the effect is apparent in the final population density). The higher the absolute values of these probabilities, the more pronounced the peaks are for f ( t). Not only that, but lower values allow the system to stabilize at lower f ( t), namely evolving towards a more viable population (compare P d = 0.9, P b = 0.1 with P d = 0.5, P b = 0.05 and P d = 0.1, P b = 0.01 in fig. 10). High relative values of P d /P b slow down evolution, while for even higher values the overall viability decreases rapidly ( f ( t) increases) and the population dies out (compare P d = 0.9, P b = 0.1 with P d = 0.9, P b = 0.05 and P d = 0.9, P b = 0.01 in fig. 10).

Conclusions

A simple model, based on physical viability, was constructed to study the evolution of biological aging. The mortality rates are shown to obey the Gompertz law. It's clear that given sufficient "room", i.e. number of individuals and amount their viability is modified, natural selection drives the system towards better phenotypes, reducing aging effects. Mutations appearing at later ages may be acted weakly upon by natural selection as suggested by the mutation accumulation theory, but beneficial mutations which extend the survival of the species will always have a tendency to increase in the population, as shown by the model.

Life exists far from thermodynamic equilibrium. Maintaining its stability could be regulated by optimally allocating the organism resources, however this does not exclude the existence of an ageless phenotype in an optimal environment. Entropy must inevitably increase within a closed system, but living beings are not closed systems. A single theory that fully decribes all aging phenomena and why aging has evolved, is still lacking. This model is not concerned with why aging evolved, but can make predictions about how various parameters, such as the age at reproductive maturity or mutation probabilities, affect the aging patterns.

More specific, in the absence of mutations lower ages at reproductive maturity R tend to reach higher values of final population density. With mutations present, the population stabilizes near the maximum density, except for species with short pre-reproductive periods, which can only tolerate low rates of deleterious mutations and thus the populations dies out. For lower values in the ratio of the probability of detrimental/beneficial mutations (P d /P b ), aging evolves almost as fast regardless of R. Higher values favor significantly species with late age of maturity, in which case there seems to be a critical value in the age at reproductive maturity, below which the population dies out and above which it survives and evolves to a stable state.

High absolute values of these probabilities produce peaks in f (t) and low values allow the system to stabilize at lower f (t). In species with low R selection pressure is more relaxed and the peaks are more pronounced. High relative values of P d /P b decrease the evolution rate, and if this ratio is further increased, f (t) increases rapidly and the population dies out. The rate of aging is controlled by the mutation probabilities in conjunction with the amount they modify each individual's phenotype, with neutral mutations allowing the system more "room" to evolve. Therefore, even in this simple model, the process of evolution of biological aging is shown to be fairly complex and very sensitive to the parameters that control it.

Figure 1 :

 1 Figure 1: N (t) vs t for different values of reproduction age R. t = 10 4 , f 0 = 0.1, no mutations.

Figure 2 :

 2 Figure 2: μ(α) vs age α for the population at t = 10 4 for various values of initial fitness f 0 . R = 8, no mutations. Inset: Slope A of μ(α) vs f 0 .

Figure 3 :

 3 Figure 3: N (t) vs t for different values of reproduction age R. t = 10 4 , f 0 = 0.1, Mutations are allowed, P b = 0.1, P d = 0.9, M = 0.01. Inset: Final population N for the stable state vs R and comparison with the case with same simulation parameters, but no mutations

Figure 4 :

 4 Figure 4: f (t) vs t for different values of reproduction age R. t = 10 4 , f 0 = 0.1, mutations are allowed, P b = 0.1, P d = 0.9, M = 0.01. Inset: Final f for the stable state vs R.

Figure 5 :

 5 Figure 5: N (t) vs t for different values of reproduction age R. t = 10 4 , f 0 = 0.1, mutations are allowed, P b = 0.01, P d = 0.1, M = 0.005. Inset: Mortality function μ(α) vs α for t = 10 4 (R = 8 only shown for clarity).

Figure 6 :

 6 Figure 6: f (t) vs t for different values of reproduction age R. t = 10 4 , f 0 = 0.1, mutations are allowed, P b = 0.01, P d = 0.1, M = 0.01.

Figure 7 :

 7 Figure 7: N (t) vs t for different values of reproduction age R. t = 10 4 , f 0 = 0.1, mutations are allowed, P b = 0.01, P d = 0.5, M = 0.005.

Figure 8 :

 8 Figure 8: f (t) vs t for different values of reproduction age R. t = 10 4 , 0 = 0.1, mutations are allowed, P b = 0.01, P d = 0.5, M = 0.005.

Figure 9 :

 9 Figure 9: f (t) vs t for different values of mutation step M . t = 10 5 , f 0 = 0.1, mutations are allowed, P b = 0.01, P d = 0.1, R = 16. Inset: Population N (t) vs t (M = 0.001 only shown for clarity).

Figure 10 :

 10 Figure 10: f (t) vs t for different values of mutation probabilities P b and P d . t = 10 4 , f 0 = 0.1, M = 0.01, R = 16. Inset: Population N (t) vs t.