Global stability of an age-structured population model on several temporally variable patches

We consider an age-structured density-dependent population model on several temporally variable patches. There are two key assumptions on which we base model setup and analysis. First, intraspecific competition is limited to competition between individuals of the same age (pure intra-cohort competition) and it affects density-dependent mortality. Second, dispersal between patches ensures that each patch can be reached from every other patch, directly or through several intermediary patches, within individual reproductive age. Using strong monotonicity we prove existence and uniqueness of solution and analyze its large-time behavior in cases of constant, periodically variable and irregularly variable environment. In analogy to the next generation operator, we introduce the net reproductive operator and the basic reproduction number \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document}R0 for time-independent and periodical models and establish the permanence dichotomy: if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0\le 1$$\end{document}R0≤1, extinction on all patches is imminent, and if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0>1$$\end{document}R0>1, permanence on all patches is guaranteed. We show that a solution for the general time-dependent problem can be bounded by above and below by solutions to the associated periodic problems. Using two-side estimates, we establish uniform boundedness and uniform persistence of a solution for the general time-dependent problem and describe its asymptotic behaviour.


Introduction
Population permanence in a patchy environment is a result of complex interactions between abiotic, biotic, and anthropogenic factors (Lewis et al. 2017). Each of these Extended author information available on the last page of the article factors have relative importance for population growth with effects that may differ for terrestrial and aquatic species, large and small populations, plants and animals, vertebrates and invertebrates, Bjørnstad and Grenfell (2001), Kareiva and Wennergren (1995) and Roughgarden (1979). Understanding how interplay of spatial heterogeneity and temporal variability of the habitat, migration patterns, density-dependent feedback and age or class distribution of a population affect its dynamics and permanence is critical for science, conservation and management of biodiversity.
Mathematical models are often used for theoretical investigation of population dynamics and establishing conditions for population permanence. Despite spatial heterogeneity, temporal variability and variation in organization and functioning of natural ecosystems affected by human actions, population models typically emphasize importance of certain factors, while neglecting others. Some of the models highlight effects of age-structure and density-or time-dependence on population dynamics, Chipot (1983), Chipot (1984), Diekmann et al. (2001), von Foerster (1959, Gurtin and Mac-Camy (1974), Iannelli (1995), Iannelli and Pugliese (2014), Kozlov et al. (2016bKozlov et al. ( , 2017 and Webb (1985). Iannelli and Milner (2017) presented a comprehensive study of age-structured population models and another significant contribution to this topic has recently been made by Inaba (2017). Age-structured models for two-sex populations have been studied in Iannelli et al. (2005). Some authors investigated N -species population with age-specific interactions and established the well-posedness and the existence of an equilibrium solution using the theory of semilinear evolution equations and derived local or asymptotic stability results, Prüß (1981) and Prüss (1983).
Spread of disease and epidemic modeling is an area where age-structured population dynamics plays a significant role. Recognition that transmission dynamics of some diseases cannot be explained by traditional models without age-structure sparked interest in analyzing the effects of age-structure on the dynamics of epidemics. These studies are sometimes limited to local results and investigating endemic equilibrium, such as Busenberg et al. (1988) and Castillo-Chavez et al. (1989), but there are also results related to global behavior of age-structured epidemic models, Busenberg et al. (1991), Feng et al. (2005), Kuniya and Iannelli (2014) and Kuniya et al. (2018).
A common feature of the above mentioned models is that they rarely include spatial heterogeneity as a factor of population growth. Models that do consider effects of habitat's spatial structure on population permanence assume either discrete or continuous space. Discrete spatial structure means that a habitat consists of several distinct patches with birth and death rate being dependent on a patch that individuals occupy. A source is a high-quality patch that yields positive population growth, while a sink is a low-quality patch and it yields negative growth rate. In isolation, a subpopulation on every patch has its own dynamics. Linking patches by dispersal leads to the source-sink dynamics, where all local subpopulations contribute to a unique global dynamics. For populations that inhabit several patches, possibility to move from one patch to another can be crucial for survival. For example, dispersal from a source to a sink can save a local sink subpopulation from extinction through the rescue effect and recolonization (Amarasekare 2004;Dias 1996;Hastings 1993). The influence of spatial heterogeneity in unstructured populations is studied in Allen (1983), Arditi et al. (2015), Cui and Chen (1998), Cui and Chen (2001), DeAngelis et al. (2016), and DeAngelis and Zhang (2014). A trade-off between competition and dispersal is investigated in Amarasekare and Nisbet (2001) and the relation between dispersal pattern and permanence is discussed in Hastings and Botsfor (2006) and Jansen and Yoshimura (1998). Age-structured models on several patches usually assume discrete age classes (immature and adults) and dispersion between a few (two or three) temporally unchangeable patches, So et al. (2001), Takeuchi (1986a), Takeuchi (1986b), Terry (2011) and Weng et al. (2010). Age-structured population growth model in continuous environment typically represent movement as diffusion (Webb 2008).
On the mathematical side, these models have been studied using several main approaches. Integrated semigroup approach is one of them, see e.g. Magal and Ruan (2018), Magal et al. (2019) and Thieme (2009). The Lyapunov function approach has been used to establish (global) stability results in some epidemic models (Chekroun et al. 2019). A universal method for systematic construction of a Lyapunov function does not exists and this method is often of little use for complex dynamical systems; moreover, some Lyapunov functions may provide better answers than others. In contrast to this, monotonicity of certain positive operators and comparison principles have been used for studying global dynamics of ordinary, delay and partial differential equations, see Hirsch and Smith (2003), Smith (1995) and Zhao (2003). Application of these methods to age-structured population models and epidemic models can be found in Busenberg et al. (1991), Diekmann et al. (2001), Kuniya and Iannelli (2014), Kuniya et al. (2018), Magal and Ruan (2018), Magal et al. (2019) and Webb (1985).
In this paper we extend analysis to an age-structured population that inhabits N temporally variable patches connected by dispersal. Two key elements of our study are density-dependent population growth represented by pure intra-cohort competition (see also Kozlov et al. 2016aKozlov et al. , 2017 and age-and time-dependent dispersal rates. From a biological point of view, intra-cohort competition occurs as a result of ageing and growth. Some species, such as certain insects, molluscs and fish, undergo metamorphosis, which is a complete change of physical appearance and structure and thereby also a change of food preferences and habitats. Individuals of different age therefore do not compete for resources, which reduces potential competitors from a whole population to a cohort. The pure intra-cohort competition is a relevant modeling approach when there are ontogenetic shifts. Age difference often correlates to body-size differences, which have effects on interactions between and within populations. Age difference is therefore viewed as a cause of diet and habitat shifts. Typically this is studied in food web theory, Petchey et al. (2008) and Zook et al. (2011), where diets are defined according to body-size, and in the theory of size-structured populations, Ebenman et al. (1996) and Narvaez et al. (2020). Pure intra-cohort competition is at one end of a continuum of how intraspecific competition may act and at the other end one finds the commonly used model of intraspecific competition that do not include structuring the population into, for example, age or body-size. In epidemiological modeling, intra-cohort transmission of a disease has been studied in Iannelli and Milner (2017, Section 1.3.5). The assumption that density-dependence has effects only on mortality rates is based on biological and modeling premises. Common individual responses to intraspecific competition and resource depletion are decreased fecundity and increased mortality. Effects of intraspecific competition on individuals may depend on individuals age, type of resource and type of competition (direct or indirect, interference or exploitation, see Gilad 2008). In some cases, lack of resources has stronger effect of mortality, while in other in affects fecundity more severely. From a modeling perspective high mortality rates of newborns can have the same effect as reduced fertility. To achieve nontrivial stability results and include these considerations, we assume that competition between individuals affects only mortality rates and that it occurs only between individuals of the same age (pure intra-cohort competition).
Some biological studies indicate that there are many different causes for dispersal, such as response to environmental conditions, prevention of inbreeding or competition for mates, and that migration can have different forms, such as one way dispersal or a round trip from a birthplace, Bowler and Benton (2005) and Dingle and Drake (2007). This may lead to differences in life-history traits, genetics and demography between dispersers and residents. Demographic studies show that dispersing females are often young individuals in their reproductive age, while old individuals usually do not engage in breeding dispersal, Gaines (1980) and Greenwood and Harvey (1982).
Taking into account these characteristics of density-dependent growth and migration, we formulate model as follow. Let n k (a, t) denote the age distribution in the population patch k at time t with the corresponding birth rate m k (a, t) and the initial distribution of population f k (a). A local subpopulation on each patch experiences intraspecific competition, which results in density-dependent mortality. The assumption that competition occurs only between members of the same age class leads to the following McKendrick-von Foerster type balance equations: in the domain subject to the birth law and the initial age distribution with n(a, t) = (n 1 (a, t), . . . , n N (a, t)) t , a, t), . . . , m N (a, t)), M(n(a, t), a, t) = diag(M 1 (n 1 (a, t), a, t), . . . , M N (n 1 (a, t), a, t)), where M k (v k , a, t) is the density-dependent mortality rate of a population on patch k and dispersion matrix D(a, t) = (D k j (a, t)) 1≤k, j≤N defines migration rates between patches and describes migration pattern.
For D(a, t) ≡ 0, there is no migration between patches and the system (1) splits into N independent balance equations. This model, under an additional assumption that M(a, t) is the logistic regulatory function (8), has recently been studied in Kozlov et al. (2017). The case D(a, t) ≡ 0 is much more challenging. The coefficients D k j (a, t), j = k, define a proportion of individuals of age a at time t on patch j that migrates to patch k, which implies that D k j (a, t) ≥ 0 for all a, t. Given that migration is costly (in terms of time and energy) and risky for migrating individuals, we define D kk (a, t) as the migration-related mortality of individuals of age a at time t on patch k that is independent of the density-dependent mortality M k (n k (a, t), a, t). The fact that individuals can disperse and move from one patch to another is important in modeling source-sink dynamics as it can explain effects of recolonization or population persistence in heterogeneous environment. It can be expected that the sign pattern and the weighted graph associated with D(a, t) affect global and asymptotic behaviour of solutions to (1)-(4).
In this paper, we provide a mathematical derivation of the results about the existence and uniqueness of a solution to the multi-patch model (1)-(4). The common point for the single-patch age-structured models is that the basic reproduction number R 0 and the characteristic equation are used to determine permanence of a population. The basic reproduction number represents expected number of offspring of an individual in constant environment. Its biological interpretation and computation in periodic environment have been discussed in Bacaër and Dads (2012) and Diekmann et al. (1990), respectively. In age-structured epidemic models, R 0 is the spectral radius of the next generation operator (Inaba 2019;Thieme 2003). The basic reproduction number is a threshold value that indicate long time behavior of a solution: if R 0 < 1, a population faces extinction, and if R 0 > 1, population persistence is granted, see Iannelli and Milner (2017), Inaba (2017, Section 1.4). Similar results hold for epidemic models, where the former condition means that the trivial disease-free equilibrium is globally asymptoticaly stable, while the latter condition implies global stability of a nontrivial equilibrium in which infected individuals persist (Chekroun et al. 2019). These results lead us to two key questions related to the multi-patch model (1)-(4): -Is it possible to define an analogue of the characteristic equation and the basic reproduction number for the multi-patch model? -If so, can they be used for the analysis of the large-time behavior of the solution and for establishing the condition for populations permanence?
The main contribution of the paper lies in a rigorous proof that both questions have affirmative answers in constant, periodic and general time-dependent cases. Similar result for time-independent case can be found in Thieme (2009, Section 6). A net reproductive dichotomy for an age-structured epidemic model in terms of disease persistence can be found in Thieme (2003, Sec.22.3). Our approach relies on the lower and upper solution technique and essentially uses monotonicity of certain integral operators associated with the balance equations. The method that we develop for the time-dependent cases allows investigation of asymptotic behaviour and global stability of the nonlinear model and considering fluctuations that are not necessarily small in amplitude. The obtained results enable discussion of conservation and management problems, such as improving survival of migrating species and pest control.
Outline A summary of the mathematical framework and our main results are presented in Sect. 2. In Sect. 3 we discuss an auxiliary model and derive some preliminary results on the corresponding lower and upper solutions. In Sect. 4 we prove the existence and uniqueness of a solution to the balance Eqs. (1)-(4) by reducing the original problem to a certain nonlinear integral equation. In Sect. 5 we define the associated characteristic equation and the maximal solution, and establish one of the key results of the paper: the basic reproduction number dichotomy. The remaining part of the paper is dedicated to the study of the asymptotic behavior and stability of the solution. We consider three cases: constant environment (i.e. the time-independent case) in Sect. 5, periodic environment in Sect. 6 and irregularly changing environment (i.e. the general time-dependent case) in Sect. 7. Notations For easy reference we fix some standard notation used throughout the paper. R N + denotes the positive cone {x ∈ R N : x i ≥ 0}. Given x, y ∈ R N we use the standard vector order relation: In particular, if D = D jk is an N × N -matrix we define D jk p for any 1 ≤ p ≤ ∞ in an obvious manner identifying D with an element of R N 2 . Given E ⊂ R N and a continuous function h : E → R, we define

The structure conditions
Before providing the main results, we give a brief summary of the structure conditions imposed on the balanced Eqs.
(1)-(4). We always assume that m(a, t) and D(a, t) are continuous 1 for (a, t) ∈B and M(v, a, t) is a continuous function of (v, a, t) ∈ R ×B. Following Iannelli (1995), we let B(t) > 0 denote the maximal length of life of individuals in population at time t ≥ 0. Then is the total population at time t. Furthermore suppose the following structure conditions hold: is a nonnegative nondecreasing function of v for v ≥ 0, and there exist real numbers μ ∞ > 0, γ > 0, and a function p(a) ≥ μ ∞ such that (H3) D C(B) < ∞ and D(a, t) is a Metzler matrix: (H5) the function f(a) is continuous and supp f ⊂ [0, B(0)).
Let us briefly explain the above conditions from the biological perspective. Concerning (H1), one usually uses a more restrictive condition that B(t) is a constant. Nevertheless, (5) is a more reasonable assumption: it means that the maximal length of life of individuals B(t) in a population may depend on t but it does not grow faster than time. Mathematically, (5) asserts that the boundary curve B(t) is transversal to the characteristics of (1).
The monotonicity assumption in (H2) ensures that increase in age-class density increases the death rate and has a negative effect on population growth. The classical example of the density independent mortality rate M k (v, a, t) = μ k (a, t) ≥ μ ∞ > 0 is compatible with γ = 0 in (H2). Another example is the logistic type model (Kozlov et al. 2017) with where L k (a, t) ∈ L ∞ (B) is the regulatory function (carrying capacity); this example fits (H2) for γ = 1. Concerning the Metzler condition in (H3), note that the dispersion coefficient D k j (a, t) expresses the proportion of population n k (a, t) that from patch j goes to patch k, which naturally yields that D k j ≥ 0. Furthermore, according the support condition in (H4), the improper integral in (3) is well-defined and actually is taken over the finite interval [a m , A m ] which lies within the domain of definition of n(a, t) for any fixed t > 0. The condition (H5) is a natural assumption that the initial distribution of population is bounded by the life length. The accessibility condition For further applications we shall also need an additional assumption on the structure of the dispersion matrix D. In order to formulate it, let us recall some relevant concepts. Given a Metzler matrix A ∈ R N ×N , one can associate a directed graph Γ (A) with nodes labeled by {1, 2, . . . , N } where an arc leads from i to j, i = j, if and only if A i j > 0. The patch j is said to be reachable from i, denoted i j, if there exists a directed path from i to j. A digraph is called connected from vertex i if i j for all j = i (Balakrishnan 1996, p. 132). A patch k is said to be accessible at age a ≥ 0 if the associated digraph Γ (D(a, t)) is connected from k for any t > 0. The accessibility condition relies on the sign pattern of the corresponding dispersion matrix and can be readily obtained by the standard tools of the nonnegative matrix theory (Minc 1988, Section 3). Now, notice that by (H4) the following value is finite: From the biological point of view,ā k is the maximal fertility age in population k. Our last condition reads as follows: (H6) For any 1 ≤ k ≤ N there exists 0 < β k <ā k such that the patch k is accessible at age β k .
In other words, (H6) asserts that for any patch k there is a moment β k > 0 such that a (composite) migration from any other patch j to k is possible within the reproductive period. This condition is based on some biological studies showing that usually young individuals in their reproductive age engage in breeding dispersal, unlike old individuals who migrate less or migrate for other reasons.

The net reproductive rate dichotomy
Let us denote by ρ(t) = n(0, t) the newborn function, i.e. a vector-function whose components denote the number of newborns on each patch. Then, the problem (1)-(4) can be reduced to the integral equation where K and F are positive nondecreasing operators with bounded ranges and Ff(t) = 0 for large t > 0. Our strategy for proving permanence results is as follows: we first establish the permanence results for time-independent and time-periodic coefficients, and then show that in the general situation, a solution of (9) can be well-controlled by these cases. If the environment is constant then the model parameters are time-independent functions. Then it is reasonable to assume that the maximal life-time is constant: Chipot (1983) and Gurtin and MacCamy (1974). Our approach relies on a fine control of large-time behaviour of an arbitrary solution to (9) by nontrivial solutions of the associated characteristic equation where the operatorK is given by the right hand side in (3) for a time-independent solution to (1) with a constant boundary condition n(0) = ρ. Clearly, ρ = 0 is a (trivial) solution of the characteristic equation.
Our goal is to establish when a nontrivial positive solution ρ 0 exists. A crucial tool here is the so-called maximal solution of the characteristic equation, i.e. a solution θ of (10) such that for an arbitrary solution ρ there holds ρ ≤ θ . In particular, θ = 0 implies that the characteristic equation has only trivial solutions. We establish the existence of the maximal solution in Sect. 5.2.
Another important ingredient is the next generation operator (see section 2.1 in Inaba 2017) We show that under conditions (H1)-(H6), R 0 : R N + → R N + is a strongly positive operator. By Perron-Frobenius theorem, its spectral radius R 0 is equal to the largest positive eigenvalue. We call this value the basic reproduction number.
To motivate the latter definition, observe that in the single-patch case, the basic reproduction number R 0 is given by It it related to the solution of the Euler-Lotka characteristic equations in the linear agestructured population model; see Iannelli and Pugliese (2014). According to Kozlov et al. (2017), R 0 is related to the solution ρ * of the characteristic equation in the nonlinear age-structured model. Namely, if R 0 ≤ 1, then ρ * = 0 and the population is going to extinction, while for R 0 > 1, we have ρ * > 0 and the population is permanent. The same is obviously valid if there are several patches without migration (i.e. D ≡ 0): every local subpopulation behaves accordingly to the value of R 0 on the respective patch.
The main contribution of this paper lies in the following dichotomy result on the long-term dynamics of populations.
Theorem A (The Net Reproductive Rate Dichotomy) If R 0 ≤ 1, then θ = 0 and the characteristic Eq. (10) has no nontrivial solutions. If R 0 > 1, then θ 0 and θ is the only nontrivial solution of the characteristic equation.
Let χ(t) be an arbitrary solution of (9). Then: Thus, the basic reproduction number R 0 effectively determines large time behavior of population on N patches in a constant environment. Here, as in the single-patch case, R 0 ≤ 1 implies extinction of a population on all patches, while R 0 > 1 grants the global permanence of a population. We see that the dichotomy result for a multipatch population is completely consistent with the single-patch case when the next generation operator R 0 coincides with the multiplication by R 0 .
It is important to emphasize that the function ϕ(a; θ) in (11) is exactly the unique equilibrium point of the problem (1), (3) provided that θ satisfies the characteristic equation. In other words, Theorem A implies the global stability result: any solution of the principal model converges at infinity to the unique equilibrium point given by the characteristic equation.
The proof of Theorem A, along with certain related results, occupies Sect. 5 and makes an essential use of the auxiliary monotonicity results collected in Sect. 3 and functional theoretic properties of the integral Eq. (9) given in Sect. 4. Our approach relies on the following steps and can be described as follows. First, we associate certain lower and upper monotone sequences to an arbitrary solution χ of (9). The existence of an upper sequence relies on the boundedness of the image of K. The construction of a lower sequence is more tricky and involves certain fine properties of the maximal solution and some previous auxiliary monotonicity results together with the accessibility condition (H6). The main problem here is to control a nonzero asymptotic behaviour of the lower approximants as t → ∞. Next, we show that the large-time behaviour of χ can be well controlled by the limits at infinity of constructed monotone approximants. Furthermore, we are able to identify the common limits as the maximal solution θ . This finally establishes that the constructed sequences converge to the equilibrium point of the original problem. Notice that the monotonicity of the lower and upper approximations is crucial because the convergence established in the first steps is valid only on any bounded interval.

Two-side estimates of R 0 and Â
A life-history trade-off between reproduction and migration has been noted for many species, including migratory birds and some insects (see for example Guerra 2011;Mole and Zera 1993;Schmidt-Wellenburg et al. 2008). This trade-off is caused by energy constraints because both reproduction and migration are energetically costly for organisms. Keeping the assumption that the environment is constant and using the specific form of the balance system, we investigate the consequences of this trade-off.
The fact that individuals do not reproduce during migration is biologically justified and mathematically stated as: The relation (12) between dispersion coefficients represents total migration from patch j toward all other patches and means that some migrants that are leaving patch j will eventually die before reaching patch k, but they will not give birth during migration. This migration related mortality is represented by the nonpositive coefficient D j j (a), such that |D j j (a)| ≥ N k=1,k = j D k j (a). Then, we establish in Sect. 5.6 the following two-side estimates for the basic reproduction number.
Theorem B Under additional assumption that (12) holds we have where m(a) is the maximal birth rate and μ(a) is the minimal death rate on all patches.
In addition, in Proposition 11 below we establish a priori estimates for the basic reproduction number and for the maximal solution θ .

Periodically and irregularly changed environment
Natural habitats are usually positively autocorrelated, see for example Steele (1985). Therefore, the assumption that the vital rates, regulating function and dispersal coefficients are changing periodically with respect to time is a reasonable approximation. In the study of the large-time behavior of a solution to Eq. (9) in a periodically changing environment, the pivotal role belongs to the characteristic equation where the operator K is given by the right hand side of (3) and n(a, t) solves (1) with a periodic boundary condition n(0, t) = ρ(t), ρ(t + T ) = ρ(t). We establish in Sect. 6 that the operator K is absolutely continuous which allows us to extend the methods of Sect. 5 to the periodic case. In particular, the corresponding next generation operator R 0 defined on space of periodic continuous functions is strictly positive and its spectral radius R 0 is equal to the largest eigenvalue. We are also able to establish the corresponding dichotomy result for a periodic environment.
If the environment is changing irregularly, the structure parameters of the principal model (1)-(4) can be estimated from above and below by nonnegative periodic functions. Using these periodic functions as structure parameters for new models, we formulate two associated periodic problems. One of them is the best-case scenario and its solution is an upper bound for the original problem. The other is the worst-case scenario and its solution is a lower bound. In other words, a solution for the general time-dependent problem can be bounded for large values of t by above and below by the solution to the associated periodic problems, as stated in Theorem 7.

Source-sink dynamics
Using the source-sink dynamics it is possible to explain permanence of a population on several patches provided that at least one patch is a source and that all patches are connected by dispersion. In Sect. 8.1 we assume that the environment is constant and consists of several patches. Then it is possible to show that survival of population on both patches is possible provided that emigration from the source is sufficiently small. Furthermore, in Sect. 8.2, we show that permanence is possible even if all patches are sinks provided that dispersion is appropriately chosen. This is especially important for migratory birds, since both of their habitats can be seen as sinks (one because of the low reproduction due to insufficient resources, and the other because of the high mortality in the winter). This example can be related to the results in Jansen and Yoshimura (1998), where a simple model is used for analysis of connection between population permanence and allocation of offspring in a population that lives on several patches. One of the results is that permanence is possible even if all patches are sinks.

Upper and lower solutions
Below we establish some auxiliary monotonicity results for lower and upper solutions to a general system of ordinary differential equations where for almost all w ∈ R N and all x ∈ [0, b). We assume additionally that F satisfies In particular, this implies that w(x) ≡ 0 is a solution of (13). We shall also exploit a weaker version of the concept of irreducibility. More precisely, let F(w, x) = (F 1 (w, x), . . . , F N (w, x)) be continuously differentiable with respect to w and let DF(w, x) := ( ∂ F k (w,x) ∂w j ) denote the corresponding Jacobi matrix.
In this paper, we are mostly interested in the particular case when In this case the condition (15) is trivially satisfied.
. The next lemmas generalize the corresponding facts for the cooperative system (cf. Smith 1995, Remark 1.2) on lower (upper) solutions of (13) with Lipschitzian F. Notice also that our proofs are somewhat different from those given in Smith (1995).
First notice that F satisfies the so-called quasimonotone condition (Hirsch and Smith 2003;Smith 1995).

Lemma 1 If F satisfies the Kamke-Müller condition then u
is absolutely continuous in [0, 1], hence applying by the fundamental theorem of calculus and (14) that as desired.

Lemma 2 Let w(x) be an upper solution of
and Lemma 1 that holds everywhere in the neighbourhood of x 0 . Thus, the claim is proved. We also claim that any upper solution to (13) . By the continuity of w(x), there exists δ such that w(x) 1 < for any |x − T | < δ. Let the set E be defined as above and x ∈ [T , T + δ). Since by (15)

The latter inequality yields
To finish the proof, let us suppose that w j (0) > 0. Since F j (y, x) is locally Lipschitz in y, for any r > 0 there exists C(r ) such that [in virtue of (15)] |F j (y, x)| ≤ C(r ) y 1 for all y ∈ R N and y ≤ r . Let 0 < β < b be chosen arbitrarily and let r = sup x∈ [0,β] where e j is the jth coordinate vector, Lemma 1 and the nonnegativity of w j (x) yield that β], and therefore in the whole interval [0, b).

Lemma 3 Let w(x) be an upper solution of (13) with w(0) > 0 and such that the k-th patch is
. Therefore we may without loss of generality assume that w k (β) = 0. Let us suppose by contradiction that there exists In particular, w k (β) = 0. Since w(0) > 0, there exists j such that w j (0) > 0 and, thus, w j (β) > 0. By the assumption, there exists a directed path k j in the graph Γ (DF(w, β)). Equivalently, there exists a sequence of pair-wise where e i denotes the ith coordinate unit vector in R N . Then Since w(x) is an upper solution of (13), we have by Lemma 1 for j 0 = k that Arguing as in (17) we find It follows from (20), the nonnegativity of (v 0 − v 1 ) i and the partial derivatives (for i = j 0 ) that all summands of the latter sum must vanish. Since the integrands are nonnegative continuous functions, they must vanish identically for t ∈ [0, 1]. In particular, (18) readily implies that (v 0 − v 1 ) j 1 = 0. Thus, w j 1 (β) = 0, and by the above we have w j 1 (β) = 0 Repeating the same argument for the pair ( j 1 , j 2 ) etc. implies w j 2 (β) = 0 etc., thus yielding that w j s (β) = w j (β) = 0, a contradiction follows.

Proposition 1 (Comparison principle) Let u(x) and v(x) be resp. upper and lower solutions to
We have for the corresponding Jacobi matrices i.e. L and L g satisfy simultaneously the Kamke-Müller condition. This readily yields the first claim of the proposition. Now let us assume that u(0) > v(0) and for some k and β

Corollary 1 Let u(x) be an lower (resp. upper) solution to
Proof Follows immediately from the fact that w(x) ≡ 0 is a solution of (13).
Proposition 2 (Existence and Uniqueness) Let (13) satisfy the Kamke-Müller condition and there exists C(F) > 0 such that Then for any Proof By the Cauchy-Peano Existence Theorem, (13) has a unique solution w( be the maximal interval of existence of the solution: b := sup{β > 0 : there exists a solution of (13) on [0, β)}.
We claim that b = b. It suffices to show that a solution w(x) is uniformly bounded on any existence interval [0, β), i.e. there exists M > 0 such that for any β < b To this end, we make a more general assumption, that w(x) is a nonnegative lower solution to (13) on [0, β) and consider In particular, H (x) is locally Lipschitz on [0, β), and thus a.e. differentiable there.
Then for any point of differentiability (21) that Integrating the latter inequality (note that H is absolutely continuous) yields This proves (22). Furthermore, since the latter upper bound is independent of β, this implies b = b, and thus the existence and the uniqueness of solution of (13) on [0, b).

Further estimates for concave F
To proceed we consider some additional assumptions on F. Namely, a vector-function A concave vector-function F is said to be strongly concave if for any α > 1 and any u ≥ 0 with u k > 0 there holds (13).

Corollary 2 Let F be a concave vector-function satisfying the Kamke-Müller condition. Let v(x) be a lower and u(x) be an upper solutions of
is a lower solution of (13).

Corollary 3 Let F be a concave vector-function satisfying the Kamke-Müller condition and
Proof By the assumptions . Then by Proposition 1, w(x) ≥ 0 for any x ∈ [0, b). Therefore by Corollary 2 w is a (nonnegative) lower solution to (13), thus by Proposition 2 we In the general case, let w(x) be the solution of (13) with the initial conditions Since which by virtue of (27) yields On the other hand, by our choice, for any k there holds that which yields (26).
Suppose F satisfy the Kamke-Müller condition and that it is concave. Then Let additionally F(w, x) be strongly concave, ξ > 0 and α > 1. If the patch k is (x, ξ). By the concavity condition, where Lu = du dx − F(u(x), x). In other words, w(x) is an upper solution with . This yields (28). Now, suppose that F(w, x) is strongly concave, ξ > 0, α > 1 and patch k is Faccessible at some β ∈ [0, b). By virtue of (28), it suffices to show that the equality On the other hand, by the assumption u(0) = ξ > 0 and Corollary 1 we have Using the strong concavity condition (25), (31) completes the proof.

The main representation
We start with an auxiliary model (36) below and then prove the existence of a unique positive solution of (1)-(4) and examine asymptotic behavior of the obtained solution.
Everywhere in this section we assume the conditions (H1)-(H4) are satisfied.

The balanced equations
Now we consider the particular case of (13) with F(w, x) given by (16). In other words, we consider the differential operator For further applications, it is useful to specify the properties of M k . Recall that in the Lotka-McKendrick-von Foester model (1) with (8) each M k (v, x) is actually an increasing linear function in v. While keeping monotonicity, we also impose some additional growth conditions on M k . Namely, we suppose that each M k (v, x) satisfies (H2), i.e. it is a nonnegative continuous function on R × [0, b), and there exist γ > 0 and μ ∞ > 0 such that where Proposition 4 Let L be given by (32) satisfying (H2) and (H3). Then for any The solution is nonnegative and bounded, and furthermore where μ(x) = min k μ k (x).
Proof Using the notation of (16), the Metzler property on D implies that F satisfies the Kamke-Müller condition in [0, b). Furthermore, since M k ≥ 0 one also has Thus, the assumptions of Proposition 2 are fulfilled. This yields the existence of the initial problem (36) and (37) which readily yields (38).

The main represenation
Lemma 4 Let B be defined by (2) and let

Then each of B − and B + is a connected open set.
Proof It suffices to prove that for any y ≥ 0, the set {s ≥ 0 : (s, y + s) ∈B} is connected. To this end let us suppose that (0, y) ∈B and let S be the closed component of {s ≥ 0 : (s, y + s) ∈B} containing (0, y). Let (s 1 , y + s 1 ) be the right endpoint of S. Then (s 1 , y + s 1 ) ∈ ∂B. We claim that (s, y + s) ∈ R 2 \B for s > s 1 . Indeed, arguing by contradiction, one concludes that there exists s 2 > s 1 such that (s 2 , y + s 2 ) ∈ ∂B. This yields B(y + s i ) = s i , i = 1, 2, thus by (5) 1 = B(y + ks 2 ) − B(y + ks 1 ) s 2 − s 1 < 1.
The contradiction yields our claim and, thus, the desired connectedness. Fig. 1 The domains B + and B + as it is shown on Figs. 1 and 2, respectively. Next, let Φ(x; ρ, y) denote respectively Ψ (x; f, y) the solutions h(x) of the initial value problems Lemma 5 Let ρ ∈ C(R + , R N + ) ∩ L ∞ (R + , R N + ) and let f ∈ C(R + , R N + ) satisfy (H5). Then Φ(x; ρ, y) (resp. Ψ (x; f, y)) is a nonnegative function non-decreasing in ρ (resp. f ). Furthermore, where ω 1 is defined by (39), and Proof It follows from Proposition 4 that (40) and (41) have a unique nonnegative solution. Next, given two arbitrary ρ and ρ * , let h(x) and h * (x) be the corresponding solutions of (40). If ρ ≥ ρ * then Proposition 1 imply h(x) ≥ h * (x) for x ≥ 0 and the monotonicity Φ(x; ρ, y) ≥ Φ(x; ρ * , y) follows. Similarly one shows the monotonicity of Ψ . Furthermore, if ρ(t) and ρ * (t) are two arbitrary nonnegative vector-functions, then Corollary 3 and Proposition 4 yield Proposition 5 and Proof First let (a, t) ∈ B and t > a. Then in the new variables (a, t) = (x, x + y) one has (x, y) ∈ B + 1 and the initial value problem (1)-(4) becomes (40) for h(x) = n(x, x + y). This yields n(x, x + y) = Φ(x; ρ, y) for each y > 0, thus, returning to the old variables yields n(a, t) = Φ(a; ρ, t − a) for any t > a > 0. This proves the first part of representation (46). The second part is similarly obtained by the change of variables (a, t) = (x + y, x). Furthermore, the continuity of n(a, t) follows from (46) and the standard facts on continuity of solutions on parameters. Finally, plugging (46) in (3) yields (47).

The integral equation
It is straightforward to see that if M, D, m and f are sufficiently smooth functions, then the function n(a, t) in (46) is a classical solution of the boundary value problem (1)-(4) in B. On the other hand, in application it is natural to assume that these functions are merely continuous (or even measurable). In that case, one can interpret the representation (46) with ρ satisfying (47) as a weak solution of (1)-(4). Furthermore, since a solution ρ(t) of the integral Eq. (47) completely determines the population dynamics n(a, t), it is natural to characterize all nonnegative solutions of (47) (with a given function f ). To this end, we observe that (47) can be thought of as an (nonlinear) operator equation on ρ: where the operators K and F are defined resp. by In this section we treat some general properties of L f . We fix some notation which will be used throughout the remaining part of the paper. Let ω 1 be defined by (39) and let where A m , a m are the constants from (H4) and m = m(a, t) = (m 1 (a, t), . . . , m N (a, t)) is the birth rate. Let us also consider the following subsets of R N + : Lemma 6 Let (H4) be satisfied. Then the operators F and K are positive on the cone of nonnegative continuous vector-functions C(R + , R N + ) and have bounded ranges: Furthermore, K is non-decreasing and Lipschitz continuous on C(R + , R N + ).
Proof It readily follows from the nonnegativity of m and Lemma 5 that K and F preserve the cone of nonnegative functions C(R + , R N + ) and non-decreasing there.
Furthermore, using (H4) we have from (43) This yields (53) and thus the boundedness of the range of K. The corresponding property for F is established similarly. Next, by (H4) m(a, t) (54). Finally, if ρ and ρ * are bounded functions then by (44), which yields that K is a Lipschitz continuous operator.
In order to establish the uniqueness we assume that ρ andρ are two solutions to (48). The tautological iterations ρ (i) := ρ andρ (i) :=ρ, i = 0, 1, 2, . . . obviously satisfy (55) which by virtue of (57) yields thus ρ(t) ≡ρ(t). Finally, by Lemma 5 Φ k and Ψ k are continuous, which yields the continuity of operators K and F, and, thus, all iterations given by (55) are continuous and so is the limit ρ. This completes the proof.

Lemma 7
Let ρ ∈ C(R + , R N + ) and ρ(t) > 0 for t ∈ [s 1 , s 2 ] ⊂ R + . Let for some k there exists β k < sup supp m k such that the patch k is accessible at β k . Then there exist a k , b k such that [a k , b k ] supp m k , β k ≤ a k , and (Kρ(t)) k > 0 for all t ∈ [s 1 + a k , s 2 + b k ].
Proof There are δ > 0 and points a k , b k , a k < a k < b k < b k , such that (i) m k (a) ≥ δ for a ∈ [a k , b k ], and (ii) the patch k is accessible at β k ≤ a k . By Lemma 3 we have hence by the maximum principle ξ(t) > 0 for any t ∈ (s 1 + a k , s 2 + b k ). This yields the desired conclusion.

Constant environment
The model (1)-(4) is more complicated for analysis under the assumption that a population lives in a temporally variable environment because the structure parameters are functions of age and time. In this section we analyze a constant environment, then in Sect. 6 we continue with a periodically changing environment, and finally in Sect. 7 we describe an irregularly changing environment. Throughout this section, we assume the conditions (H1)-(H5) are fulfilled and that the maximal life-time is constant: B(t) ≡ b. This condition is natural and is commonly used for both finite and infinite values of b, see Chipot (1983), Cushing (1984) and Gurtin and MacCamy (1974).

The characteristic equation
Under assumptions that the vital rates, carrying capacity and dispersion coefficients are time-independent functions, the system (1)-(4) becomes

), a)n(a, t) + D(a)n(a, t),
According to Proposition 6, there exists a unique solution n(a, t) of the problem (58) given by where the newborns function

