Approximate reduction of linear population models governed by stochastic differential equations: application to multiregional models

ABSTRACT

In this work we develop approximate aggregation techniques in the context of slow-fast linear population models governed by stochastic differential equations and apply the results to the treatment of populations with spatial heterogeneity. Approximate aggregation techniques allow one to transform a complex system involving many coupled variables and in which there are processes with different time scales, by a simpler reduced model with a fewer number of ‘global’ variables, in such a way that the dynamics of the former can be approximated by that of the latter. In our model we contemplate a linear fast deterministic process together with a linear slow process in which the parameters are affected by additive noise, and give conditions for the solutions corresponding to positive initial conditions to remain positive for all times. By letting the fast process reach equilibrium we build a reduced system with a lesser number of variables, and provide results relating the asymptotic behaviour of the first- and second-order moments of the population vector for the original and the reduced system. The general technique is illustrated by analysing a multiregional stochastic system in which dispersal is deterministic and the rate growth of the populations in each patch is affected by additive noise.

KEYWORDS:

• Time scales
• stochastic differential equations
• approximate aggregation
• multiregional models

2010 MATHEMATICS SUBJECT CLASSIFICATION:

• 92D25
• 92D40
• 00A69

1. Introduction

Nature offers many examples of systems with an inherent complexity whose study leads to mathematical models with a large number of state variables whose analytical study is, in most cases, not feasible. In order to be able to extract important information about the behaviour of some of these complex models, one can resort to ‘approximate aggregation methods’. These are mathematical techniques, which are usually applied in systems governed by processes with different time scales, in which appropriate approximations are introduced in order to transform the system under consideration into a reduced system with a lower number of variables, called ‘global variables’. In this way, the behaviour of the original system can be approximated, but not known with exactitude, in terms of the knowledge of the behaviour of the reduced system.

Approximate aggregation techniques in population dynamics have been widely studied in the context of deterministic systems with different time scales both in continuous and discrete time (see the review in [2] and the references therein), as well as in discrete stochastic models incorporating either environmental [19] or demographic stochasticity [21].

Stochastic differential equations (SDEs) can be thought of as resulting from the introduction of environmental stochasticity in the coefficients of deterministic ODEs. In spite of the difficulty of their analytical treatment, they have become popular as population modeling tools (see for example [6, 8] for the scalar case and [13, 24] for applications to competition models and epidemiology respectively).

The aim of this work is to formulate an approximate aggregation technique valid to reduce fast-slow linear population models governed by SDEs, to give sufficient conditions for the solutions of the models to be positive and to relate the asymptotic behaviour of the first- and second-order moments of the population vector for the original and the reduced system. The original model is built by considering a fast deterministic process that converges to an equilibrium, together with a slow process in which the parameters are affected by additive noise. In the resulting system the quotient between the rate of diffusion squared and the speed of drift is $O\left(\epsilon \right)$, $\epsilon \to 0$, where ϵ is a measure of the difference of time scales between the fast and the slow process. This is in contrast with most previous works in the field [5, 10, 11, 16] which are valid in other situations. An exception is [14] which covers our case of interest but which only deals with the relationships between the original and the reduced systems in finite time intervals.

The existing literature on the field of the approximate reduction of SDEs has been developed mainly in the context of physics and control theory. Our approach, although general in nature, is aimed to population dynamics applications. Specifically, we will employ it to study a multiregional model consisting of a single population living in a multipatch environment in such a way that dispersal is deterministic and the growth rate of the population in each patch is affected by additive noise. Models with spatial heterogeneity have been widely studied through aggregation techniques (see [3] for the continuous time case and [2] for discrete time), normally making use of the fact that in many practical situations the dispersal of individuals amongst the different patches is fast with respect to other processes like demography or competition (an exception to this is [20]). In this work we contemplate a situation that has not been dealt with so far: dispersal between some spatial patches can be fast with respect to demography whereas migration between some other patches can happen at the same time scale of demography. For example we can think of a population of birds living in different islands in such a way that inter-island movements happen at the scale of demography but, in comparison, intra-island migrations are fast.

The manuscript is organized as follows: in Section 2 we present the general formulation of linear SDE models and, in order to be able to use them in population dynamics applications, we give sufficient conditions that guarantee that if initial conditions are positive, the solution to such a system remains positive for all times. Section 3 presents the general formulation of a linear two time scale population model with stochasticity affecting the slow process. Section 4 carries out the reduction of the system, which is accomplished by letting the fast process reach equilibrium and defining and adequate set of new variables. In Section 5 we apply the general technique to the case of a stochastic multiregional model, with special attention to the case in which all migrations are fast with respect to demography, since it is the most frequent situation in practice and moreover in that case the reduced system is scalar. Section 6 presents a result that allows one to relate the asymptotic behaviour of the first and second statistical moments of the population vector for the original and the reduced system, and applies it to study the asymptotic behaviour of the multiregional models of Section 5. A brief discussion of results and the Appendix with mathematical proofs complete the manuscript.

2. Linear population models with two time scales

Throughout the paper we assume we are working in a complete probability space $(\mathrm{\Omega },\mathcal{F},\{{\mathcal{F}}_t{\}}_{t\ge 0},\mathbb{P})$ where the filtration $\{{\mathcal{F}}_t{\}}_{t\ge 0}$ satisfies the usual conditions [18]. Let us consider a structured population modelled by an linear autonomous homogeneous stochastic differential equation. The model has the form (1) $\mathrm{d}\mathbf{X}\left(t\right)=\mathbf{A}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0,$(1) where $\mathbf{X}\left(t\right)=\left(X^1\right(t),X^2(t),\dots ,X^d(t\left)\right)\in {\mathbb{R}}^d$ is the population vector, $\mathbf{A}\mathbf{=}\mathbf{[}{a}_{ij}]\in {\mathbb{R}}^{d\times d}$, ${\mathbf{B}}^h=\mathbf{[}{b}_{ij}^h]\in {\mathbb{R}}^{d\times d},$ $h=1,2,\dots ,m$ and $\mathbf{W}\left(t\right)=\left(W^1\right(t),W^2(t),\dots ,W^m(t){)}^{\mathrm{T}}$ is a m-dimensional standard Wiener process defined in the previous probability space. Moreover, we assume that ${\mathbf{X}}_0\in {\mathbb{R}}^d$ is a non-random vector. Models of the kind (1(1) $\mathrm{d}\mathbf{X}\left(t\right)=\mathbf{A}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0,$(1) ) are obtained from a linear deterministic model (2) $\frac{\mathrm{d}\mathbf{X}\left(t\right)}{\mathrm{d}t}=\mathbf{A}\mathbf{X}\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0$(2) if we add noise to the population vital rates $a_{ij}$. Indeed, system (1(1) $\mathrm{d}\mathbf{X}\left(t\right)=\mathbf{A}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0,$(1) ) can be written in the form $\mathrm{d}{X}^i\left(t\right)=\stackrel{d}{\sum _{j=1}}(a_{ij}\mathrm{d}t+b_{ij}^1\mathrm{d}{W}^1(t)+\cdots +b_{ij}^m\mathrm{d}{W}^m(t\left)\right)X^j\left(t\right),i=1,\dots ,d$ from where we see that each coefficient $b_{ij}^h$ characterizes the intensity of the noise $\mathrm{d}{W}^h$ affecting $a_{ij}$. Note that the case in which the noises $\mathrm{d}{W}^j\left(t\right)$ are correlated can be reduced to this setting through an appropriate transformation [1, p. 126]. Systems of the kind (1(1) $\mathrm{d}\mathbf{X}\left(t\right)=\mathbf{A}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0,$(1) ) can be interpreted in the sense of Ito or in the sense of Stratonovich, and the results obtained in the two cases differ [23]. However, one can choose any of the two interpretations as long as one defines appropriately the parameters of the model [7]. In this work we will make use of the Ito interpretation. It is well known [1, p. 126] that in the previous conditions there exists a unique solution to (1(1) $\mathrm{d}\mathbf{X}\left(t\right)=\mathbf{A}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0,$(1) ) which is continuous with probability one.

In order to be a valid model for a population, the solution of (1(1) $\mathrm{d}\mathbf{X}\left(t\right)=\mathbf{A}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0,$(1) ) for any positive initial condition must remain non-negative for all times, so we turn our attention to this property. Given a vector or matrix $\mathbf{x}$ we will write $\mathbf{x}\ge \mathbf{0}$ (resp. $\mathbf{x}>\mathbf{0}$) to denote that all the components of $\mathbf{x}$ are non-negative (resp. positive). Let us recall that a square matrix $\mathbf{A}$ is said to be a Metzler matrix or a essentially non-negative matrix when it has non-negative off-diagonal entries, that is, $a_{ij}\ge 0$ for $i\ne j$ [22]. It is well known that if $\mathbf{A}$ is a Metzler matrix then solutions of the deterministic system (2(2) $\frac{\mathrm{d}\mathbf{X}\left(t\right)}{\mathrm{d}t}=\mathbf{A}\mathbf{X}\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0$(2) ) meet the desired property. However, in model (1(1) $\mathrm{d}\mathbf{X}\left(t\right)=\mathbf{A}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0,$(1) ) additional requirements are needed in order to ensure positivity of solutions. The next result gives sufficient conditions:

Theorem 2.1

If matrix $\mathbf{A}$ is a Metzler matrix and matrices ${\mathbf{B}}^1,\dots ,{\mathbf{B}}^k$ are diagonal, then given $\mathbf{X}\left(0\right)={\mathbf{X}}_0>\mathbf{0}\mathbf{,}$ the solution of (1(1) $\mathrm{d}\mathbf{X}\left(t\right)=\mathbf{A}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0,$(1) ) verifies that $\mathbb{P}\left(\mathbf{X}\right(t)>\mathbf{0}$ for all $t\ge 0)=1$.

Proof.

See Appendix.

3. A model with two time scales

We suppose a stage-structured population in which individuals are classified into stages or groups attending to any characteristic of the life cycle. Moreover, each of these groups is divided into several subgroups that can correspond to different spatial patches, different individual activities or any other characteristic that could change the life cycle parameters. The model is therefore general in the sense that we do not state in detail the nature of the population or the subpopulations.

We consider the population being subdivided in q populations (or groups). Each group is subdivided in subpopulations (subgroups) in such a way that for each $i=1,2,\dots ,q$, group i has $N_i$ subgroups. Therefore, the total number of subgroups is $N:=N_1+N_2+\cdots +N_q$. We denote the fast time as τ, while the slow time will be denoted by t. In this way we have $t=\epsilon \tau$ where $\epsilon <<1$ is a small number that represents the ratio between the slow and the fast times.

Le $x^{ij}\left(\tau \right)$ be the density of subpopulation j of population i at time τ, with $i=1,2,\dots ,q$ and $j=1,2,\dots ,N_i$. In order to describe the population of group i at time τ we will use vector ${\mathbf{x}}^i\left(\tau \right)=\left(x^{i1}\right(\tau ),x^{i2}(\tau ),\dots ,x^{i{N}_i}(\tau ){)}^{\mathrm{T}}\in {\mathbb{R}}^{N_i},$ $i=1,2,\dots ,q,$ where T denotes transposition. The composition of the total population is then given by vector $\mathbf{X}\left(\tau \right)=\left({\mathbf{x}}^1\right(\tau {)}^{\mathrm{T}},{\mathbf{x}}^2(\tau {)}^{\mathrm{T}},\dots ,{\mathbf{x}}^1(\tau {)}^{\mathrm{T}}{)}^{\mathrm{T}}\in {\mathbb{R}}^N$.

In the evolution of the population we will consider two linear processes whose corresponding characteristic time scales, are very different from each other. In order to include in our model both time scales we will model these two processes, to which we will refer as the fast and the slow dynamics, by two different matrices.

In principle, we will make no special assumptions regarding the characteristics of the slow dynamics other than linearity and restrictions to guarantee positivity of solutions to the slow process. Thus, we will assume that the parameters of the slow process are defined by (3) $\mathbf{S}+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\frac{\mathrm{d}{W}^h\left(t\right)}{\mathrm{d}t},$(3) where: S.1. Matrix $\mathbf{S}\in {\mathbb{R}}^{N\times N}$ models the deterministic part of the slow process and has non-negative off-diagonal entries, that is, $\mathbf{S}$ is a Metzler matrix. We consider $\mathbf{S}$ divided into blocks ${\mathbf{S}}_{ij},1\le i,j\le q$ in such a way that $\mathbf{S}=\left[\begin{array}{ccc}{\mathbf{S}}_{11}& \cdots & {\mathbf{S}}_{1q}\\ ⋮& \ddots & ⋮\\ {\mathbf{S}}_{q1}& \cdots & {\mathbf{S}}_{qq}\end{array}\right].$ Each block ${\mathbf{S}}_{ij}=\left[s_{ij}^{\alpha \beta }\right]$ has dimensions $N_i\times N_j$ and characterizes the rates of transference of individuals from the subgroups of group j to the subgroups of group i. More specifically, for each $\alpha =1,2,\dots ,N_i$ and each $\beta =1,2,\dots ,N_j$, entry $s_{ij}^{\alpha \beta }$ represents the rate of transference of individuals from subgroup β of group j to subgroup α of group i.

