## Abstract

In this work we deal with a multiregional model in discrete time for an age-structured population which lives in an environment that changes randomly with time and is distributed in different spatial patches. In addition, and as is often the case in applications, we assume that migration is fast with respect to demography. Using approximate aggregation techniques we make use of the existence of different time scales in the model and reduce the dimension of the system obtaining a stochastic Leslie model in which the variables are the total population in each age class. Literature shows that, under reasonable conditions, the distribution of population size in matrix models with environmental stochasticity is asymptotically lognormal, and is characterized by two parameters, stochastic growth rate (s.g.r.) and scaled logarithmic variance (s.l.v.), that, in most practical cases, cannot be computed exactly. We show that the s.g.r. and the s.l.v. of the original multiregional model can be approximated by those corresponding to the reduced stochastic Leslie model, therefore simplifying its analysis. Moreover, we illustrate the usefulness of the reduction procedure by presenting some practical cases in which, although the explicit computation of the s.g.r. and the s.l.v. of the original multiregional model is not feasible, we can calculate its analogues for the reduced model.

## 1. Introduction

The so-called ‘approximate aggregation methods’ are mathematical techniques in which appropriate approximations are introduced in order to transform the system under consideration into a reduced or ‘aggregated’ system with a lesser number of variables called ‘global variables’. These techniques are usually applied in systems governed by processes with different time scales, and the approximation usually involves freezing the slow process and letting the fast process reach equilibrium. Once this is done one can build a reduced or ‘aggregated’ system. The problem then is to justify that one can extract important information about the behaviour of the original system through the study of the simpler reduced one. The reduction in the dimension of the system makes these techniques useful to study some of the complex models that we find in nature in which there are different time scales involved.

Approximate aggregation techniques have been widely studied in the context of time continuous systems with different-time scales for both linear and density-dependent models (see Auger *et al.* (2008) for a review and a list of references). We will focus our attention on the discrete time case, which has been thoroughly explored in linear, nonlinear and nonautonomous deterministic contexts (see among others Marvá *et al.* (2008); Sanz & Bravo de la Parra (1999)).

Nowadays, a great part of the research in ecological modelling is devoted to models which incorporate stochasticity as a way of taking into account different factors which either have an intrinsic random nature or are too complex to be taken into account in a deterministic way. The literature contemplates two sources of stochasticity: demographic (Caswell 2001) and environmental. The reduction of models that incorporate demographic stochasticity has been addressed in Sanz *et al.* (2003).

In this work we contemplate linear discrete models that incorporate ‘environmental stochasticity’, i.e. the randomness introduced when we consider random fluctuations in the environment and, consequently, in the vital rates which affect the population (see Caswell (2001) and Tuljapurkar & Caswell (1997, ch. 3) for an introduction; Tuljapurkar (1990) for a full discussion; and Cohen *et al.* (1983) for an application to actual fish populations). These models are analogous to the deterministic ones but in this case the matrix of vital rates in each projection interval is selected within a given set of matrices according to a certain (possibly time-varying) probability distribution. Given certain hypotheses on the pattern of temporal variation and the vital rates in each environment, the distribution of total population size is asymptotically lognormal, with an expected value and a variance dependent on two parameters: the stochastic growth rate (s.g.r.), which is the stochastic analogue of the logarithm of the dominant eigenvalue for deterministic systems, and the scaled logarithmic variance (s.l.v.), which characterizes the asymptotic deviation of the population from its mean value. In most practical situations it is not possible to compute the s.g.r. and the s.l.v. analytically and one must resort to computer simulations.

The reduction of discrete time systems for populations subjected to environmental stochasticity is carried out in Sanz & Bravo de la Parra (2000), and both this reference and Sanz & Bravo de la Parra (2007) provide certain results that relate the behaviour of the original and the aggregated system. In particular, in the case in which there are a finite number of environments, Sanz & Bravo de la Parra (2007) relates the s.g.r. of the original and the aggregated system.

The purpose of this work is to extend the abovementioned results in two directions. In the first place, and placing ourselves in a context in which the number of environmental states can be infinite, we prove that under reasonable assumptions the s.g.r. and the s.l.v. of the original system can be approximated by those of the aggregated system when the separation of time scales between demography and migration is large enough. This allows us to approximate the population distribution of the original system in terms of the s.g.r. and the s.l.v. of the reduced system.

In the second place, we show how, in some particular cases, the aggregation procedure renders a reduced system for which the s.g.r. and the s.l.v. can be computed exactly. This allows one to carry out a qualitative study of the behaviour of these models as functions of their parameters that is not possible when dealing directly with the original system. Moreover, both the reduction procedure and the results that relate the original and the aggregated systems are valid in a wider context than that of multiregional models, including any stage-structured model for a population living in a multipatch system.

The structure of the paper is as follows: §2 briefly introduces the basic form of the matrix models that consider environmental stochasticity and the main results that will be used in the rest of the paper, concerning the distribution of population size when the environmental variation is a homogeneous Markov chain. Section 3 is devoted to presenting a stochastic multiregional model in which population is structured by age and patch and migration among the different patches is fast with respect to demography. By defining the global variables as the total population in each age class, we obtain a reduced stochastic Leslie model. The model, together with the aggregation procedure, is a particular case of a technique which was first presented by the authors in Sanz & Bravo de la Parra (2000) and which is valid, for example, for any stage-structured population living in a multipatch environment.

Section 4 presents the abovementioned relationships between the original system and the reduced system. The proof of the most technical results is deferred to an appendix. In §5, where ecologists mainly interested in applications may focus their attention, we show how in some cases a non-trivial multiregional system, for which the s.g.r. and s.l.v. cannot be calculated exactly, can be aggregated to a stochastic Leslie model for which that calculation is possible. We first examine a non-structured population living in a multipatch environment and then we deal with a population with juveniles and adults living in a two-patch environment.