m(a)n(a, t) da,
satisfies the following identity: Using the notation of (49) and (50), we have

Proposition 8 Let n(a, t) be the solution of the problem (58). Then the newborns function ρ(t) satisfies the integral equation
It is natural to study stationary (i.e. time independent) solutions of (61). Indeed, since m(a) has a compact support, it follows from (60) that Ff vanishes for large enough t. This yields that any solution of (61) satisfies In particular, it is easy to see that if ρ has a limit ρ ∞ = lim t→∞ ρ(t) then ρ ∞ itself is a stationary solution of (62). In the next section we study the stationary solutions in more detail.
Corollary 4 The operatorK is nondecreasing and where Q − is defined by (52).
Proof The nondecreasing property is by Proposition 1 and (66) follows from (53).

Definition 2 The equationK
is said to be the characteristic equation for the problem (61). A nonnegative solution ρ of (67) is called a stationary solution of (61).
The set of stationary solutions is nonempty because ρ = 0 is a (trivial) stationary solution. In Sect. 5.2 we characterize all nontrivial stationary solutions.
As it was noticed before, the characteristic equation describes the possible scenario of the limit at infinity of solutions to (61). The next lemma makes this observation more precise. First let us note that the limit is well-defined for any ρ ∈ S M , where Lemma 8 For any f ∈ C(R + , R N + ), Proof It follows from (54)  Next, by virtue our choice of t we have for any 0 which yields the desired conclusions.