S.2. $W^1,\dots ,W^k$ are independent standard Wiener processes (so that by $\mathrm{d}{W}^1\left(t\right)/\mathrm{d}t,\dots ,\mathrm{d}{W}^k\left(t\right)/\mathrm{d}t$ we denote the associated weakly defined gaussian white noises) and, for each $h=1,\dots ,k$, matrix ${\mathbf{B}}^h\in {\mathbb{R}}^{N\times N}$ models the contribution of noise $\mathrm{d}{W}^h\left(t\right)$ to the dynamics of the slow process. In order to be able to guarantee positivity of solutions (Theorem 2.1), we assume that matrices ${\mathbf{B}}^h$ are diagonal, that is, ${\mathbf{B}}^h=\mathrm{diag}({\mathbf{B}}_1^h,\dots ,{\mathbf{B}}_q^h),h=1,\dots ,k\mathrm{with}{\mathbf{B}}_i^h=\mathrm{diag}(b_{i1}^h,\dots ,b_{i{N}_i}^h),i=1,\dots ,q.$ Note from Equation (3(3) $\mathbf{S}+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\frac{\mathrm{d}{W}^h\left(t\right)}{\mathrm{d}t},$(3) ) that $b_{i\alpha }^h$ models the intensity of the noise $\mathrm{d}{W}^h\left(t\right)$ affecting coefficient $s_{ii}^{\alpha \alpha }$.

In this context, we say that an eigenvalue λ of a certain square matrix $\mathbf{A}$ is strictly dominant when the real part of λ is strictly larger than the real part of the rest of the eigenvalues of $\mathbf{A}$.

As far as the behaviour of the fast dynamics is concerned, we will make the following three assumptions:

F.1. The fast process is deterministic.

F.2. The fast dynamics is an internal process for each group, that is, there is no transference of individuals from one group to a different one. Therefore, for each $i=1,\dots ,q$, the fast dynamics of group i will be represented by a Metzler matrix ${\mathbf{F}}_i$ of dimensions $N_i\times N_i$. We will assume that ${\mathbf{F}}_i$ is irreducible in the sense that there exists r>0 such that $r\mathbf{I}+{\mathbf{F}}_i$ is a primitive non-negative matrix [22]. Therefore, the matrix that governs the fast dynamics for the whole population is (4) $\mathbf{F}=\mathrm{diag}({\mathbf{F}}_1,{\mathbf{F}}_2,\dots ,{\mathbf{F}}_q).$(4)

F.3. The fast process in each group has a non-trivial equilibrium point which is asymptotically stable. More specifically, for each $i=1,\dots ,q$ matrix ${\mathbf{F}}_i$ has eigenvalue 0 and it is strictly dominant so that the rest of the eigenvalues of ${\mathbf{F}}_i$ have negative real parts. Since ${\mathbf{F}}_i$ is a Metzler irreducible matrix, eigenvalue 0 is simple [22, Theorem 2.6] and moreover, there exist positive right and left eigenvectors ${\mathbf{v}}_i$ and ${\mathbf{u}}_i$ of ${\mathbf{F}}_i$ associated to eigenvalue 0 for which we choose the following normalization conditions (5) $\begin{array}{rl}{\mathbf{F}}_i{\mathbf{v}}_i& =\mathbf{0};{\mathbf{u}}_i^{\mathrm{T}}{\mathbf{F}}_i={\mathbf{0}}^{\mathrm{T}}\\ \parallel {\mathbf{v}}_i{\parallel }_1& =1;{\mathbf{u}}_i^{\mathrm{T}}{\mathbf{v}}_i=1,\end{array}$(5) where $\parallel \ast {\parallel }_1$ denotes the 1-vector norm.

In order to incorporate both time scales in our model, we will make use of parameter $\epsilon =t/\tau$. The model, to which we will refer as original system, has the following form in the slow time (6) $\mathrm{d}\mathbf{X}\left(t\right)=\frac{1}{\epsilon }\mathbf{F}\mathbf{X}\left(t\right)\mathrm{d}t+\mathbf{S}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right).$(6) Alternatively, using the fast time τ we have (7) $\mathrm{d}\mathbf{X}\left(\tau \right)=\mathbf{F}\mathbf{X}\left(\tau \right)\mathrm{d}\tau +\epsilon \mathbf{S}\mathbf{X}\left(\tau \right)\mathrm{d}\tau +\sqrt{\epsilon }\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(\tau \right)\mathrm{d}{W}^h\left(\tau \right),$(7) where we have made an abuse of notation and kept the same notation for $\mathbf{X}$ and the $W^h$ in the new time, and we have used [1, p. 47] that if $\mathbf{W}\left(t\right)$ is a standard Wiener process then so is $\sqrt{\epsilon }\mathbf{W}\left(t/\epsilon \right)$. We stress that the quotient between the rate of diffusion squared and the speed of drift is $O\left(\epsilon \right)$, $\epsilon \to 0$, which is in contrast with most approaches in the analysis of two-time scales systems [5, 10, 11, 16] which are valid in other situations.

4. Approximate reduction of the model

In order to reduce the original system (7(7) $\mathrm{d}\mathbf{X}\left(\tau \right)=\mathbf{F}\mathbf{X}\left(\tau \right)\mathrm{d}\tau +\epsilon \mathbf{S}\mathbf{X}\left(\tau \right)\mathrm{d}\tau +\sqrt{\epsilon }\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(\tau \right)\mathrm{d}{W}^h\left(\tau \right),$(7) ), we will use the fact that the fast process has an asymptotic stable equilibrium, and we will approximate this system by another one in which the fast process has reached equilibrium.

Let $i=1,\dots ,q$ be fixed and let ${\hat{\mathbf{F}}}_i:={\mathbf{v}}_i{\mathbf{u}}_i^{\mathrm{T}}$. From Equation (5(5) $\begin{array}{rl}{\mathbf{F}}_i{\mathbf{v}}_i& =\mathbf{0};{\mathbf{u}}_i^{\mathrm{T}}{\mathbf{F}}_i={\mathbf{0}}^{\mathrm{T}}\\ \parallel {\mathbf{v}}_i{\parallel }_1& =1;{\mathbf{u}}_i^{\mathrm{T}}{\mathbf{v}}_i=1,\end{array}$(5) ), we have that if the system were governed by the fast process exclusively, for any initial condition ${\mathbf{x}}^i\in {\mathbb{R}}_{+}^{N_i}$ the population of each group i would tend to vector ${\hat{\mathbf{x}}}^i:={\hat{\mathbf{F}}}_i{\mathbf{x}}^i=\left({\mathbf{u}}_i^{\mathrm{T}}{\mathbf{x}}^i\right){\mathbf{v}}_i$. From this expression we note that vector ${\mathbf{v}}_i$ defines the equilibrium population structure for group i and ${\mathbf{u}}_i$ is a vector of reproductive values, that is, the larger $u_i^j$, the higher the contribution of the j-th subgroup of the ith group to the equilibrium population. Therefore ${\mathbf{u}}_i$ characterizes the size of the equilibrium population.

We define matrix $\hat{\mathbf{F}}$ in the following way, $\hat{\mathbf{F}}=\mathrm{diag}({\hat{\mathbf{F}}}_1,{\hat{\mathbf{F}}}_2,\dots ,{\hat{\mathbf{F}}}_q),$ and then for any initial condition $\mathbf{X}\in {\mathbb{R}}_{+}^N$, the equilibrium population for the fast dynamics in the whole system is given by $\hat{\mathbf{X}}=\hat{\mathbf{F}}\mathbf{X}$. Let us define the non-negative matrices $\mathbf{V}\mathbf{:}=\mathrm{diag}({\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_q)\in {\mathbb{R}}^{N\times q},\mathbf{U}\mathbf{:}=\mathrm{diag}({\mathbf{u}}_1^{\mathrm{T}},{\mathbf{u}}_2^{\mathrm{T}},\dots ,{\mathbf{u}}_q^{\mathrm{T}})\in {\mathbb{R}}^{q\times N}$ whose interpretation is immediate bearing in mind what we pointed out about ${\mathbf{v}}_i$ and ${\mathbf{u}}_i$.

Some of the properties of these matrices are gathered in the following lemma, whose proof is straightforward:

Lemma 4.1

Matrices $\mathbf{F},\hat{\mathbf{F}},\mathbf{V}$ and $\mathbf{U}$ verify:

1. $\hat{\mathbf{F}}={\mathbf{V}\mathbf{U}\mathbf{,}\mathbf{U}\mathbf{V}\mathbf{=}\mathbf{I}}_q$ and the columns of $\mathbf{V}$ are independent and constitute a basis of $\mathrm{Im}\hat{\mathbf{F}}$.

2. $\mathbf{F}\hat{\mathbf{F}}=\mathbf{0}$.

Now, from Equation (6(6) $\mathrm{d}\mathbf{X}\left(t\right)=\frac{1}{\epsilon }\mathbf{F}\mathbf{X}\left(t\right)\mathrm{d}t+\mathbf{S}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right).$(6) ) we will build an auxiliary system replacing the state variables in the right side by its equilibrium values for the fast process and use that $\mathbf{F}\hat{\mathbf{F}}=\mathbf{0}$, (8) $\begin{array}{rl}\mathrm{d}\mathbf{X}\left(t\right)& =\frac{1}{\epsilon }\mathbf{F}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right)\\ & =\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right).\end{array}$(8) Now we define the vector of global variables as (9) $\mathbf{Y}=(y^1,\dots ,y^q{)}^{\mathrm{T}}\mathbf{:}=\mathbf{U}\mathbf{X}\in {\mathbb{R}}^q$(9) and, multiplying Equation (8(8) $\begin{array}{rl}\mathrm{d}\mathbf{X}\left(t\right)& =\frac{1}{\epsilon }\mathbf{F}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right)\\ & =\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right).\end{array}$(8) ) on the left by $\mathbf{U}$ and using that $\hat{\mathbf{F}}=\mathbf{V}\mathbf{U}$ and Equation (9(9) $\mathbf{Y}=(y^1,\dots ,y^q{)}^{\mathrm{T}}\mathbf{:}=\mathbf{U}\mathbf{X}\in {\mathbb{R}}^q$(9) ) we obtain the aggregated system (10) $\mathrm{d}\mathbf{Y}\left(t\right)=\overline{\mathbf{S}}\mathbf{Y}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\overline{\mathbf{B}}}^h\mathbf{Y}\left(t\right)\mathrm{d}{W}^h\left(t\right),$(10) where we have defined (11) $\begin{array}{rl}\overline{\mathbf{S}}& =\left[{\overline{s}}_{ij}\right]:=\mathbf{U}\mathbf{S}\mathbf{V}=\left[\begin{array}{ccc}{\mathbf{u}}_1^{\mathrm{T}}{\mathbf{S}}_{11}{\mathbf{v}}_1& \cdots & {\mathbf{u}}_1^{\mathrm{T}}{\mathbf{S}}_{1q}{\mathbf{v}}_q\\ ⋮& \ddots & ⋮\\ {\mathbf{u}}_q^{\mathrm{T}}{\mathbf{S}}_{q1}{\mathbf{v}}_1& \cdots & {\mathbf{u}}_q^{\mathrm{T}}{\mathbf{S}}_{qq}{\mathbf{v}}_q\end{array}\right]\in {\mathbb{R}}^{q\times q}\end{array}$(11) (12) $\begin{array}{rl}{\overline{\mathbf{B}}}^h& :={\mathbf{U}\mathbf{B}}^h\mathbf{V}=\mathrm{diag}\left(\stackrel{N_1}{\sum _{\alpha =1}}{v}_1^{\alpha }{u}_1^{\alpha }{b}_{1\alpha }^h,\dots ,\stackrel{N_q}{\sum _{\alpha =1}}{v}_q^{\alpha }{u}_q^{\alpha }{b}_{q\alpha }^h\right)\in {\mathbb{R}}^q,h=1,\dots ,k.\end{array}$(12) Note that the global variables $\mathbf{Y}\left(t\right)$ defined by Equation (9(9) $\mathbf{Y}=(y^1,\dots ,y^q{)}^{\mathrm{T}}\mathbf{:}=\mathbf{U}\mathbf{X}\in {\mathbb{R}}^q$(9) ) have the following expression in terms of the variables $\mathbf{X}\left(t\right)$ of the auxiliary system: $y^i\left(t\right)={\mathbf{u}}_i^{\mathrm{T}}{\mathbf{x}}^i\left(t\right)=u_i^1{x}^{i1}\left(t\right)+u_i^2{x}^{i2}\left(t\right)+\cdots +u_i^{N_i}{x}^{i{N}_i}\left(t\right),i=1,\dots ,q.$ Note that:

1. $y^i\left(t\right)$ is a linear combination of the variables corresponding to group i, being the coefficients of the combination the components of vector ${\mathbf{u}}_i$. Recall that ${\mathbf{u}}_i$ is a vector of reproductive values for the fast process in group i. Therefore, for each $j=1,\dots ,N_i$, variable $x^{ij}\left(t\right)$ has a relative weight in $y^i\left(t\right)$ which is proportional to $u_i^j$, that is, proportional to the contribution to the total equilibrium population that an individual initially present in group i and subgroup j would have in the case that the system were governed by the fast process exclusively.