## 2. Matrix models with environmental stochasticity

This section presents the basic form of the matrix models that consider environmental stochasticity. We assume then that the population lives in an ambient in which there are different environmental states that we suppose indexed with , where is a (possibly infinite, numerable or denumerable) set of indexes. The vital rates corresponding to each one of these environments are given by the non-negative matrices , in such a way that, for each *ζ*, **A**_{ζ} represents the vital rates of the population in environment *ζ* and is the set of different vital rates.

The environmental variation is characterized by a sequence of random variables *τ*_{n}, *n*=0,1,2,… defined on a certain probability space and with state space . For each realization *ω*∈Ω of the process, the population is subjected to environmental conditions *τ*_{n+1}(*ω*) between times *n* and *n*+1.

Thus, the model reads:
2.1where for each is a vector random variable in which represents the population vector at time *n*. Throughout we will assume that **Z**_{0} is a fixed (non-random) non-zero vector **z**_{ 0}≥**0**. We will also consider the structure of the population vector
2.2where ∥*∥ denotes (unless specifically stated) the 1-norm in , i.e. . Let be the set of possible population structures (column vectors with non-negative components and unitary in the 1-norm) and let be the *σ*-algebra of Borel sets of .

Following a great part of the literature on the field, we will deal exclusively with the case in which the environmental variation *τ*_{n} is a homogeneous Markov chain (which includes as a particular case the i.i.d. setting). Let us define the *n*-step transition probabilities

Then it is immediate to check that the bivariate process is a homogeneous Markov chain with state space . Following Cohen (1977) we define the transition probabilities of this chain as 2.3for , , , .

Let us recall some features regarding the asymptotic behaviour of general matrix models with environmental stochasticity like equation (2.1). In order to do so, we introduce the following concept.

## Definition 2.1.

A set of *N* × *N* non-negative matrices is said to be ergodic if there exists a positive integer g such that any product of g matrices (with repetitions allowed) drawn from is a positive matrix (i.e. its components are all positive) and moreover there exist constants *α*,*β*>0 such that for all
where is the maximum of the elements of **A**_{ζ} and min^{+}(**A**_{ζ}) is the minimum of the positive elements of **A**_{ζ}.

In particular, any finite set of *N*×*N* non-negative matrices which share a common primitive incidence matrix is an ergodic set. More generally, if the condition **C**≥**A**_{ζ}≥**B**, with **B** and **C** primitive matrices holds, then is an ergodic set. Note that if the set is ergodic and **z**_{ 0}≥**0**, **z**_{ 0}≠**0** then **Z**_{n}≠**0** for all *n* with probability 1 and so the population structure (2.2) is well defined.

In what follows we will assume that the vital rates and the pattern of environmental variation for system (2.1) meet the following hypotheses.

(H1) The chain

*τ*_{n}has a unique stationary probability distribution and moreover 2.4where , 0≤*ρ*<1 are constants.(H2) The set of matrices of vital rates is ergodic.

If there is a finite number of environmental conditions, i.e. if *τ*_{n} takes only a finite number of values then hypothesis (H1) is met if the chain is irreducible and aperiodic, i.e. if its transition probabilities matrix is primitive (Doob 1953).

Under hypotheses (H1) and (H2), from Cohen (1977) we have that system (2.1) is stochastically strongly ergodic in the sense that the *n*-step transition probabilities *G*_{n} for system (2.1) converge to a stationary probability distribution *F* uniformly in the initial state of the chain. Specifically,

## Theorem 2.2.

(Cohen 1977, theorem 3.) *Let us assume that system (2.1 ) meets hypotheses (H1) and (H2). Then there exists a stationary probability distribution F for* **D**_{n} *such that, given ε>0,* *and* *such that* *,*

Now let us turn to the study of the total population size ∥**Z**_{n}∥. In this regard we have the following theorem.

## Theorem 2.3.

(Tuljapurkar & Orzack 1980, essentially lemma 4.) *Under hypotheses (H1) and (H2)*,

*We can define the s.g.r. of system (2.1 ) through**with probability 1. Moreover, a is finite, and is independent of the initial probabilities of the chain and of the initial (non-zero) population vector***z**_{ 0}≥2.5**0**. Parameter*a*can also be computed through*where the averages with the F subscript are taken with respect to the stationary distribution F and assuming that***z**_{ 0}*is positive*.*We can also define the s.l.v. of the system as**where**denotes variance and where σ*^{2}*is independent of the initial probabilities of the chain and of the initial (non-zero) population vector***z**_{ 0}≥**0**. Moreover, σ^{2}*is finite and**where*2.6*If σ*^{2}>0*the population size is asymptotically lognormal in the sense that**where**denotes a normal distribution of zero mean and unit variance and**denotes convergence in distribution*.

## Proof.

Essentially this result is lemma 4 in Tuljapurkar & Orzack (1980) but in this reference the authors impose an additional condition (expression (11’) in Tuljapurkar & Orzack (1980) or, equivalently, condition (3.13) in Furstenberg & Kesten (1960)) that, we will show, is not necessary given (H1) and (H2). Indeed, let us define variables

Lemma 4 of Tuljapurkar & Orzack (1980) hinges on theorem 3 in Furstenberg & Kesten (1960), in which the authors apply the central limit theorem to the variables *ξ*_{n} (which in that reference are not necessarily bounded) and they need an additional condition to guarantee that the moments of *ξ*_{n} are finite. However, using (H2) it is straightforward to check that there exist constants *C* and *D* such that,
with *C* and *D* depending only on the constants *g*,*α*,*β* corresponding to the ergodic set . Therefore, |*ξ*_{n}| is bounded for each *n* and therefore its moments are finite. The rest of the proof is identical to that of Tuljapurkar & Orzack (1980). ▪