Lemma 9
The set of lower solutions of (67) is bounded:

Proposition 9
For any ρ + ∈ Q + the limit exists and θ is a solution of the characteristic equation. Furthermore, (i) θ does not depend on a particular choice of ρ + ∈ Q + ; (ii) if ρ is an arbitrary lower solution of (67) then ρ ≤ θ .
Iterating the latter inequality yields ρ ≤K i ρ ≤K i ρ + , and passing to the limit as i → ∞ we get ρ ≤ ρ + (θ ). This proves (ii). Now suppose that ρ + 1 ∈ Q + . Then θ(ρ + 1 ) is a solution of the characteristic equation, hence by (ii) which, by symmetry, yields the equality in the latter inequality. This establishes the independence of θ(ρ + ) on a choice of ρ + , implying (i). (69) is called the maximal solution of the characteristic equation.

Definition 3 The unique θ defined by
Note that the maximal solution θ does not depend on the initial population distribution f(a) and it is essentially determined by the maternity function m(a). As we shall see, the maximal solution plays a distinguished role in the asymptotic analysis.

The basic reproduction number dichotomy
Throughout this section we assume additionally that the condition (H6) is also fulfilled. Let us consider the scaled version ofK by Equivalently, we have component-wise where Y(a; x, λ) = 1 λ ϕ(a; λx), x ∈ R N + .
Thus, the existence of a nontrivial solution to the characteristic Eq. (67) is equivalent to the existence of a pair (e, λ) , where a unit vector (a direction) e ∈ R N + , e = 1 and a scalar λ > 0 are such that e = R λ e.
The next lemma establishes that for each direction e ∈ R N + there is at most one such pair. We denote by C , C up and C low the classes of solutions, upper and lower solutions of (67), respectively.