2. The global variables are conservative for the fast process. Indeed, suppose that the fast process is the only one acting in the system. Then we would have $\dot{\mathbf{X}}\left(t\right)=\mathbf{F}\mathbf{X}\left(t\right)/\epsilon$ and using that $\mathbf{U}\mathbf{F}=\mathbf{0}$, $\dot{\mathbf{Y}}\left(t\right)=\mathbf{U}\dot{\mathbf{X}}\left(t\right)=\mathbf{U}\mathbf{F}\mathbf{X}\left(t\right)/\epsilon =\mathbf{0}$.

3. The components of the matrices representing the drift and the diffusion for the reduced system are certain linear combinations of their analogues for the original system, where the coefficients of the combination are determined by the equilibrium characteristics of the fast process.

The next result together with Theorem 2.1 guarantees that the original and the aggregated systems have positive solutions for any positive initial conditions.

Proposition 4.2

$\mathbf{F}+\epsilon \mathbf{S}$ and $\overline{\mathbf{S}}$ are Metzler matrices and matrices ${\overline{\mathbf{B}}}^h,$ $h=1,\dots ,k$ are diagonal. Therefore, according to Theorem 2.1 both the original system (7(7) $\mathrm{d}\mathbf{X}\left(\tau \right)=\mathbf{F}\mathbf{X}\left(\tau \right)\mathrm{d}\tau +\epsilon \mathbf{S}\mathbf{X}\left(\tau \right)\mathrm{d}\tau +\sqrt{\epsilon }\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(\tau \right)\mathrm{d}{W}^h\left(\tau \right),$(7) ) and the aggregated system (10(10) $\mathrm{d}\mathbf{Y}\left(t\right)=\overline{\mathbf{S}}\mathbf{Y}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\overline{\mathbf{B}}}^h\mathbf{Y}\left(t\right)\mathrm{d}{W}^h\left(t\right),$(10) ) verify that for any positive initial condition the solution remains positive for all t>0 with probability one.

Proof.

$\mathbf{F}+\epsilon \mathbf{S}$ is clearly a Metzler matrix for it is the sum of Metzler matrices. Now let $i\ne j$. Since $\mathbf{S}$ is a Metzler matrix we have that ${\mathbf{S}}_{ij}\ge \mathbf{0}$ and using the fact that ${\mathbf{u}}_i$ and ${\mathbf{v}}_j$ are positive vectors, from Equation (11(11) $\begin{array}{rl}\overline{\mathbf{S}}& =\left[{\overline{s}}_{ij}\right]:=\mathbf{U}\mathbf{S}\mathbf{V}=\left[\begin{array}{ccc}{\mathbf{u}}_1^{\mathrm{T}}{\mathbf{S}}_{11}{\mathbf{v}}_1& \cdots & {\mathbf{u}}_1^{\mathrm{T}}{\mathbf{S}}_{1q}{\mathbf{v}}_q\\ ⋮& \ddots & ⋮\\ {\mathbf{u}}_q^{\mathrm{T}}{\mathbf{S}}_{q1}{\mathbf{v}}_1& \cdots & {\mathbf{u}}_q^{\mathrm{T}}{\mathbf{S}}_{qq}{\mathbf{v}}_q\end{array}\right]\in {\mathbb{R}}^{q\times q}\end{array}$(11) ) it follows that ${\overline{s}}_{ij}={\mathbf{u}}_i^{\mathrm{T}}{\mathbf{S}}_{ij}{\mathbf{v}}_j\ge 0$ and so $\overline{\mathbf{S}}$ is a Metzler matrix. Moreover, from Equation (12(12) $\begin{array}{rl}{\overline{\mathbf{B}}}^h& :={\mathbf{U}\mathbf{B}}^h\mathbf{V}=\mathrm{diag}\left(\stackrel{N_1}{\sum _{\alpha =1}}{v}_1^{\alpha }{u}_1^{\alpha }{b}_{1\alpha }^h,\dots ,\stackrel{N_q}{\sum _{\alpha =1}}{v}_q^{\alpha }{u}_q^{\alpha }{b}_{q\alpha }^h\right)\in {\mathbb{R}}^q,h=1,\dots ,k.\end{array}$(12) ) it is clear that ${\overline{\mathbf{B}}}^h$ is a diagonal matrix for each $h=1,\dots ,k$.

5. Multiregional models with two time scales

In this section we will illustrate the reduction technique by applying it to a multiregional model.

5.1. Model setting

We consider a population living in a multipatch system. We assume that there are a number N of different patches among which the individuals can migrate. The growth of the population in each patch is linear and is affected by stochasticity. Migration among patches is assumed to be deterministic.

We number the patches in the form $(i,\alpha )$, i=1, 2, $\alpha =1,\dots ,N_i$, where $N:=N_1+N_2$. Coming back to the terminology of Section 3, the first index, i, defines the ‘group’ of patches and the second, α, the ‘subgroup’ within that group. Note that in our setting we have chosen that the number q of groups is 2 just for the sake of simplicity in the expression of the matrices involved, but the generalization of the model to an arbitrary number q of groups is straightforward. Let $x^{i\alpha }\left(t\right)$ denote the population in group i and subgroup j at time t and let $\mathbf{X}=(x^{11},\dots ,x^{1N_1},x^{21},\dots ,x^{2N_2}{)}^{\mathrm{T}}\in {\mathbb{R}}^N$ be the population vector.

We assume that migration is fast within each group of patches, that is, from any patch $(i,\alpha )$ to any other patch of the form $(i,\beta )$, $\beta \ne \alpha$. Migration is assumed to be slow between patches belonging to different groups, that is, from any patch $(i,\alpha )$ to any other patch of the form $(j,\beta )$, with $j\ne i$. We can think of each group of patches as spatial regions located close to each other so that migration between them is easy, whereas different groups of patches correspond to regions amongst which migration is more difficult. For example, in a population of birds each group of patches can correspond to an island and the subgroups can correspond to the different spatial locations within an island, so that intra-island migrations are fast with respect to inter-island movements.

Let τ and t denote, respectively, the times corresponding to the fast and the slow migrations and let $\epsilon :=\tau /t$ be the ratio between both. We assume that the growth of the population takes place in the slow time scale. For each pair $(i,\alpha )$, i=1, 2, $\alpha =1,\dots ,N_i$, let $r_i^{\alpha }$ be the deterministic population growth rate in patch $(i,\alpha )$ and let us assume that this growth rate is affected by a noise defined by a certain linear combination ${\sigma }_{i\alpha }^1\mathrm{d}{W}^1\left(t\right)/\mathrm{d}t+\cdots +{\sigma }_{i\alpha }^k\mathrm{d}{W}^k\left(t\right)/\mathrm{d}t$ of (weakly defined) independent white noise processes, where ${\sigma }_{i\alpha }^h\ge 0$ for each $h=1,\dots ,k$. Now we define matrices $\begin{array}{rl}{\mathbf{R}}_i& :=\mathrm{diag}(r_i^1,\dots ,r_i^{N_i})\in {\mathbb{R}}^{N_i\times N_i},i=1,\dots ,q\\ \mathbf{R}& :=\mathrm{diag}({\mathbf{R}}_1,{\mathbf{R}}_2)\in {\mathbb{R}}^{N\times N}\\ {\mathbf{G}}_i^h& :=\mathrm{diag}({\sigma }_{i1}^h,\dots ,{\sigma }_{i{N}_i}^h)\in {\mathbb{R}}^{N_i\times N_i},i=1,\dots ,q,h=1,\dots ,k\\ {\mathbf{G}}^h& :=\mathrm{diag}({\mathbf{G}}_1^h,{\mathbf{G}}_2^h)\in {\mathbb{R}}^{N\times N},h=1,\dots ,k.\end{array}$ Regarding the slow migration between different groups of patches, for each $i\ne j$ and each $\alpha =1,\dots ,N_i$, $\beta =1,\dots ,N_j$, we define $l_{ij}^{a\beta }\ge 0$ as the (slow) migration coefficient from patch $(j,\beta )$ to patch $(i,\alpha )$. Similarly we define $l_{ii}^{a\beta }=0$ for each $\alpha ,\beta =1,\dots ,N_i$ with $\alpha \ne \beta$ (as there is no slow migration within a group of patches) and $l_{ii}^{\alpha \alpha }:=-\stackrel{q}{\sum _{j=1,j\ne i}}\stackrel{N_j}{\sum _{\beta =1,\beta \ne \alpha }}{l}_{ji}^{\beta \alpha },i=1,2,\alpha =1,\dots ,N_i.$ Let us define matrices ${\mathbf{L}}_{ij}:=\left[l_{ij}^{\alpha \beta }\right]\in {\mathbb{R}}^{N_i\times N_j}$, $i,j=1,\dots ,q$ and $\mathbf{L}\mathbf{:}=\left[\begin{array}{cc}{\mathbf{L}}_{11}& {\mathbf{L}}_{12}\\ {\mathbf{L}}_{21}& {\mathbf{L}}_{22}\end{array}\right].$ Then the slow process, that is, the joint effect of the growth process and of the slow migration between patches of different groups, can be modelled by the following system of SDEs $\mathrm{d}\mathbf{X}\left(t\right)=\left(\mathbf{L}\mathbf{+}\mathbf{R}\right)\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{G}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0.$ Regarding the fast migration between patches of the same group, for each i=1, 2 and $\alpha ,\beta =1,\dots ,N_i$ with $\alpha \ne \beta ,$ let $m_i^{\alpha \beta }\ge 0$ be the (fast) migration rate from patch $(i,\beta )$ to patch $(i,\alpha )$. For $\beta =\alpha$ we define (13) $m_i^{\alpha \alpha }:=-\stackrel{N_i}{\sum _{\beta =1,\beta \ne \alpha }}{m}_i^{\beta \alpha }.$(13) Now let ${\mathbf{M}}_i:=\left[m_i^{\alpha \beta }\right]\in {\mathbb{R}}^{N_i\times N_i}$, i=1, 2. We assume that the $m_i^{\alpha \beta }$ are such that ${\mathbf{M}}_i$ is irreducible for each i=1, 2. From Equation (13(13) $m_i^{\alpha \alpha }:=-\stackrel{N_i}{\sum _{\beta =1,\beta \ne \alpha }}{m}_i^{\beta \alpha }.$(13) ) we have that the columns of each ${\mathbf{M}}_i$ add up to zero and so matrix ${\mathbf{M}}_i+\mathbf{I}$ is a non-negative primitive column stochastic matrix. Therefore 0 is the (strictly) dominant eigenvalue of ${\mathbf{M}}_i$ and moreover it is simple. Let ${\mathbf{u}}_i:={\mathbf{1}}_i=(1,\dots ,1{)}^{\mathrm{T}}\in {\mathbb{R}}^{N_i}$ and ${\mathbf{v}}_i=(v_i^1,\dots ,v_i^{N_i}{)}^{\mathrm{T}}\in {\mathbb{R}}^{N_i}$ be its associated positive left and right eigenvectors, where we assume the normalization condition ${\mathbf{1}}_i^{\mathrm{T}}{\mathbf{v}}_i=1$. Note that vector ${\mathbf{v}}_i$ defines the equilibrium distribution between the different patches of group i when we consider the fast migration as the only process acting on the system. Now we define matrices $\begin{array}{rl}\mathbf{\text{M}}& :=\mathrm{diag}({\mathbf{M}}_1,{\mathbf{M}}_2)\in {\mathbb{R}}^{N\times N}\\ \mathbf{V}& :=\mathrm{diag}({\mathbf{v}}_1,{\mathbf{v}}_2)\in {\mathbb{R}}^{N\times 2},\mathbf{U}:=\mathrm{diag}({\mathbf{1}}_1^{\mathrm{T}},{\mathbf{1}}_2^{\mathrm{T}})\in {\mathbb{R}}^{2\times N}\\ {\hat{\mathbf{M}}}_i& :={\mathbf{v}}_i{\mathbf{1}}_i^{\mathrm{T}}\in {\mathbb{R}}^{N_i\times N_i},i=1,2,\hat{\mathbf{M}}:=\mathrm{diag}({\hat{\mathbf{M}}}_1,{\hat{\mathbf{M}}}_2)\in {\mathbb{R}}^{N\times N}.\end{array}$ Then, the complete model that takes into account the joint effect of the slow population growth in each patch and the slow and fast migrations between patches is given by $\mathrm{d}\mathbf{X}\left(t\right)=\left(\frac{1}{\epsilon }\mathbf{M}\mathbf{+}\mathbf{L}\mathbf{+}\mathbf{R}\right)\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{G}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right)$ or, using the fast time τ, (14) $\mathrm{d}\mathbf{X}\left(\tau \right)=\left(\mathbf{M}\mathbf{+}\epsilon \right(\mathbf{L}\mathbf{+}\mathbf{R}\mathbf{)}\left)\mathbf{X}\right(\tau )\mathrm{d}\tau +\sqrt{\epsilon }\stackrel{k}{\sum _{h=1}}{\mathbf{G}}^h\mathbf{X}(\tau \left)\mathrm{d}{W}^h\right(\tau ),\mathbf{X}(0)={\mathbf{X}}_0,$(14) which constitutes a system of N linear SDEs.