The convergence to the lognormal distribution for the population size has also been established under similar but weaker hypotheses in Heyde & Cohen (1985).

Note that if *a*>0 then ∥**Z**_{n}∥ grows exponentially for every realization and if *a*<0 then the population goes extinct with probability 1.

It is important to stress that, in most practical cases, the exact computation of *a* and *σ*^{2} is not feasible due to the fact that in general one cannot obtain an analytical expression for the stationary probability distribution *F* that intervenes in equations (2.5) and (2.6). In §5 we make use of the aggregation procedure to be able to obtain exact results for certain multiregional models with fast migration.

## 3. Reduction of a stochastic multiregional model with fast migration

In this section we introduce the model and the reduction procedure for a certain kind of stochastic multiregional models in which migration is fast with respect to demography. The procedure is a particular case of that introduced in Sanz & Bravo de la Parra (2000) to reduce general models which incorporate environmental stochasticity and in which there are two time scales.

Multiregional models consider the dynamics of a possibly age-structured population distributed among different spatial patches among which they can migrate. These models have been used for the study of human populations (Rogers 1995) and in ecology (Caswell 2001). The usual approach has been deterministic, but the stochastic setting has also been used (Cohen 1982).

The above references do not explicitly consider the existence of different time scales in the multiregional system. The fact that in many situations migration and demography take place with different time scales (usually migration is fast with respect to demography although the contrary may also happen; see Liaw (1980)) has been exploited in several works in order to reduce the complexity of the multiregional model. Some examples in a linear deterministic setting are that of Sanz & Bravo de la Parra (1999) (autonomous case) and Sanz & Bravo de la Parra (2001) (non-autonomous case), and in the nonlinear case we refer the reader to Sanz *et al.* (2008). The use of multiregional models with different time scales to study specific instances of fish populations can be found in Charles *et al.* (2000).

We think of a population living in a habitat in which there are different environmental conditions, which we suppose indexed by an index that influences the vital rates. The population is structured in *q* age classes (corresponding to groups) and spread out in different spatial patches (subgroups) among which they may migrate. We assume that individuals in each age class *i* may migrate among *r* spatial patches. Therefore, the total number of subgroups is *N*:=*q r* and the composition of the population is then given by vector , where is the number of individuals in age class *i* living in patch *j* at time *n* and *T* denotes transposition.

Demography and migration are responsible for the transference of individuals among the different stages, and we suppose that migration is a fast process in comparison with demography. Moreover, we choose as the time step Δ_{n}=(*n*,*n*+1) for the model, the duration of each age class.

For each *i* and each environment , migration for individuals of age *i* is modelled by a matrix that, since migration is a conservative process for the total number of individuals, is column stochastic. Besides, we suppose that **P**_{i}(*ζ*) is primitive for each *i* and each *ζ*. This is the case, for example, if the fast process in each environment verifies: (i) transition from any patch to any other, in a sufficient number of steps, is allowed and (ii) individuals of at least one patch are allowed to stay in that patch. Finally, let **P**_{ζ}=diag(**P**_{1}(*ζ*),…,**P**_{q}(*ζ*)) be the matrix that models migration for the population and let the eigenvalues of **P**_{ζ} (i.e. the union of the eigenvalues of the different **P**_{i}(*ζ*)) ordered by decreasing modulus be

Now we assume that there exists *δ*<1 such that
3.1

The demographic process in each environment *ζ* is defined by the following generalized Leslie matrix
3.2where
and the coefficients and denote, respectively, the fertility and survival rates for individuals of age *i* in patch *j* in environment *ζ*.

In order to approximate the effect of migration over the time step of the model, which is much longer than its corresponding projection interval, we assume that if the population is subjected to environment *ζ* during Δ_{n} matrix **P**_{σ} operates a number *k* of times, where *k* (which we assume to be an integer) can be interpreted as the ratio between the projection intervals corresponding to demography and migration. Thus, the set of vital rates for our system in the different environments is
3.3

The pattern of environmental variation is defined by a sequence of random variables *τ*_{n},*n*=1,2,… which take values in the set of environmental states . Therefore, the proposed model, which we will refer to as the ‘original system’, consists in the following system of *N* random difference equations:
3.4with the initial condition **X**_{0}=**x**_{0} with probability 1, where **x**_{0}≥0. We will denote the variables of the original system as whenever we want to stress that they depend on the *k* chosen.

Now, following Sanz & Bravo de la Parra (2000), we make use of approximate aggregation to reduce the original system (3.4), consisting of *N* variables (microvariables) associated with the different subgroups, by an aggregated system of *q* variables (global variables), each of them associated with one group.

Let **v**_{i}(*ζ*) be the right Perron eigenvector of **P**_{i}(*ζ*), i.e. the eigenvector uniquely defined by the conditions **P**_{i}(*ζ*)**v**_{i}(*ζ*)=**v**_{i}(*ζ*),**v**_{i}(*ζ*)>0,**1**^{T}**v**_{i}(*ζ*)=1 where . **v**_{i}(*ζ*) would characterize the equilibrium structure of migration for individuals in age class *i* if the population were constantly subjected to environment *ζ*.

In this way, the matrix that characterizes the conditions for migration in environment *ζ* for group *i* is
3.5and, for the total population, we have matrix
Note that, due to equation (3.1), we have that
3.6Now we define matrices
and introduce the ‘auxiliary system’ as the stochastic model defined by
3.7with the initial condition , and let
be the set of matrices of vital rates for this system. Note that the auxiliary system can be interpreted as the result of letting migration reach equilibrium in the original system.

For each age class *i* we define a global variable corresponding to the total population of the auxiliary system with age *i*, i.e.
and then we define the vector **Y**_{n} of global variables
3.8

From equation (3.5) it follows that for all . Now we multiply both sides of equation (3.7) by **U** and obtain
which can be written in terms of the global variables exclusively, i.e. we have an aggregated system defined by
3.9with the initial condition **Y**_{0}=**Ux**_{0} where, for each *n*, is given by