Lemma 10
The operator R λ is decreasing with respect to λ: In particular, given an arbitrary direction e ∈ R N + , e = 1, where card is the cardinality of the corresponding set.
In the course of the proof of the lemma we have established the following property.
The limit case λ = 0 plays a distinguished role in the further analysis. Notice that Y k (a; x, λ) is non-decreasing in λ > 0 and by (42) Y k (a; x, λ) ≤ e N D b , where the constant b is from (H1). This implies that the limit with μ k (x) is defined by (35). Since m k ≥ 0, the limit is well defined for each x ∈ R N + . To proceed, we recall some standard concepts of the nonnegative matrix theory. A matrix A is called reducible (Minc 1988) if for some permutation matrix P where A 11 , A 22 are square matrices, otherwise A is called irreducible. There is the following combinatorial characterization of the irreducibility, see Berman and Plemmons (1979, p. 27), Meyer (2000, p. 671): the condition that a nonnegative matrix A of order n ≥ 2 is irreducible is equivalent to any of the following conditions: (a) no nonnegative eigenvector of A has a zero coordinate; (b) A has exactly one (up to scalar multiplication) nonnegative eigenvector, and this vector is positive; (c) αx ≥ Ax and x > 0 implies x 0; (d) the associated graph Γ (A) is strongly connected.