Note that under the above conditions we are in the general setting of Section 3 by taking $\mathbf{S}=\mathbf{L}\mathbf{+}\mathbf{R},\mathbf{F}\mathbf{=}\mathbf{M},\hat{\mathbf{F}}=\hat{\mathbf{M}},{\mathbf{B}}^h={\mathbf{G}}^h,h=1,\dots ,k$ and Hypotheses S1, S2, F1, F2 and F3 for the slow and fast processes are met.

5.2. Model reduction

Therefore we can proceed to the aggregation of the original system following the procedure developed in Section 4. We define the global variables $\mathbf{Y}:=\mathbf{U}\mathbf{X}$, that is, $y^i\left(t\right)={\mathbf{1}}_i^{\mathrm{T}}{\mathbf{x}}^i\left(t\right)=x^{i1}\left(t\right)+x^{i2}\left(t\right)+\cdots +x^{i{N}_i}\left(t\right),i=1,2,$ so that $y^i$ is the total population in all patches of group i, and the reduced aggregated system is the two dimensional SDE (15) $\mathrm{d}\mathbf{Y}\left(t\right)=\overline{\mathbf{S}}\mathbf{Y}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\overline{\mathbf{G}}}^h\mathbf{Y}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{Y}\left(0\right)={\mathbf{U}\mathbf{X}}_0,$(15) where $\begin{array}{rl}\overline{\mathbf{S}}& :=\mathbf{U}(\mathbf{L}\mathbf{+}\mathbf{R}\mathbf{)}\mathbf{V}=\left(\begin{array}{cc}{\mathbf{1}}_1^{\mathrm{T}}({\mathbf{L}}_{11}+{\mathbf{R}}_1){\mathbf{v}}_1& {\mathbf{1}}_1^{\mathrm{T}}{\mathbf{L}}_{12}{\mathbf{v}}_2\\ {\mathbf{1}}_2^{\mathrm{T}}{\mathbf{L}}_{21}{\mathbf{v}}_1& {\mathbf{1}}_1^{\mathrm{T}}({\mathbf{L}}_{22}+{\mathbf{R}}_2){\mathbf{v}}_2\end{array}\right)\in {\mathbb{R}}^{2\times 2}\\ {\overline{\mathbf{G}}}^h& :={\mathbf{U}\mathbf{G}}^h\mathbf{V}=\mathrm{diag}({\mathbf{1}}_1^{\mathrm{T}}{\mathbf{G}}_1^h{\mathbf{v}}_1,{\mathbf{1}}_2^{\mathrm{T}}{\mathbf{G}}_2^h{\mathbf{v}}_2),h=1,\dots ,k.\end{array}$ In order to be more specific we will consider the particular case in which $N_1=N_2=2$, that is, there are 4 patches $(1,1)$, $(1,2)$, $(2,1)$ and $(2,2),$ and the population vector is $\mathbf{X}=(x^{11},x^{12},x^{21},x^{22}{)}^{\mathrm{T}}$. Dispersal is fast between patches $(1,1)$ and $(1,2)$ and between $(2,1)$ and $(2,2)$, and is slow in the rest of cases.

Let i=1, 2 be fixed. We will keep the notation introduced above except in the case of the fast migrations, where in order to simplify it we will denote by $p_i$ and $q_i$ the (fast) migration rates from $(i,1)$ to $(i,2)$ and from $(i,2)$ to $(i,1)$ respectively. Therefore ${\mathbf{M}}_i=\left(\begin{array}{cc}-p_i& q_i\\ p_i& -q_i\end{array}\right),i=1,2.$ We assume $0<p_i,q_i<1$ so that matrix ${\mathbf{M}}_i$ is irreducible and vector ${\mathbf{v}}_i$ has the form ${\mathbf{v}}_i=(q_i,p_i{)}^{\mathrm{T}}/(p_i+q_i)$. Then the complete model has the form (14(14) $\mathrm{d}\mathbf{X}\left(\tau \right)=\left(\mathbf{M}\mathbf{+}\epsilon \right(\mathbf{L}\mathbf{+}\mathbf{R}\mathbf{)}\left)\mathbf{X}\right(\tau )\mathrm{d}\tau +\sqrt{\epsilon }\stackrel{k}{\sum _{h=1}}{\mathbf{G}}^h\mathbf{X}(\tau \left)\mathrm{d}{W}^h\right(\tau ),\mathbf{X}(0)={\mathbf{X}}_0,$(14) ) with ${\mathbf{G}}^h=\mathrm{diag}({\sigma }_{11}^h,{\sigma }_{12}^h,{\sigma }_{21}^h,{\sigma }_{22}^h)$, $h=1,\dots ,k$ and $\mathbf{M}+\epsilon (\mathbf{L}\mathbf{+}\mathbf{R}\mathbf{)}$ given by $\left(\begin{array}{cccc}-p_1-\epsilon (l_{21}^{11}+l_{21}^{21}-r_1^1)& q_1& \epsilon l_{12}^{11}& \epsilon l_{12}^{12}\\ p_1& -q_1-\epsilon (l_{21}^{12}+l_{21}^{22}-r_1^2)& \epsilon l_{12}^{21}& \epsilon l_{12}^{22}\\ \epsilon l_{21}^{11}& \epsilon l_{21}^{12}& -p_2-\epsilon (l_{12}^{11}+l_{12}^{21}-r_2^1)& q_2\\ \epsilon l_{21}^{21}& \epsilon l_{21}^{22}& p_2& -q_2-\epsilon (l_{12}^{12}+l_{12}^{22}-r_2^2)\end{array}\right).$ The global variables are $\mathbf{Y}=(y^1,y^2{)}^{\mathrm{T}}=(x^{11}+x^{12},x^{21}+x^{22}{)}^{\mathrm{T}}$ and the reduced system has the form (15(15) $\mathrm{d}\mathbf{Y}\left(t\right)=\overline{\mathbf{S}}\mathbf{Y}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\overline{\mathbf{G}}}^h\mathbf{Y}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{Y}\left(0\right)={\mathbf{U}\mathbf{X}}_0,$(15) ) with (16) $\begin{array}{rl}\overline{\mathbf{S}}& =\left(\begin{array}{cc}{\overline{s}}_{11}& {\overline{s}}_{12}\\ {\overline{s}}_{21}& {\overline{s}}_{22}\end{array}\right)\\ & =\left(\begin{array}{cc}\frac{q_1(-l_{21}^{11}-l_{21}^{21}+r_1^1)+p_1(-l_{21}^{12}-l_{21}^{22}+r_1^2)}{p_1+q_1}& \frac{q_1(l_{12}^{11}+l_{12}^{21})+p_1(l_{12}^{12}+l_{12}^{22})}{p_2+q_2}\\ \frac{q_1(l_{21}^{11}+l_{21}^{21})+p_1(l_{21}^{12}+l_{21}^{22})}{p_1+q_1}& \frac{q_2(-l_{12}^{11}-l_{12}^{21}+r_2^1)+p_2(-l_{12}^{12}-l_{12}^{22}+r_2^2)}{p_2+q_2}\end{array}\right)\\ {\overline{\mathbf{G}}}^h& =\mathrm{diag}({\overline{\sigma }}_1^h,{\overline{\sigma }}_2^h)=\mathrm{diag}\left(\frac{q_1{\sigma }_{11}^h+p_1{\sigma }_{12}^h}{p_1+q_1},\frac{q_2{\sigma }_{21}^h+p_2{\sigma }_{22}^h}{p_2+q_2}\right),h=1,\dots ,k.\end{array}$(16)

5.3. Case in which migration is fast

Let us now consider the particular but relevant case in which there is only one group of patches, and therefore all the migrations among the patches are fast with respect to demography. This has been the usual setting in the literature of approximate aggregation techniques [2, 3]. We can consider this case as a particular instance of the setting in Section 5.1 when we take q equal to 1 instead of equal to 2. In order to simplify the notation, we omit the index 1 regarding the group of patches. Therefore, we denote the patches as $j=1,\dots ,N$, and the vectors and matrices have the form $\mathbf{X}=(x^1,\dots ,x^N{)}^{\mathrm{T}}$, $\mathbf{R}:=\mathrm{diag}(r^1,\dots ,r^N)$, $\mathbf{L}=\mathbf{0}$ (there are no slow migrations), $\mathbf{M}:=\left[m^{\alpha \beta }\right]$ and ${\mathbf{G}}^h:=\mathrm{diag}({\sigma }_1^h,\dots ,{\sigma }_N^h)$. Therefore, the original system (14(14) $\mathrm{d}\mathbf{X}\left(\tau \right)=\left(\mathbf{M}\mathbf{+}\epsilon \right(\mathbf{L}\mathbf{+}\mathbf{R}\mathbf{)}\left)\mathbf{X}\right(\tau )\mathrm{d}\tau +\sqrt{\epsilon }\stackrel{k}{\sum _{h=1}}{\mathbf{G}}^h\mathbf{X}(\tau \left)\mathrm{d}{W}^h\right(\tau ),\mathbf{X}(0)={\mathbf{X}}_0,$(14) ) takes the form $\mathrm{d}\mathbf{X}\left(\tau \right)=\left(\mathbf{M}\mathbf{+}\epsilon \mathbf{R}\right)\mathbf{X}\left(\tau \right)\mathrm{d}\tau +\sqrt{\epsilon }\stackrel{k}{\sum _{h=1}}{\mathbf{G}}^h\mathbf{X}\left(\tau \right)\mathrm{d}{W}^h\left(\tau \right).$ In this case there is only one global variable $y\left(t\right)=x^1\left(t\right)+x^2\left(t\right)+\cdots +x^N\left(t\right)$ and if $\mathbf{v}=(v^1,\dots ,v^N{)}^{\mathrm{T}}$ and $\mathbf{1}=(1,\dots ,1{)}^{\mathrm{T}}$ are the positive right and left eigenvectors of $\mathbf{M}$ associated to eigenvalue 0 scaled so that ${\mathbf{1}}^{\mathrm{T}}\mathbf{v}=1$, we have $\mathbf{V}=\mathbf{v},$ $\mathbf{U}={\mathbf{1}}^{\mathrm{T}}$. Therefore, the reduced system is the scalar SDE $\mathrm{d}y\left(t\right)=\overline{s}y\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\overline{g}}^{h}y\left(t\right)\mathrm{d}{W}^h\left(t\right),y\left(0\right)=\stackrel{N}{\sum _{j=1}}{x}^j\left(0\right),$ where (17) $\overline{s}:={\mathbf{1}}^{\mathrm{T}}\mathbf{R}\mathbf{v}\mathbf{=}\stackrel{N}{\sum _{j=1}}{v}^j{r}^j,{\overline{g}}^h={\mathbf{1}}^{\mathrm{T}}{\mathbf{G}}^h\mathbf{v}=\stackrel{N}{\sum _{j=1}}{v}^j{\sigma }_j^h.$(17)

6. Relationships between the original and the reduced systems

The aim of this Section is to relate the asymptotic behaviour of the first- and second-order moments of the solution of the original and the reduced systems (8(8) $\begin{array}{rl}\mathrm{d}\mathbf{X}\left(t\right)& =\frac{1}{\epsilon }\mathbf{F}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right)\\ & =\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right).\end{array}$(8) ) and (10(10) $\mathrm{d}\mathbf{Y}\left(t\right)=\overline{\mathbf{S}}\mathbf{Y}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\overline{\mathbf{B}}}^h\mathbf{Y}\left(t\right)\mathrm{d}{W}^h\left(t\right),$(10) ) introduced in Section 3. The two systems fit in the framework of [14] in which the quotient between the rate of diffusion squared and the speed of drift is $O\left({\epsilon }^p\right)$, $\epsilon \to 0$, with $p>1/2$, but the results obtained in that reference, essentially a stochastic extension of the Tikhonov theory of singular perturbations, are valid only for finite time intervals.

