François Massol
email: francois.massol@cefe.cnrs.fr

Vincent Calcagno
email: calcagno@isem.univ-montp2.fr

Julien Massol
email: julien.massol@ens-lyon.fr

The metapopulation fitness criterion: proof and perspectives

Keywords: dispersal, fitness, invasibility, metapopulation, persistence criterion, viability

Metapopulation theory has recently been stirred by the development of a metapopulation persistence criterion (R m ) quantifying the "lifetime dispersal success" of a newly colonized deme in a sparsely occupied metapopulation. No rigorous proof of this criterion in continuous time was available so far. Here, we show that this criterion can be mathematically justified from standard Jacobian arguments. A summary of the key elements of this proof emphasizes the assumptions under which this criterion is valid. Examples illustrate how to generally compute metapopulation fitness in continuous time models. The underlying assumptions behind the criterion are

Introduction

A topic of primary interest in population biology is the possibility of persistence [START_REF] May | Unanswered questions in ecology[END_REF]. A biological system is said persistent when all its components are protected from deterministic extinction. This means that the boundaries at which one component gets extinct are repellent, or equivalently that each component of the system has positive growth rate when rare (for more details on definitions of persistence, see [START_REF] Jansen | Shaken not stirred: On permanence in ecological communities[END_REF].

This "protection" criterion is used extensively in population genetics, population dynamics and ecology [START_REF] Roughgarden | Theory of Population Genetics and Evolutionary Ecology: An Introduction[END_REF][START_REF] Charlesworth | Evolution in Age-Structured Populations, 2nd Edition[END_REF][START_REF] Caswell | Matrix Population Models: Construction, Analysis, and Interpretation[END_REF][START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF].

For single-component systems, such as a single-species population, it is usually called a viability criterion. The viability of a population is commonly assessed by computing the basic reproduction ratio (R 0 ) of the species when it is rare [START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio r 0 in models for infectious diseases in heterogeneous populations[END_REF][START_REF] Metz | How should we define fitness for general ecological scenarios[END_REF][START_REF] Mylius | On evolutionarily stable life histories, optimization and the need to be specific about density dependence[END_REF][START_REF] Caswell | Matrix Population Models: Construction, Analysis, and Interpretation[END_REF]. In physiologically structured (e.g. age-structured) populations, this computation requires knowledge of the structure of the population (e.g. the age pyramid) and its fecundity function (the expected rate of offspring sired by one adult in each category). R 0 is the fecundity averaged over all categories, weighted according to the population structure [START_REF] Charlesworth | Evolution in Age-Structured Populations, 2nd Edition[END_REF][START_REF] Caswell | Matrix Population Models: Construction, Analysis, and Interpretation[END_REF].

In multiple-component systems, e.g. communities of several species, persistence can be studied for each species individually, with the background environment (parameter values in the model) being set by the other species in the community. R 0 can be calculated as before but assuming that the focal species has no impact on its environment because it is rare [START_REF] Diekmann | On the formulation and analysis of general deterministic structured population models -i. linear theory[END_REF][START_REF] Diekmann | On the formulation and analysis of general deterministic structured population models ii. nonlinear theory[END_REF][START_REF] Diekmann | Steady-state analysis of structured population models[END_REF][START_REF] Chesson | Mechanisms of maintenance of species diversity[END_REF]. In this scenario, the persistence criterion is called an invasibility criterion, since it indicates whether the focal species (the invader) is able to invade the other species (the residents). This criterion constitutes a useful fitness measure to predict evolutionary trajectories under recurrent substitutions of mutants [START_REF] Metz | How should we define fitness for general ecological scenarios[END_REF]. Persistence of the whole community is achieved when each species can invade the others (a condition sometimes called "mutual invasibility").

A difficulty arises when one wants to apply persistence criteria to spatially structured systems, such as metapopulations [START_REF] Levins | Some demographic and genetic consequences of environmental heterogeneity for biological control[END_REF][START_REF] Hanski | Metapopulation Biology: Ecology, Genetics, and Evolution[END_REF] and metacommunities [START_REF] Holt | From metapopulation dynamics to community structure: some consequences of spatial heterogeneity[END_REF][START_REF] Leibold | The metacommunity concept: a framework for multiscale community ecology[END_REF][START_REF] Holyoak | Metacommunities: Spatial Dynamics and Ecological Communities[END_REF]. Such systems are made of local subsystems (demes or patches) connected by dispersal. In this case, a raw use of R 0 as a persistence criterion leads to tedious computations and often lacks simplicity (e.g. [START_REF] Lebreton | Demographic models for subdivided populations: The renewal equation approach[END_REF].

Dispersal multiplies the number of categories to keep a census of in order to compute the fecundity function of individuals. In addition, qualifying the rarity of a species is more difficult in a spatially structured system because a species can be rare in different ways (for instance it can have very few individuals in multiple demes or have a little more of them in a unique deme). As changes of deme state are not unidirectional (a local population can increase or decrease in abundance regardless of whether it increased or decreased in the recent past), a stable deme state distribution is a bit trickier to compute than a stable age distribution in a closed population (but see [START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio r 0 in models for infectious diseases in heterogeneous populations[END_REF] for complicated, though spatially unstructured, computations of R 0 ).

The quest for a tractable persistence criterion and its application to metapopulations and metacommunities has been a recurrent subject of theoretical research (e.g. [START_REF] Chesson | Persistence of a markovian population in a patchy environment[END_REF][START_REF] Lebreton | Demographic models for subdivided populations: The renewal equation approach[END_REF][START_REF] Casagrandi | A persistence criterion for metapopulations[END_REF][START_REF] Hastings | Persistence of spatial populations depends on returning home[END_REF]. One solution is to study the dynamics of state probabilities (e.g. [START_REF] Chesson | Persistence of a markovian population in a patchy environment[END_REF][START_REF] Metz | How should we define fitness in structured metapopulation models? including an application to the calculation of evolutionarily stable dispersal strategies[END_REF][START_REF] Casagrandi | A persistence criterion for metapopulations[END_REF] rather than the deterministic dynamics of each deme (as when one uses a logistic model of population growth to predict the trajectory of the mean population abundance in each deme, e.g. [START_REF] Hastings | Persistence of spatial populations depends on returning home[END_REF]. This kind of formalism is often employed in statistical physics [START_REF] Van Kampen | Stochastic Processes in Physics and Chemistry, 3rd Edition[END_REF]. Focusing on state probabilities without monitoring the state of each deme is justified as long as the number of demes is large (effectively infinite).

Under this assumption, mean values of stochastic variables equal their average values over all the demes studied.

In two recent papers, Metz and Gyllenberg introduced the metapopulation closest equivalent to the single-population R 0 [START_REF] Gyllenberg | On fitness in structured metapopulations[END_REF][START_REF] Metz | How should we define fitness in structured metapopulation models? including an application to the calculation of evolutionarily stable dispersal strategies[END_REF]. Their fitness measure, called R m , measures the lifetime dispersal success of occupied demes in a metapopulation: R m is computed as the average number of dispersers produced by a typical colonized deme, from its colonization by a unique immigrant up to its extinction. The Whereas [START_REF] Chesson | Persistence of a markovian population in a patchy environment[END_REF] proved the validity of the R m criterion for discretetime Markovian metapopulation dynamics (under specific assumptions), no such proof has been provided for continuous-time metapopulation models. [START_REF] Gyllenberg | On fitness in structured metapopulations[END_REF] seemed to imply that the R m criterion is always true in continuous-time metapopulation models, and only mentioned that "in most cases of interest, positivity arguments" should suffice to ensure this (Gyllenberg and Metz, 2001). However, this has not been checked in the papers presenting the basic metapopulation model [START_REF] Gyllenberg | On fitness in structured metapopulations[END_REF][START_REF] Metz | How should we define fitness in structured metapopulation models? including an application to the calculation of evolutionarily stable dispersal strategies[END_REF]. And, accordingly, positivity of a matrix is a mathematical property with little biological meaning. This article has two primary goals: (i) remedying the absence of proof for the continuous-time version of the criterion and (ii) assessing whether the R m criterion of Metz and Gyllenberg depends on specific assumptions.

We use a general but sensible model of metapopulation with continuous-time dynamics (similar to the one in [START_REF] Gyllenberg | On fitness in structured metapopulations[END_REF][START_REF] Metz | How should we define fitness in structured metapopulation models? including an application to the calculation of evolutionarily stable dispersal strategies[END_REF][START_REF] Casagrandi | A persistence criterion for metapopulations[END_REF]. Classically, persistence analysis of a population or metapopulation studied in continuous time (i.e. through differential equations) is carried out by linearization of the dynamical equations. Persistence holds as soon as the Jacobian matrix qualifying the trivial (extinction) equilibrium admits at least one eigenvalue with positive real part. This method of investigation has solid mathematical background (e.g. [START_REF] Horn | Topics in Matrix Analysis[END_REF].

Using these standard Jacobian methods, we find necessary and sufficient conditions for persistence. We discuss their biological meaning and their connection to Metz and Gyllenberg's R m criterion. Our main finding is that the latter verbal criterion is acceptable in all continuous-time Markov process metapop-ulation models. We provide a recipe to obtain analytical expressions of the criterion for a broad class of models, and illustrate it with examples. Combined with [START_REF] Chesson | Persistence of a markovian population in a patchy environment[END_REF]'s proof for discrete-time Poisson processes, our results suggest that the R m criterion can be applied to typical "memoryless" metapopulations (i.e. models with no time lag between population abundance and recruitment).

Given the similarity between the viability criterion in a metapopulation and the invasion criterion in a metacommunity, we present results for single-species systems (viability) only. The extension of our results to invasion are straightforward, as discussed in the last section. We first summarize the issue and model chosen by [START_REF] Metz | How should we define fitness in structured metapopulation models? including an application to the calculation of evolutionarily stable dispersal strategies[END_REF], using a very simple structured metapopulation.

2 Rationale of the R m criterion

A simple toy model

Consider a closed spatially implicit metapopulation, formed by an infinity of demes, in which every deme has two potential microsites, and thus a population size of 0, 1 or 2 individuals. At birth, individuals may choose to disperse and, if they do so, they enter a disperser pool. Only individuals that have settled in microsites can reproduce (asexually) and there are no explicit life stages except the potential disperser stage between birth and settlement. Let µ i , λ i and γ i be the rates of per capita mortality, per capita birth, and deme extinction respectively, in a deme containing i individuals (i equals 1 or 2). Newborns in a deme with i settled individuals choose to disperse with probability d i . Let δ be the density of dispersers in the disperser pool, and assume that this density is homogeneous in space and can be described using a deterministic equation in continuous time. Let α and µ D be the rates at which dispersers try to settle in a deme and the rate at which they die, respectively.

Dispersers arriving in a deme with i settled individuals stay with probability s i . A disperser that does not stay in a deme goes back to the disperser pool.

Intuitively, s 2 must be zero so that fully occupied demes cannot receive new migrants, and d 2 must be one, so that newborns in fully occupied demes always disperse. This model is summarized in Fig. 1.

Let p k (t) be the proportion of demes that are occupied by k individuals at time t. The four quantities, p 0 , p 1 , p 2 and δ obey the following master equation:

dp 0 dt = (γ 1 + µ 1 )p 1 + γ 2 p 2 -αδs 0 p 0 dp 1 dt = αδs 0 p 0 + 2µ 2 p 2 -(γ 1 + µ 1 + λ 1 (1 -d 1 ) + αδs 1 )p 1 dp 2 dt = (αδs 1 + λ 1 (1 -d 1 ))p 1 -(γ 2 + 2µ 2 )p 2 dδ dt = λ 1 d 1 p 1 + 2λ 2 p 2 -(µ D + αs 0 p 0 + αs 1 p 1 )δ (1)
or, in matrix form:

dP dt = G(P )P (2)
where P = [p 0 , p 1 , p 2 , δ] and G(P ) is the following transition matrix function:

8

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

G(P ) =                     -αδs 0 γ 1 + µ 1 γ 2 0 αδs 0 -αδs 1 -γ 1 -µ 1 -λ 1 (1 -d 1 ) 2µ 2 0 0 αδs 1 + λ 1 (1 -d 1 ) -γ 2 -2µ 2 0 0 λ 1 d 1 2λ 2 -(α (p 0 s 0 + p 1 s 1 ) + µ D )                    
(3)

Applying the R m criterion

We want to assess whether the metapopulation is viable, i.e. grows when individuals are scarce. Intuitively, when individuals are scarce, a good criterion for viability is obtained by looking at the "disperser output" yielded by a unique colonized deme, from its colonization up to its eventual extinction. The R m criterion (or metapopulation fitness criterion) can be formulated verbally as [START_REF] Gyllenberg | On fitness in structured metapopulations[END_REF][START_REF] Metz | How should we define fitness in structured metapopulation models? including an application to the calculation of evolutionarily stable dispersal strategies[END_REF][START_REF] Parvinen | A novel fitness proxy in structured locally finite metapopulations with diploid genetics, with an apllication to dispersal evolution[END_REF]:

Criterion 1 A metapopulation is viable if, and only if, the average number of dispersers (emigrants) produced by a population founded by a typical rare immigrant is greater than 1.

Here, we develop this intuitive criterion. Suppose that the metapopulation and the disperser pool are empty. Let us impose an extrinsic flow of immigrants, say at rate ǫ. Each immigrant has probability αs 0 µ D +αs 0 of surviving dispersal and landing into an empty deme, founding a deme with one individual. Demes occupied by one individual are thus created at rate ǫ αs 0 µ D +αs 0 . We want to mon-itor the flow of dispersers that will be produced by settled individuals, under the sole effect of the extrinsic flow of immigrants (the density of dispersers is thus kept nil). Let z 1 and z 2 be the probabilities that a deme contains 1 and 2 individuals, respectively. Equ. 1 now applies to the system at the viability boundary:

dz 1 dt = 2µ 2 z 2 -(γ 1 + µ 1 + λ 1 (1 -d 1 ))z 1 + ǫ αs 0 µ D + αs 0 dz 2 dt = λ 1 (1 -d 1 )z 1 -(γ 2 + 2µ 2 )z 2 (4)
Eq. 4 has an invariant measure Z(ǫ) given by:

z 1 (ǫ) = ǫ αs 0 µ D + αs 0 γ 2 + 2µ 2 (γ 1 + µ 1 )(γ 2 + 2µ 2 ) + γ 2 λ 1 (1 -d 1 ) z 2 (ǫ) = ǫ αs 0 µ D + αs 0 λ 1 (1 -d 1 ) (γ 1 + µ 1 )(γ 2 + 2µ 2 ) + γ 2 λ 1 (1 -d 1 ) z 0 (ǫ) = 1 -z 1 (ǫ) -z 2 (ǫ) (5)
The rate at which dispersers will be intrinsically produced from settled indi-

viduals is thus [λ 1 d 1 , 2λ 2 ].[z 1 (ǫ), z 2 (ǫ)].
For the metapopulation to be viable (following our heuristic criterion), this intrinsic rate (emigrant flow) must exceed the imposed rate of immigrants (immigrant flow), ǫ, i.e.:

[λ 1 d 1 , 2λ 2 ].[z 1 (ǫ), z 2 (ǫ)] > ǫ (6)
Since the left-hand side of Eq. 6 is linearly dependent on ǫ (cf. Eq. 5), the choice of ǫ has no incidence on the viability property. We can rephrase this criterion using symbols from [START_REF] Metz | How should we define fitness in structured metapopulation models? including an application to the calculation of evolutionarily stable dispersal strategies[END_REF]

A = [0, λ 1 d 1 , 2λ 2 ] (7)
and vector Z as Z = Z(ǫ) ǫ , then the heuristic R m viability criterion is

R m = A.Z > 1 (8)
which develops as:

R m = s 0 αλ 1 (2λ 2 (1 -d 1 ) + d 1 (γ 2 + 2µ 2 )) (s 0 α + µ D )(γ 2 (γ 1 + µ 1 + λ 1 (1 -d 1 )) + 2µ 2 (γ 1 + µ 1 )) (9)
Vector Z actually describes the "typical rarity" of the species, i.e. it contains the weigths that should be given to every patch state. These weights are analogous to the "l x " coefficient used in stage-structured populations [START_REF] Caswell | Matrix Population Models: Construction, Analysis, and Interpretation[END_REF], i.e. the proportion of individuals that survive from birth to stage x.

Here, z x yields the relative weight of state x in the computation of the R m , i.e. a measure of the total time passed in state x in the course of a patch lifetime.

If Y is the probabilistic state of a patch that has just received a rare immigrant

(in our present example Y = [0, 1, 0]), then Z = -G -1 0 .Y , with G 0 being G
computed without the immigration pressure of the rare species. In terms of dimensions, Y is dimensionless and G is a matrix of (typically negative) rates, so that G -1 0 and Z are times. R m is dimensionless since it is the product of vector A (a vector of rates, i.e. inverse of times) and vector Z (eq. 8).

11

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

Proving the criterion mathematically

The criterion given in section 2.2 is formulated in biological terms: it addresses the problem of metapopulation viability using a quantity (the "lifetime production of emigrants" of a patch) with a tangible ecological meaning. The intuitive explanation of this criterion has already been given by [START_REF] Metz | How should we define fitness in structured metapopulation models? including an application to the calculation of evolutionarily stable dispersal strategies[END_REF]. This does not constitute a proof of its validity, and the simple formulation retained may not be suitable for the whole class of related metapopulation models.

Assessing whether a metapopulation described by equation 1 is viable is equivalent to assessing whether the linearized dynamical system implied by equation 1, when the metapopulation is almost empty, is unstable. Instability may be looked for by studying the eigenvalues of the linearized system: when all eigenvalues have negative real parts, the system is stable; when at least one eigenvalue has positive real part, the system is unstable. The R m criterion should thus be shown to summarize the behavior of the dominant eigenvalue of the linearized system, using a single number.

That's what we establish in the next section, providing a more general formulation of the R m criterion, and translating it into biological terms.

12

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

3 Mathematical proof of the criterion: a digest

The complete proof of the R m criterion can be found in Appendix A. The Appendix is self-sufficient, so that interested readers may read it without coming back incessantly to the main text. Here, we only give a sketch of the proof, highlighting its key issues and assumptions.

General model and assumptions

We consider a closed spatially implicit metapopulation occupied by individuals from one species. We note P the vector that describes the state of the metapopulation, i.e. given first by all the probabilities that a sample patch is found containing 0, 1, 2,...,N individuals (possibly structured by classes, sex, etc.), and then by all densities of "free" individuals that do not reside in a patch (e.g. the density of dispersers, if the model uses a disperser pool as in the toy model above). The size of vector P is at least N + 1.

We note G(P ) the "transition matrix", i.e. the matrix function that determines the master equation of P :

dP dt = G(P )P (10)
Elements of G correspond to transition rates, i.e. the rates at which populations change in state (upper diagonal block of G), the rates of free individual production (lower left-hand block), and the mortality and flows among categories of free individuals (lower diagonal block).

More specifically, the first column of G (g .,1 ) represents the rates (i) of transition from an empty population to a population containing at least one indi- and are functions of P . We assume that there are q such colonization rates,

noted (M k ) k∈[1;q] .
We are interested in determining whether eq. 1 displays an unstable equilib-

rium at P * = [1, 0, 0, ..., 0]
. This is achieved through the knowledge of the eigenvalues of equation 10, linearized near P * .

Linearization of the equation, and definition of the initial states

Colonization rates M k are assumed to obey two simple rules (see Appendix A):

(1) colonization rates are not directly affected by the proportion of empty patches (i.e. ∂M k /∂p 0 = 0 for all k);

(2) when the metapopulation is almost empty, colonization rates increase with an increase of any population state probability (except empty population state) and any density of free individuals (i.e. for all k, the vector ∂ P M k taken at P = P * has only non-negative components).

These two conditions allow the following transformation of equation 10 near the equilibrium P * (section A.3):

dP dt ≈ G(P * ).P + q k=1 (A k .P )Y k (11)
where the vectors Y k correspond to the different possible initial states of recently colonized patches (formally,

Y k = ∂ M k g .,1
), and the vectors A k are the first-order production rates of emigrants that will found populations with ini-

tial state Y k (formally, A k = ∂ P M k ).
Following equation 11, the viability problem is equivalent to finding the dominant eigenvalue of matrix G + Y k A T k at the equilibrium P * . The eigenvalue corresponding to its first row and first column is always 0 (the empty state is neutrally stable), which simplifies the issue.

From now on, tilde quantities ( X and so forth) indicate vectors and matrices whose first row/column have been removed. The final eigenvalue problem is reduced to finding the dominant eigenvalue of matrix J defined by:

J = G + q k=1 Y k A k T (12)
(for clarity, the dependence of G on P * will now be omitted).

Simplification of the eigenvalue calculation

Finding the dominant eigenvalue of J is made easier by following two simplication steps. First, it can be proved that matrix J admits one real dominant eigenvalue (section A.4).

Second, some manipulations (sections A.6 and A.7) show that if we consider 15

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

the q × q matrix R(λ), defined by its elements r ij (λ):

r ij (λ) = -A i . G -λI -1 Y j (13)
then the problem of finding the greatest real eigenvalue of J, λ, is turned into finding a matrix R(λ) that admits 1 as an eigenvalue.

Matrix R(0) contains the flows of migrant individuals produced, weighted by population lifetimes. Indeed, Y j is the j th possible initial state of a newly colonized patch.

-G -1 Y j corresponds to the quasi-equilibrium state probabilities for a population that began in state Y j (cf. the Z vector of section 2). On average, a population that began in state Y j spends a proportion -G -1 Y j k of its lifetime (i.e. up to its eventual extinction) in state k.

The product -A i . G -1 Y j thus gives the amount of emigrants that will found new populations starting in state Y i , produced by populations that initially started in state Y j . For some positive λ, the matrix R(λ) corresponds to the same quantities, but with an artificial increase in the extinction rates of the populations, quantified by λ.

Final steps and the R m criterion

The R m criterion is found by focusing on the dominant eigenvalue m(λ) of R(λ)

(see sections A.7 and A.8). The dominant eigenvalue of R is a non-negative continuous decreasing functions of λ, that converges towards 0 (since all r ij have these properties). Therefore, by checking that the dominant eigenvalue of R( 0) is strictly greater than 1, one effectively checks that there exists some positive λ such that the dominant eigenvalue of R(λ) is 1. Conversely, if such a λ exists, the dominant eigenvalue of R(0) must be greater than 1, given that 16

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

it is a decreasing function of λ.

The mathematical formulation of the R m criterion is thus: The dominant eigenvalue of this matrix, m(0), can be interpreted as the "metapopulation basic reproduction ratio", by analogy with the R 0 criterion, the basic reproduction ratio of stage-structured populations. If greater than one, the population of patches initially grows (viability), otherwise it will decline.

m(0) > 1 ⇔ Sp( J) ∩ R * + = ∅ (
When newly colonized patches have only one initial state ( Y , with the corresponding emigrant production vector, A), R reduces to a 1 × 1 matrix with only one eigenvalue equal to -A. G -1 Y (the expression given by [START_REF] Metz | How should we define fitness in structured metapopulation models? including an application to the calculation of evolutionarily stable dispersal strategies[END_REF], just like a matrix population model reduces to a scalar population model when individuals are not structured. When there are several population initial states, the computation of the R m criterion becomes more complicated, paralleling the complexification of R 0 in stage-structured populations [START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio r 0 in models for infectious diseases in heterogeneous populations[END_REF].

4 A recipe for the analytical computation of R m , with examples

The recipe

Based on our proof of the R m criterion, we give a simple five-point recipe that helps make no mistake in computing the R m value for a theoretical metapopulation model:

(1) clearly state the transition matrix G(P ), extract its first column g .,1 when P = P * , and find the colonization rates M k ;

(2) find the corresponding initial population states Y k = ∂ M k g .,1 ;

(3) compute the emigrant production vectors, A k = ∂ P M k ;

(4) compute the inverse of the reduced transition matrix, G -1 ;

(5) collect the r ij (0) = -A i . G -1 Y j in matrix R(0) (tilde versions of matrices and vectors are obtained by removing the first row/column) and compute

R m = ρ[R(0)], the dominant eigenvalue of R(0).
Obviously this recipe is only efficient to get analytical expressions for R m . For numerical calculations, it is probably better to use an approximated procedure, such as the one provided by [START_REF] Metz | How should we define fitness in structured metapopulation models? including an application to the calculation of evolutionarily stable dispersal strategies[END_REF], instead of computing the inverse of the reduced transition matrix.

We now illustrate our recipe with three examples. We first provide two examples in which the R m criterion is very simple, and then one example in which the R m criterion is more complicated to compute. All these models are simple extensions of Metz and Gyllenberg's model presented in section 2.1, and they retain the same notations and deme size N = 2. These models have not been explicitly designed to answer biological questions (even if they convey some interesting results, assumptions of these models are quite unrealistic), but rather to illustrate how to apply the R m fitness concept. Each example is presented briefly, then the recipe given in this section is applied to the model, and finally we perform a simplification of the example to get some tractable results.

Example 1: an animal model with mobile adults

The model

Consider a metapopulation as in section 2.1, but in which all individuals (adults and juveniles) can disperse (e.g. to escape crowded demes) and enter the disperser pool. Assume for simplicity that after dispersal, adults behave just as juveniles migrants do. As before, individuals that were not accepted in a deme get back to the disperser pool, hence the qualifier of active disperser pool. We slightly change the parameterization and state that e i is the per capita rate of emigration from a deme with i individuals. The model is summarized in Fig. 2.

Step 1: transition matrix and colonization rates

The G transition matrix (when N = 2) is rewritten as (cf. Casagrandi and 356 [START_REF] Casagrandi | A persistence criterion for metapopulations[END_REF][START_REF] Parvinen | Evolution of dispersal in metapopulations with local density dependence and demographic stochasticity[END_REF]:

Gatto

G(P ) =             -αδs 0 γ 1 + µ 1 + e 1 γ 2 0 αδs 0 -αδs 1 -γ 1 -µ 1 -λ -e 1 2(µ 2 + e 2 ) 0 0 αδs 1 + λ -γ 2 -2(µ 2 + e 2 ) 0 0 d 1 2d 2 -(α s i + µ D )             (15)
In this example, the first column of G only depends on one colonization rate:

M = αs 0 δ (16) 4.2.3
Step 2: initial population states 358

The initial state vector Y = ∂ M g .,1 is given by:

Y = [-1, 1, 0, 0] (17)

Step 3: emigrant production vectors

The emigrant production vector A = ∂ P M is given by:

A = [0, 0, 0, αs 0 ] (18

Step 4: inverse of the transition matrix

The reduced transition matrix is:

G =               -γ 1 -µ 1 -λ -e 1 2(µ 2 + e 2 ) 0 λ -γ 2 -2(µ 2 + e 2 ) 0 e 1 2e 2 -(αs 0 + µ D )               ( 19
)
Its inverse is:

G -1 = 1 det G         a 2 (αs 0 + µ D ) 2(µ 2 + e 2 ) (αs 0 + µ D ) 0 λ (αs 0 + µ D ) a 1 (αs 0 + µ D ) 0 a 2 e 1 + 2λe 2 2µ 2 e 1 + 2(a 1 + e 1 )e 2 2(µ 2 + e 2 )(µ 1 + γ 1 + e 1 ) + γ 2 a 1         (20)
where

a 1 = γ 1 + µ 1 + λ + e 1 (21)
a 2 = γ 2 + 2(µ 2 + e 2 ) ( 22
) det G = -[2(µ 2 + e 2 )(µ 1 + γ 1 + e 1 ) + γ 2 a 1 ] (αs 0 + µ D ) (23) 4.2.6
Step 5: R matrix and R m calculation

The value of R m is obtained from the four previous steps:

R m = s 0 α(2e 2 (e 1 + λ) + e 1 (γ 2 + 2µ 2 )) (s 0 α + µ D )(γ 2 λ + (γ 1 + µ 1 + e 1 )(γ 2 + 2(e 2 + µ 2 ))) (24)
The condition R m = 1 in the three-parameter space (e 1 , e 2 , λ) is represented in Fig. 3. This condition is a surface which separates the region of viability (above the surface) from the region of extinction (under the surface we observe is that viability is easier when e 1 is low, λ is high and e 2 is not too low (Fig. 3).

Application

We note c = µ D αs 0 +µ D the cost of dispersal. When the species is "environmentally blind" (i.e. e cannot be adjusted according to local population size) and parameters are locally population size-independent (i.e. µ 1 = µ 2 = µ and

γ 1 = γ 2 = γ), we obtain: R m = (1 -c) [2e(e + λ) + e(γ + 2µ)] γλ + (γ + µ + e) [γ + 2(e + µ)] (25)
so that the condition R m > 1 requires:

(1) λ > µ + γ, i.e. births more than compensate deaths and extinctions;

(2) a low value of c, i.e. the reward of dispersing to other patches is not too low (the maximal value of c is a complicated function of parameters;

when γ → 0, this can be approximated by c < λ-µ λ+µ 2 );

(3) an intermediate value of e (the maximal and minimal values of e are complicated functions, but when γ → 0, this simplifies to

λ-µ-c(λ+µ)- √ (1-c)[(λ-µ) 2 -c(λ+µ) 2 ] 2c < e < λ-µ-c(λ+µ)+ √ (1-c)[(λ-µ) 2 -c(λ+µ) 2 ] 2c
). This condition is intuitive since low e prevents the rescue of extinct patches (due to stochastic deaths or extinction) by occupied patches, and high e diminishes the average occupancy of patches, and thus lowers the metapopulation reproduction rate.

We can visually check these analytical conditions of viability by looking at the intersection of the plane e 1 = e 2 = e with the surface depicted in Fig. 3.

Viability is easiest to achieve when λ is high and e is intermediate (with the parameters of Fig. 3, this occurs for e ≈ 1.7).

The model

Consider now a metapopulation of plants. In this case, adults are sessile (as in section 2.1) and dispersal is passive (propagules that fail to settle in a deme just die, and do not get back to the dispersal pool: there is no "second chance", as in seed bank models). A potential colonizer arriving in a deme with i individuals has probability 1s i to die. Let us assume that propagules immediately try to settle in a deme and have vanishingly short lifespan (you may think of vegetative propagules with little nutritive and protective tissues).

This amounts to saying that the dynamics of the disperser pool is very fast compared to the dynamics of demes, so that the disperser pool always attains its equilibrium. In mathematical terms:

dδ dt = N k=0 kλ k d k p k -(α + µ D )δ = 0 (26)
which yields:

(α + µ D ) δ = iλ i d i ( 27
)
with iλ i d i = λ 1 d 1 p 1 + 2λ 2 p 2 when N = 2. We note c = µ D α+µ D the potential cost of dispersal. The realized cost of dispersal c ′ (measuring the loss in probability of settlement of a propagule relative to local reproduction) is given as the probability of dying during dispersal, plus the probability of surviving dispersal but not being able to settle in a patch, i.e.:

c ′ = c + (1 -c) (1 -s i ) (28)
The realized cost of dispersal depends on the state of the metapopulation since there is a greater cost associated to dispersal in a more crowded metapopulation. The model is summarized in Fig. 4. We need not bother about the δ column and row of the transition matrix any more since immigration pressures can be directly expressed as functions of iλ i d i (cf. equation 27).

Step 1: transition matrix and colonization rates

The G transition matrix (when N = 2) is rewritten as (cf. [START_REF] Casagrandi | A persistence criterion for metapopulations[END_REF][START_REF] Parvinen | Evolution of dispersal in metapopulations with local density dependence and demographic stochasticity[END_REF]:

G(P ) =             -iλ i d i s 0 (1 -c) γ 1 + µ 1 γ 2 iλ i d i s 0 (1 -c) -iλ i d i s 1 (1 -c) -γ 1 -µ 1 -λ 1 (1 -d 1 ) 2µ 2 0 iλ i d i s 1 (1 -c) + λ 1 (1 -d 1 ) -γ 2 -2µ 2             (29)
In this example, the first column of G only depends on one colonization rate:

M = s 0 (1 -c) iλ i d i (30)

Step 2: initial population states

The initial state vector Y = ∂ M g .,1 is given by:

Y = [-1, 1, 0] (31)

Step 3: emigrant production vectors

The emigrant production vector A = ∂ P M is given by:

A = s 0 (1 -c)[0, λ 1 d 1 , 2λ 2 ] (32) 24
Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

Step 4: inverse of the transition matrix

The reduced transition matrix is:

G =         -γ 1 -µ 1 -λ 1 (1 -d 1 ) 2µ 2 λ 1 (1 -d 1 ) -γ 2 -2µ 2         ( 33
)
Its inverse is:

G -1 = 1 det G         -γ 2 -2µ 2 -2µ 2 -λ 1 (1 -d 1 ) -γ 1 -µ 1 -λ 1 (1 -d 1 )         (34)
where

det G = [γ 1 + µ 1 + λ 1 (1 -d 1 )] [γ 2 + 2µ 2 ] -2µ 2 λ 1 (1 -d 1 ) (35)

4.3.6

Step 5: R matrix and R m calculation

The value of R m is obtained from the four previous steps:

R m = s 0 (1 -c)λ 1 (2λ 2 (1 -d 1 ) + d 1 (γ 2 + 2µ 2 )) γ 2 (γ 1 + µ 1 + λ 1 (1 -d 1 )) + 2µ 2 (γ 1 + µ 1 ) (36)
Let us compare this result to the expression of R m that we found for the opening model of section 2.1 (Eq. 9). Two major modifications were introduced in we replace s 0 (1-c) by (1-c * ) in Eq. 36, we exactly recover the R m expression of section 2.1 (see Eq. 9). This means that the second modification (instantaneous dispersal) did not have any effect, except on the experienced cost of dispersal. It is remarkable that taking into account the delay due to dispersal has no effect on the formulation of the R m criterion.

Application

When parameters are locally population size-independent (i.e.

λ 1 = λ 2 = λ, µ 1 = µ 2 = µ and γ 1 = γ 2 = γ)
, we obtain (with the notation d = d 1 ):

R m = s 0 (1 -c)λ [2λ(1 -d) + d(γ + 2µ)] γ [γ + µ + λ(1 -d)] + 2µ(γ + µ) (37)
The condition R m > 1 thus requires:

(1) λ > γ + µ (as in the previous example);

(2) one of the three following conditions:

(a) c < 1-µ+γ λ , (γ+2µ)(γ+µ+λ)-2λµ 2λ 2 (1-c) < s < µ+γ λ(1-c) and d < γ(λ+3µ+γ)+2µ 2 -2(1-c)sλ 2 λ[(1-c)s[2(λ-µ)-γ]-γ] ; (b) c < 1 -µ+γ λ and s > µ+γ λ(1-c) ; (c) 1 -µ+γ λ < c < [λ-µ-γ][γ+2(µ+λ)] 2λ 2 , s > (γ+2µ)(γ+µ+λ)-2λµ 2λ 2 (1-c) and d < γ(λ+3µ+γ)+2µ 2 -2(1-c)sλ 2 λ[(1-c)s[2(λ-µ)-γ]-γ] .
Conditions on c, s and d are intuitive: either the dispersal cost is low and the acceptance rate is high (case b), or dispersal must be limited to counteract the effect of a low acceptance rate (case a) or of a high dispersal cost (case c). A strictly positive d is not mandatory since completely occupied patches provide emigrants anyway (i.e. d 2 = 1). When dispersal cost is too high, or both acceptance rate is low and dispersal cost is high, viability is not possible, whatever the value of d.

Example 3: a model with hermaphrodites and males

The model

We now elaborate a more complex model based on the example given in section 2.1. We assume, as before, that each deme contains from 0 to 2 individuals, and we consider parameters describing the emigration rates of adults (e, as in section 4.2) and the dispersal probabilities of newborn offspring (d, as in section 4.3). We consider that the metapopulation is made up of two classes of individuals: hermaphrodites (all parameters given by index H) and males (all parameters given by index M ).

There are two mortality rates (µ H and µ M ), two emigration rates (e H and e M ) and two dispersal probabilities for offspring born in not fully occupied demes (d H and d M ; offspring born in fully occupied demes always disperse) corresponding to these two classes. All migration pressures are computed directly (i.e. there is no density of dispersers, like in section 4.3). These assumptions lead to 6 state probabilities, noted as: state 0 (empty patch), state 1 (one hermaphrodite), state 2 (2 hermaphrodites), state 3 (one male), state 4 (two males), and state 5 (one hermaphrodite and one male). In this model, hermaphrodite individuals can reproduce in two ways: they either self-fertilize (always the case in state 1, and it happens with probability σ in states 2 and 5), or outcross with another individual (probability 1σ in states 2 and 5). Self-fertilized offspring suffers from inbreeding depression, and a fraction D of these offspring die. The per capita fecundities of hermaphrodites in states 1, 2 and 5 are thus modelled as (λ is the base fecundity):

λ 1 = (1 -D)λ (38) λ 2 = λ 5 = (1 -σD)λ (39)
At birth, a proportion f of offspring are hermaphrodites. Catastrophes happen 444 at a rate γ, the cost of dispersal is c, and the two colonization rates are noted M H and M M .

Step 1: transition matrix and colonization rates

The transition matrix is quite cumbersome. Instead we give its (transposed) first column:

g T .,1 = [-(1 -c)s 0 (m H + m M ), (1 -c)s 0 m H , 0, (1 -c)s 0 m M , 0, 0] (40)
with m H and m M the raw immigration pressures defined by:

m H = e H (p 1 + 2p 2 + p 5 ) + f [λ 1 d H p 1 + 2λ 2 p 2 + λ 5 p 5 ] (41)
m M = e M (p 3 + 2p 4 + p 5 ) + (1 -f ) [λ 1 d M p 1 + 2λ 2 p 2 + λ 5 p 5 ] (42)
m H and m M acount for (i) mobile adults that emigrate (first term) and (ii) newborns that disperse (second term). The colonization rates M H and M M are given by:

M H = (1 -c)s 0 m H (43) M M = (1 -c)s 0 m M (44)
After accounting for the effect of inbreeding depression (equations 38 and 39), colonization rates are:

M H = (1-c)s 0 {e H (p 1 + 2p 2 + p 5 ) + λf [(1 -D)d H p 1 + (1 -σD) (2p 2 + p 5 )]} (45)
M M = (1-c)s 0 {e M (p 3 + 2p 4 + p 5 ) + λ(1 -f ) [(1 -D)d M p 1 + (1 -σD) (2p 2 + p 5 )]} (
Y H = [1, 0, 0, 0, 0] (47)
Y M = [0, 0, 1, 0, 0] (48)
i.e. rare hermaphrodites begin in state 1 and rare males in state 3. Note that the first element of the non-reduced Y vectors is always -1, but we directly presented reduced vectors in this example.

Step 3: emigrant production vectors

Emigrant production rates are given by equations 45 and 46:

A H = (1 -c)s 0 [e H + λf (1 -D)d H , 2e H + 2λf (1 -σD), 0, 0, e H + λf (1 -σD)] (49)
A M = (1-c)s 0 [λ(1-f )(1-D)d M , 2λ(1-f )(1-σD), e M , 2e M , e M +λ(1-f )(1-σD)] (50

Step 4: inverse of the transition matrix

The calculation of G -1 is tedious and its writing uneasy, so we will stick here to the expressions of Z H and Z M . The quasi-equilibrium distribution of states after initial hermaphrodite colonization is given by:

Z H = (γ + 2a H )(γ + a H + a M ), c H (γ + a H + a M ), c M a H (γ+2a H ) γ+a M , 0, c M (γ + 2a H ) (γ + a H + a M ) [(γ + a H + b)(γ + 2a H ) -2c H a H ] -c M a M (γ + 2a H ) (51)
where

a H = e H + µ H (52) a M = e M + µ M (53) b = (1 -D)(1 -(1 -f )d M -d H f )λ (54) c H = (1 -D)(1 -d H )λf (55)
c M = (1 -D)(1 -d M )(1 -f )λ (56)
The only component of Z H equal to 0 is the one corresponding to state 4 (two males in a patch), as this state can only be reached through the immigration of a male in a patch already occupied by one male (state 3) and hence cannot be attained in the absence of immigrants.

The quasi-equilibrium distribution of states after initial male colonization is given by:

Z M = [0, 0, 1 γ + e M + µ M , 0, 0] (57)
As males do not reproduce on their own, they either die in the patch they found or emigrate towards some other patch. Thus, the average time passed by a male in a patch (state 3) is the inverse of the sum of death, catastrophe and emigration rates.

4.4.6

Step 5: R matrix and R m calculation

The term r HM = A H .Z M equals 0, so that eigenvalues of R are easily obtained (they are the diagonal elements of the matrix). The eigenvalue associated with the male patch cycle is:

r M M = (1 -c)s 0 e M γ + e M + µ M (58)
which is always inferior to 1. The eigenvalue associated with the hermaphrodite patch cycle (r HH ) yields the expression of R m :

R m = (1 -c)s 0 {(γ + 2a H ) [(e H + (1 -D)d H f λ)a H + c M (e H + x)] + 2c H (γ + a H + a M )(e H + x)} (γ + a H + a M ) [(γ + a H + b)(γ + 2a H ) -2c H a H ] -c M a M (γ + 2a H ) (59)
where

x = λf (1 -σD) (60)

Application

The formula for R m is quite complicated to analyze in detail. In Fig. 5, regions of viability are presented using projections on two-parameter planes. The different panels of Fig. 5 suggest the existence of constraints on the covariations of c and e (for instance). These constraints seem rather strong (except perhaps for constraints due to inbreeding depression on d, e, and σ).

We can remove the dimension of λ from the expression of

R m . Let ζ = γ/λ, q = b/λ, y = x/λ, f H = e H /λ, k H = a H /λ, k M = a M /λ, l H = c H /λ, and
l M = c M /λ.
The new expression for R m is independent of λ and given by:

R m = (1 -c)s 0 {(ζ + 2k H ) [(f H + (1 -D)d H f )k H + l M (f H + y)] + 2l H (ζ + k H + k M )(f H + y)} (ζ + k H + k M ) [(ζ + k H + q)(ζ + 2k H ) -2l H k H ] -l M k M (ζ + 2k H ) (61)
At the limit f ≈ 1, the value of R m can be computed as the sum of the pure 31

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

hermaphrodite fitness, R H m and an infinitesimal term (1-f )r M that represents the contribution of males to the metapopulation fitness.

R H m = (1 -c)s 0 {f H [ζ + 2(1 -D + k H )] + 2l H y + d H (1 -D)(ζ + 2(k H -f H ))} (ζ + 2k H )(ζ + k H ) + l H ζ (62)
5 From viability to invasibility

Add-ons to the recipe

The R m criterion formulated so far is valid only to infer the viability of a metapopulation, i.e. to assess whether a single-species almost empty metapopulation can grow. A germane problem is the invasion of a species in a metacommunity. A species can invade a metacommunity if it can grow in abundance when it is very rare. The proof given in appendix A is also true for the invasibility problem (the problem only needs to be rewritten in a more general fashion, but the core arguments hold).

For invasion problems, the methodology we described has to be adapted as follows:

(1) for invasibility problems, the G matrix is much larger: there is a patch state for each combination of individuals from each species (in general, if the number of potential species is S and each patch can host at most N individuals, there are N +S N = (N +S)! N !S! patch states; when S ≫ N , the number of patch states is equivalent to

-N +S-1 N = N +S-1 N -1
dimensions, following Pascal's rule);

(2) the equilibrium state P * is not [1, 0, ..., 0] in a metacommunity , but rather correponds to the ecological equilibrium obtained when all other species are present and the focal species is absent;

(3) consequently, the initial state vectors Y k are found by taking the columns of G that correspond to non-null elements of P * , differentiating each of these columns with respect to M k and averaging the resulting vectors following the weights found in P * . Note that an initial state corresponds to a given M k colonization rate, not to a particular patch state (i.e. it is not necessary to separate average initial states into a list of elementary initial states that are obtained with a certain probability, following P * );

(4) the R m criterion is then found as in the viability case, i.e. through deriving the R(0) matrix and looking for its dominant eigenvalue. Obviously, the larger size of the matrix will make it more difficult to get analytical expressions of R m .

5.2 Does R m = 1 hold when invading identical species?

For the models studied here, R m is both an invasibility and a single-species persistence criterion. Persistence of the whole community is achieved when all species have R m > 1. However, we have left unaddressed an important issue that only arises for invasibility and not for viability: is the R m invasibility 33

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005 criterion equal to one when a species tries to invade a community of like species? Like species, i.e. species that are identical in all relevant parameters, should behave as "neutral" species.

To simplify the problem, suppose that only two species are considered, a resident species with phenotype X, and a potential invader with phenotype X ′ . Invasion of the community by the potential invader happens when R m (X ′ , X) > 1 (the first variable indicates the invader's phenotype, the second, the resident's). The problem of neutrality is: does R m (X, X) = 1 for all phenotypes X? The property R m (X, X) = 1 is highly desirable as we expect R m to behave as other fitness proxies (e.g. R 0 , r, etc.). [START_REF] Metz | How should we define fitness in structured metapopulation models? including an application to the calculation of evolutionarily stable dispersal strategies[END_REF] hypothesized that this property holds in ecologically sensible models, but no real proof of it exists. Moreover, we have found what looks like a counter-example: the metapopulation model developed by [START_REF] Slatkin | Competition and regional coexistence[END_REF] to explain the regional coexistence of two species. In this model, exactly similar species (i.e. species for which all parameters in the equations are given the same value) can generally invade one another (R m (X, X) > 1).

This odd result means mutual invasion of neutral species. Such "neutral invasions" do not occur in related models developed to address the effect of competition-colonization trade-off on species coexistence at the metacommunity level [START_REF] Hastings | Disturbance, coexistence, history, and competition for space[END_REF][START_REF] Calcagno | Coexistence in a metacommunity: the competition-colonization trade-off is not dead[END_REF].

The key to this paradoxical situation is that Slatkin's model is built under the assumption that the carrying capacity of a patch varies with the composition of the community inhabiting the patch. Specifically, the model of Slatkin describes a situation in which the carrying capacity of patch is positively correlated to the number of species in the patch. As the number of species has nothing to do with the traits of this species, two "neutral" species still benefit from increased carrying capacity. Actually, there can be "mixed" patches in which all species have identical parameters.

The opposite situation can be conceived, for instance for organisms with a strong altruistic behavior directed only towards conspecifics. In such a model, the community carrying capacity would be negatively correlated to community diversity, yielding R m (X, X) < 1.

These models both make the assumption that there is some way to discriminate species, even when these species are identical in all aspects (from the model's point of view). Hence, there can be an advantage (as in Slatkin's model) or a disadvantage (as in our hypothetical model) of being a rare species, irrespective of the traits considered in the model. This violates the symmetry in individual interactions required for neutrality, since individuals discriminate heterospecifics from conspecifics [START_REF] Chesson | Resolving the biodiversity paradox -commentary[END_REF].

Hence, R m (X, X) = 1, but this is not a paradox. Still, it makes R m unsuitable as a fitness measure. It is likely that building a "sensible" model in which R m (X, X) = 1 to address an evolutionary question is always possible (e.g. by adding a dummy trait addressing "species membership", i.e. a "green beard"), but this remains to be clarified.

35

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

Essential results

We have proved that the R m criterion in continuous time [START_REF] Gyllenberg | On fitness in structured metapopulations[END_REF][START_REF] Metz | How should we define fitness in structured metapopulation models? including an application to the calculation of evolutionarily stable dispersal strategies[END_REF]) is equivalent to the general Jacobian criterion for instability of the trivial equilibrium. We have proved that the simple R m formulation of Metz and Gyllenberg is valid only when colonization rates can be summarized using only one function of metapopulation state variables (e.g. disperser density, or the average of local population sizes, cf.

the first two examples). When colonization rates depend on several functions, we derived a more general mathematical formulation of the R m criterion. Although the original R m criterion had not been presented for models with more than one class of initial states, we have shown it still applies (see also [START_REF] Parvinen | A novel fitness proxy in structured locally finite metapopulations with diploid genetics, with an apllication to dispersal evolution[END_REF] for an application to diploid genetics).

Our proof of the validity of the R m criterion for one-species metapopulations (viability) can be generalized directly to invasibility in n-species metacommunities. In metacommunities, the species R m procedure consists in computing the mean number of new settlements due to each deme occupied by the focal species when it is rare, at a resident metacommunity steady state computed for the n -1 other species.

Assumptions required for validity

The R m criterion has been shown to hold for at least two types of models:

(i) discrete-time Markov models with Poisson-distributed immigration, and binomial emigration [START_REF] Chesson | Persistence of a markovian population in a patchy environment[END_REF], and (ii) continuous-time Markov models (present study). The fact that this criterion holds in different contexts calls for a careful examination of the assumptions underlying its validity. We first delineate the key assumptions we made to prove the criterion, and then we compare the assumptions of this general class of models with those of models studied by [START_REF] Chesson | Persistence of a markovian population in a patchy environment[END_REF].

To prove the R m criterion for continuous-time Markov metapopulation models, we have made six major assumptions:

(1) colonization rates do not depend on the proportion of empty sites, and depend positively on all other metapopulation variables when the metapopulation is almost empty;

(2) metapopulation dynamics can be described by a matrix-based master equation (such as equation 10), i.e. there is an upper limit to the number of individuals present in a site.

(3) the metapopulation consists of a large (nearly infinite) number of patches (in other words, we neglected demographic stochasticity at the metapopulation level);

(4) spatial structure was implicit (or equivalently, dispersal was uniform among demes);

(5) metapopulation dynamics are continuous, rather than discrete (i.e. there is no fixed time lag between events, events can happen any time);

(6) dynamics of the metapopulation are memoryless, i.e. the state of the metapopulation at time t only depends on its state during the immediately previous time interval of length dt.

Two questions are essential here: (i) are these assumptions really necessary for the R m criterion to hold? (ii) are these assumptions "natural" and/or "realistic"?

6.2.1 Assumption 1: colonization rates depend positively on all variables except empty patch frequency

This assumption is (intuitively) necessary for the criterion to hold because we have explicitly used it (to prove that the r ij are non-negative and decreasing in section A.7). In fact, it is possible that the R m criterion still holds when this assumption is not true, but the corresponding proof is bound to be considerably more complicated than the present one. Anyway, this assumption appears quite realistic and natural. Its contrary would seem rather odd.

6.2.2 Assumption 2: an upper limit to the number of individuals present in a site

This may not be necessary since the Perron-Frobenius theorem, which is one of the main arguments behind the proof of the criterion, also holds for linear operator in infinite dimensions. However, some points of the proof may necessitate some changes to accomodate for this change in model assumptions.

Besides, it seems quite natural to work with bounded local populations, just because the conservation of mass applies and restricts potential birth rate at the local scale when resources are limited (which is always the case).

Assumption 3: a metapopulation consisting of a very large number of patches

This assumption actually implies three things at the same time:

(1) averages equal means, i.e. theoretical expectations based on patch state probabilities are equal to empirical expectations based on actual patch states;

(2) there is no demographic stochasticity at the metapopulation level (the state of the metapopulation never deviates from its mean-field expectation, even though the state of each patch deviates stochastically from its mean-field expectation);

(3) there is no risk of stochastic absorption in the empty metapopulation state (i.e. the metapopulation never gets empty when it started from a steady state non-empty situation if the species has R m > 1).

Relaxing this assumption would mean letting different effects affect the model.

First, the absence of relation between average and mean implies that we cannot analytically resolve immigration pressures from the theoretical distribution of patch states. Instead, we have to focus on the empirical distribution of patch states, which may deviate significantly from its theoretical expectation when the number of patches is low. From a mathematical viewpoint, accepting this first effect involves at least switching from ordinary differential equations to stochastic differential equations. However, a more rigorous approach to this effect would necessitate accounting for the absence of continuity in the value of patch state probabilities (i.e. when the number of patch is K, the elementary increment of patch state probabilities is 1/K). Second, the existence of important deviations of the metapopulation state around its expected mean implies a risk of accidental extinction [START_REF] Nisbet | Modelling Fluctuating Populations[END_REF][START_REF] Otto | A Biologist's Guide to Mathematical Modeling in Ecology and Evolution[END_REF]. This risk is maximal when there is only one population and decreases monotonically as the number of patches increases. The most simple and most realistic way to model a metapopulation that displays such an effect 39 Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www. elsevier.com -doi:10.1016/j.tpb.2009.02.005 is to account for the state of each and every population, and to count the number of individuals in each population (as in [START_REF] Lebreton | Demographic models for subdivided populations: The renewal equation approach[END_REF]. However, generalizations based on such models do not look easy.

Is this assumption realistic? In real systems, the actual number of populations that are part of the same connected system can be very important. However, when it is necessary to keep the assumption of uniform dispersal (see next section), the size of the metapopulation plays the role of the mean dispersal area and can be quite low. It is thus very helpful to simultaneously state assumptions on patch numbers and dispersal range distribution.

Assumption 4: dispersal is uniform

Assuming that the metapopulation does not represent the range of populations over which dispersal can be approximated as uniformly distributed (see previous section), but rather represents the total number of populations concerned by the model, it seems questionable to assume a uniform dispersal of offspring among patches. Indeed, this assumption is only realistic for organisms that undergo passive isotropic dispersal and have numerous fast-moving offspring to attain remote patches in a small increment of time. Moreover, another issue is linked to the existence of uniform dispersal, i.e. the fact that dispersal cost is uniform for all pairs of patches. Intuitively, long distance dispersal is expected to be more costly than dispersing to a neighboring patch.

A mathematical model that incorporates non-uniform dispersal is necessarily based on explicit spatial structure (see e.g. [START_REF] Neubert | Dispersal and pattern formation in a discrete-time predator-prey model[END_REF]. When immigration in a certain patch only depends on the state of a few other patches, pair approximations can be used to lower the complexity of the problem and make 40 Although the R m criterion may apply to such a model, its mathematical formulation is not clear since the equivalent of the transition matrix is a linear operator on real-valued functions defined on space. Such a linear operator has an infinite number of dimensions, so that its eigenvalues may be more difficult to obtain.

In the present paper, reproduction, mortality and migration have been much more detailed than the ecological interactions underlying these rates. Indeed, nothing is known about the competitive interactions that imply density-dependence and competition among individuals. Implicitly, competition is assumed to happen only among individuals sharing the same patch. However, it may be that the scale at which individuals reproduce (i.e. the population scale) is different from the range at which competitive effects occur. In other words, offspring may be born in the same population as their parent(s), but individuals may compete for resources over several populations (thus implying that carrying capacities and density-dependence are shared over these populations). Taking into account these other scales (competition, predation, etc.) would complicate the model and make it deviate from its original "metapopulation" spirit. However, it is worth mentioning that dispersal is only one of the processes whose spatial scale needs to be addressed in metapopulation (or metacommunity) models.

Assumption 5: dynamics are continuous

Mathematical models aimed at describing ecological situations are classically described as either continuous-time or discrete-time models. Here, we focused on certain continuous-time models to provide a proof of the R m criterion.

However, such a proof exists for some discrete-time models [START_REF] Chesson | Persistence of a markovian population in a patchy environment[END_REF].

Thus, it does not seem that this assumption is necessary for the criterion to hold. As far as realism is concerned, discrete and continuous models can be argued for or against as easily. For instance, populations with pulsed birth events may be more realistically represented by discrete-time models, whereas populations of organisms with no stage structure and non-pulsed births are conveniently described by continuous-time models.

Relaxing the assumption of continuous dynamics affects the way immigration, emigration, birth and death are accounted for. In classical continuous-time models, these processes are described by their rates (this is the case in the present study), so that each event concerns one individual (one birth, one migrant or one death). In discrete-time models, these processes occur at fixed time steps. Thus modulating the strength of the different processes implies that an event must concern a variable number of individuals. For instance, in [START_REF] Chesson | Persistence of a markovian population in a patchy environment[END_REF], emigration of newborn offspring is binomial, i.e. each offspring has a certain probability of emigrating at birth; in the same model, immigration follows a Poisson distribution, which is the natural limit for a binomial distribution (n trials, p probability of success) when the product n.p (here, the ratio of immigrants divided by the number of patches) has a finit limit when n goes towards infinity.

It is worth remarking that events described in continuous time using a constant rate or in discrete time using a Poisson distribution yield exactly the same dynamics. Indeed, in the first case, the time between events is exponetially distributed, so that averaging over a fixed time period, the number of events that occured follows a Poisson distribution (the terms "Poisson process", "Markov arrival process", etc. actually stem from this equivalency). Hence, [START_REF] Chesson | Persistence of a markovian population in a patchy environment[END_REF] and our proofs actually show that any metapopulation model in which immigration is a Poisson process, be it in discrete or continuous time, obeys the R m criterion.

Assumption 6: dynamics are memoryless

As a consequence of using a very simple master equation to describe the dynamics of the metapopulation, the present class of models does not account for delayed effects or lags between migrant departure and arrival times. However, such a model can be envisaged as an extension of equation 10. If G τ describes the marginal effects (generalized "transition" matrix) imputable to the state of the metapopulation at a time tτ on the metapopulation at time t, the master equation could be written as:

dP dt (t) = ∞ 0 G τ [P (t -τ )]P (t -τ )dτ (63)
This equation is already a simplification of a more general problem since we have assumed that marginal effects only depend on time lag (τ ) through the state of the metapopulation (P (tτ )), and not on the actual time (tτ ).

Equation 10 is recovered by letting G τ (P ) = G(P )δ 0 (τ ) where δ 0 is Dirac's delta. Intuitively, a metapopulation described by equation 63 responds differently to the R m criterion than more simple, memoryless models:

(1) it looks quite likely that the R m computed using matrix G 0 describes ity, or stage-structure of individuals within populations, could also be incorporated in the model. Another formalism, developed by O. Diekmann and colleagues [START_REF] Diekmann | On the formulation and analysis of general deterministic structured population models -i. linear theory[END_REF][START_REF] Diekmann | On the formulation and analysis of general deterministic structured population models ii. nonlinear theory[END_REF][START_REF] Diekmann | Steady-state analysis of structured population models[END_REF], and used in the more mathematically-oriented of Metz and Gyllenberg's papers on the R m criterion [START_REF] Gyllenberg | On fitness in structured metapopulations[END_REF]) could be used to construct the proof on fewer mathematical assumptions. It could also quite simply account for stage structure. Most interestingly, a proof built within this framework would be independent of the differential equation approach we used here.

In the two existing proofs of the R m criterion's validity, viability (growing when rare) and persistence (being protected from total extinction) are equivalent.

However, several complications of the model can hinder this equivalence (e.g.

incorporating demographic stochasticity, letting the dynamics depend on the state of the metapopulation in the distant past, etc.). Clearly, the definition of the R m criterion is more suited for defining viability than persistence, i.e. a model can be conceived in which a rare species produces more emigrants than its immigrant input, and yet does not resist total extinction.

Once viability and persistence are clearly separated, the effect of demographic stochasticity on these two properties has to be assessed. It is quite possible that a metapopulation with R m > 1 is not viable due to the effect of demographic stochasticity. A well-known example of such an effect is given by studies on stochastic differential equations (Oksendal, 2007, chap. 5): when the growth rate of a population is the sum of a constant deterministic part r and of a random part αW t (where W t is a white-noise process), the population grows exponentially only when r > α 2 /2 (and not r > 0 as predicted by purely deterministic models). When 0 < r < α 2 /2, demographic stochasticity is strong enough to drive an exponentially growing population towards sure extinction. An interesting improvement to the existing R m criterion would be to account for this effect at the metapopulation scale.

PROOF. The exponential of a block diagonal matrix is still a block diagonal matrix. If we decompose M using Jordan normal form, its decomposition is block diagonal and defines its invariant subspaces [START_REF] Horn | Topics in Matrix Analysis[END_REF].

Thus, invariant subspaces of e tM are subspaces of the invariant subspaces of M . From proposition 18, the ⊃ inclusion of relation B.5 is known to be true.

The reverse inclusion comes from the fact that the exponential of a Jordan block has diagonal terms equal to e λ i t where λ i 's are eigenvalues of M (see e.g. [START_REF] Horn | Topics in Matrix Analysis[END_REF], and thus, using the characteristic polynomial of e tM , we can prove the lemma. 2

B.4 Stability of Ω through G -λI -1
Proposition 20 If λ ≥ 0 and W ∈ Ω, (i.e. ∀i, z i ≥ 0), the solution vector W ′ to the equation

W ′T = -W T G -λI -1 (B.6)
is in Ω.

PROOF. Let w ′ i be the elements of vector W ′ and g ij the elements of G.

W ′T G -λI = -W T can be decomposed for each column j:

(g jj -λ)w ′ j + i =j g ij w ′ i = -w j (B.7)
Given that all diagonal elements g jj of G are strictly negative (due to extinction and mortality rates), λg jj > 0 and thus Eq. B.7 can be further modified:

w ′ j - i =j g ij λ -g jj w ′ i = w j λ -g jj (B.8)
Now let χ j = i =j g ij w ′ i . Multiplying Eq. B.8 by g jk and summing over j, we get for all k: χ kj =k g jk λg jj χ j = j =k g jk w j λg jj ≥ 0 (B.9) Let χ j = min k (χ k ). Then, Eq. B.9 implies that: Proposition 21 G -λI -1 → 0 when λ ∈ R + + iR and ℜ(λ) → ∞.

χ j   1 - j =k g jk λ -g jj

994

PROOF. First, assume that G is diagonalizable in C. Let Q be the invertible complex matrix which changes G into its Jordan normal form, i.e.

Q -1 GQ =               λ 1 0 0 0 ... 0 0 0 λ n               (B.11)
where λ i 's all have negative real parts.

62

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

Let λ be a complex number with positive real part. From Eq. B.11, we obtain:

G -λI = Q               λ 1 -λ 0 0 0 ... 0 0 0 λ n -λ               Q -1 (B.12)
and thus

G -λI -1 = Q -1               1 λ 1 -λ 0 0 0 ... 0 0 0 1 λn-λ               Q (B.13)
From Eq. B.13 and using the definition of the spectral norm (which is consistent with the natural scalar product on vectors), we immediately deduce that

G -λI -1 ≤ Q -1 max 1 λ i -λ 2 Q (B.14)
Now we observe that

1 λ i -λ = ℜ(λ i ) -ℜ(λ) -i (ℑ(λ i ) -ℑ(λ)) (ℜ(λ i ) -ℜ(λ)) 2 + (ℑ(λ i ) -ℑ(λ)) 2 ∼ -1 ℜ(λ) (B.

  interpretation of the criterion Interpreting the R m criterion in biological terms requires understanding the signification of m(0). As already stated, the matrix R(0) contains the flows of emigrant productions weighted by population lifetimes. Each column of this matrix corresponds to a colonization pathway, i.e. an initial state Y k . The life cycles of patches are thus distinguished in different classes, depending on the initial state in which they began. Each class of life cycle initiates other classes of life cycles (the rows of R(0)) throughout its life. R(0) can thus be seen as a transition matrix between classes of patch life cycles, just as a Lefkovitch matrix is a transition matrix between classes of individuals.

  the present model: passive colonization and instantaneous dispersal. Passive colonization introduces a difference in the way dispersers that failed to settle are treated: they die instead of getting back to the disperser pool. This basically changes the cost of dispersal. With active colonization, the physiological cost of dispersal would be c(P ) = µ D α s i +µ D Unsurprisingly, it would be reduced (as dispersers have "second chances") and it would depend on the metapopulation state (as dispersers are more likely to leave the disperser pool in less crowded metapopulations). The value of this cost in a nearly empty metapopulation (as is relevant for viability analysis) would be c * = µ D αs 0 +µ D If

  446

  sum of each column of G is negative, and thus χ j ≥ 0. If we inject this result into Eq. B.8, the proposition is then finally proved. 2 992 B.5 Convergence of G -λI -1

Figure 1 .Figure 2 .Figure 3 .

 123 Figure 1. Schematic representation of Metz and Gyllenberg's model when N = 2. Solid arrows represent flows proportional to the source, dotted arrows represent flows proportional to the source times the density of dispersers. Symbols are as in the text.

20

  )Massol F.,Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

  Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

). What 21

  Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

386

23

  Massol F., Calcagno V.,Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005 vectors and matrices (i.e. X) are obtained by removing all rows and columns that correspond to a patch state empty of the focal species and densities of free individuals from other species (so, for patch states, tilde vectors have N +S N

	S N N ! as S → ∞). For instance, if
	only 2 species are present and the maximum number of individuals per
	patch is fixed at 3, there are 10 possible patch states. Moreover, reduced
	32

  Massol F., Calcagno V.,Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005 it analytically tractable. When populations are small and far apart, the model can be approximated by a point process in continuous space. In that context, a good approximation is to consider diffusion equations associated with the probability density p k (x, t) of observing a population with k individuals at point x at time t. Such a probability density function obeys a Fokker-Planck equation that can be analytically studied and/or numerically solved for steady states.

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

Acknowledgements We thank P. David, A. Duputié, J.-D. Lebreton, and L. Mercier for advice, P. Chesson for acess to his early papers, and J. Chave, P. Chesson, P. Jarne, M. Leibold, K. Parvinen, O. Ronce, F. Rousset and an anonymous reviewer for useful comments on the manuscript. This work was funded by the French Ministry of Agriculture to F. M., by the French Ministry of Research to V. C., and by the French Ministry of Education to J. M.

well the initial growth of a sparsely occupied metapopulation: when the metapopulation receives its first immigrant, the long-term past is unimportant;

(2) on the other hand, when this R m value is superior to 1, the metapopulation may eventually ends in the extinct state due to delayed negative effects. For instance, if a metapopulation can overshoot its steady state mean density (due to, e.g., accelerating birth rates) and that overshooting this steady state mean density at time t causes deterministic extinction of the whole metapopulation at time t + T , then the metapopulation will periodically starts from nothing (R m > 1) and then fall to an empty state.

As the empty metapopulation state is an absorbing state when there is no extraneous immigration, the fitness of an initially rare species that obeys the above example should be 0. However, a very small extraneous immigration pressure creates a time-averaged metapopulation state quite different from the empty state. To resolve such a paradoxical situation, we need to take into account the existence of infinitesimal immigration from outside of the metapopulation, and to define weak persistence as the fact that an initially rare species becomes more common (on average over time). To assess this property we have two options:

(1) compute a long-term average fitness proxy for the metapopulation. For population models, this is usually done by computing the geometric average of the R 0 value over an infinite number of generations. However, such a long-term average requires that the fitness proxy can be used far from the trivial empty equilibrium. In the case of metapopulations, this translates to the simple question: can metapopulation growth be represented 44 Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www. elsevier.com -doi:10.1016/j.tpb.2009.02.005 by the R m when the metapopulation is already in a non-empty state? If this proves true, then the basic reproduction ratio of a metapopulation over a time period of length τ would be obtained as 1 τ exp τ 0 log R m (t)dt.

(2) compute the growth rate of the mean density of individuals in populations, N , using the fact that d N /dt = [0, 1, 2, ...].(dP/dt) (i.e. increase in the mean density can be computed using equation 63). This fitness proxy is usually less convenient than the former one because in most case this computation involves the knowledge of higher moments (variance, skewness, ...) of the individual density distribution.

In any case, the two roles that have been played by the R m criterion in simple models, i.e. (i) assessing whether a lonely individual can invade a landscape and (ii) assessing whether a metapopulation is protected from deterministic extinction, are dissociated in memoryless models. Thus, a metapopulation might initially grow when rare but also be subject deterministic extinction once well developed, or, on the other hand, a metapopulation might be incapable of starting from scratch but displays a steady state protected from deterministic extinctions (similarly to a population with an Allee effect).

Challenges and perspectives

The R m criterion has not been proved in spatially explicit metapopulation model. A conceptual and mathematical puzzle now lies in the way spatially explicit metapopulations can be described in order to make them realistic while keeping the model analytically tractable.

Other complications, such as temporal variability, environmental heterogene-A Appendix: proof of the criterion

A.1 Notations

The mathematical notations we use are presented in Table A.1.

With these notations, the invariant subspace of matrix M associated with eigenvalue λ is ker(M -λI) ν (λ) . By Jordan decomposition theorem [START_REF] Horn | Matrix Analysis[END_REF], λ) .

A matrix is non-negative if all its elements are positive or zero.

47

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives The ith row / the jth column vector of M X T Transposed of X, be it a vector or a matrix X.Y Scalar product of vectors X and Y

The set of eigenvalues of matrix M Vect(X) Vector subspace induced by vector X

ker(M )

The kernel of matrix M (subspace of vectors X such that M.X = 0)

The spectral radius of matrix M (highest modulus of its eigenvalues)

The index of eigenvalue λ of a matrix

The space of square real matrices of dimension n f • g Composition of functions f and g We consider a closed spatially implicit metapopulation occupied by individuals from one species. We note P the vector that describes the state of the metapopulation, i.e. given first by all the probabilities that a sample patch is found containing 0, 1, 2,... individuals (possibly structured by classes, sex, etc.), and then by all densities of "free" individuals that do not reside in a patch (e.g. the density of dispersers if the model makes use of a disperser pool). The size of vector P is noted N + 1. We note G(P ) the "transition matrix", i.e. the matrix function that determines the master equation of P :

Matrix G(P ) contains the rates of birth, dispersal, death and extinction in the metapopulation. As any transition matrix, G has negative diagonal elements.

The upper diagonal block of G determines the transition matrix of patch state probabilities. Thus, it is a conservative transition matrix, i.e. the sum of each of its column equals 0 (because the sum of probabilities remains equal to 1). The upper right-hand block of G is full of 0 (by construction: we model only flows among patch states mediated by densities of free individuals, not flows from densities of free individuals to patch states mediated by other patch states), its lower diagonal block describes the dynamics of the densities of free individuals (possibly in different classes), and its lower left-hand block describes the rates of production of free individuals in the different patch states. In general, the following relationships hold: Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

where g ij is the element of G on the i th row and j th column. As G depends on P , each element g ij is a function g ij (P ). We assume that these functions are sufficiently smooth and derivable to allow the required differential expressions.

In particular, the first column of G, g .,1 , gives the transition rates from an empty patch state to other states. In the most simple models, only g 2,1 ≥ 0 (i.e. the only possible transition for an empty patch is to be colonized by one individual). More complicated models can be envisaged in which individuals can belong to different classes (e.g. wingless and winged individuals in aphids), so that more than one component of g .,1 may be positive. In any case, positive components of g .,1 correspond to colonization rates.

We assume that q of these colonization rates are positive, and these rates,

, are functions of P . We consider that the components of P are independent in the determination of M k (i.e. the fact that the sum of state probabilities equates 1 is not used to simplify the expression for M k ), so that ∂M k /∂p 0 = 0 for all k (i.e. colonization rates do not depend directly on the proportions of empty patches). Moreover, we assume that all colonization rates respond positively to all increases of non-empty state probabilities or of densities of free individuals when the actual metapopulation is completely empty (state P * = [1, 0, 0, ..., 0]), i.e. for all k, the vector ∂ P M k taken at P = P * has only non-negative components.

A.3 Linearization of the master equation

Around an equilibrium P * of equation A.1, the dynamics of small deviations of P can be linearized using the differential form ∂ P G of G:

The viability of the metapopulation can be studied through the instability of Eq. A.2 near the special equilibrium in which all demes are empty (P * = [1, 0, 0, ..., 0]). In this case, the only components of ∂ P G(P * )(P -P * ) that will be used to compute dP dt are the ones associated with the first column of G(P ), as it can be seen in the following development:

As we have assumed that g .,1 depends on only q variables (the M k ), we can further develop ∂ P G(P * )dP.P * :

Because g .,1 is a vector function of q real variables, partial differentials, ∂ M k g .,1 (P * ), are vectors. ∂ P M k is a vector, noted A k (which has only non-negative components, cf. section A.2). Let Y 1 , Y 2 , ..., Y q be the vectors defined by:

Eq. A.2 can be rewritten as:

Instability of Eq. A.6 is achieved if and only if the Jacobian matrix J(P * ) defined by

has at least one eigenvalue with positive real part.

As we consider a closed metapopulation, colonization rates equal 0 when all patches are empty. In mathematical terms, this statement amounts to g .,1 (P * ) = 0. Because the first column of Y k A T k is a null vector (the first components of A k is equal to ∂M k /∂p 0 , which equals 0), it follows that matrix J(P * ) is block upper triangular. The first diagonal element of G(P * ), which corresponds to the first eigenvalue of J(P * ), is 0. Thus, proving the strict instability of Eq. A.6 is equivalent to studying whether the lower diagonal block of J(P * ) has an eigenvalue with positive real part.

To that end, we introduce the tilde operator, which is defined as follows: if M is a N + 1 square matrix, M is defined as the lower N × N diagonal block of M (obtained by removing its first row and first column). Similarly, if V is a N + 1 vector, V is defined as the vector obtained from V by removing its first element. Matrix J (= J(P * )) verifies the following relation:

Using these notations, viability of the metapopulation is achieved if, and only if, matrix J has at least one eigenvalue with positive real part. The reduced version of Eq. A.6 is:

We shall now derive necessary and sufficient conditions for viability of the 52

Remark 2 We note that the vectors Y k are independent. Indeed, each Y k corresponds to the derivative of g .,1 by one of the M k , which in turn applies to a specific colonization transition.

Remark 3 We note that Y k and A k are in Ω (i.e. have only non-negative components). This is a result of the construction of Y k (equation A.5) and a "natural" constraint on the relationship between state probabilities, free individual densities and migration pressure at the equilibrium (section A.2).

Remark 4 If X ∈ Ω and Y ∈ Ω, X.Y ≥ 0. This is straightforward: all vectors in Ω have only non-negative components, so that scalar products between two elements of Ω are non-negative.

A.4 Uniqueness of the dominant eigenvalue

Differential equation A.9 has all its solutions of the form t → e t J .V . In appendix B, we prove that all these solutions stay in Ω if they begin in it. The following proposition, based on results given in Appendix B reduces the problem of looking for eigenvalues of J with positive real parts.

Proposition 5 J only admits one eigenvalue with highest real part, and this eigenvalue is real.

PROOF. From proposition 17 given in Appendix B applied to t = 1, there is an eigenvalue of e J equal to ρ(e t J ). From lemma 19 given in Appendix B, there is an eigenvalue λ of J such that ρ(e t J ) = e λ . For all other eigenvalues µ of J, e µ is an eigenvalue of e J , and thus, from proposition 17, e ℜ(µ) < e ℜ(λ) which proves that λ has the highest real part of all elements of Sp( J). The same line of reasoning leads to conclude that for all non-negative t, ρ(e t J ) = e λt . From proposition 17, this means that e λt is real for all t ≥ 0, or equivalently that ℑ(λt) ≡ 0 [2π] for all t ≥ 0. This last statement implies that ℑ(λ) = 0, and λ must be real. Uniqueness of λ is also achieved since all candidates λ must be real, so that if their real parts are equal, eigenvalues are equal. 2

A.5 Invertibility of the reduced transition matrix

Proposition 6 If λ ≥ 0, G -λI is invertible.

PROOF. We remark that G is invertible because:

(i) its upper diagonal block (patch state probabilities) obeys Hadamard's theorem (its diagonal elements are negative and greater, in absolute value, than the sum of all non-diagonal elements in each column),

(ii) its lower diagonal block (densities of free individuals) has only strictly negative components on its diagonal (the effect of mortality and outward flow rates).

Using the same argument, we conclude that G-λI is invertible for all positive λ. 2

54

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

Corollary 7 G has no positive eigenvalue.

914

A.6 Transformation of the eigenvalue problem

Consider the eigenvalue problem defined by J.V = λV . This problem is equivalent to finding λ and V such that ( J -λI).V = 0, i.e. whether 0 ∈ Sp( J -λI).

In the previous section, we have shown that G -λI is invertible when λ ≥ 0.

For non-negative λ, let Z k (λ) be defined as:

Remark 8 As a corollary of proposition 20 and remark 3, we note that all 916 scalar products A i .Z j (λ) are non-negative (cf. remark 4)..

The following identity: .11) insures that, for all non-negative λ:

We now define matrix H(λ) (or more simply H):

Equation A.12 can be reformulated as:

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

A.7 Reduction of the H(λ) matrix

Let Ξ be the vector subspace induced by the family of Z k (λ) vectors , i.e. Ξ = Vect(Z 1 , ..., Z q ). Following remark 2, the family of Z k (λ) vectors is independent, and thus dim Ξ = q. We remark that Im (H) ⊂ Ξ, so that H.Ξ ⊂ Ξ.

Thus, we may define the reduction H Ξ of H on Ξ. This also implies that all eigenvectors of H associated with non-null eigenvalues are elements of Ξ. We note R(λ) (or more simply R) the q × q matrix representing H Ξ on the basis formed by vectors Z 1 , Z 2 , ..., Z q . Elements r ij (λ) of R(λ) are given by the scalar product:

Proposition 9 For all i and j, and as long as λ remains non-negative, r ij is a non-negative decreasing function of λ. Besides, r ij (λ) → 0 when λ → ∞.

PROOF. All r ij are non-negative (cf. 8). Now

real, non-negative and converges to 0 as λ → ∞. It is a continuous function 56 Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005 of λ.

PROOF. From proposition 9, R(λ) is non-negative, so that the Perron-Frobenius theorem applies, i.e. m(λ) is real and

which proves almost all results (cf. proposition 9). m is clearly continuous since all r ij are continuous functions of λ. 2

In order to prove that m is a decreasing function of λ, we need the following lemma:

Lemma 11 Let Φ be the set of decreasing functions of R + → R + . Let f and g be elements of Φ. The following propositions hold:

(i) f + g ∈ Φ;

(ii) f g ∈ Φ;

(iii) ∀α > 0, f α ∈ Φ.

PROOF. The proof of it is straightforward. 2

Corollary 12 For all k ≥ 1, all elements of R(λ) k are non-negative decreasing functions of λ.

We can now prove that:

Proposition 13 m is a decreasing function of λ.

57

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Let R m be defined as:

The "R m " theorem can be stated as:

PROOF. If R m > 1, there exists λ > 0 such that m(λ) = 1 (because m is decreasing, continuous and converges towards 0 as its variable tends to infinity).

B.1 Boundedness of Eq. A.9 solutions

Differential equation A.9 has all its solutions of the form t → e t J .V . To show that all these solutions stay in Ω if they begin in it, we show the following proposition:

Proposition 15 For all i ∈ [1; N ] and for all V ∈ Ω, if E i .V = 0, then

PROOF. Suppose that V = (v i ) is a vector on the boundary of Ω such that

JV can be developed:

so that

The second series of terms in the above sum is positive since the A k and Y k are in Ω (remark 4). The first term develops as follows:

Because G is pseudo-stochastic, we have g ii ≤ 0, and ∀j = i, g ij ≥ 0. As v i = 0 and ∀j, v j ≥ 0, it follows that E i . GV ≥ 0. 2

B.2 Positivity property

Proposition 16 For all non-negative t, matrix e t J is non-negative.

59

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005

PROOF. Proposition 15 means that for all V on the boundary of Ω, the speed vector of differential equation A.9 points inwards. This entails that for all V in Ω and for all non-negative t, e t J V stays in Ω. Taking the canonical basis (E i ), which is in Ω, it follows that all columns of e t J are in Ω, and thus e t J is non-negative. 2

Proposition 17 For all non-negative t, ρ(e t J ) is a non-negative real eigenvalue of e t J . In particular, for all λ in Sp(e t J ) \ {ρ(e t J )}, |λ| < ρ(e t J ).

PROOF. The proof results from the application of the Perron-Frobenius theorem to e t J , which is non-negative. 2

B.3 A spectral mapping lemma

Proposition 18 If λ ∈ C is an eigenvalue of matrix M and V ∈ C n an associated eigenvector, V is also an eigenvector of matrix e M associated with eigenvalue e λ .

PROOF.

72

Massol F., Calcagno V., Massol J. (2009) The metapopulation fitness criterion: proof and perspectives, Theoretical Population Biology, 75, 2-3. author-produced version of the final draft post-refeering the original publication is available at http://www.elsevier.com -doi:10.1016/j.tpb.2009.02.005