Lemma 11
The map R 0 : R N → R N defined by (75) is linear and strongly positive, i.e. x > 0 implies R 0 x 0. In particular, R 0 is an irreducible matrix. Furthermore, Proof Indeed, the linearity follows immediately by (75) and (74). Since the matrix M(0, a) is diagonal, the associated digraphs of the matrices D(a) and D(a) − M(0, a) are equal. Therefore, using (H6) readily yields that Y k (a; x, 0) > 0 for any a > β k . Hence, repeating the argument of Lemma 10 we have from (75) and (H4) that (R 0 x) k > 0 for any k. This proves R 0 x 0. Suppose by contradiction that R 0 is reducible. Then for some permutation matrix P where A 11 , A 22 are square matrices. Let x > 0 be a vector in R N + with all first m coordinates zero, where m is the order of A 11 . By (77) PR 0 P t x has the same property, i.e. the vector R 0 P t x has at least m zero coordinates which contradicts to the fact that R 0 P t x 0. This proves the irreducibility. Finally, (76) follows from (73).

Corollary 6
If R 0 e ≤ e for any e ∈ R N + , e = 1, then the characteristic Eq. (67) admits only trivial solutions.
Proof Indeed, if ρ = 0 is a nontrivial solution of (67) then by (70) e = ρ/ ρ is a solution of R λ e = e for λ = ρ . On the other hand, using the assumption and (76) we obtain a contradiction follows.
Let us denote by R 0 the spectral radius of the linear map R 0 . Combining the irreducibility of R 0 with the Perron-Frobenius theorem (Berman and Plemmons 1979, Theorem 1.3.26) implies the following important observation.