The first- and second-order moments of the solution of general linear stochastic equations with a deterministic initial condition are finite and can be calculated as the solution of certain linear ordinary differential equations. Specifically [1], for system (1(1) $\mathrm{d}\mathbf{X}\left(t\right)=\mathbf{A}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0,$(1) ) $\mathbb{E}\mathbf{X}\left(t\right)$ is the unique solution to the equation $\frac{\mathrm{d}}{\mathrm{d}t}\mathbb{E}\mathbf{X}\left(t\right)=\mathbf{A}\mathbb{E}\mathbf{X}\left(t\right),\mathbb{E}\mathbf{X}\left(0\right)={\mathbf{X}}_0,$ whereas the matrix of second-order moments $\mathbf{H}\left(t\right):=\mathbb{E}\left(\mathbf{X}\right(t\left)\mathbf{X}\right(t{)}^{\mathrm{T}})\in {\mathbb{R}}^{d\times d}$ is the unique non-negative-definite symmetric solution of the equation (18) $\frac{\mathrm{d}\mathbf{H}\left(t\right)}{\mathrm{d}t}=\mathbf{A}\mathbf{H}\left(t\right)+\mathbf{H}\left(t\right){\mathbf{A}}^{\mathrm{T}}+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{H}\left(t\right)({\mathbf{B}}^h{)}^{\mathrm{T}},\mathbf{H}(0)={\mathbf{X}}_0{\mathbf{X}}_0^{\mathrm{T}}.$(18) Regarding Equation (18(18) $\frac{\mathrm{d}\mathbf{H}\left(t\right)}{\mathrm{d}t}=\mathbf{A}\mathbf{H}\left(t\right)+\mathbf{H}\left(t\right){\mathbf{A}}^{\mathrm{T}}+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{H}\left(t\right)({\mathbf{B}}^h{)}^{\mathrm{T}},\mathbf{H}(0)={\mathbf{X}}_0{\mathbf{X}}_0^{\mathrm{T}}.$(18) ), in order to work with a more tractable expression, we will make use of the Kronecker matrix product and the ‘vec’ operator (see [9] for details). The Kronecker matrix product for two matrices $\mathbf{A}=\left[a_{ij}\right]\in {\mathbb{C}}^{m\times n}$ and $\mathbf{B}=\left[b_{ij}\right]\in {\mathbb{C}}^{r\times s}$ is defined as a matrix of size $mr\times ns$ with mn blocks in which the block in position $(i,j)$ has the form $a_{ij}\mathbf{B}$. For any matrix $\mathbf{A}\in {\mathbb{C}}^{m\times n}$, vec$\mathbf{A}$ is defined as the column vector that contains, in order, the columns of $\mathbf{A}$. For any matrices $\mathbf{A}$, $\mathbf{Y}$, and $\mathbf{B}$ for which the product $\mathbf{A}\mathbf{Y}\mathbf{B}$ makes sense, one has vec$\left(\mathbf{A}\mathbf{Y}\mathbf{B}\right)=({\mathbf{B}}^{\mathrm{T}}\otimes \mathbf{A})$vec$\mathbf{Y}$ [9, Theorem 6.8], and so applying the ‘vec’ operator to both sides of Equation (18(18) $\frac{\mathrm{d}\mathbf{H}\left(t\right)}{\mathrm{d}t}=\mathbf{A}\mathbf{H}\left(t\right)+\mathbf{H}\left(t\right){\mathbf{A}}^{\mathrm{T}}+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{H}\left(t\right)({\mathbf{B}}^h{)}^{\mathrm{T}},\mathbf{H}(0)={\mathbf{X}}_0{\mathbf{X}}_0^{\mathrm{T}}.$(18) ) we obtain that the second-order moments verify $\frac{\mathrm{d}}{\mathrm{d}t}\mathrm{vec}\mathbf{H}\left(t\right)={\mathbf{A}}_2\mathrm{vec}\mathbf{H}\left(t\right),\mathrm{vec}\mathbf{H}\left(0\right)=\mathrm{vec}\left({\mathbf{X}}_0{\mathbf{X}}_0^{\mathrm{T}}\right),$

where ${\mathbf{A}}_2{\mathbf{=}\mathbf{I}}_N\otimes \mathbf{A}\mathbf{+}\mathbf{A}\otimes {\mathbf{I}}_N+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\otimes {\mathbf{B}}^h\in {\mathbb{R}}^{d^2\times d^2}$ and ${\mathbf{I}}_N$ denotes the $N\times N$ identity matrix.

Therefore, the asymptotic behaviour of the first- and second-order moments of system (1(1) $\mathrm{d}\mathbf{X}\left(t\right)=\mathbf{A}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),\mathbf{X}\left(0\right)={\mathbf{X}}_0,$(1) ) can be characterized in terms of the dominant eigenvalue and eigenvectors of matrices $\mathbf{A}$ and ${\mathbf{A}}_2$. We will now use this fact to study the moments of systems (8(8) $\begin{array}{rl}\mathrm{d}\mathbf{X}\left(t\right)& =\frac{1}{\epsilon }\mathbf{F}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right)\\ & =\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right).\end{array}$(8) ) and (10(10) $\mathrm{d}\mathbf{Y}\left(t\right)=\overline{\mathbf{S}}\mathbf{Y}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\overline{\mathbf{B}}}^h\mathbf{Y}\left(t\right)\mathrm{d}{W}^h\left(t\right),$(10) ).

In the case of the original system (8(8) $\begin{array}{rl}\mathrm{d}\mathbf{X}\left(t\right)& =\frac{1}{\epsilon }\mathbf{F}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right)\\ & =\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right).\end{array}$(8) ), the analogous to matrices $\mathbf{A}$ and ${\mathbf{A}}_2$ in the previous reasonings are matrices $\mathbf{C}\in {\mathbb{R}}^{N\times N}$ and ${\mathbf{C}}_2\in {\mathbb{R}}^{N^2\times N^2}$ given by (19) $\mathbf{C}\left(\epsilon \right):=\mathbf{F}+\epsilon \mathbf{S},{\mathbf{C}}_2\left(\epsilon \right):={\mathbf{F}}_2+\epsilon {\mathbf{S}}_2,$(19) where we have defined ${\mathbf{F}}_2:={\mathbf{I}}_N\otimes {\mathbf{F}\mathbf{+}\mathbf{F}\otimes \mathbf{I}}_N{\mathbf{,}\mathbf{S}}_2:={\mathbf{I}}_N\otimes {\mathbf{S}\mathbf{+}\mathbf{S}\otimes \mathbf{I}}_N+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\otimes {\mathbf{B}}^h.$ In the case of the aggregated system (10(10) $\mathrm{d}\mathbf{Y}\left(t\right)=\overline{\mathbf{S}}\mathbf{Y}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\overline{\mathbf{B}}}^h\mathbf{Y}\left(t\right)\mathrm{d}{W}^h\left(t\right),$(10) ) the corresponding matrices are $\overline{\mathbf{C}}\left(\epsilon \right)=\epsilon \overline{\mathbf{S}}\in {\mathbb{R}}^{q\times q},{\overline{\mathbf{C}}}_2\left(\epsilon \right)=\epsilon {\overline{\mathbf{S}}}_2\in {\mathbb{R}}^{q^2\times q^2}$ where ${\overline{\mathbf{S}}}_2:={\mathbf{I}}_q\otimes \overline{\mathbf{S}}+\overline{\mathbf{S}}\otimes {\mathbf{I}}_q+\stackrel{k}{\sum _{h=1}}{\overline{\mathbf{B}}}^h\otimes {\overline{\mathbf{B}}}^h.$ We introduce the following two hypotheses:

Hypothesis 6.1

Matrix $\overline{\mathbf{S}}$ has a simple and strictly dominant real eigenvalue μ associated to non-negative right and left eigenvectors $\mathbf{r}$ and $\mathbf{l}$, respectively.

Hypothesis 6.2

Matrix ${\overline{\mathbf{S}}}_2$ has a simple and strictly dominant real eigenvalue ${\mu }_2$ associated to non-negative right and left eigenvectors ${\mathbf{r}}_2$ and ${\mathbf{l}}_2$, respectively.

From Hypothesis 6.1 we have that the first-order moments of the reduced system have the following asymptotic behaviour $\underset{t\to \mathrm{\infty }}{lim}\frac{\mathbb{E}\mathbf{Y}\left(t\right)}{\exp \left(\mu t\right)}=\frac{⟨\mathbf{l},{\mathbf{y}}_0⟩}{⟨\mathbf{l},\mathbf{r}⟩}\mathbf{r},$ where ${\mathbf{y}}_0:={\mathbf{U}\mathbf{X}}_0$ and $⟨\ast ,\ast ⟩$ denotes the standard vector scalar product. Analogously, from Hypothesis 6.2 we have, for the second-order moments, $\underset{t\to \mathrm{\infty }}{lim}\frac{\mathrm{vec}\mathbb{E}\left(\mathbf{Y}\right(t\left)\mathbf{Y}\right(t{)}^{\mathrm{T}})}{\exp \left({\mu }_{2}t\right)}=\frac{⟨{\mathbf{l}}_2,{\mathbf{y}}_0\otimes {\mathbf{y}}_0⟩}{⟨{\mathbf{l}}_2,{\mathbf{r}}_2⟩}{\mathbf{r}}_2.$ The next result relates the asymptotic behaviour of the first- and second-order moments of the original system (8(8) $\begin{array}{rl}\mathrm{d}\mathbf{X}\left(t\right)& =\frac{1}{\epsilon }\mathbf{F}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right)\\ & =\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right).\end{array}$(8) ) with that of the reduced system when ϵ is small, that is, when the separation between the time scales of migration and of demography is large enough:

Theorem 6.1

(a) Let Hypothesis 6.1 hold. Then for small enough $\epsilon >0,$ matrix $\mathbf{C}\left(\epsilon \right)=\mathbf{F}+\epsilon \mathbf{S}$ has a simple and strictly dominant eigenvalue $\lambda \left(\epsilon \right)$ that can be written in the form $\lambda \left(\epsilon \right)=\epsilon \mu +O\left({\epsilon }^2\right)$ with associated right and left non-negative eigenvectors $r\left(\epsilon \right)$ and $l\left(\epsilon \right)$ such that $\mathbf{r}\left(\epsilon \right)=\mathbf{V}\mathbf{r}+\mathbf{O}\left(\epsilon \right);\mathbf{l}\left(\epsilon \right)={\mathbf{U}}^{\mathrm{T}}\mathbf{l}+\mathbf{O}\left(\epsilon \right).$ Consequently, for the original system (8(8) $\begin{array}{rl}\mathrm{d}\mathbf{X}\left(t\right)& =\frac{1}{\epsilon }\mathbf{F}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right)\\ & =\mathbf{S}\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}t+\stackrel{k}{\sum _{h=1}}{\mathbf{B}}^h\hat{\mathbf{F}}\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right).\end{array}$(8) ) we have: (20) $\underset{\tau \to \mathrm{\infty }}{lim}\frac{\mathbb{E}\mathbf{X}\left(\tau \right)}{\exp (\epsilon \mu +O({\epsilon }^2\left)\right)}=\frac{⟨\mathbf{l},{\mathbf{U}\mathbf{X}}_0⟩}{⟨\mathbf{l},\mathbf{r}⟩}\mathbf{V}\mathbf{r}+\mathbf{O}\left(\epsilon \right).$(20) (b) Let Hypothesis 6.2 hold. Then for small enough $\epsilon >0,$ matrix ${\mathbf{C}}_2\left(\epsilon \right)={\mathbf{F}}_2+\epsilon {\mathbf{S}}_2$ has a simple and strictly dominant eigenvalue ${\lambda }_2\left(\epsilon \right)$ that can be written in the form (21) ${\lambda }_2\left(\epsilon \right)=\epsilon {\mu }_2+O\left({\epsilon }^2\right)$(21) with associated right and left non-negative eigenvectors ${\mathbf{r}}_2\left(\epsilon \right)$ and ${\mathbf{l}}_2\left(\epsilon \right)$ such that (22) ${\mathbf{r}}_2\left(\epsilon \right)=(\mathbf{V}\otimes \mathbf{V}){\mathbf{r}}_2+\mathbf{O}\left(\epsilon \right);{\mathbf{l}}_2\left(\epsilon \right)=(\mathbf{U}\otimes \mathbf{U}{)}^{\mathrm{T}}{\mathbf{l}}_2+\mathbf{O}(\epsilon ).$(22) Consequently, the asymptotic behaviour of the second-order moment of the original system is characterized by: (23) $\underset{\tau \to \mathrm{\infty }}{lim}\frac{\mathrm{vec}\mathbb{E}\left(\mathbf{X}\right(t\left)\mathbf{X}\right(t{)}^{\mathrm{T}})}{\exp (\epsilon {\mu }_2+O({\epsilon }^2\left)\right)}=\frac{⟨{\mathbf{l}}_2,\left(\mathbf{U}\otimes \mathbf{U}\right)\mathrm{vec}\left({\mathbf{X}}_0{\mathbf{X}}_0^{\mathrm{T}}\right)⟩}{⟨{\mathbf{l}}_2,{\mathbf{r}}_2⟩}\left(\mathbf{V}\otimes \mathbf{V}\right){\mathbf{r}}_2+\mathbf{O}\left(\epsilon \right).$(23)

Proof.

See Appendix.

From the last theorem it follows in particular that if $\mu <0$ (resp. $\mu >0$), then $\lambda \left(\epsilon \right)<0$ (resp. $\lambda \left(\epsilon \right)>0$) for small enough ϵ, so that if the expected value of the population vector tends to zero (resp. infinity) in the reduced system the same happens for the original one when ϵ is small. Something analogous happens for the second-order moments regarding ${\mu }_2$ and ${\lambda }_2\left(\epsilon \right)$.