The aggregated system (3.9) can be interpreted as a stochastic model in which the pattern of environmental variation coincides with that of the original system and in which the matrix of vital rates in each environment *ζ* is
in such a way that the set of matrices for the different environmental conditions is

Note that, for each *ζ*, is a classical Leslie matrix given by
where the vital rates have the form
i.e. each fertility rate *f*_{i}(*ζ*) in the aggregated system is a weighted linear combination of the fertility rates in the general system corresponding to individuals of age class *i* in environment *ζ*, being the weights of the coefficients of the equilibrium spatial distribution for migration in environment *ζ*. Something analogous holds for the survival rates.

The original multiregional model has been transformed into a reduced system in which the spatial distribution has been averaged in a certain way and the population appears structured only by age.

As we mentioned above, both the reduction procedure we presented and the relationships between the original and the reduced systems are valid for a more general kind of model subjected to the effects of environmental stochasticity. Specifically, the results are valid for a whole class of models governed by different time scales presented in Sanz & Bravo de la Parra (2000), which in particular include the following settings:

The population is divided into

*q*groups (age classes in model (3.4)) attending to any characteristic of the life cycle, and each group*i*=1,2,…,*q*is itself split into*N*_{i}subgroups (spatial patches in model (3.4)). Then, is the number of individuals in group*i*and subgroup*j*at time*n*and the population vector is where*N*:=*N*_{1}+ ⋯ +*N*_{q}.The projection interval of the model is that corresponding to the slow dynamics, on which we impose no special assumptions. Therefore, for any environmental condition

*ζ*, the slow process will be modelled by a non-negative projection matrix , which can be considered as divided into blocks**M**_{i j}(*ζ*), 1≤*i*,*j*≤*q*of dimensions*N*_{i}×*N*_{j}that characterize the rates of transference of individuals from the subgroups of group*j*to the subgroups of group*i*in environment*ζ*.Regarding the fast process,

— the fast dynamics is an internal process for each group, i.e. there is no transference of individuals from one group to a different one, and

— for each

*i*and*ζ*, matrix is a column stochastic primitive matrix; in particular, the fast process is conservative of the total number of individuals in each group.

This setting allows one to model a population structured by any characteristic of its life cycle (stage structured models, see Caswell (2001) for comments and references) and living in a multipatch environment with migration among the different patches. A particular case of this is the case of spatially heterogeneous size-structured models.

We have to stress that both the aggregation procedure and the results obtained in the next section remain valid under the more general assumptions above.

## 4. Relationship between the distribution of population size in the original system and the aggregated system

Although as we have seen that the aggregation procedure is valid independently of the pattern of environmental variation for system (3.4), in order to study the relationships between the behaviour of the original and the aggregated system, we impose conditions under which the total population size is asymptotically lognormal. More specifically we will assume,

*A*1. *The chain τ*_{n}*verifies condition (H1 ) in §2, i.e. it has a unique stationary probability distribution**that verifies* (*2.4*).

Let us now investigate the relationships between the ergodicity of the sets and . In order to do so let us introduce the following starting hypothesis:

*A*2. *The non-zero entries of matrices***M**_{ζ}*and***P**_{ζ}*are bounded away from zero and infinity, i.e. there exist positive constants ε*, *K*, *such that for all**we have*
4.1

The next result shows that, under condition *A*2, the positive entries corresponding to the original and aggregated systems are bounded away from zero and infinity and, moreover, that given a very weak condition if is an ergodic set then so are and for large enough *k*. First we introduce the following definition: a non-negative matrix **A** is said to be row-allowable if it has at least a non-zero element in each one of its rows.

## Proposition 4.1.

*(i) Assume A*2 *holds. Then there exist positive constants α,β such that for all*

*and all*

*we have*4.2

*i.e. positive entries of the matrices of the sets*

*,*

*and*

*can be bounded away from zero and infinity and, moreover, the bounds can be chosen independent of k.*

*(ii) If any product of g matrices drawn from* *is positive and matrices* **M**_{ζ} *are row-allowable, then any product of any g* + 1 *matrices drawn from* *is positive and there exists* *such that for all k≥k*_{0} *the same thing happens for any product of g* +1 *matrices drawn from**.*

## Proof.

(i) Since the entries of matrices **P**_{i}(*ζ*) are non-negative and its columns add up to 1, the entries of **P**_{i}(*ζ*) must be bounded from above for all *i* and *ζ*. In part (a) of proposition 3 in Sanz & Bravo de la Parra (2001) (where the role of the subindex *n* there is played by *ζ* here) it is proved that if *A*2 holds then the normalized right eigenvectors **v**_{i}(*ζ*) for matrices **P**_{i}(*ζ*) have entries bounded away from zero and infinity, i.e. for all , for certain *γ* and . Now, it is clear that if and are non-negative matrices and **AB**≠**0**, then min^{+}(**AB**)≥min^{+}(**A**)*m i n*^{+}(**B**) and . Therefore, from *A*2 and the definition of and it follows that there exist constants , , such that for all , ; ; ; . Now using this fact and the convergence of to uniform in , it follows that there exist and such that for all , ; . The result follows taking and .

(ii) This result is part (a) of proposition 8 in Sanz & Bravo de la Parra (2007). The fact that here the set of matrices can be infinite does not introduce any changes in the proof. ▪

As a consequence of the previous result we have,

## Corollary 4.2.

*Assume A1 and A2 holds. If the set* *is ergodic and the matrices* *are row-allowable, then the set* *is ergodic for all k≥k*_{0} *(where k*_{0} *is given by proposition 4.1). Moreover, the constants that bound above and below the positive entries of matrices of* *can be chosen independent of k.*