Corollary 7 The spectral radius R 0 > 0 and it is a simple eigenvalue of
Definition 4 The linear map R 0 is called the net reproductive map associated to the problem (58). Its spectral radius R 0 is called the basic reproduction number.
The latter definition can be motivated as folows. For a single patch model, i.e. N = 1, the linear system (74) becomes a single equation where The quantity R 0 is well-established and is known as the (inherent) basic reproduction number in the linear time-independent model on a single patch (Iannelli 1995;Cushing 1998); see also Kozlov et al. (2016b) or Kozlov et al. (2017). Note that in this case, is the survival probability, i.e. the probability for an individual to survive to age v. Then R 0 is the expected number of offsprings per individual per lifetime. Recall that in the one-dimensional case, R 0 is related to the intrinsic growth rate of population by the characteristic equation. Namely, when R 0 > 1 population is growing, while for R 0 ≤ 1 population is declining. The next result extends this dichotomy onto the general multipatch case. Proof First let us assume that R 0 ≤ 1 and suppose by contradiction thatKρ = ρ for some ρ > 0. Let λ = ρ and e = ρ/λ, then by (70) and (76), The latter easily implies that there exists t > 1 such that R 0 e ≥ te. On iterating the obtained inequality yields R k 0 e ≥ t k e, thus a contradiction. Now suppose that R 0 > 1. By Corollary 7, there exists a positive eigenvector e 0 0 of R 0 . Since e 0 0 there exists λ > 0 such that λe 0 ≥ ω 2 1 N , where ω 2 is defined by (51). By (52), ρ + := λe 0 ∈ Q + , hence Lemma 9 implies that is a solution to (67). On the other hand, since R 0 > 1 we have hence, by the continuity argument for some λ > 0 small enough there holds Therefore, setting ρ − := λe 0 we obtain i.e. ρ − is an lower solution of (67). In other words, ρ − ∈ C low , thus (ii) of Proposition 9 yields thus θ is a nontrivial solution.

Asymptotic behaviour of a general solution of (61)
Let us return to the general Eq. (61). If the initial distribution of population vanishes: n(a, 0) = f(a) = 0, the uniqueness of solution of (58) immediately implies that the population density n(a, t) ≡ 0 for all a, t ≥ 0. This conclusion also holds true even under a weaker assumption that Ff ≡ 0. The latter is evident from the biological point of view: the population disappears if its initial distribution is older that the maternity period. Taking into account these observations, it is naturally to assume that The main result of this section states that under this assumption, any solution of (61) behaves asymptotically as the maximal solution.
Theorem 2 Let χ be the solution to (61) satisfying (81). Then We start with two results describing the upper and lower solutions to Eq. (61).

Lemma 12 Let χ be a solution to (61). Then
where the latter inequality should be understood component-wise.
Notice that that the class of stationary upper solutions is nonempty. Indeed, it follows from (53) that, for example, 2(ω 2 + )1 is such a an upper solution for any > 0. Now, let us define the iterative sequence by Then applying the argument of the proof of Proposition 9 yields that {ρ (i) } is nonincreasing: Also, since ρ (i) is a constant vector function, it follows by Lemma 8 that L f ρ (i) ∈ S A m and also that We claim that for any j ≥ 0 The proof is by induction. Notice that (b) holds trivially for j = 0, and by the assumption (84) which yields (a) for j = 0. Let the claims (a)-(b) hold true for some j ≥ 1. Then (a) follows from the monotonicity of L f : Furthermore by the assumption χ ( j) ∈ S j A m and χ ( j) ∞ = ρ ( j) . Hence Lemma 8 yields which yields (b) for j + 1. Next, it follows from (a) and the boundedness of the image of L that {χ ( j) } is nonincreasing and bounded from below, thus has a limit which obviously is a solution of (61). By the uniqueness, lim j→∞ χ ( j) (t) = χ(t). Now, let 1 ≤ k ≤ N . Then the sequence of the coordinate functions χ ( j) k (t) is non-increasing with respect to j and lim j→∞ χ ( j) k (t) ≤ θ k + for all j ≥ j 0 and t ≥ j A m . Passing to the limit j → ∞ we obtain χ k (t) ≤ θ k + for t ≥ j A m which easily implies (83).

Proof of Theorem 2
If R 0 ≤ 1, then Theorem 1 yields θ ≡ 0, then (83) immediately yields (82). Therefore we shall suppose that R 0 > 1. Let χ be the unique solution to (61) and let θ 0 be the unique maximal solution of (67). By Lemma 13, it suffices to show that there exists a lower solution ρ − to (60) such that ρ − ∈ S M for some M ≥ 0 and ρ − ∞ = 0. In the remained part of the proof we shall construct such a solution. Let us consider an auxiliary sequence of iterations We claim that the new function ρ − (t) defined by is a lower solution to Eq. (61) for certain M > A m , sufficiently large j ≥ 1 and sufficiently small λ > 0 to be specified later. To this end, first notice that hence using an induction by j ≥ 1, one gets i.e. the sequence ρ ( j) is non-decreasing in j. It also follows from the latter inequality that ρ ( j) ≤ L f ρ ( j) , i.e. ρ ( j) is a lower solution to (61). Hence, ρ − (t) defined by (86) is a lower solution to (60) in the interval t ∈ [0, M]. In particular, Next, we assume that t ∈ [M, M + A m ]. By the assumption M > A m , hence one has (Ff)(t) = 0 and L f ρ − = Kρ − . We have by (86) and condition (H4) that  (65) that We claim that the integrals (87) and (88) are nonnegative. The first integral is nonnegative by virtue of Corollary 5. To show that (88) is nonnegative, let us estimate where a k ≥ β k and b k are the same as in Lemma 7. Since F(f) is not identically zero, there exists an interval [s 1 , s − 2], where this function is positive. Applying Lemma 7 for ρ = ρ (1) = F(f), we get that Therefore and, in particular this is true for k = 1. Repeating this argument yields This implies that Φ k (a; ρ ( j) , t − a) > 0 for a ≥ β k and t − a ∈ [s 1 + ( j − 1)a 1 , s 2 + ( j − 1)b 1 ].