Let us apply this result to relate the behaviour of the multiregional models of Section 5. In the first place we will consider the model with q=2, $N_1=N_2=2$. Note from (16(16) $\begin{array}{rl}\overline{\mathbf{S}}& =\left(\begin{array}{cc}{\overline{s}}_{11}& {\overline{s}}_{12}\\ {\overline{s}}_{21}& {\overline{s}}_{22}\end{array}\right)\\ & =\left(\begin{array}{cc}\frac{q_1(-l_{21}^{11}-l_{21}^{21}+r_1^1)+p_1(-l_{21}^{12}-l_{21}^{22}+r_1^2)}{p_1+q_1}& \frac{q_1(l_{12}^{11}+l_{12}^{21})+p_1(l_{12}^{12}+l_{12}^{22})}{p_2+q_2}\\ \frac{q_1(l_{21}^{11}+l_{21}^{21})+p_1(l_{21}^{12}+l_{21}^{22})}{p_1+q_1}& \frac{q_2(-l_{12}^{11}-l_{12}^{21}+r_2^1)+p_2(-l_{12}^{12}-l_{12}^{22}+r_2^2)}{p_2+q_2}\end{array}\right)\\ {\overline{\mathbf{G}}}^h& =\mathrm{diag}({\overline{\sigma }}_1^h,{\overline{\sigma }}_2^h)=\mathrm{diag}\left(\frac{q_1{\sigma }_{11}^h+p_1{\sigma }_{12}^h}{p_1+q_1},\frac{q_2{\sigma }_{21}^h+p_2{\sigma }_{22}^h}{p_2+q_2}\right),h=1,\dots ,k.\end{array}$(16) ) that if there is at least a non-zero coefficient $l_{ij}^{\alpha \beta }$, matrix $\overline{\mathbf{S}}$ is irreducible and then Hypothesis 6.1 holds. Its dominant eigenvalue is given by $\mu =\frac{{\overline{s}}_{11}+{\overline{s}}_{22}+\sqrt{({\overline{s}}_{11}-{\overline{s}}_{22}{)}^2+4{\overline{s}}_{12}{\overline{s}}_{21}}}{2}$ and so the rate of growth of the first-order moments for the original system is $\epsilon \mu +O\left({\epsilon }^2\right)$ for small enough ϵ. This expression allows one to study how the different combinations of the migration and growth parameters affect the asymptotic behaviour of the expected value of the population. In particular, $\mu <0$ if and only if $\mathrm{tr}\overline{S}<0$ and $det\overline{S}>0$. Since $\overline{\mathbf{S}}$ has order two we can calculate right and left eigenvectors $\mathbf{r}$ and $\mathbf{l}$ associated to μ and then apply Equations (20(20) $\underset{\tau \to \mathrm{\infty }}{lim}\frac{\mathbb{E}\mathbf{X}\left(\tau \right)}{\exp (\epsilon \mu +O({\epsilon }^2\left)\right)}=\frac{⟨\mathbf{l},{\mathbf{U}\mathbf{X}}_0⟩}{⟨\mathbf{l},\mathbf{r}⟩}\mathbf{V}\mathbf{r}+\mathbf{O}\left(\epsilon \right).$(20) ) and (23(23) $\underset{\tau \to \mathrm{\infty }}{lim}\frac{\mathrm{vec}\mathbb{E}\left(\mathbf{X}\right(t\left)\mathbf{X}\right(t{)}^{\mathrm{T}})}{\exp (\epsilon {\mu }_2+O({\epsilon }^2\left)\right)}=\frac{⟨{\mathbf{l}}_2,\left(\mathbf{U}\otimes \mathbf{U}\right)\mathrm{vec}\left({\mathbf{X}}_0{\mathbf{X}}_0^{\mathrm{T}}\right)⟩}{⟨{\mathbf{l}}_2,{\mathbf{r}}_2⟩}\left(\mathbf{V}\otimes \mathbf{V}\right){\mathbf{r}}_2+\mathbf{O}\left(\epsilon \right).$(23) ) to obtain the full information about the asymptotic behaviour of the first-order moments of the system.

We could argue analogously for the second-order moments: in this case matrix ${\overline{\mathbf{S}}}_2$ is of order four and so an analytical computation of ${\mu }_2$, ${\mathbf{r}}_2$ and ${\mathbf{l}}_2$ is non-feasible, but for instance we can use the Routh–Hurwitz criterion [17] to obtain conditions for ${\mu }_2$ to be negative, so that $\epsilon {\mu }_2+O\left({\epsilon }^2\right)$ will also be negative for small ϵ.

In the case of the system of Section 5.3, the analysis is simpler cause the aggregated system is scalar and so are $\overline{\mathbf{S}}$ and ${\overline{\mathbf{S}}}_2$, so that Hypotheses 6.1 and 6.2 hold trivially. Moreover, $\mu =\overline{s}$ and ${\mu }_2=2\overline{s}+\sum _{h=1}^k({\overline{g}}^h{)}^2$ where $\overline{s}$ and the ${\overline{g}}^h$ are given by Equation (17(17) $\overline{s}:={\mathbf{1}}^{\mathrm{T}}\mathbf{R}\mathbf{v}\mathbf{=}\stackrel{N}{\sum _{j=1}}{v}^j{r}^j,{\overline{g}}^h={\mathbf{1}}^{\mathrm{T}}{\mathbf{G}}^h\mathbf{v}=\stackrel{N}{\sum _{j=1}}{v}^j{\sigma }_j^h.$(17) ) and we can take $\mathbf{r}=\mathbf{l}={\mathbf{r}}_2={\mathbf{l}}_2=1$. Therefore from Equations (20(20) $\underset{\tau \to \mathrm{\infty }}{lim}\frac{\mathbb{E}\mathbf{X}\left(\tau \right)}{\exp (\epsilon \mu +O({\epsilon }^2\left)\right)}=\frac{⟨\mathbf{l},{\mathbf{U}\mathbf{X}}_0⟩}{⟨\mathbf{l},\mathbf{r}⟩}\mathbf{V}\mathbf{r}+\mathbf{O}\left(\epsilon \right).$(20) ) and (23(23) $\underset{\tau \to \mathrm{\infty }}{lim}\frac{\mathrm{vec}\mathbb{E}\left(\mathbf{X}\right(t\left)\mathbf{X}\right(t{)}^{\mathrm{T}})}{\exp (\epsilon {\mu }_2+O({\epsilon }^2\left)\right)}=\frac{⟨{\mathbf{l}}_2,\left(\mathbf{U}\otimes \mathbf{U}\right)\mathrm{vec}\left({\mathbf{X}}_0{\mathbf{X}}_0^{\mathrm{T}}\right)⟩}{⟨{\mathbf{l}}_2,{\mathbf{r}}_2⟩}\left(\mathbf{V}\otimes \mathbf{V}\right){\mathbf{r}}_2+\mathbf{O}\left(\epsilon \right).$(23) ) we have $\begin{array}{rl}\underset{\tau \to \mathrm{\infty }}{lim}\frac{\mathbb{E}\mathbf{X}\left(\tau \right)}{\exp (\epsilon \mu +O({\epsilon }^2\left)\right)}& =y\left(0\right)\mathbf{v}+\mathbf{O}\left(\epsilon \right),\\ \underset{\tau \to \mathrm{\infty }}{lim}\frac{\mathrm{vec}\mathbb{E}\left(\mathbf{X}\right(t\left)\mathbf{X}\right(t{)}^{\mathrm{T}})}{\exp (\epsilon {\mu }_2+O({\epsilon }^2\left)\right)}& =y\left(0{)}^2\right(\mathbf{v}\otimes \mathbf{v})+\mathbf{O}(\epsilon ).\end{array}$ Note in particular that, when ϵ is small enough, the expected value of the population vector tends to zero if $\overline{s}<0$, and the second-order moments also tend to zero if $\overline{s}<-\sum _{h=1}^k({\overline{g}}^h{)}^2/2$.

7. Discussion

In this work we have presented a technique to carry out the reduction of linear two-scale population models governed by SDEs, therefore allowing one to simplify the treatment of these models. The reduction of the original model with N variables is carried out by letting the fast process reach equilibrium and defining an appropriate set of q global variables, q<N, which are linear combinations of the state variables and are conservative for the fast process. Moreover, we have obtained conditions that guarantee the positivity of solutions to the model.

We have also presented a result that allows one to know the the asymptotic behaviour of the first- and second-order statistical moments of the population vector of the original system through the computation of the dominant eigenvalues and eigenvectors of certain matrices of dimensions $q\times q$ and $q^2\times q^2$ associated to the reduced system.

The aggregation technique has been applied to study stochastic multiregional models in which migration between some patches can be fast with respect to demography whereas dispersal between other patches happens at the scale of demography. In the simpler case in which all migrations are fast, the reduced system is a scalar SDE whose analysis is straightforward, therefore allowing one to easily characterize the asymptotic behaviour of those moments for the original multiregional model.

This work suggests possible lines of future development: on the first hand, and still within the linear setting, work needs to be done to try to relate the stochastic stability [15] of the origin in the original and the reduced system, for this will provide more information regarding the persistence-extinction of the population. Secondly, in order to be able to study more realistic population models, the technique should be extended to nonlinear settings.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Luis Sanz http://orcid.org/0000-0002-1054-4568

Additional information

Funding

Authors are supported by MINECO (Ministerio de Economía, Industria y Competitividad), Spain, project MTM2014-56022-C2-1-P.

Appendix

Proof

Proof of Theorem 2.1.

Let us fix ${\mathbf{X}}_0>\mathbf{0}$ and let $k_0$ be sufficiently large so that for every $i=1,\dots ,d$, the ith component of ${\mathbf{X}}_0$ verifies ${\mathbf{X}}_0^i>1/k_0$. For each integer $k\ge k_0$ we define the stopping time ${\tau }_k:=inf\{t\ge 0:X^i(t)\le 1/k\text{ for some }i=1,\dots ,d\}$. Clearly ${\tau }_k$ is increasing as $k\to \mathrm{\infty }$. Let us set ${\tau }_{\mathrm{\infty }}:=\underset{k\to \mathrm{\infty }}{lim}{\tau }_k$. We want to show that ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$ with probability one, and to do so we will proceed by contradiction. Let us assume that the previous statement is false. Then there exists T>0 and $\epsilon \in (0,1)$ such that $\mathbb{P}({\tau }_{\mathrm{\infty }}\le T)>\epsilon$, and therefore there exists an integer $k_1\ge k_0$ such that (A1) $\mathbb{P}({\tau }_k\le T)\ge \epsilon \mathrm{forall}k\ge k_1.$(A1) Let ${\mathring{\mathbf{R}}}_{+}^d=(0,\mathrm{\infty })\times \stackrel{\left(d\right)}{\cdots }\times (0,\mathrm{\infty })$, let $v:{\mathring{\mathbf{R}}}_{+}^d\to \mathbf{R}$ be any $C^2$ function and let $V\left(t\right):=v\left(\mathit{X}\right(t\left)\right)$. From Ito's theorem [1] we have (A2) $\mathrm{d}V\left(t\right)=Lv\left(\mathbf{X}\right(t\left)\right)\mathrm{d}t+\frac{\mathrm{\partial }v}{\mathrm{\partial }x}\left(\mathbf{X}\right(t\left)\right)\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),$(A2) where $\begin{array}{rl}Lv\left(\mathbf{x}\right)& =\stackrel{d}{\sum _{i,j=1}}\frac{\mathrm{\partial }v}{\mathrm{\partial }{x}_i}{a}_{ij}{x}_j+\frac{1}{2}\stackrel{m}{\sum _{h=1}}\stackrel{d}{\sum _{i,j=1}}\left({\mathbf{B}}^h{\mathbf{x}\mathbf{x}}^{\mathrm{T}}\right({\mathbf{B}}^h{)}^{\mathrm{T}}{)}_{ij}{x}_i{x}_j\frac{{\mathrm{\partial }}^{2}v}{\mathrm{\partial }{x}_i\mathrm{\partial }{x}_j}\\ & =\stackrel{d}{\sum _{i,j=1}}\frac{\mathrm{\partial }v}{\mathrm{\partial }{x}_i}{a}_{ij}{x}_j+\frac{1}{2}\stackrel{m}{\sum _{h=1}}\stackrel{d}{\sum _{i=1}}{b}_i^h{b}_j^h{x}_i{x}_j\frac{{\mathrm{\partial }}^{2}v}{\mathrm{\partial }{x}_i\mathrm{\partial }{x}_j}\end{array}$ and we have used that the ${\mathbf{B}}^h$ are diagonal, ${\mathbf{B}}^h=\mathrm{diag}(b_1^h,\dots ,b_d^h)$. Let us choose $v\left(\mathbf{x}\right):=-\sum _{i=1}^d\log x_i/(1+x_i)$, which is positive and $C^2\left({\mathring{\mathbf{R}}}_{+}^d\right)$. Straightforward calculations show that (A3) $Lv\left(\mathbf{x}\right)=-\stackrel{d}{\sum _{i=1}}\frac{a_{ii}}{1+x_i}-\stackrel{d}{\sum _{i=1}}\sum _{j\ne i}\frac{a_{ij}{x}_j}{1+x_i}+\frac{1}{2}\stackrel{m}{\sum _{h=1}}\stackrel{d}{\sum _{i=1}}(b_i^h{)}^2\frac{1+2x_i}{(1+x_i{)}^2}.$(A3) Note that, if $x_i>0$ for all $i=1,\dots ,d$, then the first and the third term in Equation (A3(A3) $Lv\left(\mathbf{x}\right)=-\stackrel{d}{\sum _{i=1}}\frac{a_{ii}}{1+x_i}-\stackrel{d}{\sum _{i=1}}\sum _{j\ne i}\frac{a_{ij}{x}_j}{1+x_i}+\frac{1}{2}\stackrel{m}{\sum _{h=1}}\stackrel{d}{\sum _{i=1}}(b_i^h{)}^2\frac{1+2x_i}{(1+x_i{)}^2}.$(A3) ) are bounded from above in ${\mathring{\mathbf{R}}}_{+}^d$, and moreover, since $\mathbf{A}$ is a Metzler matrix, the second term is negative in ${\mathring{\mathbf{R}}}_{+}^d$. Therefore there exists $K\ge 0$ such that $Lv\left(\mathbf{x}\right)\le K$ in ${\mathring{\mathbf{R}}}_{+}^d$. Clearly $\mathbf{X}(t\wedge {\tau }_k)$, where we define $a\wedge b:=min\{a,b\}$, belongs to ${\mathring{\mathbf{R}}}_{+}^d$ for all $t\ge 0,$ and therefore, using the previous bound we have from Equation (A2(A2) $\mathrm{d}V\left(t\right)=Lv\left(\mathbf{X}\right(t\left)\right)\mathrm{d}t+\frac{\mathrm{\partial }v}{\mathrm{\partial }x}\left(\mathbf{X}\right(t\left)\right)\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right),$(A2) ) $V\left(\mathbf{X}\right(T\wedge {\tau }_k\left)\right)\le V\left({\mathbf{X}}_0\right)+K(T\wedge {\tau }_k)+{\int }_0^{T\wedge {\tau }_k}\frac{\mathrm{\partial }v}{\mathrm{\partial }x}\left(\mathbf{X}\right(t\left)\right)\stackrel{m}{\sum _{h=1}}{\mathbf{B}}^h\mathbf{X}\left(t\right)\mathrm{d}{W}^h\left(t\right)$ and therefore, using that the expected value of the stochastic integral is zero (A4) $\mathbb{E}V\left(\mathbf{X}\right(T\wedge {\tau }_k\left)\right)\le V\left({\mathbf{X}}_0\right)+K\mathbb{E}(T\wedge {\tau }_k)\le V\left(\mathbf{X}\right(0\left)\right)+KT.$(A4) On the other hand, let us define ${\mathrm{\Omega }}_k=\{{\tau }_k\le T\}$ for $k\ge k_1$, so that Equation (A1) means that $\mathbb{P}\left({\mathrm{\Omega }}_k\right)\ge \epsilon$. Now for all $\omega \in {\mathrm{\Omega }}_k$ there exists some $i=1,\dots ,d$ such that $X^i(t,\omega )=1/k$ and, since v is a sum of positive summands, $V\left(\mathbf{X}\right({\tau }_k,\omega \left)\right)\ge -\log \left(1/k\right)/(1+1/k)$. Since V is non-negative, if we denote by $1_{{\mathrm{\Omega }}_k}$ the indicator function of ${\mathrm{\Omega }}_k$ we can write $\mathbb{E}V\left(\mathbf{X}\right(T\wedge {\tau }_k\left)\right)\ge \mathbb{E}\left[1_{{\mathrm{\Omega }}_k}V\right(\mathbf{X}\left({\tau }_k\right)\left)\right]\ge -\epsilon \log \left(1/k\right)/(1+1/k)$ and taking the limit $k\to \mathrm{\infty }$ we obtain $\underset{k\to \mathrm{\infty }}{lim}\mathbb{E}V\left(\mathbf{X}\right(T\wedge {\tau }_k\left)\right)=\mathrm{\infty }$ which clearly contradicts (A4(A4) $\mathbb{E}V\left(\mathbf{X}\right(T\wedge {\tau }_k\left)\right)\le V\left({\mathbf{X}}_0\right)+K\mathbb{E}(T\wedge {\tau }_k)\le V\left(\mathbf{X}\right(0\left)\right)+KT.$(A4) ) and therefore it must be ${\tau }_{\mathrm{\infty }}=\mathrm{\infty }$ a.s. as we wanted to prove.