*Therefore, if the aggregated system meets the sufficient conditions of theorems 2.2 and 2.3 for the existence of stochastic strong ergodicity, a s.g.r. and a s.l.v., and matrices* **M**_{ζ}*,* *are row-allowable, then those sufficient conditions are also met by the auxiliary system and by the original system for k≥k*_{0}*.*

Let us introduce,

*A*3. *The set**is ergodic and the matrices***M**_{ζ}, *are row-allowable.*

Note that matrices (3.2) are row-allowable iff all are not zero and for every environment *ζ* and every patch *j* there is at least one age class *i* in which . Moreover, if for all , with **A** and **B** primitive matrices, then is ergodic.

Therefore, under conditions *A*1–*A*3, we can define the parameters

We are now interested in finding a relationship between the parameters *a*^{(k)} and *a*, on the one hand, and between *σ*^{(k)} and *σ*, on the other hand. Specifically, we would like to be able to approximate *a*^{(k)} and *σ*^{(k)} through *a* and *σ*, easier to compute when *k* is large. Theorem 9 in Sanz & Bravo de la Parra (2007) guarantees that we can carry out this approximation for the s.g.r. in the case in which the number of environmental states is finite, but the technique used in the proof does not allow one to establish a similar result for the s.l.v. Our main result guarantees, in this more general framework with possibly infinite environmental states, that both *a*^{(k)} and *σ*^{(k)} can be approximated by *a* and *σ*.

## Theorem 4.3.

*Let us assume conditions A1, A2 and A3. Then*
4.3
4.4*and so the lognormal asymptotic distribution for the original system can be approximated through the parameters a and σ that correspond to the reduced system.*

## Proof.

See the appendix. ▪

As a consequence of this result we have that the lognormal asymptotic distribution for the original system can be approximated through that corresponding to the reduced system.

## 5. Applications: exact results for the reduced system

As we previously pointed out, in most ecological models the exact computation of the s.g.r. and the s.l.v. is not feasible, due to the fact that there is not an explicit expression for the stationary distribution of age structure. Therefore, one must resort to computer simulations to approximate those parameters (Tuljapurkar & Caswell 1997). The results of the last section allow one to carry out these simulations using the simpler aggregated system. Furthermore, in some cases our reduction procedure allows transformation of a non-trivial model into a simpler reduced one in which we can compute the s.g.r. and the s.l.v. exactly. Below we introduce two such cases.

### (a) Non-structured population living in a multipatch environment

Let us consider the multiregional model of §3 in the case in which the population is not structured by age (*q*=1), is distributed among *r* spatial patches and there are a finite number *s* of environmental conditions that vary according to a Markov chain *τ*_{n} with values in the set

The population vector is , where is the population in patch *i*. Moreover we have that demography in environment *ζ* is modelled by
where *m*_{i}(*ζ*), *i*=1,…,*r*; *ζ*=1,…,*s* is the growth rate in patch *i* and environment *ζ*. Migration among the different patches in each environment *ζ* is modelled by a column-stochastic primitive matrix . If we assume that all the growth rates are positive *m*_{i}(*ζ*)>0, *i*=1,…,*r*; *ζ*=1,…,*s*, then condition *A*3 of §4 is met.

The Markov chain *τ*_{n} is defined by a certain vector of initial probabilities and a row-stochastic matrix of transition probabilities such that , 1≤*u*,*v*≤*s*. We assume that *τ*_{n} is irreducible and ergodic, i.e. that **Q** is a primitive matrix, and so *τ*_{n} has a stationary probability distribution, which we denote ** π**=(

*π*

_{1},

*π*

_{2},…,

*π*

_{s})

^{T}, and condition

*A*1 of §4 is met.

Let **v**_{ζ}>0 be the equilibrium distribution for migration in environment *ζ* normalized so that . Then the reduced system (3.9) takes the form of a scalar stochastic system,
5.1where is the total population and

Clearly the set is ergodic and matrices **M**_{ζ} are row-allowable, so both the reduced system (5.1) and the original system (for large enough *k*) meet the sufficient hypotheses for the existence of an asymptotic lognormal distribution. We can explicitly compute the s.g.r. and the s.l.v. that characterize the asymptotic distribution for system (5.1), the reason being that in this case the normalized population *y*_{n}/∥*y*_{n}∥ is trivial (equal to 1 with probability 1), and so the stationary distribution of the chain (*τ*_{n},*y*_{n}/∥*y*_{n}∥) is simply that of *τ*_{n}.

In order to proceed we assume that **Q** is diagonalizable so that **Q**=**RDR**^{−1}, with **D**=diag(1,*λ*_{2},…,*λ*_{s}) and **R**=(**1**,**r**_{2},…,**r**_{s}), where 1,*λ*_{2},…,*λ*_{s} are the eigenvalues of **Q** and **1**,**r**_{2},…,**r**_{s} corresponding right eigenvectors, **Qr**_{j}=*λ*_{j}**r**_{j}, 2≤*j*≤*s* and **Q1**=**1**.

Let us define the stochastic variables . Then, according to theorem 2.3 the s.g.r. of system (5.1) is 5.2

Let us now define the variables (with zero mean with respect to the stationary distribution) , which take values , 1≤*ζ*≤*s*. Then the variables *c*_{r} defined in theorem 2.3 are

The *r*-step transition probabilities are defined by matrix **Q**^{r} so that . A straightforward calculation shows that
where the column vectors **z** and **z**_{π} are defined as

Taking into account that **R**=(**1**,**r**_{2},**r**_{3},…,**r**_{s}) and the fact that we have
so that the s.l.v. *σ*^{2} of the reduced system can be written as
5.3

Summing up, if the separation of time scales is large enough, the aggregation procedure allows one to approximate the s.g.r. and the s.l.v. of the original system, whose explicit computation is not feasible, through the explicit expressions (5.2) and (5.3).