Now we choose the index j and the number M to satisfy
Therefore, for such a and t if λ is sufficiently small positive number. This gives positivity of (88) for t ≥ M + a k . If t ≤ M + a k then the first integral in (88) This proves that the function ρ − (t) defined by (86) is a lower solution to Eq. (67), therefore by Lemma 13 we have the desired convergence that completes the proof.

Asymptotics of total population
According to the assumption made in the beginning of this section, the maximal length of life is constant: Then the total (multipatch) population P(t) at time t is the vector-function Then we have the following result.
Theorem 3 Let n(a, t) be the solution of (58) and let the condition (81) hold. Then the following dichotomy holds: if R 0 ≤ 1 then P(t) → 0 as t → ∞, and if R 0 > 1 then where θ is the maximal solution to the characteristic equation.
Proof Denote by ρ(a) the newborns function determined by f(a) by virtue of (61). We have for general t > 0 On the other hand, by (H5) supp f ⊂ [0, b], hence using (59) we have for any t > b that Next, by Theorem 1 and Theorem 2 we have lim t→∞ ρ(t) = θ and furthermore by By continuity of solutions (92) with respect to a parameter and (64), we have for any fixed a > 0 that lim t→∞ Φ(a; ρ, t − a) = ϕ(a; θ).

Estimates for the basic reproduction number and for the maximal solution
In this section we shall assume that the condition (12) hold, i.e.
The biological meaning of the latter inequality is that individuals do not reproduce during migration (but can die). This condition immediately implies that D kk (a) ≤ 0.
Throughout this section, we use the following notation: Proposition 10 Under the made assumptions, Proof By Corollary 7 there exists an eigenvector ρ 0 of R 0 corresponding the maximal eigenvalue R 0 , i.e. R 0 ρ = R 0 ρ. Let us consider the problem (74) with the initial condition x = ρ. Using the assumption (12) and summing up the Eq. (74) for Then by (75) Since the sum N k=1 ρ k > 0 we arrive at the right hand side of (93). Now, in order to prove the left hand side inequality in (93), notice that in the made notation by virtue of D k j (a) ≥ 0 for k = j and Y j (a, ρ) ≥ 0 for all admissible a we have Combining this with (75) we obtain thus implying (93) by virtue of ρ k > 0.

Remark 1
The estimates (93) are optimal. Indeed, if D k j ≡ 0, the system (74) splits into separate equations implying that each e k is an eigenvector of R 0 with eigenvalue In order to establish the corresponding estimates for the maximal solution θ we consider an auxiliary function where the minimum is taken over the simplex Lemma 14 In the above notation, M(t, a) is nondecreasing in t > 0 and Furthermore, where p(a) is the function from (H2).
is the minimum point of (94) then ξ = λξ ∈ S(t ), hence using the monotonicity condition in (H2) and ξ ≥ ξ we obtain which yields the nondecreasing monotonicity. In particular the limit in (95) Passing to the limit as t → +0 in the latter inequality yields μ ≤ μ(a), thus implying (95). Finally, assume again that ξ ∈ S(t) is the minimum point of (94) for t > 0. Then using (H2) and the Hölder inequality we obtain which yields (96).
Remark 2 Let us comment on (98) from the biological point of view. Notice by Theorem 2 that N k=1 θ k is the asymptotical value of the total number of newborns on all patches. By the dichotomy, R 0 ≤ 1 implies θ = 0, thus the total asymptotical number of newborns is zero. On the other hand, in the nontrivial case R 0 > 1, hence by (93) ∞ 0 m(a)e − a 0 μ(v)dv da > 1, which easily implies that (97) has a positive solution.
The next proposition provides a lower estimate for the maximal solution.
Proposition 12 Let there exist a function q(a) > 0 such that If for some k where θ − k is the unique solution to equation and Proof First notice that (101) implies by (93) that R 0 > 1, thus θ 0. Since D k j (a) ≥ 0 for j = k, the k-th equation in (64) yields ϕ k (a, θ), a) − D kk (a))ϕ k (a; θ), hence using (100) we obtain by virtue of M k (0, a) = μ k (a) that Arguing similar to the proof of Proposition 11 we get from ϕ k (0; θ) = θ k that Since θ 0, one has θ k > 0, hence Then I (t) is decreasing, I (θ k ) ≤ 1 and by (101) I (0) > 1, thus there exists (a unique) solution θ − k of (103) such that θ k ≥ θ − k .

Periodically varying environment
We consider an important particular case of the main problem (1)-(4) when the environment is periodically changing. In this section and in the rest of the paper, it is assumed that the vital rates, regulating function and dispersion coefficients are timedependent and periodic with a period T > 0. The boundary-initial value problem (1)- (4) for any 1 ≤ k, j ≤ N . Throughout this section, we assume that the conditions (H1)-(H5) are satisfied.
Notice that the existence and uniqueness of a solution n(a, t) to the periodic problem follows from the general result given by Proposition 6 and it is given explicitly by (46). Note also that n(a, t) need not to be periodic in t but it is natural to expect that n(a, t) converges to a T -periodic function ρ(t) for t sufficient large, where ρ(t) solves the associated characteristic equation Here the operator K is defined by where the initial condition We shall assume that the nonnegative cone C T (R + , R N + ) is equipped with the supremum norm ρ(t) C([0,T ]) . It follows from the uniqueness results of Sect. 4.2 that the function Φ(x; ρ, y) is T -periodic in y.
A function ρ ∈ C T (R + , R N + ) is said to be an upper (resp. lower) solution to (106) if ρ ≥ Kρ (resp. ρ ≤ Kρ). It follows from Lemma 5 and condition (H4), it follows that K has a bounded range: In particular, any solution of the characteristic Eq. (106) is bounded by ω 2 .
Recall that a (nonlinear) operator is called absolutely continuous if it is continuous and maps bounded sets into relatively compact sets.
is an absolutely continuous operator.
Proof By the Arzela-Ascoli theorem it suffices to show that the family of functions is uniformly bounded and equicontinuous for any R > 0. The first property is by (108).
In order to prove that the family is equicontinuous, we estimate | Kρ(t 1 ) − Kρ(t 2 )| for |t 1 − t 2 | < δ 1 and for any ρ To this end, we assume that τ := t 2 − t 1 > 0 is such that where a m , A m and b 1 are the structure constants in (H1) and (H4) . Rewriting and using the property that m(a, t i ) = 0 for any a outside [a m , A m ] for i = 1, 2 we obtain component-wise estimates We have by (43) that for any τ > 0 and where C 1 depends only on the structural constants. Next, since m k (a, t) is a T -periodic in t, by (H4) m k is uniformly continuous on the strip [a m , A m ] × R. Since supp m k ⊂ [a m , A m ] × R, there exists δ 2 > 0 such that for any 1 ≤ k ≤ N , a ∈ [0, A m ] and |τ | < δ 2 one has the inequality This yields I 1 < /2. In order to estimate I 2 , we notice that Φ k (x) := Φ k (x, t 1 −a; ρ) is the solution of the initial problem (107). Notice that by (42) where the latter equality is by the periodicity. Therefore, applying the mean value theorem to (107) we obtain for any 0 ≤ x 1 < x 2 < A m + δ 1 and for some ξ ∈ (x 1 , x 2 ) that where C 4 depends only on the structure conditions and R. This readily implies Choosing δ 1 small enough, yields the desired conclusion.
Proposition 13 For any ρ + (t) such that ρ + (t) ≥ ω 2 · 1 N , where ω 2 is defined by (51), the limit exists and is a solution to the characteristic Eq. (106). Furthermore, the limit θ(t) does not depend on a particular choice of ρ + (t) and it is the maximal solution to Eq. (106) in the sense that if ρ(t) is any solution to the characteristic Eq. (106) then Proof Since Kρ + (t) ≤ ω 2 · 1 N ≤ ρ + (t) and by the monotonicity of K we get: which implies that {ρ ( j) (t)} is a non-increasing sequence. The sequence is bounded from below because K j ρ + ≥ 0, therefore there exists a pointwise lim j→∞ K j ρ + (t) =: θ(t). The sequence {ρ ( j) (t)} is uniformly bounded by the constant ω 2 . Applying Lemma 15 to family {ρ ( j) k (t)} implies that the convergence is in fact uniform on each compact subset of R. Thus θ is a nonnegative continuous T -periodic solution of (106). The rest of the proof is analogous to the proof of Proposition 9.
In the remaining part of this section we additionally assume that additionally condition (H6) holds. In that case, due to the periodicity, the infimum in (H6) can be replaced by the minimum. Then arguing similarly to Lemma 10, one can verify that for any ρ(t) ∈ C T (R + , R N + ) and 0 < λ 1 < λ 2 , hence the corresponding next generation operator is well-defined defined by where Y(x, y; ρ) is the solution of the linear system Let R 0 denote the largest eigenvalue of R 0 and let θ = θ(t) ∈ C T (R + , R N + ) be the maximal solution of Eq. (106). Then the following results are established similarly to Theorems 1, 2 and 3 respectively.
where θ is the maximal solution to the characteristic Eq. (106).