We state as a Theorem a collection of results from Baumgärtel [4] about perturbation of semisimple eigenvalues of matrices that we will use in the proof of Theorem 6.1:

Theorem A.1

Essentially Theorems 2 and 3 in [4], pp. 267–269

Let E be a complex linear space of finite dimension N and let $\mathbf{A}\left(\epsilon \right)$ be a linear operator defined on E that admits the following holomorphic expansion $\mathbf{A}\left(\epsilon \right)={\mathbf{A}}_0+\epsilon {\mathbf{A}}_1+{\epsilon }^2{\mathbf{A}}_2+\cdots ,\left|\epsilon \right|<R.$ Let ${\lambda }_0$ be a semi-simple eigenvalue of ${\mathbf{A}}_0$ and let $\mathbf{P}$ be its associated eigenprojection.

1. Let $\lambda \left(\epsilon \right)$ be a perturbed eigenvalue such that $\lambda \left(\epsilon \right)\to {\lambda }_0$ when $\epsilon \to 0$. Then $\lambda \left(\epsilon \right)$ has the form $\lambda \left(\epsilon \right)={\lambda }_0+{\mu }_j\epsilon +o\left(\epsilon \right),\epsilon \to 0,$ for some $j=1,\dots ,q,$ where ${\mu }_j,j=1,\dots ,q$ are the eigenvalues of the restriction $\mathbf{T}{\mid }_{\mathrm{Im}\mathbf{P}}$ of operator $\mathbf{T}:={\mathbf{P}\mathbf{A}}_1\mathbf{P}$ to $\mathrm{Im}\mathbf{P}$. Moreover, associated to $\lambda \left(\epsilon \right)$ we can choose a $($right$)$ eigenvector $\mathbf{x}\left(\epsilon \right)$ that is continuous at $\epsilon =0$ and such that ${\mathbf{x}}_0:=\mathbf{x}\left(0\right)\ne 0$ verifies ${\mathbf{T}\mathbf{x}}_0={\mu }_j{\mathbf{x}}_0$.

2. In the case that μ is a simple eigenvalue of $\mathbf{T},$ $\lambda \left(\epsilon \right)$ and $\mathbf{x}\left(\epsilon \right)$ are holomorphic functions of ϵ.

Proof

Proof of Theorem 6.1

For small ϵ we will consider matrix $\mathbf{C}\left(\epsilon \right)=\mathbf{F}+\epsilon \mathbf{S}$ as a perturbation of matrix $\mathbf{F}$. From Equation (4(4) $\mathbf{F}=\mathrm{diag}({\mathbf{F}}_1,{\mathbf{F}}_2,\dots ,{\mathbf{F}}_q).$(4) ) and the properties of each matrix ${\mathbf{F}}_i$, we have that 0 is a semisimple eigenvalue of $\mathbf{F}$ and that it is strictly dominant. Therefore (A5) $\mathbf{F}=(\mathbf{V}\mid {\mathbf{V}}^{\mathrm{\prime }})\mathrm{diag}({\mathbf{0}}_q,{\mathbf{\Sigma }}^{\mathrm{\prime }})\left(\frac{\mathbf{U}}{{\mathbf{U}}^{\mathrm{\prime }}}\right)$(A5) is a Jordan canonical decomposition of $\mathbf{F}$ [12] where ${\mathbf{V}}^{\mathrm{\prime }}\in {\mathbb{R}}^{N\times (N-q)}$, ${\mathbf{U}}^{\mathrm{\prime }}\in {\mathbb{R}}^{N\times (N-q)}$ and ${\mathbf{\Sigma }}^{\mathrm{\prime }}\in {\mathbb{R}}^{(N-q)\times (N-q)}$ are certain matrices and ${\mathbf{\Sigma }}^{\mathrm{\prime }}$ is associated to eigenvalues with strictly negative real part. From Equation (A5(A5) $\mathbf{F}=(\mathbf{V}\mid {\mathbf{V}}^{\mathrm{\prime }})\mathrm{diag}({\mathbf{0}}_q,{\mathbf{\Sigma }}^{\mathrm{\prime }})\left(\frac{\mathbf{U}}{{\mathbf{U}}^{\mathrm{\prime }}}\right)$(A5) ) and Lemma 4.1, we have that $\hat{\mathbf{F}}=\mathbf{V}\mathbf{U}$ is the eigenprojection matrix of $\mathbf{F}$ associated to eigenvalue 0.

Using the continuity of eigenvalues on matrix entries, the dominant eigenvalue of $\mathbf{C}\left(\epsilon \right)$ for small enough ϵ necessarily must belong to the set of eigenvalues of $\mathbf{C}\left(\epsilon \right)$ that tend to 0 when $\epsilon \to 0$. Using part 1 of Theorem A.1 applied to $\mathbf{C}\left(\epsilon \right)=\mathbf{F}+\epsilon \mathbf{S}$ with ${\mathbf{A}}_0=\mathbf{F}$, ${\mathbf{A}}_1=\mathbf{S},$ ${\lambda }_0=0$ and $\mathbf{P}=\hat{\mathbf{F}}$, we conclude that this set of eigenvalues has the form (A6) $0+\epsilon {\mu }_j+o\left(\epsilon \right),j=1,\dots ,q,$(A6) where the ${\mu }_j$ are the eigenvalues of the restriction $\mathbf{T}{\mid }_{\mathrm{Im}\hat{\mathbf{F}}}$ of $\mathbf{T}:=\hat{\mathbf{F}}\mathbf{S}\hat{\mathbf{F}}$ to $\mathrm{Im}\hat{\mathbf{F}}$. Clearly $\mathrm{Im}\hat{\mathbf{F}}$ is invariant for $\mathbf{T}$ and, using Equation (11(11) $\begin{array}{rl}\overline{\mathbf{S}}& =\left[{\overline{s}}_{ij}\right]:=\mathbf{U}\mathbf{S}\mathbf{V}=\left[\begin{array}{ccc}{\mathbf{u}}_1^{\mathrm{T}}{\mathbf{S}}_{11}{\mathbf{v}}_1& \cdots & {\mathbf{u}}_1^{\mathrm{T}}{\mathbf{S}}_{1q}{\mathbf{v}}_q\\ ⋮& \ddots & ⋮\\ {\mathbf{u}}_q^{\mathrm{T}}{\mathbf{S}}_{q1}{\mathbf{v}}_1& \cdots & {\mathbf{u}}_q^{\mathrm{T}}{\mathbf{S}}_{qq}{\mathbf{v}}_q\end{array}\right]\in {\mathbb{R}}^{q\times q}\end{array}$(11) ) and Lemma 4.1, we have $\mathbf{T}=\hat{\mathbf{F}}\mathbf{S}\hat{\mathbf{F}}=\mathbf{V}\mathbf{U}\mathbf{S}\mathbf{V}\mathbf{U}=\mathbf{V}\overline{\mathbf{S}}\mathbf{U}$, so that $\overline{\mathbf{S}}$ is the matrix of the restriction $\mathbf{T}{\mid }_{\mathrm{Im}\hat{\mathbf{F}}}$ expressed in the basis of $\mathrm{Im}\hat{\mathbf{F}}$ defined by the columns of $\mathbf{V}$. Therefore, $\mathbf{T}{\mid }_{\mathrm{Im}\hat{\mathbf{F}}}$ and $\overline{\mathbf{S}}$ have the same eigenvalues and, moreover, if $\mathbf{x}$ and $\mathbf{y}$ are, respectively, right and left eigenvectors of $\overline{\mathbf{S}}$ associated to a certain eigenvalue λ, then $\mathbf{V}\mathbf{x}\ne \mathbf{0}$ and ${\mathbf{U}}^{\mathrm{T}}\mathbf{y}\ne \mathbf{0}$ are right and left eigenvectors of $\mathbf{T}{\mid }_{\mathrm{Im}\hat{\mathbf{F}}}$ associated to eigenvalue λ.

Therefore, using Hypothesis 6.1 it follows that μ is a simple eigenvalue for $\mathbf{T}{\mid }_{\mathrm{Im}\hat{\mathbf{F}}}$ and is associated to right and left eigenvectors $\mathbf{V}\mathbf{r}\ne \mathbf{0}$ and ${\mathbf{U}}^{\mathrm{T}}\mathbf{l}\ne \mathbf{0}$, respectively. Moreover, μ is strictly dominant for $\mathbf{T}{\mid }_{\mathrm{Im}\hat{\mathbf{F}}}$ so that, using Equation (A6(A6) $0+\epsilon {\mu }_j+o\left(\epsilon \right),j=1,\dots ,q,$(A6) ), $\mu \left(\epsilon \right):=\epsilon \mu +o\left(\epsilon \right)$ is a simple eigenvalue for $\mathbf{C}\left(\epsilon \right)$ and strictly dominant when ϵ is small enough.

Now, using part 2 of Theorem A.1 and the fact that μ is simple for $\mathbf{T}{\mid }_{\mathrm{Im}\hat{\mathbf{F}}}$, it follows that $\mu \left(\epsilon \right)$ has the form $\mu \left(\epsilon \right):=\epsilon \mu +O\left({\epsilon }^2\right)$ and that is associated to right and left eigenvectors which can be written in the form $\mathbf{V}\mathbf{r}+\mathbf{O}\left(\epsilon \right)$ and ${\mathbf{U}}^{\mathrm{T}}\mathbf{l}+\mathbf{O}\left(\epsilon \right)$ respectively. Then, for small enough ϵ $\underset{\tau \to \mathrm{\infty }}{lim}\frac{\mathbb{E}\mathbf{X}\left(\tau \right)}{\exp (\epsilon \mu +O({\epsilon }^2\left)\right)}=\frac{({\mathbf{l}}^{\mathrm{T}}\mathbf{U}+\mathbf{O}(\epsilon \left)\right){\mathbf{X}}_0}{({\mathbf{l}}^{\mathrm{T}}\mathbf{U}+\mathbf{O}(\epsilon \left)\right)({\mathbf{V}\mathbf{r}}^{\mathrm{T}}+\mathbf{O}(\epsilon \left)\right)}(\mathbf{V}\mathbf{r}+\mathbf{O}(\epsilon \left)\right)=\frac{⟨\mathbf{l},{\mathbf{U}\mathbf{X}}_0⟩}{⟨\mathbf{l},\mathbf{r}⟩}\mathbf{V}\mathbf{r}+\mathbf{O}\left(\epsilon \right)$ and so Equation (20(20) $\underset{\tau \to \mathrm{\infty }}{lim}\frac{\mathbb{E}\mathbf{X}\left(\tau \right)}{\exp (\epsilon \mu +O({\epsilon }^2\left)\right)}=\frac{⟨\mathbf{l},{\mathbf{U}\mathbf{X}}_0⟩}{⟨\mathbf{l},\mathbf{r}⟩}\mathbf{V}\mathbf{r}+\mathbf{O}\left(\epsilon \right).$(20) ) is proved.