For example, expression (5.2) makes it possible to study how changes in the parameters affect the long run growth (*a*>0) or extinction (*a*<0) of the population.

### (b) A population with juveniles and adults living in a two-patch environment

#### (i) Model and reduction

Let us consider the setting of §3 in the case *q*=2, and *r*=2, i.e. the population has two age classes, juveniles and adults, and is living in a habitat with two patches. Then the population vector in the original system is and
where *S*^{j}(*ζ*) stands for the survival rate of juveniles in patch *j*, *p*_{i}(*ζ*) denotes the migration rates in environment *ζ* of individuals of age *i* from patch 1 to patch 2 and *t*_{i}(*ζ*) has an analogous meaning for migration from patch 2 to patch 1. If we assume that, for all , *p*_{1}(*ζ*), *p*_{2}(*ζ*), *t*_{1}(*ζ*) and *t*_{2}(*ζ*) are different from 0 or 1 then matrices **P**_{1}(*ζ*) and **P**_{2}(*ζ*) are primitive for all , and so we can reduce system (3.4).

Vectors **v**_{i}(*ζ*) have the form **v**_{i}(*ζ*)=(*t*_{i}(*ζ*)/(*p*_{i}(*ζ*)+*t*_{i}(*ζ*)),*p*_{i}(*ζ*)/(*p*_{i}(*ζ*)+*t*_{i}(*ζ*)))^{T}, *i*=1,2, and therefore the reduced system (3.9) has the form
5.4with and

We proceed to analyse the reduced system (5.4) in three cases regarding the kind of environmental variation for which we can obtain exact results.

#### (ii) Case 1

Let us assume that the parameters that govern demography are independent of the environment, i.e. , *S*^{j}(*ζ*)=*S*^{j}, where , *S*^{j}<1. Moreover, let
where *p*_{i},*t*_{i}∈(0,1) are non-random and *η*_{n} is a homogeneous Markov chain with values in [*ε*,1], where 0<*ε*<1, that meets hypothesis *A*1 in §4 (with *η*_{n} playing the role of *τ*_{n}). Note that hypothesis *A*2 is also met.

This model can be interpreted by saying that reproduction and survival are not affected by the environment and that migration rates depend on the environment through a common multiplicative random variable. A good/bad environment increases/decreases all migration rates by the same factor.

Then we have
where
i.e. the reduced system is the deterministic system . Since is a primitive matrix hypothesis *A*3 of §4 is met and so the results of that section hold.

Note that the s.g.r. and the s.l.v. of the reduced system are

#### (iii) Case 2

Let us assume that migration is deterministic, i.e. *p*_{i}(*ζ*)=*p*_{i}, *t*_{i}(*ζ*)=*t*_{i} (with *p*_{i},*t*_{i}∈(0,1)) and that
where we assume that are non-random, positive and *s*^{j}<1. *η*_{τn} is a scalar homogeneous Markov chain, which we assume irreducible and ergodic that takes values in the set {*δ*_{1},*δ*_{2},…,*δ*_{s}}, where *δ*_{i}∈(0,1]. We denote its stationary probability distribution as ** π**.

Note that this case can be interpreted by saying that all fertility and survival rates depend on the environment through a common multiplicative random variable. A good/bad environment increases/decreases all vital rates by the same factor.

Then the aggregated system is where , **A** being the non-random matrix

The variables *η*_{n} are bounded away from zero and infinity and matrix **A** is primitive, so is an ergodic set, and hypotheses *A*1–*A*3 of §4 are met.

Now we proceed to calculate the s.g.r. and the s.l.v. for the reduced system (5.4). In order to do so we will use the fact that the evolution of the population structure is a deterministic one and only the total number of individuals evolves in a non-deterministic way, and we will reduce the problem to the scalar case, which we have already solved in a previous example. Specifically, let
where *G*_{τn}:=*η*_{τn}*η*_{τn−1}⋯*η*_{τ1}. Let *μ*_{1}≥*μ*_{2} be the eigenvalues of **A** and let **u**_{1},**u**_{2} be associated right eigenvectors. Since **A** is primitive then *μ*_{1}>0 and *μ*_{1}>*μ*_{2} and so **A** can be diagonalized as **A**=**VDV**^{−1} with **D**=diag(*μ*_{1},*μ*_{2}), **V=**(**u**_{1},**u**_{2}). Defining new stochastic variables , we have that **S**_{n}=*G*_{τn}**D**^{n}**S**_{0} and, coming back to the **Y**_{n},
As we mentioned before, the population structure **Y**_{n}/∥**Y**_{n}∥ has a deterministic evolution and converges to **u**_{1}/∥**u**_{1}∥. The logarithm of the total population ∥**Y**_{n}∥ is
where the second term in the right-hand side is deterministic. Now we apply to the results that we obtained for the scalar system (5.1) so that , where and *σ*^{2} is given by (5.3) with . Taking into account that
we finally conclude that where .

#### (iv) Case 3

Let us assume that migration is deterministic, i.e. *p*_{i}(*ζ*)=*p*_{i}, *t*_{i}(*ζ*)=*t*_{i} (with *p*_{i},*t*_{i}∈(0,1)) and that
where we assume that are non-random, positive and *s*^{j}<1. In this model survival is not affected by the environment, whereas fertility rates are. A good/bad environment increases/decreases fertility rates by the same multiplicative constant leaving the survival rates unaltered.