Irregularly varying environment
In order to study asymptotic behavior of the solution to the model (1)-(4) in the case when temporal variation is irregular, we assume that the vital rates, regulating function and dispersion coefficients are bounded from below and above by equiperiodic functions for large t. These periodic functions define two auxiliary periodic problems, whose solutions provide upper and lower bounds to a solution of the original problem. This leads us to two-side estimates of a solution to the original problem for large t.
More precisely, throughout this section we shall suppose that there exists T 1 ≥ 0 and T -periodic functions m ± k , M ± k and D ± k j such that for any a ≥ 0 and t ≥ T 1 m − (a, t) ≤ m(a, t) ≤ m + (a, t), As in Sect. 6, one can consider the corresponding characteristic equations where ν denote − or +, and the operators K ν are defined component-wise by and Φ ν (x, y; ρ) is the unique solution of the system with ρ ∈ C T (R + , R N + ). Then by Proposition 6 Also let us denote by R ± 0 and R ± 0 the corresponding next generation operators and basic reproduction numbers. The main result of this section states that a solution of the population problem in an irregularly changing environment can be estimated by the corresponding solutions of the associated periodically varying population problems.
Proof Without loss of generality T 1 ≥ B(0). Let R > 2ω 2 and let us define {χ ( j) (t)} j≥0 and {χ where the operators K and F are defined by (49) and (50) respectively. Arguing as in the proof of Proposition 7, we obtain the existence of the limit lim j→∞ χ ( j) (t) = χ(t), where χ(t) is a solution to (48). Also by Proposition 13, lim j→∞ χ ( j) is the maximal solution to (110). We will prove by induction that for any j ≥ 0 there holds For j = 0 the claim follows from χ (0) (t) = χ (0) + , y) = R for any y ≥ 0 and the structure parameters are estimated by (109), one easily deduces from the definition of Φ ν (x, y; ρ) that Φ(a; χ (0) , t − a) ≤ Φ + (a; χ (0) + , t − a) for a ≥ 0 and t − a ≥ T 1 . Since and t − a > T 1 for all a ∈ [a m , A m ] and t > T 1 + A m we obtain This proves the induction assumption for j = 1. Now suppose that the induction claim holds for some j ≥ 1. Arguing similarly, we obtain for any t > T 1 + ( j + 1)A m that which proves (113). Therefore, passing to the limit we obtain If R + 0 ≤ 1, then by Theorem 5 lim t→∞ χ + (t) = 0, hence (114) implies (i).
To proceed with (ii) notice that (114) already yields the upper estimate in (112). It remains to show that there exists a lower solution χ − (t) to χ(t) = L − f χ(t). We use auxiliary sequence {χ where ρ − is a solution to the characteristic equation K − ρ − (t) = ρ − (t) and λ > 0 is sufficiently small.

Applications
In this section we consider two simple applications of our approach showing how dispersion promotes survival of a population on sink patches. In the usual situation, a habitat is a mixture of sources and sinks. Our first example shows that permanence on all patches is possible if the patches are connected and if emigration from sources is sufficiently small and does not cause extinction of a local subpopulation. Some researchers indicate that survival of migrating species is possible even if all occupied patches are sinks, see Jansen and Yoshimura (1998). Taking migratory birds as an example, we demonstrate that this is possible under certain conditions.

A single source and multiple sinks
In order to demonstrate the influence of dispersion on persistence of population, we compare a system with N isolated patches with the corresponding system with dispersion. Recall that in the isolated case, D(a, t) ≡ 0 implying by (79) that the basic reproduction number of the kth patch is given by where Π k (v) = e − a 0 μ k (v) dv is the survival probability. In this case the spectrum of the next generation operator is We assume that R (1) 0 > 1 and R (k) 0 ≤ 1, for k ≥ 2. In the biological terms, this is equivalent to saying that the first patch is a source and all other patches are sinks. Without migration, the population will persist on the first patch and become extinct on all other patches. For details about the age-structured logistic model that we used to describe isolated patches, we refer readers to Kozlov et al. (2017). Under the made assumptions, is given by Therefore, the next generation operator takes the form Then latter relation yields Now, recall that if A is a symmetric matrix and x is an eigenvector with a simple eigenvalue λ then the corresponding perturbed eigenvalue of A + B (B may not be symmetric) is given by For ε = 0, the largest eigenvalue is R (1) 0 with the eigenvector e 1 = (1, 0, . . . , 0). The perturbed eigenvalue, which will be the basic reproduction number for the next generation operator R 0 , is thus R 0 > 1 for small enough ε > 0 provided that tyhe function B 11 (a) ≤ 0 everywhere and strictly negative in at least one point of the support of m 1 . Thus shows that survival on all patches is possible if emigration from the source is sufficiently small.

Multiple sinks, without a source
We consider an extreme situation when a population inhabits two patches and the basic reproduction number on each patch is less or equal to one. A realistic example for this kind of situation is a population of migratory birds. Their habitats consists of two patches: breeding range (characterized by the high birth rate in summer and high death rate in winter) and non-breeding range (low birth and death rates). Thus, the breeding range is a sink because of the winter conditions, and the non-breeding range is a sink because of too few births. We will demonstrate that, even in this case, there is a chance of survival if the structure parameters are suitably chosen.
In the matrix form this becomes Thus, to show that R 0 > 1, it is sufficient to show that Pρ > 0 for some choice of parameters and certain vector ρ.
It follows that Pρ > 0 and hence R 0 ρ > ρ. The latter implies that that R 0 > 1, thus R 0 has an eigenvalue greater than one, which proves the permanence of population on both patches.

Discussion
Our work advances understanding of dynamics of age-structured, density-dependent populations in patchy, temporally-variable habitats. Intra-specific competition can act in different ways and have a whole spectrum of different forms, where the pure intraspecific competition and unstructured population-wide competition are on the opposite ends. Population growth models are often based on the assumption that competition is unstructured. In contrast to this, we focus on pure intra-cohort competition as we find it more relevant for exploring ecological questions related to, for example, ontogenetic shifts.
Our results show that the use of the basic reproduction number R 0 for determining permanence criteria for age structured populations can be expanded from a single-patch models into multi-patch models. The main result, Theorem A (the Net Reproductive Rate Dichotomy), is applicable if such populations are regulated by density dependent mortality and live in temporally unchanging habitats or if they face periodic environmental variation as well as fluctuations of large amplitudes. This opens up the opportunity to analyse a range of possible outcomes given specific settings relevant for ecological understanding or given specific parameter sets that define the ecological question.
We use the results to further explore the dynamics of a population that inhabits both source and sinks patches. We also study if a metapopulation can persist if all patches are sinks. Our results for the first question show that permanence is possible if all patches are connected and emigration from source patches is sufficiently small and does not cause extinction of a local subpopulation. Analysis of the second question shows that permanence of a population on pure sink patches is possible for certain patterns of dispersion.
The methodology used in the proofs, for example the storng monotonicity of lower and upper bounds, is possible to expand into other areas of population dynamics. We believe that it can be applied in analysis of multispecies dynamics, such as predatorprey, or epidemiological dynamics under assumption that model parameters are timedependent.