(b) The proof is analogous to that of part (a). A notable property of the Kronecker product that we will use in the sequel is [9, Theorem 6.1] $({\mathbf{A}}_1\otimes {\mathbf{A}}_2\otimes \cdots \otimes {\mathbf{A}}_m)({\mathbf{B}}_1\otimes {\mathbf{B}}_2\otimes \cdots \otimes {\mathbf{B}}_m)=\left({\mathbf{A}}_1{\mathbf{B}}_1\right)\otimes \left({\mathbf{A}}_2{\mathbf{B}}_2\right)\otimes \cdots \otimes \left({\mathbf{A}}_m{\mathbf{B}}_m\right).$ We will use Equation (19(19) $\mathbf{C}\left(\epsilon \right):=\mathbf{F}+\epsilon \mathbf{S},{\mathbf{C}}_2\left(\epsilon \right):={\mathbf{F}}_2+\epsilon {\mathbf{S}}_2,$(19) ) and for small ϵ will consider matrix $\mathbf{C}\left(\epsilon \right)$ as a perturbation of matrix ${\mathbf{F}}_2={\mathbf{I}}_N\otimes \mathbf{F}\mathbf{+}\mathbf{F}\otimes {\mathbf{I}}_{\mathbf{N}}$. First, let us study the spectral properties of ${\mathbf{F}}_2$. Let ${\hat{\mathbf{v}}}_i$, $i=1,\dots ,N$, denote the ith. column of matrix $(\mathbf{V}\mid {\mathbf{V}}^{\mathrm{\prime }})$, that is, ${\hat{\mathbf{v}}}_i$ is an eigenvector of $\mathbf{F}$ associated to eigenvalue 0 if $1\le i\le q$ or to an eigenvalue with strictly negative real part if $q+1\le i\le N$.

From [9, Theorem 6.5] we know that the eigenvalues of ${\mathbf{I}}_N\otimes {\mathbf{F}\mathbf{+}\mathbf{F}\otimes \mathbf{I}}_N$ are the set $\{{\lambda }_i+{\lambda }_j:{\lambda }_i\in \sigma (\mathbf{F}),i=1,\dots ,N\}$ and that, associated to each eigenvalue ${\lambda }_i+{\lambda }_j$ we have right and left eigenvectors ${\hat{\mathbf{v}}}_i\otimes {\hat{\mathbf{v}}}_j$ and ${\hat{\mathbf{u}}}_i\otimes {\hat{\mathbf{u}}}_j$ where ${\hat{\mathbf{v}}}_i$ and ${\hat{\mathbf{u}}}_i$ are, respectively, right and left eigenvectors of $\mathbf{F}$ associated to eigenvalue ${\lambda }_i$, $i,j=1,\dots ,N$. Taking into account the spectrum of $\mathbf{F}$ we conclude that $\sigma \left({\mathbf{F}}_2\right)$ consists of eigenvalue 0 with multiplicity $q^2$ associated to right (resp. left) left eigenvectors ${\hat{\mathbf{v}}}_i\otimes {\hat{\mathbf{v}}}_j$, (resp. ${\hat{\mathbf{u}}}_i\otimes {\hat{\mathbf{u}}}_j$) $i,j=1,\dots ,q$, and $N^2-q^2$ eigenvalues with strictly negative real parts associated to right (resp. left) eigenvectors ${\hat{\mathbf{v}}}_i\otimes {\hat{\mathbf{v}}}_j$ (resp. ${\hat{\mathbf{u}}}_i\otimes {\hat{\mathbf{u}}}_j$), where either i or j belong to the set $\{q+1,\dots ,N\}$. Clearly, eigenvalue $\lambda =0$ is semisimple and is strictly dominant for ${\mathbf{F}}_2$.

It is easy to check that if $\mathbf{A}$ and $\mathbf{B}$ are any matrices with the block structure $\mathbf{A}=({\mathbf{a}}_1\mid \cdots \mid {\mathbf{a}}_k)$, $\mathbf{B}=({\mathbf{B}}_1\mid \cdots \mid {\mathbf{B}}_s)$ where the ${\mathbf{a}}_i$ are column vectors, then $\mathbf{A}\otimes \mathbf{B}=({\mathbf{a}}_1\otimes {\mathbf{B}}_1\mid \cdots \mid {\mathbf{a}}_1\otimes {\mathbf{B}}_s\mid \cdots \mid {\mathbf{a}}_k\otimes {\mathbf{B}}_1\mid \cdots \mid {\mathbf{a}}_k\otimes {\mathbf{B}}_s).$ Applying this result we can write $\begin{array}{rl}({\hat{\mathbf{v}}}_1\otimes {\hat{\mathbf{v}}}_1\mid \cdots \mid {\hat{\mathbf{v}}}_1\otimes {\hat{\mathbf{v}}}_q\mid \cdots \mid {\hat{\mathbf{v}}}_q\otimes {\hat{\mathbf{v}}}_1\mid \cdots \mid {\hat{\mathbf{v}}}_q\otimes {\hat{\mathbf{v}}}_q)& =\mathbf{V}\otimes \mathbf{V}\\ ({\hat{\mathbf{u}}}_1\otimes {\hat{\mathbf{u}}}_1\mid \cdots \mid {\hat{\mathbf{u}}}_1\otimes {\hat{\mathbf{u}}}_q\mid \cdots \mid {\hat{\mathbf{u}}}_q\otimes {\hat{\mathbf{u}}}_1\mid \cdots \mid {\hat{\mathbf{u}}}_q\otimes {\hat{\mathbf{u}}}_q)& ={\mathbf{U}}^{\mathrm{T}}\otimes {\mathbf{U}}^{\mathrm{T}}\end{array}$ so that we can write a Jordan canonical decomposition of ${\mathbf{F}}_2$ in the form ${\mathbf{F}}_2=(\mathbf{V}\otimes \mathbf{V}\mid {\mathbf{V}}_2^{\mathrm{\prime }})\mathrm{diag}({\mathbf{0}}_{q^2},{\mathbf{\Sigma }}_2^{\mathrm{\prime }})\left(\frac{\mathbf{U}\otimes \mathbf{U}}{{\mathbf{U}}_2^{\mathrm{\prime }}}\right),$ where ${\mathbf{\Sigma }}_2^{\mathrm{\prime }}$, ${\mathbf{V}}_2^{\mathrm{\prime }}$ and ${\mathbf{U}}_2^{\mathrm{\prime }}$ are appropriate matrices. From this expression we see that ${\hat{\mathbf{F}}}_2:=(\mathbf{V}\otimes \mathbf{V})(\mathbf{U}\otimes \mathbf{U})=\left(\mathbf{V}\mathbf{U}\right)\otimes \left(\mathbf{V}\mathbf{U}\right)=\hat{\mathbf{F}}\otimes \hat{\mathbf{F}}$ is the eigenprojection of ${\mathbf{F}}_2$ corresponding to eigenvalue 0.

Like in part (a), we know that, for small enough ϵ, the dominant eigenvalue of ${\mathbf{C}}_2\left(\epsilon \right)$ necessarily must belong to the set of eigenvalues of ${\mathbf{C}}_2\left(\epsilon \right)$ that tend to 0 when $\epsilon \to 0$. Let us consider the restriction ${\mathbf{T}}_2{\mid }_{\mathrm{Im}{\hat{\mathbf{F}}}_2}$of matrix ${\mathbf{T}}_2:={\hat{\mathbf{F}}}_2{\mathbf{S}}_2{\hat{\mathbf{F}}}_2$ to $\mathrm{Im}{\hat{\mathbf{F}}}_2$. Clearly $\mathrm{Im}{\hat{\mathbf{F}}}_2$ is invariant for ${\mathbf{T}}_2$ and, using Equation (11(11) $\begin{array}{rl}\overline{\mathbf{S}}& =\left[{\overline{s}}_{ij}\right]:=\mathbf{U}\mathbf{S}\mathbf{V}=\left[\begin{array}{ccc}{\mathbf{u}}_1^{\mathrm{T}}{\mathbf{S}}_{11}{\mathbf{v}}_1& \cdots & {\mathbf{u}}_1^{\mathrm{T}}{\mathbf{S}}_{1q}{\mathbf{v}}_q\\ ⋮& \ddots & ⋮\\ {\mathbf{u}}_q^{\mathrm{T}}{\mathbf{S}}_{q1}{\mathbf{v}}_1& \cdots & {\mathbf{u}}_q^{\mathrm{T}}{\mathbf{S}}_{qq}{\mathbf{v}}_q\end{array}\right]\in {\mathbb{R}}^{q\times q}\end{array}$(11) ), Lemma 4.1 and the properties of the Kronecker product we have ${\mathbf{T}}_2={\hat{\mathbf{F}}}_2{\mathbf{S}}_2{\hat{\mathbf{F}}}_2=(\mathbf{V}\otimes \mathbf{V}){\overline{\mathbf{S}}}_2(\mathbf{U}\otimes \mathbf{U})$ so that ${\overline{\mathbf{S}}}_2$ is the matrix of ${\mathbf{T}}_2{\mid }_{\mathrm{Im}{\hat{\mathbf{F}}}_2}$ expressed in the basis of $\mathrm{Im}{\hat{\mathbf{F}}}_2$ defined by the columns of $\mathbf{V}\otimes \mathbf{V}$. Therefore, ${\mathbf{T}}_2{\mid }_{\mathrm{Im}{\hat{\mathbf{F}}}_2}$ and ${\overline{\mathbf{S}}}_2$ have the same eigenvalues and, moreover, if $\mathbf{x}$ and $\mathbf{y}$ are, respectively, right and left eigenvectors of ${\overline{\mathbf{S}}}_2$ associated to a certain eigenvalue λ, then $(\mathbf{V}\otimes \mathbf{V})\mathbf{x}\ne \mathbf{0}$ and $(\mathbf{U}\otimes \mathbf{U}{)}^{\mathrm{T}}\mathbf{y}\ne \mathbf{0}$ are right and left eigenvectors of ${\mathbf{T}}_2{\mid }_{\mathrm{Im}{\hat{\mathbf{F}}}_2}$ associated to eigenvalue λ.

Now let us assume Hypothesis 6.2. Then we can reason exactly as we did in part (a), substituting $\mathbf{C}\left(\epsilon \right)$, $\mathbf{F}$, $\mathbf{S}$, $\overline{\mathbf{S}}$, $\hat{\mathbf{F}}$, $\mathbf{V}$ and $\mathbf{U}$ for ${\mathbf{C}}_2\left(\epsilon \right)$, ${\mathbf{F}}_2$, ${\mathbf{S}}_2$, ${\overline{\mathbf{S}}}_2$, ${\hat{\mathbf{F}}}_2$, $\mathbf{V}\otimes \mathbf{V}$ and $\mathbf{U}\otimes \mathbf{U}\mathbf{,}$ respectively, and obtain that for small enough ϵ, ${\mathbf{C}}_2\left(\epsilon \right)$ has a simple and strictly dominant eigenvalue that can be written in the form (21(21) ${\lambda }_2\left(\epsilon \right)=\epsilon {\mu }_2+O\left({\epsilon }^2\right)$(21) ) and associated to eigenvectors with the form (22(22) ${\mathbf{r}}_2\left(\epsilon \right)=(\mathbf{V}\otimes \mathbf{V}){\mathbf{r}}_2+\mathbf{O}\left(\epsilon \right);{\mathbf{l}}_2\left(\epsilon \right)=(\mathbf{U}\otimes \mathbf{U}{)}^{\mathrm{T}}{\mathbf{l}}_2+\mathbf{O}(\epsilon ).$(22) ), from where the asymptotic behaviour (23(23) $\underset{\tau \to \mathrm{\infty }}{lim}\frac{\mathrm{vec}\mathbb{E}\left(\mathbf{X}\right(t\left)\mathbf{X}\right(t{)}^{\mathrm{T}})}{\exp (\epsilon {\mu }_2+O({\epsilon }^2\left)\right)}=\frac{⟨{\mathbf{l}}_2,\left(\mathbf{U}\otimes \mathbf{U}\right)\mathrm{vec}\left({\mathbf{X}}_0{\mathbf{X}}_0^{\mathrm{T}}\right)⟩}{⟨{\mathbf{l}}_2,{\mathbf{r}}_2⟩}\left(\mathbf{V}\otimes \mathbf{V}\right){\mathbf{r}}_2+\mathbf{O}\left(\epsilon \right).$(23) ) follows.