The matrices of vital rates in the reduced system are
and we assume that *η*_{n} is a sequence of i.i.d. random variables such that 1/*η*_{n} has a probability density function *g*(*x*) given by
where *ε* and *K* are arbitrary positive constants. Were it not for the fact that *g*(*x*) has been made 0 for *x* values smaller than *ε* and larger than *K*, *g*(*x*) would be the probability density function of a gamma distribution used in Tuljapurkar (1990, p. 66) in the context of a two-age class model for which an exact distribution of age structure can be obtained. We use *ε* and *C* as cut-offs that make the support of *g* compact so that the elements of the **M**_{ζ} (and with them those of ) are bounded away from zero and infinity, and so the set is ergodic. Since *ε* and *K* can be made, respectively, as small and large as we want, results in Tuljapurkar (1990) are, from a practical point of view, as good an approximation to our case as desired.

Thus, following Tuljapurkar (1990), we have that
where *C* is the probability density function defined by
with *κ* being a normalizing constant.

## Acknowledgements

This work was supported by Proyecto de Investigación MTM2008-06462 (MCYT).

## Appendix

In the first place we will establish some results for the generic model (2.1) that will be useful further on. If we define
we can write (2.1) as **Z**_{n}=**T**_{n}**z**_{0}. Note that each realization of **T**_{n} is a product of *n* matrices belonging to the set that we assume to be ergodic. The literature shows,

### Lemma A.1.

(Equation in Tuljapurkar & Orzack (1980), lemma 3 in Furstenberg & Kesten (1960).) *If* *is an ergodic set with constants g*, *α and β* (*see definition 2.1*) *then, for any realization of the process* **T**_{n}*, we have that, for 1 ≤ i,j,l ≤ N,*
A1

*where η <*1

*is a constant that can be written as*A2

*and*

*is a constant that depends only on N,g,*

*α*and*β*.Note that, if a set of matrices **A**_{ζ} is ergodic with constants *g*,*α* and *β*, so is the set of transposed matrices so that equation (A1) is also valid if we interchange the subindexes of all the matrices **T**_{n+1} that appear in the inequality. Therefore, from now on we will refer to equation (A1) to denote any one of these inequalities.

Using equation (A1) we will prove the so-called ‘weak ergodic theorem in demography’ (Cohen 1977, 1979; Tuljapurkar & Orzack 1980; Seneta 1981) for products of matrices belonging to an ergodic set, without resorting to Hilbert’s pseudometric. Moreover, we show that the speed with which vectors of population structure corresponding to two initial populations converge to each other depends only on the constants of the ergodic set. This fact will be necessary later on.

### Theorem A.2.

*Let us assume that the set* *is ergodic with constants g, α and β. Given δ>0 there exists n*

_{0}

*such, that for all n≥n*

_{0}

*, for all initial states*

*and for any realization of the process we have*A3

*where, given δ, the value of n*

_{0}

*depends only on N and on the constants of the ergodic set.*

### Proof.

Given any **T**_{n},**w**_{0} and we define for all 0≤*i*≤*N* and we have
A4with the definitions
A5
A6

Now we will prove that reaches its maximum when **h**_{0} and are vectors belonging to the canonical basis of and therefore, according to (A4)–(A6), the first side of A3 also reaches its maximum when belong to the canonical basis.

Denoting by **t**_{1},…,**t**_{N} the columns of and taking into account that we have, whenever , that

The norm in the right-hand side of this equation does not depend on the value of and only on the direction of the unitary vector with components with 2≤*j*≤*N*, from where it follows that necessarily the maximum of is reached either when or when . Reasoning along the same lines for the rest of the components we conclude that the maximum of is reached when , i.e. **h**_{0} and must be two vectors of the canonical basis and the same happens to **w**_{0} and .

Now let **w**_{0}=**e**^{j} and be two canonical vectors . Then using equation (A1) we have that for all
where *η* is given by equation (A2). Given *δ*>0 the result follows taking *n*_{0} such that *δ*=*N D*_{2}*η*^{n0−1}. □

Now let us consider the process and its transition probabilities *G*_{n} defined by equation (2.3). We will denote *G*:=*G*_{1}. Under hypotheses (H1) and (H2), theorem 2.2 guarantees that the transition probabilities *G*_{n} converge to a probability distribution *F* uniformly in the initial state of the chain. The proof of theorem 2.2 hinges on a series of lemmas, among which we have the following ones, which will be relevant thereafter.

### Lemma A.3.

(Lemma 5 in Cohen 1977.) *Let us assume hypotheses (H1) and (H2). Given ε>0 and* *such that for all* *such that* *and for all* *we have*

### Lemma A.4.

(Lemma 6 in Cohen 1977.) *Let us assume hypotheses (H1) and (H2). Given ε>* 0 *and* *such that* *and* *we have*
*where* *are obtained from* *through equation (2.2 ).*

Once we have established some results regarding general matrix models with environmental stochasticity, we will consider the original system (3.4) and the auxiliary system (3.7), so that with and we will introduce several results that relate these two systems.

As a direct consequence of proposition 4.1 we have, using theorem A.2,

### Proposition A.5.

*Let us assume A2 and A3. Then, given δ>0 there exists* *such that ∀n≥n*_{0}*,k≥k*_{0} *(k*_{0} *defined in proposition 4.1 ), ω∈Ω and for any initial states*
*Moreover, given δ the value of n*_{0} *depends only on N and on the constants g, α and β (defined in proposition 4.1).*

Let us also consider the bivariate chains and , where and and let their transition probabilities be defined through with and . We will denote and .

If we assume hypotheses *A*1–*A*3, then, as a consequence of corollary 4.2, for all *k*≥*k*_{0} the chains and meet the hypotheses of theorem 2.2, and so there exist stationary probability distributions and *F*^{(k)} for and , respectively. Specifically we have,

### Lemma A.6.

*Given hypotheses A*1–*A*3 *we have*,

*Given ε*>0,*and*,,*such that**For all k*≥*k*_{0}*we have that, given ε*> 0,*and**there exists*(*possibly depending on k*)*such that*∀*n*≥*n*_{0}, ,

Our goal now is to prove that *F*^{(k)} converges to when . In the first place we establish a number of relationships between the systems (3.4) and (3.7) when *n* is fixed and .

### Lemma A.7.

*Let us assume hypothesis A*1–*A*3. *Then*,

*Given**uniformly in ω ∈ Ω*.*Given**uniformly in**and ω ∈ Ω*.*For all**and*,*uniformly in**and*

### Proof.

(i) It suffices to prove the result for *n*=2. Let ∥*∥ be any matrix norm induced by a vector norm. Let us denote and . Then
where we have used (3.6) and the bounds (4.1) and (4.2).

(ii) Let ∥*∥ denote, depending on the context, the 1-vector norm and its associated matrix norm in . Since is an ergodic set, its matrices cannot have any zero columns, so that, using equation (4.2), and moreover . Therefore, A7

Now we can write where we have taken into account equation (A7). The result now follows from (i). (iii) Immediate taking into account (ii), the definition of and and the fact that convergence with probability 1 implies convergence in distribution. ▪

The following two lemmas extend lemmas A.2 and A.3 (lemmas 5 and 6 in Cohen (1977)) to the system (3.4), the results being uniform in *k*. Specifically,

### Lemma A.8.

*Let us assume hypothesis A1–A3. Given ε>0 and* *, ∃δ,k*_{1}*>0 such that if* *verify that* *and* *,* *, then*

### Proof.

and the result follows applying part (iii) of lemma A.7 with *n*=1 to the first and third terms in the right-hand side and lemma A.3 to the second term. ▪

### Lemma A.9.

*Let us assume hypotheses A1–A3. Given ε>0 and* *,* *, such that* *and* *, we have*
A8

### Proof.

From lemma A.8 we have that such that , (A8) holds. Moreover, from proposition A.5 we have that given *δ*>0, such that , *l*≥*l*_{0}. The result follows by taking . ▪

From lemmas A.8 and A.9, and proceeding analogously to the proof of theorems 2 and 3 of Cohen (1977) one can prove the following two results, equivalent to theorems 2 and 3 of Cohen (1977) for the chains uniformly in *k*:

### Proposition A.10.

*Let us assume hypotheses A1–A3. Given ε>0,* *and* *,* *and k*_{2} *given by lemma A.9 such that ∀k≥k*_{2}*,∀n≥n*_{0}*,* *one has*

### Proposition A.11.

*Let us assume hypotheses A1–A3. For all k≥k*_{2} *with k*_{2} *given by lemma A.9 there exists a stationary probability distribution F*^{(k)} *such that given ε>0,* *and* *,* *such that* *,*

As a consequence of the previous results, we finally have the result we were aiming at.

### Theorem A.12.

*Let us assume hypotheses A1–A3. Given* *and* *we have*

### Proof.

Taking any , , , we can write Now the result follows if we first take and use part (iii) of lemma A.7, and then take and use part (ii) of lemma A.6 and proposition A.11. ▪

### Proof of theorem 4.3.

Applying theorem 2.3, the hypotheses guarantee that we can define a s.g.r. and a s.l.v. for the auxiliary system (3.7) through and . The definition of the global variables in (3.8) implies that and therefore we have that and . Consequently, all we have to do is to prove that and .

Applying theorem 2.3 to the auxiliary and the original (for *k*≥*k*_{0}) systems, we know that
A9and
A10Moreover, we can define, for *r*=1,2,…
where **x**_{0}>0 and we have
A11and
A12Denoting the expression in brackets in (A9) and the expression in brackets in (A10) we have
A13
Now we will show that both terms in the right-hand side of (A13) tend to zero when .

In the first place we will prove that . Let us denote , , , so that and . Using the bounds (4.2) we have that, for all *ω*∈Ω,
A14

Moreover, from part (i) of lemma A.7 it follows that uniformly in *ω*∈Ω and uniformly in *ω*∈Ω. Now we make use of the result by which if two sequences of functions *M*^{(k)}(*ω*) and *S*^{(k)}(*ω*) converge when to *M*(*ω*) and *S*(*ω*) uniformly in *ω*∈Ω and, moreover, for all *ω* |*M*(*ω*)|≤*C* and |*S*(*ω*)|≥*δ* for certain positive constants *C* and *δ*, then uniformly on *ω*∈Ω (*). As a consequence we have that uniformly on *ω*∈Ω. Now, using (A14) we have that there exist positive constants 0<*ε*<*K* such that for all *ω*∈Ω. Since the logarithm is a uniform continuous function in [*ε*,*K*] it follows, using (*), that as we wanted to show.

From the discussion above it follows that for a certain constant *L*>0. Now using this fact, theorems 5 and 25.8 in Billingsley (1986), it follows that the second term in the right-hand side of equation (A13) tends to zero when and so (4.3) holds.

Now let be fixed. Proceeding in the same way as above it follows that A15

We will show that the series in the right-hand side of equation (A11) converges uniformly in *k* for *k*≥*k*_{0}. Expressions 3.26 and 3.27 in Furstenberg & Kesten (1960) applied to the auxiliary system show that with , where *C*>1 is defined in equation (3.2) of that reference and *λ*_{1}<1 in the hypothesis *A I I* (p. 464). Therefore, the series that defines converges.

Reasoning in the same way for the original system we have that for all *k*≥*k*_{0}
where *η* is given in lemma A.1 when applied to system (3.4) (recall that, owing to proposition 4.1, the constants of the ergodic set can be chosen independent of *k* and therefore *η* is independent of *k*) and *ρ* in (2.4). Therefore, the series in the right-hand side of (A12) converges uniformly in *k* for *k*≥*k*_{0}. Therefore taking the limit when in (A12), interchanging the series with the limit due to the uniform convergence of the series, and taking into account (A15) we obtain (4.4) as desired. ▪

## Footnotes

One contribution of 17 to a Theme Issue ‘From biological and clinical experiments to mathematical models’.

- © 2009 The Royal Society