On the Multispecies Delayed Gurtin-MacCamy Model

Abstract

The paper deals with the description of multispecies model with delayed dependence on the size of population. It is based on the Gurtin and MacCamy model. The existence and uniqueness of the solution for the new problem of 𝑛 populations dynamics are proved, as well as the asymptotical stability of the equilibrium age distribution.

1. Introduction

In this paper we consider mutual influence of 𝑛 populations. We assume that each population develops differently, affecting each other. Populations do not destroy each other. However, they share the same natural resources and the space. In our model 𝑢 denotes the density of the ecosystem consisting of 𝑛 different populations. Therefore 𝑢=(𝑢¹,…,𝑢^𝑛) is the vector function. Each component 𝑢^𝑖 for 𝑖=1,…,𝑛 denotes the 𝑖th population density. The development of 𝑖th population can be expressed by the following system of the equations: 𝐷𝑢^𝑖(𝑎,𝑡)=−𝜆^𝑖𝑎,𝑧_𝑡𝑢^𝑖(),𝑢𝑎,𝑡_𝑖(0,𝑡)=^∞_0[−𝑟,0]𝛽^𝑖(𝑎,𝑧(𝑡+𝑠))𝑢^𝑖𝑎,𝑡+𝑠^𝑖𝑧𝑑𝑠𝑑𝑎,_𝑖(𝑡)=^∞₀𝑢^𝑖(𝑎,𝑡)𝑑𝑎 (1)(2)(3) with the initial condition 𝑢^𝑖(𝑎,𝑠)=𝜑^𝑖(𝑎,𝑠)for𝑠∈−𝑟^𝑖,0, (4) where [−𝑟,0]=−𝑟¹,0×−𝑟²,0×⋅⋅⋅×−𝑟^𝑛,𝑧()=𝑧,0𝑡+𝑠¹𝑡+𝑠¹,𝑧²𝑡+𝑠²,…,𝑧^𝑛𝑡+𝑠^𝑛,𝑧^𝑖_𝑡∶−𝑟^𝑖,0→ℝ₊,𝑟^𝑖⩾0,ℝ₊=[),𝑧0,∞_𝑖𝑡(𝑠)=𝑧^𝑖(𝑡+𝑠). (5) The rate 𝐷𝑢^𝑖(𝑎,𝑡)=lim_ℎ→0𝑢^𝑖(𝑎+ℎ,𝑡+ℎ)−𝑢^𝑖(𝑎,𝑡)ℎ (6) denotes the intensity of changing of the 𝑖th population in time. In particular, if 𝑢^𝑖 is differentiable then 𝐷𝑢^𝑖=𝜕𝑢^𝑖/𝜕𝑎+𝜕𝑢^𝑖/𝜕𝑡. The quantity 𝑧^𝑖(𝑡) is the total 𝑖th population at time 𝑡. We can express total population of the whole ecosystem at time 𝑡 by 𝑧(𝑡)=𝑧¹(𝑡),…,𝑧^𝑛(𝑡). (7) Thus 𝑧_𝑡=𝑧¹_𝑡,…,𝑧^𝑛_𝑡. (8) We assume that 𝑧^𝑖₀(0)=^∞₀𝜑^𝑖(𝑎,0)𝑑𝑎. (9) The birth and death processes of the 𝑖th population are described by coefficients 𝜆^𝑖 and 𝛽^𝑖. We assume that both processes depend on the population size not only at the moment 𝑡, but also at any preceding period of time (so-called history segment). Introducing the delay parameter to the model has deep biological approach (see, for instance, [1]). All natural processes occur with some delay with respect to the moment of their initiation. We can take into consideration, for example, the period of pregnancy (constant delay) or morbidity (variable delay). In our model, considered processes occur with 𝑛 various delays, typical for 𝑛 populations. It is a novel approach to the population dynamics research. Moreover the system (1)–(4) consists of the general description of 𝑛 populations dynamics. If populations develop independently, then the functions 𝜆^𝑖 and 𝛽^𝑖 depend on 𝑖th coordinate of the function 𝑧_𝑡 only and the system consists of 𝑛 independent von Foerster equations. However, if there are populations relationships, we can consider three cases. (1)  Competition for food: in this case 𝛽^𝑖 is decreasing function of variables describing the numbers of species competing for food with the 𝑖th one. For example, 𝛽^𝑖(𝑎,𝑧)=max{𝑎(𝐹−∑_𝑗∈𝐼𝑘_𝑗𝑧_𝑗),0}, where 𝐹 denotes total amount of food and 𝐼 is the set of species competing with the 𝑖th one. (2)  Schema predator prey: if 𝑖th species feeds on 𝑗th species individuals, then 𝜆^𝑗 and 𝛽^𝑖 are increasing functions of variables 𝑧^𝑗 and 𝑧^𝑖, respectively. The classical models show that balance state between two species competing for food is not possible. In this case the model does not have nonzero equilibrium age distribution. However, if we additionally consider predator feeding on one competing species, then we get balance between three considered species (see, for instance, [2–5]). (3)  Symbiosis: if 𝑖th and 𝑗th species live in symbiosis with each other, then 𝛽^𝑖 can be an increasing function of 𝑧^𝑗 (when the presence of 𝑗th species individuals is the reproduction of 𝑖th species favour) or 𝜆^𝑗 is a decreasing function of 𝑧^𝑖 (when the presence of 𝑖th species protects 𝑗th species individuals from death). In this paper we develop the idea of Gurtin and MacCamy model [6] for one age-dependent population ()((𝑡𝑢(),𝑢()=𝐷𝑢𝑎,𝑡=−𝜆𝑎,𝑧))𝑎,𝑡0,𝑡^∞₀𝛽((𝑡𝑢()𝑧(𝑡)=𝑎,𝑧))𝑎,𝑡𝑑𝑎,_∞0𝑢()𝑢()(𝑎).𝑎,𝑡𝑑𝑎,𝑎,0=𝜑 (10) This classical model has many generalizations [7–16]. Considerations of the dynamics of age-structured populations have received substantial treatment on various fields [17–22]. In particular our latest extension of this theme was presented in the article [10], where we considered an age-dependent population dynamics with a delayed dependence on the structure. The right-hand side of the equation in our model is not in the form 𝜆𝑢 but Λ𝑢. Here Λ is an operator that does not apply to the value of the function 𝑢, as in the classical model, but to its restriction to some particular space: 𝐷𝑢(𝑎,𝑡)=−Λ𝑎,𝑧_𝑡𝑢(),𝑢()=𝑎,𝑡+⋅0,𝑡^∞₀⁰_−𝑟𝛽((𝑢()𝑧(𝑡)=𝑎,𝑧𝑡+𝑠))𝑎,𝑡+𝑠𝑑𝑠𝑑𝑎,_∞0𝑢()[],𝑢()()[].𝑎,𝑡𝑑𝑎for𝑡∈−𝑟,𝑇𝑎,𝑠=𝜑𝑎,𝑠for𝑠∈−𝑟,0 (11) In our paper [10] we present the proof of the existence and uniqueness of the new problem solution. Conditions of the exponential asymptotical stability for this model still remain to be formulated. On account of the significant difficulties with formulating conditions of the stability we consider slightly simpler version of the model than the one described in [10]. We resign from the delay in the structure during a deliberation on mortality. However, the delayed dependence of the structure still remains the element of the birth process description. In this paper we aim to state stability conditions for the ecosystem of 𝑛 populations. The plan of the paper is as follows. In Section 2 we formulate the problem in the terms of operator equations. Section 3 contains the proof of the existence and uniqueness of the solution. We also study equilibrium age distributions, that is, solutions to the problem which are independent of time. In Section 4 exponential asymptotic stability of the equilibrium age distributions is analyzed. We use generalized Laplace transform to study the stability.

We consider our model (1)–(4) under the following assumptions. (𝐻₁)𝜑=(𝜑¹,…,𝜑^𝑛) where 𝜑^𝑖∈𝐶(ℝ₊×[−𝑟^𝑖,0]) for 𝑖=1,…,𝑛. (𝐻₂)Φ^𝑖=∫^∞₀sup_{𝑠∈[−𝑟𝑖,0]}𝜑^𝑖(𝑎,𝑠)𝑑𝑎<∞ for 𝑖=1,…,𝑛. (𝐻₃)The function [−𝑟^𝑖,0]∋𝑡↦∫^∞₀𝜑^𝑖(𝑎,𝑡)𝑑𝑎 is continuous. (𝐻₄)𝜆=(𝜆¹,…,𝜆^𝑛), 𝛽=(𝛽¹,…,𝛽^𝑛) where 𝜆^𝑖∈𝐶(ℝ₊×⨉^𝑛_𝑖=1𝐶([−𝑟^𝑖,0])), 𝛽^𝑖∈𝐶(ℝ₊^𝑛+1); the Fréchet derivatives 𝐷𝜆^𝑖 of 𝜆^𝑖(𝑎,𝜓) with respect to 𝜓 exist for all 𝑎⩾0 and 𝜓⩾0. (𝐻₅)The components of the function 𝜆(⋅,𝜓) =(𝜆¹(⋅,𝜓),…,𝜆^𝑛(⋅,𝜓)) belong to 𝐶(𝐶([−𝑟^𝑖,0]);  𝐿^∞(ℝ₊)), respectively, for 𝑖=1,…,𝑛. (𝐻₆)The components of Fréchet derivative 𝐷_𝜓0𝜆=(𝐷_𝜓0𝜆¹,…,𝐷_𝜓0𝜆^𝑛) in the point 𝜓₀ as functions of 𝜓₀ belong to 𝐶(𝐶([−𝑟^𝑖,0]);ℒ(𝐶([−𝑟^𝑖,0]),𝐿^∞(ℝ₊))), 𝑖=1,…,𝑛. Here ℒ(𝑋,𝑌) denotes the Banach space of all bounded linear operators from 𝑋 to 𝑌. (𝐻₇)𝜑^𝑖⩾0,𝜆^𝑖⩾0,𝛽^𝑖⩾0 for 𝑖=1,…,𝑛. (𝐻₈)The functions 𝜆 and 𝛽 are bounded; that is, ‖𝜆(𝑎,𝜓)‖⩽𝜆₀ and ‖𝛽(𝑎,𝑧)‖⩽𝛽₀ where 𝜆₀ and 𝛽₀ are finite quantities. (𝐻₉)There exist 𝜆₁ and 𝛽₁ such that 𝜆𝑎,𝜓₁−𝜆𝑎,𝜓₂⩽𝜆₁𝜓₁−𝜓₂ (12)      for every 𝑎,𝜓₁,𝜓₂ and ^𝑛𝑛_𝑖=1𝑗=1𝜕𝛽^𝑖(𝑎,𝑧)𝜕𝑧^𝑗⩽𝛽₁ (13)       for every 𝑎,𝑧.

2. Equivalent Formulation of the Problem

In this section we will formulate problem (1)–(4) in terms of operator equations. Thanks to this we will prove local and global existences of the presented problem solution. The next theorem is analogous to these well-known results, for example, von Foerster [23] or Gurtin and MacCamy models.

Theorem 1. Nonnegative continuous functions 𝑧_𝑡=(𝑧¹_𝑡,…,𝑧^𝑛_𝑡) and 𝐵=(𝐵¹,…,𝐵^𝑛), where 𝐵^𝑖(𝑡)=_[−𝑟,0]^{𝑡+𝑠𝑖}₀𝛽^𝑖𝑡+𝑠^𝑖−𝑎,𝑧(𝑡+𝑠)×𝐵^𝑖(𝑎)𝑒^{−∫𝑖𝑎𝑡+𝑠𝜆𝑖(𝜏−𝑎,𝑧𝜏)𝑑𝜏}+𝑑𝑎𝑑𝑠_{[−𝑟,0]∞0}𝛽^𝑖𝑎+𝑡+𝑠^𝑖,𝑧(𝑡+𝑠)×𝜑^𝑖(𝑎)𝑒^{−∫𝑖0𝑡+𝑠𝜆𝑖(𝑎+𝜏,𝑧𝜏)𝑑𝜏}𝑧𝑑𝑎𝑑𝑠,_𝑖𝑡𝑠^𝑖=^{𝑡+𝑠𝑖}₀𝐵^𝑖(𝑎)𝑒^{−∫𝑖𝑎𝑡+𝑠𝜆𝑖(𝜏−𝑎,𝑧𝜏)𝑑𝜏}+𝑑𝑎_∞0𝜑^𝑖(𝑎,0)𝑒^{−∫𝑖0𝑡+𝑠𝜆𝑖(𝑎+𝜏,𝑧𝜏)𝑑𝜏}𝑑𝑎for𝑡+𝑠^𝑖𝑧⩾0,_𝑖𝑡𝑠^𝑖=𝑧^𝑖₀𝑡+𝑠^𝑖for𝑡+𝑠^𝑖<0, (14)(15) for 𝑖=1,…,𝑛, are the solutions of the problems (2) and (3) up to time 𝑇>0 if and only if the function 𝑢=(𝑢¹,…,𝑢^𝑛) is the solution of the age-dependent 𝑛 populations problem (1)–(4) on [0,𝑇] and 𝑢 is defined by the formula 𝑢^𝑖()=𝜑𝑎,𝑡^𝑖(𝑎−𝑡,0)𝑒^{−∫𝑡0𝜆𝑖(𝑎−𝑡+𝜏,𝑧𝜏)𝑑𝜏}𝐵for𝑎⩾𝑡,_𝑖(𝑡−𝑎)𝑒^{−∫𝑎0𝜆𝑖(𝛼,𝑧𝑡−𝑎+𝛼)𝑑𝛼}for𝑡>𝑎, (16) where 𝐵^𝑖(𝑡)=𝑢^𝑖(0,𝑡) for each 𝑖=1,…,𝑛.

Proof. The idea of the proof is analogous to the result included in the paper [9], treating the theme of the age-dependent population problem for one species. Let 𝑢=(𝑢¹,…,𝑢^𝑛) be a solution of the problem up to time 𝑇. Let 𝑢^𝑖(ℎ)=𝑢^𝑖(𝑎₀+ℎ,𝑡₀+ℎ), 𝜆^𝑖(ℎ)=𝜆^𝑖(𝑎₀+ℎ,𝑧_𝑡0+ℎ) then we can rewrite (1) as the equation 𝑑𝑢^𝑖/𝑑ℎ+𝜆^𝑖(ℎ)𝑢^𝑖=0 with the unique solution 𝑢^𝑖𝑎₀+ℎ,𝑡₀+ℎ=𝑢^𝑖𝑎₀,𝑡₀𝑒^{−∫ℎ0𝜆𝑖(𝜂)𝑑𝜂}. (17) Substituting (𝑎₀,𝑡₀)=(𝑎−𝑡,0), ℎ=𝑡 and (𝑎₀,𝑡₀)=(0,𝑡−𝑎), ℎ=𝑎 into (17) yields the formula (16). Applying (16) to (2) and (3) we obtain the operator equations (14) and (15).
To prove the second part of the theorem we should assume that 𝑧^𝑖_𝑡⩾0 and 𝐵^𝑖⩾0 for 𝑖=1,…,𝑛 are continuous functions on the interval [0,𝑇] fulfilling conditions (14) and (15). Let 𝑢 for each component 𝑢^𝑖 be defined on ℝ₊×[0,𝑇] by the formula (16). The function 𝑢 is nonnegative because of (𝐻₇). An easy computation shows that (4) holds, and 𝑢(0,𝑡)=𝐵(𝑡) for 𝑡>0. 𝑢∈𝐿¹(ℝ₊) because 𝜆, 𝛽, and 𝑧_𝑡 are continuous and 𝜑∈𝐿¹(ℝ). It follows from (14)–(16) that (2) and (3) are satisfied. To complete the proof, let us notice that (6) and (16) imply existing 𝐷𝑢 on ℝ₊×[0,𝑇] and the equality (1) holds.

Let us make some estimation for the necessity of the next section. We turn back to the operator equation (14). By the assumptions (𝐻₂) and (𝐻₈) we get 𝐵^𝑖(𝑡)⩽𝛽₀𝑟¹⋯𝑟^𝑛𝑡₀sup_{𝑠∈[−𝑟𝑖,0]}𝐵^𝑖(𝑎+𝑠)𝑑𝑎+𝛽₀𝑟¹⋯𝑟^𝑛Φ^𝑖. (18) Denoting ℬ^𝑖(𝑡)=sup_{𝜏∈[−𝑟𝑖,𝑡]}𝐵^𝑖(𝜏), (19) we obtain ℬ^𝑖(𝑡)⩽𝛽₀𝑟¹⋯𝑟^𝑛𝑡₀ℬ^𝑖(𝑎)𝑑𝑎+𝛽₀𝑟¹⋯𝑟^𝑛Φ^𝑖, (20) and by Gronwall's inequality we get ℬ^𝑖(𝑡)⩽𝛽₀𝑟¹⋯𝑟^𝑛Φ^𝑖𝑒^{𝛽0𝑟1⋯𝑟𝑛𝑡}. (21)

3. Existence and Uniqueness

According to the theorem in Section 2 to solve the population problem up to time 𝑇 it is sufficient to find functions 𝑧=(𝑧¹,…,𝑧^𝑛)∈⨉^𝑛_𝑖=1𝐶⁺([−𝑟^𝑖,𝑇]) and 𝐵=(𝐵¹,…,𝐵^𝑛)∈⨉^𝑛_𝑖=1𝐶⁺([−𝑟^𝑖,𝑇]) satisfying the system of operator equations (14) and (15). We can notice that (14) is a Volterra equation of 𝐵^𝑖 with a unique solution ℬ^𝑖_𝑇(𝑧). Let us define on ⨉^𝑛_𝑖=1𝐶⁺([−𝑟^𝑖,𝑇]) the new operator 𝒵_𝑇=(𝒵¹_𝑇,…,𝒵^𝑛_𝑇) where 𝒵^𝑖_𝑇(𝑧)(𝑡)=^𝑡₀ℬ^𝑖_𝑇(𝑧)(𝑎)𝑒^{−∫𝑡𝑎𝜆𝑖(𝜏−𝑎,𝑧𝜏)𝑑𝜏}+𝑑𝑎_∞0𝜑^𝑖(𝑎,0)𝑒^{−∫𝑡0𝜆𝑖(𝑎+𝜏,𝑧𝜏)𝑑𝜏}𝑑𝑎. (22) We define the previous operator using the system of (15) for 𝑖=1,…,𝑛 with 𝐵=(𝐵¹,…,𝐵^𝑛) replaced by ℬ_𝑇=(ℬ¹_𝑇,…,ℬ^𝑛_𝑇).

Theorem 2. The operator 𝒵_𝑇∶⨉^𝑛_𝑖=1𝐶⁺[−𝑟^𝑖,𝑇]→⨉^𝑛_𝑖=1𝐶⁺[−𝑟^𝑖,𝑇] defined by (22) has a unique fixed point for any 𝑇>0.

Proof. We prove that the operator 𝒵_𝑇 is contracting, and then the assertion will be a consequence of the Banach fixed point theorem. Let us consider the Banach space 𝐶[−𝑟^𝑖,𝑇] with the Bielecki norm ‖𝑓^𝑖‖_𝑇=sup_{𝑡∈[−𝑟𝑖,𝑇]}𝑒^−𝐾𝑡|𝑓^𝑖(𝑡)| for any 𝐾>0. Such norm is equivalent to classical supremum norm ‖⋅‖ in 𝐶[−𝑟^𝑖,𝑇]. Choose 𝑧,𝑧∈⨉^𝑛_𝑖=1𝐶[−𝑟^𝑖,𝑇]. Applying the definition (22) we obtain 𝒵^𝑖_𝑇(𝑧)−𝒵^𝑖_𝑇(𝑧)_𝑇⩽sup_{𝑡∈[−𝑟𝑖,𝑇]}𝑒^−𝐾𝑡×^𝑡₀𝑒^{−∫𝑡𝑎𝜆𝑖(𝜏−𝑎,𝑧𝜏)𝑑𝜏}−𝑒^{−∫𝑡𝑎𝜆𝑖(𝜏−𝑎,}𝑧_{𝜏)𝑑𝜏}ℬ^𝑖_𝑇(𝑧)(𝑎)𝑑𝑎+sup_{𝑡∈[−𝑟𝑖,𝑇]}𝑒^−𝐾𝑡×^𝑡₀𝑒^{−∫𝑡𝑎𝜆𝑖(𝜏−𝑎,}𝑧_{𝜏)𝑑𝜏}ℬ^𝑖_𝑇(𝑧)(𝑎)−ℬ^𝑖_𝑇(𝑧)(𝑎)𝑑𝑎+sup_{𝑡∈[−𝑟𝑖,𝑇]}𝑒^−𝐾𝑡×^∞₀𝑒^{−∫𝑡0𝜆𝑖(𝑎+𝜏,𝑧𝜏)𝑑𝜏}−𝑒^{−∫𝑡0𝜆𝑖(𝑎+𝜏,}𝑧_𝜏^)𝑑𝜏𝜑^𝑖(𝑎,0)𝑑𝑎=𝐼₁+𝐼₂+𝐼₃. (23) Estimating the quantities 𝐼₁,  𝐼₂, and 𝐼₃ we use the assumptions (𝐻₂), (𝐻₇), and (𝐻₈), inequalities (21) and |𝑒^𝑘−1|⩽|𝑘|𝑒^|𝑘|. Thus 𝐼₁⩽sup_{𝑡∈[−𝑟𝑖,𝑇]}𝑒^−𝐾𝑡×^𝑡₀1−𝑒^{−∫𝑡𝑎[𝜆𝑖(𝜏−𝑎,𝑧𝜏)−𝜆𝑖(𝜏−𝑎,}𝑧_{𝜏)]𝑑𝜏}ℬ^𝑖_𝑇(𝑧)(𝑎)𝑑𝑎⩽𝜆₁‖𝑧−𝑧‖_𝑇sup_{𝑡∈[−𝑟𝑖,𝑇]}𝑒^−𝐾𝑡×^𝑡₀^𝑡_𝑎𝑒^{2𝜆0(𝑡−𝑎)}𝑒^𝐾𝜏ℬ^𝑖_𝑇(𝑧𝑎)⩽𝜆)(𝑑𝜏𝑑𝑎₁𝐾𝛽₀𝑟¹⋯𝑟^𝑛Φ^𝑖𝛽₀𝑟¹⋯𝑟^𝑛−2𝜆₀𝑒^{𝛽0𝑟1⋯𝑟𝑛𝑇}‖𝑧−𝑧‖_𝑇,𝐼₃⩽sup_{𝑡∈[−𝑟𝑖,𝑇]}𝑒^−𝐾𝑡×^∞₀^𝑡_𝑎𝑒^{2𝜆0(𝑡−𝑎)}𝑒^𝐾𝜏𝜆₁‖𝑧−𝑧‖_𝑇𝜑^𝑖()⩽𝜆𝑎,0𝑑𝜏𝑑𝑎₁𝐾Φ^𝑖𝑒^2𝜆0𝑇‖𝑧−𝑧‖_𝑇. (24) From (14) and the definition of ℬ_𝑇 we have ℬ^𝑖_𝑇(𝑧)(𝑡)−ℬ^𝑖_𝑇(𝑧𝑡)=)(_{[−𝑟,0]𝑡+𝑠𝑖0}𝛽^𝑖𝑡+𝑠^𝑖−𝑎,𝑧(𝑡+𝑠)⋅𝑒^{−∫𝑖𝑎𝑡+𝑠𝜆𝑖(𝜏−𝑎,𝑧𝜏)𝑑𝜏}×ℬ^𝑖_𝑇(𝑧)(𝑎)−ℬ^𝑖_𝑇(𝑧𝑎)+)(𝑑𝑎𝑑𝑠_{[−𝑟,0]𝑡+𝑠𝑖0}ℬ^𝑖_𝑇(𝑧𝑎)×𝛽)(_𝑖𝑡+𝑠^𝑖−𝑎,𝑧(𝑡+𝑠)×𝑒^{−∫𝑖𝑎𝑡+𝑠𝜆𝑖(𝜏−𝑎,𝑧𝜏)𝑑𝜏}−𝛽^𝑖𝑡+𝑠^𝑖−𝑎,𝑧(𝑡+𝑠)×𝑒^{−∫𝑖𝑎𝑡+𝑠𝜆𝑖(𝜏−𝑎,}𝑧_{𝜏)𝑑𝜏}𝑑𝑎𝑑𝑠+_[−𝑟,0]^∞₀𝜑^𝑖()×𝛽𝑎,0_𝑖𝑎+𝑡+𝑠^𝑖,𝑧(𝑡+𝑠)×𝑒^{−∫𝑖0𝑡+𝑠𝜆𝑖(𝑎+𝜏,𝑧𝜏)𝑑𝜏}−𝛽^𝑖𝑎+𝑡+𝑠^𝑖,𝑧(𝑡+𝑠)×𝑒^{−∫𝑖0𝑡+𝑠𝜆𝑖(𝑎+𝜏,}𝑧_{𝜏)𝑑𝜏}𝑑𝑎𝑑𝑠. (25)(26) Let us denote (26) by 𝑓^𝑖(𝑡); then ℬ^𝑖_𝑇(𝑧)(𝑡)−ℬ^𝑖_𝑇(𝑧)(𝑡)⩽𝛽₀𝑟¹⋯𝑟^𝑛𝑡₀ℬ^𝑖_𝑇(𝑧)(𝑎)−ℬ^𝑖_𝑇(𝑧𝑎)+𝑓)(𝑑𝑎_𝑖(𝑡), (27) and hence, by Gronwall's inequality we have ℬ^𝑖_𝑇(𝑧)(𝑡)−ℬ^𝑖_𝑇(𝑧𝑡)⩽𝑓)(_𝑖(𝑡)+𝛽₀𝑟¹⋯𝑟^𝑛𝑡₀𝑓^𝑖(𝑎)𝑒^{𝛽0𝑟1⋯𝑟𝑛(𝑡−𝑎)}𝑑𝑎. (28) Let us estimate |𝑓^𝑖(𝑡)|𝑓^𝑖(𝑡⩽)_{[−𝑟,0]𝑡+𝑠𝑖0}𝐶₁ℬ^𝑖_𝑇(𝑧𝑎)×𝛽)(_𝑖𝑡+𝑠^𝑖−𝑎,𝑧(𝑡+𝑠)−𝛽^𝑖𝑡+𝑠^𝑖−𝑎,𝑧(𝑡+𝑠)𝑑𝑎𝑑𝑠+_[−𝑟,0]^∞₀𝐶₂𝜑^𝑖()×𝛽𝑎,0_𝑖𝑎+𝑡+𝑠^𝑖,𝑧(𝑡+𝑠)−𝛽^𝑖𝑎+𝑡+𝑠^𝑖,𝑧()⩽𝑡+𝑠𝑑𝑎𝑑𝑠𝐶𝛽₁𝑟¹⋯𝑟^𝑛𝑒^𝐾𝑡‖𝑧−𝑧‖_𝑇×Φ^𝑖𝑡₀𝛽₀𝑟¹⋯𝑟^𝑛𝑒^{𝛽0𝑟1⋯𝑟𝑛𝑎}𝑑𝑎+^∞₀𝜑^𝑖()⩽𝑎,0𝑑𝑎𝐶𝛽₁𝑟¹⋯𝑟^𝑛Φ^𝑖𝑒^𝐾𝑡𝑒^{𝛽0𝑟1⋯𝑟𝑛𝑡}‖𝑧−𝑧‖_𝑇, (29) where 𝐶₁ and 𝐶₂ are positive constants and 𝐶=max{𝐶₁,𝐶₂}. Finally, we have 𝐼₂⩽sup_{𝑡∈[−𝑟𝑖,𝑇]}𝑒^−𝐾𝑡×^𝑡₀𝑓^𝑖(𝑎)+𝛽₀𝑟¹⋯𝑟^𝑛𝑎₀𝑓^𝑖(𝜉)𝑒^{𝛽0𝑟1⋯𝑟𝑛(𝑎−𝜉)}⩽𝑑𝜉𝑑𝑎𝐶𝛽₁𝑟¹⋯𝑟^𝑛Φ^𝑖𝑒^𝐾𝑡‖𝑧−𝑧‖_𝑇×sup_{𝑡∈[−𝑟𝑖,𝑇]}𝑒^{−𝐾𝑡𝑡}₀𝑒^{(𝐾+𝛽0𝑟1⋯𝑟𝑛)𝑎}𝑑𝑎+𝛽₀𝑟¹⋯𝑟^𝑛𝑡₀^𝑎₀𝑒^𝐾𝜉𝑒^{𝛽0𝑟1⋯𝑟𝑛𝑎}⩽𝑑𝜉𝑑𝑎𝐶𝛽₁𝑟¹⋯𝑟^𝑛Φ^𝑖𝐾+𝛽₀𝑟¹⋯𝑟^𝑛‖𝑧−𝑧‖_𝑇𝑒^{𝛽0𝑟1⋯𝑟𝑛𝑇}𝛽1+⁰𝑟¹⋯𝑟^𝑛𝐾. (30) Therefore, we can choose sufficiently large constant 𝐾 for fixed 𝑇 that 𝐼₁, 𝐼₂, and 𝐼₃ are less than 𝐶‖𝑧−𝑧‖_𝑇 with the constant 𝐶∈[0,(1/3𝑛)) (independent of 𝑧 and 𝑧). This shows the contraction of 𝒵_𝑇 and completes the proof.

According to the previous mentioned there exists the exact one solution for any interval [0,𝑇] so the solutions defined on two different intervals coincide on their intersection. By the extension property we have the existence and uniqueness of the 𝑛 populations problem solution for all times.

4. Stability of Equilibrium Age Distribution

A stationary solution 𝑢(𝑎)=(𝑢¹(𝑎),…,𝑢^𝑛(𝑎)) of the model (1)–(4) satisfies the following system of the equations: 𝑢^𝑖₀^′(𝑎)+𝜆^𝑖𝑎,𝑧₀𝑢^𝑖₀𝑧=0,_𝑖0=^∞₀𝑢^𝑖₀(𝑎)𝑢𝑑𝑎,_𝑖0(0)=𝑟¹⋯𝑟^𝑛∞₀𝛽^𝑖𝑎,𝑧₀𝑢^𝑖₀(𝑎)𝑑𝑎 (31) for𝑖=1,…,𝑛.

The population of the whole ecosystem 𝑧₀=(𝑧¹₀,…,𝑧^𝑛₀) and its birthrate 𝐵₀=𝑢₀(0)=(𝑢¹₀(0),…,𝑢^𝑛₀(0)) are constants. The quantity 𝑢₀∈𝐶¹(ℝ₊) is the solution of the system (31), and it will be referred to the equilibrium age distribution. The probability that an individual of 𝑖th population survives to age 𝑎 if the population of the ecosystem is on the constant size level 𝑧₀ can be expressed by 𝜋^𝑖₀(𝑎)=𝑒^{−∫𝑎0𝜆𝑖(𝛼,𝑧0)𝑑𝛼}. (32) The quantity 𝑅^𝑖𝑧₀=𝑟¹⋯𝑟^𝑛∞₀𝛽^𝑖𝑎,𝑧₀𝜋^𝑖₀(𝑎)𝑑𝑎 (33) is the number of offsprings expected to be born to an individual of 𝑖th population when the population of the whole ecosystem equals 𝑧₀. We can formulate the following theorem describing the connection between these three quantities.

Theorem 3. Let 𝑧₀=(𝑧¹₀,…,𝑧^𝑛₀) and 𝑧^𝑖₀>0 for 𝑖=1,…,𝑛, and assume that 𝜋^𝑖₀(⋅),  𝛽^𝑖(⋅,𝑧₀)𝜋^𝑖₀(⋅)∈𝐿¹(ℝ₊) for each 𝑖=1,…,𝑛. Then 𝑅𝑧₀=𝑅¹𝑧₀,…,𝑅^𝑛𝑧₀ (34) with 𝑅^𝑖𝑧₀=1,𝑖=1,…,𝑛 (35) is a necessary and sufficient condition that an equilibrium age distribution exists. The unique equilibrium age distribution 𝑢₀=(𝑢¹₀,…,𝑢^𝑛₀) corresponding to 𝑧₀=(𝑧¹₀,…,𝑧^𝑛₀) is given by 𝑢^𝑖₀(𝑎)=𝐵^𝑖₀𝜋^𝑖₀(𝑎),𝑖=1,…,𝑛, (36) where 𝐵^𝑖₀=𝑧^𝑖₀∫^∞₀𝜋^𝑖₀(𝑎).𝑑𝑎 (37)

Proof. The function (36) is the unique solution of (31)₁ with the initial condition 𝑢^𝑖₀(0)=𝐵^𝑖₀. By (31)₂ we obtain the formula (37) for 𝐵^𝑖₀. An easy computation shows that (35) is equivalent to (31)₃ for each 𝑖=1,…,𝑛.

We now turn to the problem of the equilibrium age distribution stability. We consider “perturbations” 𝜉^𝑖(𝑎,𝑡) and 𝑝^𝑖(𝑡). Let us write 𝑢(𝑎,𝑡)=𝑢₀(𝑎)()=𝑢+𝜉𝑎,𝑡₁₀(𝑎)+𝜉¹(𝑎,𝑡),𝑢²₀(𝑎)+𝜉²(𝑎,𝑡),…,𝑢^𝑛₀(𝑎)+𝜉^𝑛()𝑢𝑎,𝑡_𝑖(𝑎,𝑡)=𝑢^𝑖₀(𝑎)+𝜉^𝑖(𝑎,𝑡), (38)(39) where 𝜉^𝑖(𝑎,𝑡)=𝜑^𝑖(𝑎,𝑡)−𝑢^𝑖₀(𝑎) for 𝑡∈[−𝑟^𝑖,0]. And 𝑧(𝑡)=𝑧₀(𝑡)=𝑧+𝑝₁₀+𝑝¹(𝑡),𝑧²₀+𝑝²(𝑡),…,𝑧^𝑛₀+𝑝^𝑛(𝑡),𝑧_𝑡=𝑧₀+𝑝_𝑡=𝑧¹₀+𝑝¹_𝑡,𝑧²₀+𝑝²_𝑡,…,𝑧^𝑛₀+𝑝^𝑛_𝑡,𝑧^𝑖(𝑡)=𝑧^𝑖₀+𝑝^𝑖(𝑡),𝑧^𝑖_𝑡=𝑧^𝑖₀+𝑝^𝑖_𝑡, (40)(41) where 𝑝^𝑖_𝑡(𝑠)=𝑝^𝑖(𝑡+𝑠) for 𝑠∈[−𝑟^𝑖,0]. Our goal is to formulate relations for 𝜉^𝑖 and 𝑝^𝑖 which guarantee that 𝑢^𝑖 and 𝑧^𝑖₀ obey the basic equations (1)–(4). From (1) we obtain 𝐷𝜉^𝑖(𝑎,𝑡)+𝜆^𝑖₀(𝑎)𝜉^𝑖(𝑎,𝑡)+𝜔^𝑖(𝑎)𝑝^𝑖_𝑡=𝑥^𝑖(𝑎,𝑡) (42) for 𝑥^𝑖(𝑎,𝑡)=−𝐷_𝑧𝑖𝜆^𝑖𝑎,𝑧₀𝑝^𝑖_𝑡𝜉^𝑖(𝑎,𝑡)−Λ^𝑖𝑎,𝑝_𝑡𝐵^𝑖₀𝜋^𝑖₀(𝑎)+𝜉^𝑖(,𝜆𝑎,𝑡)_𝑖0(𝑎)=𝜆^𝑖𝑎,𝑧₀, (43) where 𝜔^𝑖∶ℝ₊×𝐶([−𝑟^𝑖,0])→ℝ is the functional 𝜔^𝑖(𝑎)𝑝^𝑖_𝑡=𝐷_𝑧𝑖𝜆^𝑖𝑎,𝑧₀𝑝^𝑖_𝑡𝐵^𝑖₀𝜋^𝑖₀(𝑎),Λ^𝑖𝑎,𝑝_𝑡=𝜆^𝑖𝑎,𝑧₀+𝑝_𝑡−𝜆^𝑖₀(𝑎)−𝐷_𝑧𝑖𝜆^𝑖𝑎,𝑧₀𝑝^𝑖_𝑡. (44) We conclude from (2) that 𝜉^𝑖(0,𝑡)=^∞_0[−𝑟,0]𝛽^𝑖₀(𝑎)𝜉^𝑖𝑎,𝑡+𝑠^𝑖𝑑𝑠𝑑𝑎+𝜅^𝑖𝑝^𝑖_𝑡+𝜓^𝑖(𝑡), (45) where 𝛽^𝑖𝑎,𝑧₀=𝛽^𝑖₀(𝑎),𝜓^𝑖(𝑡)=^∞_0[−𝑟,0]𝐷_𝑧𝑖𝛽^𝑖𝑎,𝑧₀𝑝^𝑖𝑡+𝑠^𝑖𝜉^𝑖𝑎,𝑡+𝑠^𝑖+𝑑𝑠𝑑𝑎_{∞0[−𝑟,0]}Ω^𝑖()×𝐵𝑎,𝑝𝑡+𝑠_𝑖0𝜋^𝑖₀(𝑎)+𝜉^𝑖𝑎,𝑡+𝑠^𝑖𝜅𝑑𝑠𝑑𝑎_𝑖𝑝^𝑖_𝑡=𝐵^𝑖₀^∞_0[−𝑟,0]𝐷_𝑧𝑖𝛽^𝑖𝑎,𝑧₀𝑝^𝑖_𝑡𝜋^𝑖₀(𝑎)𝑑𝑠𝑑𝑎. (46) Here 𝜅^𝑖∶𝐶([−𝑟^𝑖,0])→ℝ is the functional. Moreover, we have Ω^𝑖𝑎,𝑝(𝑡+𝑠)=𝛽^𝑖𝑎,𝑧₀+𝑝(𝑡+𝑠)−𝛽^𝑖₀(𝑎)−𝐷_𝑧𝑖𝛽^𝑖𝑎,𝑧₀𝑝^𝑖𝑡+𝑠^𝑖. (47) A trivial verification shows that (3) gives the condition for 𝑝 in the form 𝑝^𝑖(𝑡)=^∞₀𝜉^𝑖(𝑎,𝑡)𝑑𝑎. (48) We consider the system (42)–(48) with the initial condition 𝜉^𝑖(𝑎,0)=𝜂^𝑖(𝑎), (49) where 𝜂^𝑖(𝑎)=𝜑^𝑖(𝑎,0)−𝐵^𝑖₀. (50) We first express the solution of (42)–(48) in the form of the matrix equation. Let 𝜉^𝑖 be the solution of (42) up to time 𝑇. Let (𝑎₀,𝑡₀)∈ℝ₊×[0,𝑇], 𝜉^𝑖(ℎ)=𝜉^𝑖(𝑎₀+ℎ,𝑡₀+ℎ), 𝜆^𝑖₀(ℎ)=𝜆^𝑖₀(𝑎₀+ℎ), 𝜔^𝑖(ℎ)=𝜔^𝑖(𝑎₀+ℎ), 𝑝^𝑖_ℎ=𝑝^𝑖_𝑡+ℎ, and 𝑥^𝑖(ℎ)=𝑥^𝑖(𝑎₀+ℎ,𝑡₀+ℎ). Then we can rewrite (42) as the equation 𝑑𝑑ℎ𝜉^𝑖(ℎ)+𝜆^𝑖₀(ℎ)𝜉^𝑖(ℎ)=𝑥^𝑖(ℎ)−𝜔^𝑖(ℎ)𝑝^𝑖_ℎ (51) with the unique solution 𝜉^𝑖𝑎₀+ℎ,𝑡₀+ℎ=𝜉^𝑖𝑎₀,𝑡₀𝑒^{−∫ℎ0𝜆𝑖0(𝑎0+𝜏)𝑑𝜏}+𝑒^{−∫ℎ0𝜆𝑖0(𝑎0+𝜏)𝑑𝜏}×^ℎ₀𝑥^𝑖𝑎₀+𝜏,𝑡₀+𝜏−𝜔^𝑖𝑎₀+𝜏𝑝^𝑖_𝑡0+𝜏×𝑒^{∫𝜏0𝜆𝑖0(𝑎0+𝛼)𝑑𝛼}𝑑𝜏. (52) Substituting (𝑎₀,𝑡₀)=(𝑎−𝑡,0), ℎ=𝑡 for 𝑎⩾𝑡 and (𝑎₀,𝑡₀)=(0,𝑡−𝑎), ℎ=𝑎 for 𝑎<𝑡 yields the formulas 𝜉^𝑖(𝑎,𝑡)=𝜂^𝑖(𝑎−𝑡)𝜋^𝑖()+𝑎−𝑡,𝑎_𝑡0𝑥^𝑖(𝑎−𝑡+𝜏,𝜏)−𝜔^𝑖(𝑎−𝑡+𝜏)𝑝^𝑖_𝜏×𝜋^𝑖₀()𝜉𝑎−𝑡+𝜏,𝑎𝑑𝜏for𝑎⩾𝑡,_𝑖(𝑎,𝑡)=𝑏^𝑖(𝑡−𝑎)𝜋^𝑖₀(𝑎)+^𝑡_𝑡−𝑎𝑥^𝑖(𝜏+𝑎−𝑡,𝜏)−𝜔^𝑖(𝜏+𝑎−𝑡)𝑝^𝑖_𝜏×𝜋^𝑖₀(𝜏+𝑎−𝑡,𝑎)𝑑𝜏for𝑎<𝑡, (53) where 𝜋^𝑖₀()=𝜋𝑎−𝑡,𝑎^𝑖0(𝑎)𝜋^𝑖₀(),𝑏𝑎−𝑡_𝑖(𝑡−𝑎)=𝜉^𝑖(0,𝑡−𝑎). (54)(55) Furthermore, (45) implies that 𝑝^𝑖(𝑡)=^𝑡₀𝜉^𝑖(𝑎,𝑡)𝑑𝑎+^∞_𝑡𝜉^𝑖(𝑎,𝑡)𝑑𝑎. (56) It follows that 𝑝^𝑖(𝑡)+^𝑡₀^∞₀𝜔^𝑖(𝑎)𝑝^𝑖_𝜏𝜋^𝑖₀()−𝑎,𝑎+𝑡−𝜏𝑑𝑎𝑑𝜏_𝑡0𝑏^𝑖(𝜏)𝜋^𝑖₀()=𝑡−𝜏𝑑𝜏_{𝑡0∞𝑡−𝜏}𝑥^𝑖(𝜏+𝑎−𝑡,𝜏)𝜋^𝑖₀()+𝜏+𝑎−𝑡,𝑎𝑑𝑎𝑑𝜏_∞0𝜂^𝑖(𝜏)𝜋^𝑖₀(𝜏,𝜏+𝑡)𝑑𝜏. (57) We conclude from (48) that 𝑏^𝑖(𝑡)=_[−𝑟,0]^{𝑡+𝑠𝑖}₀𝛽^𝑖₀(𝑎)𝜉^𝑖𝑎,𝑡+𝑠^𝑖+𝑑𝑎𝑑𝑠_{[−𝑟,0]∞𝑡+𝑠𝑖}𝛽^𝑖₀(𝑎)𝜉^𝑖𝑎,𝑡+𝑠^𝑖𝑑𝑎𝑑𝑠+𝜅^𝑖(𝑡)𝑝^𝑖_𝑡+𝜓^𝑖(𝑡), (58) and hence −𝜅^𝑖𝑝^𝑖_𝑡+𝑏^𝑖(𝑡)+^𝑡_0−[𝑟,0]^∞₀𝛽^𝑖₀𝑎−𝜏+𝑡+𝑠^𝑖𝜔^𝑖(𝑎)𝑝^𝑖_𝜏×𝜋^𝑖₀𝑎,𝑎−𝜏+𝑡+𝑠^𝑖−𝑑𝑎𝑑𝑠𝑑𝜏_{𝑡0[−𝑟,0]}𝛽^𝑖₀𝑡+𝑠^𝑖−𝜏𝑏^𝑖(𝜏)𝜋^𝑖₀𝑡+𝑠^𝑖=−𝜏𝑑𝑠𝑑𝜏_{[−𝑟,0]𝑡+𝑠𝑖0∞𝑡+𝑠𝑖−𝜏}𝛽^𝑖₀(𝑎)𝑥^𝑖𝜏+𝑎−𝑡−𝑠^𝑖×𝜋,𝜏_𝑖0𝜏+𝑎−𝑡−𝑠^𝑖+,𝑎𝑑𝑎𝑑𝜏𝑑𝑠_{[−𝑟,0]∞0}𝛽^𝑖₀𝜏+𝑡+𝑠^𝑖𝜂^𝑖(𝜏)×𝜋^𝑖₀𝜏,𝜏+𝑡+𝑠^𝑖𝑑𝜏𝑑𝑠+𝜓^𝑖(𝑡). (59) Conditions (57) and (59) imply that 𝐲(𝑡)=𝐲¹(𝑡)⋮𝐲^𝑛(𝑡), (60) where 𝐲^𝑖(𝑡)=𝑝^𝑖(𝑡)𝑏^𝑖(𝑡)𝐲_𝑖𝑡,(𝑠)=𝐲^𝑖(𝑡+𝑠)for𝑠∈−𝑟^𝑖,0,𝑖=1,2,…,𝑛, (61)(62) satisfy the matrix equation 𝐀𝐲_𝑡+^𝑡₀𝐊(𝑡−𝜏)𝐲_𝜏𝑑𝜏=𝐟(𝑡), (63) where 𝐀 and 𝐊 are block diagonal matrices. We have 𝐀𝐀=¹𝐎…𝐎𝐎𝐀²…𝐎⋮⋮⋱⋮𝐎𝐎…𝐀^𝑛𝐊,𝐊=¹𝐎…𝐎𝐎𝐊²…𝐎⋮⋮⋱⋮𝐎𝐎…𝐊^𝑛 (64) for zero matrix 𝐎, 𝐀^𝑖=¹⁰_−𝜅𝑖1, 𝐀^𝑖𝐲_𝑡=^{𝑝𝑖𝑡(0)0}_{−𝜅𝑖𝑝𝑖𝑡𝑏𝑖𝑡(0)}=^{𝑝𝑖(𝑡)0}_{−𝜅𝑖𝑝𝑖𝑡𝑏𝑖(𝑡)} and𝐊^𝑖(𝑡)=∫^∞₀𝜋^𝑖₀(𝑎,𝑎+𝑡)𝜔^𝑖(𝑎)𝑑𝑎−𝜋^𝑖₀(𝑡)∫⁰_−𝑟𝑖∫^∞₀𝛽^𝑖₀(𝑎+𝑡+𝑠)𝜋^𝑖₀(𝑎,𝑎+𝑡+𝑠)𝜔^𝑖(𝑎)𝑑𝑎𝑑𝑠−∫⁰_−𝑟𝑖𝛽^𝑖₀(𝑡+𝑠)𝜋^𝑖₀(𝑡+𝑠)𝑑𝑠. (65)

Since 𝜔(𝑎) is the functional so in the formula of 𝐊(𝑡) we consider 𝜔 as the function with functional values. Furthermore, in the notation (63) looks like the ordinary matrix equation, but in fact it is the operator one. Moreover, 𝐟(𝑡)=𝐟¹(𝑡)⋮𝐟^𝑛(𝑡) (66) for 𝐟^𝑖(𝑡)=𝑓^𝑖₁(𝑡)𝑓^𝑖₂(𝑡)−𝐀^𝑖𝐲^𝑖_1𝑡−^𝑡₀𝐊^𝑖(𝑡−𝜏)𝐲^𝑖_1𝜏𝑓𝑑𝜏_𝑖1(𝑡)=^𝑡₀^∞_𝑡−𝜏𝑥^𝑖(𝜏+𝑎−𝑡,𝜏)𝜋^𝑖₀(𝜏+𝑎−𝑡,𝑎)𝑑𝑎𝑑𝜏+^∞₀𝜂^𝑖(𝜏)𝜋^𝑖₀()𝑓𝜏,𝜏+𝑡𝑑𝜏,_𝑖2(𝑡)=_[−𝑟,0]^{𝑡+𝑠𝑖}₀^∞_{𝑡+𝑠𝑖−𝜏}𝛽^𝑖₀(𝑎)𝑥^𝑖𝜏+𝑎−𝑡−𝑠^𝑖×𝜋,𝜏_𝑖0𝜏+𝑎−𝑡−𝑠^𝑖+,𝑎𝑑𝑎𝑑𝜏𝑑𝑠_{[−𝑟,0]∞0}𝛽^𝑖₀𝜏+𝑡+𝑠^𝑖𝜂^𝑖(𝜏)×𝜋^𝑖₀𝜏,𝜏+𝑡+𝑠^𝑖𝑑𝜏𝑑𝑠+𝜓^𝑖(𝑡) (67) and 𝐲^𝑖₁(𝑡)=^𝑝𝑖₀^(𝑡)|_{[−𝑟𝑖,0]}, 𝐲^𝑖_1𝑡(𝑠)=𝐲^𝑖₁(𝑡+𝑠) for 𝑠∈[−𝑟^𝑖,0].

We can notice that the last two elements of the expression defining 𝐟^𝑖(𝑡) equal zero for 𝑡>𝑟.

We return to the previous deliberation in the next theorem. We take some additional assumptions.(𝐻₁₀)The Fréchet derivatives 𝐷_𝑧𝑖𝜆^𝑖(𝑎,𝑧₀) and the derivation 𝐷_𝑧𝑖𝛽^𝑖(𝑎,𝑧₀) for 𝑖=1,2…,𝑛 as functions of the variable 𝑎 belong to 𝐿^∞(ℝ₊). (𝐻₁₁)(Λ^𝑖(𝑎,𝜌)/‖𝜌‖₀) and (Ω^𝑖(𝑎,𝜌)/‖𝜌‖₀) tend to zero as ‖𝜌‖₀→0 uniformly for 𝑎→0; here ‖⋅‖₀ denotes the supremum norm in the Banach space 𝐶([−𝑟^𝑖,0]). (𝐻₁₂)𝜆_∗=inf_𝑎⩾0{𝜆^𝑖(𝑎,𝑧₀),𝑖=1,2,…,𝑛}>0. Denote 𝑒_𝜃∶(−∞,0]→ℝ by the formula 𝑒_𝜃(𝑠)=𝑒^𝜃𝑠. The exponential asymptotical stability of the model is established by the following theorem.

Theorem 4. Let 𝜋^𝑖₁(𝑎,𝑎+𝑡)=∫_[−𝑟,0]𝛽^𝑖₀(𝑎+𝑡+𝑠^𝑖)𝜋^𝑖₀(𝑎,𝑎+𝑡+𝑠^𝑖)𝑑𝑠, and let 𝜋^𝑖₁(𝑡)=∫_[−𝑟,0]𝛽^𝑖₀(𝑡  +  𝑠^𝑖)  𝜋^𝑖₀(𝑡+𝑠^𝑖)𝑑𝑠. Let us assume that there exists some 𝜇>0 that the equation 1=^∞₀𝑒^−𝑡𝛾𝜋^𝑖₁(𝑡)−𝑑𝑡_∞0𝑒^{−𝑡𝛾∞}₀𝜋^𝑖₀(𝑎,𝑎+𝑡)𝜔^𝑖(𝑎)𝑒_𝛾⋅𝑑𝑎𝑑𝑡1−^∞₀𝑒^−𝑡𝛾𝜋^𝑖₁(𝑡)−𝑑𝑡_∞0𝑒^−𝑡𝛾𝜋^𝑖₀(𝑡)⋅𝑑𝑡_∞0𝑒^{−𝑡𝛾∞}₀𝜋^𝑖₁(𝑎,𝑎+𝑡)𝜔^𝑖(𝑎)𝑒_𝛾𝑑𝑎𝑑𝑡−𝜅^𝑖𝑒_𝛾 (68) has no solution 𝛾 with Re(𝛾)⩾−𝜇. Then there exist real numbers 𝛿>0 and 𝜇>0 such that for any initial data 𝜑 with ‖𝜑−𝑢₀‖_𝐿¹<𝛿, the corresponding solution of the population problem (1)–(4), if it exists for 𝑡>0, satisfies 𝑧_𝑡−𝑧₀=𝑂𝑒^−𝜇𝑡𝑢(),𝑎,𝑡−𝑢₀(𝑎)=𝑂𝑒^−𝜇𝑡foreach𝑎as𝑡→∞. (69)(70)

Proof. Let 𝑢 be the solution of the population problem (1)–(4) for 𝑡>0. In the proof we will use the properties of the generalized Laplace transform. Let 𝜚∶[0,∞)→𝐶([−𝑟^𝑖,0])^∗ be the function with measure value. Define 𝜚̂(𝜃)=^∞₀𝜚(𝑡)𝑒_𝜃𝑒^−𝜃𝑡𝑑𝑡. (71) If we define the convolution by the formula 𝜚∗𝑣_𝜏(𝑡)=^𝑡₀𝜚(𝑡−𝜏)𝑣(𝜏+⋅)𝑑𝜏, (72) where 𝑣(𝑡)=0 for 𝑡<0, then we have the equality analogous to the property of classical Laplace transform, that is, 𝜚∗𝑣_𝜏(𝜃)=𝜚̂(𝜃)𝑣(𝜃). (73) Using generalized Laplace transform to (63) we get 𝐀(𝜃)𝐲(𝜃)+𝐊(𝜃)𝐲(𝜃)=𝐟̂(𝜃), (74) where 𝐀(𝜃)=𝐀¹(𝜃)𝐎𝐀𝐎⋯𝐎₂(𝜃)𝐀⋯𝐎⋮⋮⋱⋮𝐎𝐎⋯^𝑛(𝜃),𝐀^𝑖(𝜃)=10−𝜅^𝑖𝑒_𝜃1,𝐊(𝜃)=𝐊¹(𝜃)𝐎𝐊𝐎⋯𝐎₂(𝜃)𝐊⋯𝐎⋮⋮⋱⋮𝐎𝐎⋯^𝑛(𝜃),𝐊^𝑖(𝜃)=∫^∞₀𝑒^−𝑡𝜃∫^∞₀𝜋^𝑖₀(𝑎,𝑎+𝑡)𝜔^𝑖(𝑎)𝑒_𝜃𝑑𝑎𝑑𝑡−𝜋^𝑖₀(𝜃)∫^∞₀𝑒^−𝑡𝜃∫^∞₀𝜋^𝑖₁(𝑎,𝑎+𝑡)𝜔^𝑖(𝑎)𝑒_𝜃𝑑𝑎𝑑𝑡−𝜋^𝑖₁(𝜃) (75) for 𝜋^𝑖₁(𝑎,𝑎+𝑡)=∫_[−𝑟,0]𝛽^𝑖₀(𝑎+𝑡+𝑠^𝑖)𝜋^𝑖₀(𝑎,𝑎+𝑡+𝑠^𝑖)𝑑𝑠 and 𝜋^𝑖₁(𝑡)=∫_[−𝑟,0]𝛽^𝑖₀(𝑡+𝑠^𝑖)𝜋^𝑖₀(𝑡+𝑠^𝑖)𝑑𝑠. From the equality (68) we conclude that the matrix 𝐀+𝐊 has an analytic inverse for Re(𝜃)⩾−𝜇. Moreover, we have 𝐀⁻¹(𝜃)=𝐀¹⁻¹(𝜃)𝐎𝐀𝐎⋯𝐎₂₋₁(𝜃)𝐀⋯𝐎⋮⋮⋱⋮𝐎𝐎⋯^𝑛−1(𝜃),𝐀^𝑖−1(𝜃)=𝜅10_𝑖𝑒_𝜃1, (76) so 𝐀 is analytic for 𝜃∈ℂ.
The solution of the matrix equation (63) exists and is given by 𝐲(𝑡)=𝐀⁻¹𝐟_𝑡+^𝑡₀𝐉(𝑡−𝜏)𝐟_𝜏𝑑𝜏, (77) where 𝐟_𝜏(𝑠)=^𝐟1⋮_𝐟^(𝜏+𝑠)_{𝑛(𝜏+𝑠)}, 𝐟^𝑖_𝜏(𝑠)=𝐟^𝑖(𝜏+𝑠) for 𝑠∈[−𝑟^𝑖,0] and 𝐉(𝑡)=1𝑒2𝜋^{−𝜇𝑡∞}_−∞𝑒^𝑖𝜁𝑡𝐉−𝜇+𝑖𝜁𝑑𝜁. (78)𝐉=(𝐀(𝜃)+𝐊(𝜃))⁻¹−𝐀⁻¹(𝜃) is the Laplace transform of the function 𝐉 and |𝐉(𝑡)|⩽𝐶₁𝑒^−𝜇𝑡, (79) where |𝐉| denotes the sum of moduli of the elements of the matrix 𝐉. Here and in the whole proof 𝐶₁, 𝐶₂,  𝐶₃,… denote positive constants.
By (𝐻₁₂), (32), and (54) 𝜋^𝑖₀(𝑎)⩽𝑒^{−𝜆∗𝑎},𝜋^𝑖₀𝑎^′,𝑎⩽𝑒^{−𝜆∗(𝑎−𝑎′)}𝑎^′⩽𝑎. (80)
For 𝜆_∗>𝜇 and ‖𝛽₀‖_𝐿^∞<∞ we can estimate that ‖𝐟(𝑡)‖⩽𝐶₂^𝑡₀𝑒^{−𝜆∗(𝑡−𝜏)}‖𝑥(⋅,𝜏)‖_𝐿¹+𝜂𝑑𝜏_𝐿1𝑒^{−𝜆∗𝑡}+𝜓(𝑡)_𝐿1,‖𝐲(𝑡)‖⩽𝐶₃𝜂_𝐿1𝑒^−𝜇𝑡+𝜓(𝑡)_𝐿1+^𝑡₀𝑒^{−𝜇(𝑡−𝜏)}𝜓(𝜏)_𝐿1+‖𝑥(⋅,𝜏)‖_𝐿¹𝑑𝜏, (81)(82) where ‖𝐲(𝑡)‖=(∑^𝑛_𝑖=1(𝑝^𝑖)²+(𝑏^𝑖)²)^1/2, 𝜉(⋅,𝑡)_𝐿1⩽^𝑡₀𝑒^{−𝜆∗𝜏}‖𝑏(𝑡−𝜏)‖𝑑𝜏+𝜂_𝐿1𝑒^{−𝜆∗𝑡}+^𝑡₀𝑒^{−𝜆∗(𝑡−𝜏)}‖𝑥(⋅,𝜏)‖_𝐿¹+‖𝜔‖_𝐿¹𝑝_𝜏0𝑑𝜏⩽𝐶₄𝜂_𝐿1𝑒^−𝜇𝑡+^𝑡₀𝑒^{−𝜇(𝑡−𝜏)}×𝜓(𝜏)_𝐿1+‖𝑥(⋅,𝜏)‖_𝐿¹𝑑𝜏. (83) Let 𝜀>0. By (𝐻₁₁) there exists 𝛿=𝛿(𝜀) such that Λ𝑎,𝑝_𝑡⩽𝜀𝑝_𝑡0,Ω𝑎,𝑝_𝑡⩽𝜀𝑝_𝑡0 (84) for ‖𝑝_𝑡‖₀<𝛿. By the previous inequalities and the fact that ‖𝐷_𝑧𝜆(𝑎,𝑧₀)‖_𝐿∞, ‖𝐷_𝑧𝛽(𝑎,𝑧₀)‖_𝐿∞, and ‖𝑢₀‖_𝐿¹ are finite, we can estimate that ‖𝑥(⋅,𝑡)‖_𝐿¹,𝜓(𝑡)_𝐿1⩽𝐶₅𝑝_𝑡0⋅𝜉(⋅,𝑡)_𝐿1+𝜀𝑝_𝑡0. (85) We require that 𝜀<1 and 𝛿<𝜀. Hence ‖𝑝_𝑡‖₀<𝜀. According to the previous remarks, we have ‖𝑥(⋅,𝑡)‖_𝐿¹,𝜓(𝑡)_𝐿1⩽𝐶₅(𝑡),𝜎(𝑡)=𝜉(𝜀𝜎⋅,𝑡)_𝐿1+𝑝_𝑡0. (86) Thus (82)–(86) imply that 𝜎(𝑡)⩽𝑀𝜂_𝐿1𝑒^−𝜇𝑡+2𝜀^𝑡₀𝑒^{−𝜇(𝑡−𝜏)}𝜎(𝜏)𝑑𝜏, (87) where 𝑀>0. Gronwall's inequality yields 𝜎(𝑡)⩽𝑀𝜂_𝐿1𝑒^{(−𝜇+2𝑀𝜀)𝑡}. (88) Let us choose 𝛿 and 𝜀>0 such that 𝜇=𝜇−2𝑀𝜀>0 and 𝛿<min(𝛿(𝜀),𝛿(𝜀)/𝑀). Assume that ‖𝜑−𝑢₀‖_𝐿¹=‖𝜂‖_𝐿¹<𝛿. In that case, from what has already been proved, it follows that ‖𝑥(⋅,𝑡)‖_𝐿¹, ‖𝜓(𝑡)‖_𝐿¹, 𝜎(𝑡) and ‖𝑏(𝑡)‖ are 𝑂(𝑒^−𝜇𝑡). Therefore (41), (86)₂ imply (69) and (53)₂, (80), and (39) imply (70).

Example 5. Let 𝑛=2. Let us consider 𝜆^𝑖(𝑎,𝑧)=𝜆^𝑖𝑎,𝑧¹₁,𝑧¹₂,𝑧²₁,𝑧²₂,𝛽^𝑖(𝑎,𝑧)=𝛽^𝑖𝑎,𝑧¹,𝑧², (89) where 𝜆^𝑖(𝑎,𝑧)=𝜆^𝑖𝑧¹(0),𝑧¹−𝑟¹,𝑧²(0),𝑧²−𝑟²,𝛽^𝑖(𝑎,𝑧)=𝛽^𝑖𝑧¹,𝑧²𝑒^{−𝛼𝑖𝑎} (90) with 𝜆^𝑖>0,  𝛽^𝑖>0,  𝑟^𝑖>0,  𝛼^𝑖>0 for 𝑖=1,2. Let us define 𝜆^𝑖₀=𝜆^𝑖𝑧¹₀,𝑧¹₀,𝑧²₀,𝑧²₀,𝛽^𝑖₀=𝛽^𝑖𝑧¹₀,𝑧²₀𝑒^{−𝛼𝑖𝑎} (91) for an arbitrary 𝑧₀=(𝑧¹₀,𝑧²₀)∈ℝ². From (32) and (33) we conclude that 𝜋^𝑖₀(𝑎)=𝑒^{−𝑎𝜆𝑖0},𝑅^𝑖𝑧₀=𝑟¹𝑟²𝛽^𝑖₀𝛼^𝑖+𝜆^𝑖₀. (92) By Theorem 3 there exists an equilibrium age distribution 𝑢^𝑖₀(𝑎)=𝑧^𝑖₀𝜆^𝑖₀𝑒^{−𝜆𝑖0𝑎} where 𝐵₀=𝑧^𝑖₀𝜆^𝑖₀ if and only if 𝑟¹𝑟²𝛽^𝑖₀/(𝛼^𝑖+𝜆^𝑖₀)=1. To investigate the stability of the equilibrium age distribution we consider (68). Let 𝑧₀=(𝑧¹₀,𝑧²₀) be the solution of the system of the equations 𝑅^𝑖𝑧₀=1,𝑖=1,2. (93) First, let us notice that 𝜔^𝑖𝑒_𝛾=𝜕𝜆^𝑖𝜕𝑧^𝑖₁+𝜕𝜆^𝑖𝜕𝑧^𝑖₂𝑒^{−𝛾𝑟𝑖}𝜆^𝑖₀𝑧^𝑖₀𝑒^{−𝑎𝜆𝑖0} (94) with the derivatives in the point 𝑧₀. It is easy to notice that |𝜕𝜆^𝑖/𝜕𝑧^𝑖₁+(𝜕𝜆^𝑖/𝜕𝑧^𝑖₂)𝑒^{−𝛾𝑟𝑖}| is bounded on half-plane Re(𝛾)>𝜇 for every real 𝜇 (also negative). Analogously 𝜅^𝑖𝑒_𝛾=𝜆^𝑖₀𝛽^𝑖₀⋅𝑧^𝑖₀𝑟^𝑖⋅𝜕𝛽^𝑖𝜕𝑧^𝑖𝑧₀⁰_−𝑟𝑖𝑒^𝛾𝑠𝑑𝑠. (95) It is also possible to prove that for every real 𝜇 the function |𝜅^𝑖(𝑒_𝛾)| is bounded on half-plane Re(𝛾)>𝜇. Consider the right-hand side of (68). All components can be presented in the form 𝐴,𝐵+𝛾 (96) where 𝐴 is bounded and 𝐵⩾max{𝛼^𝑖,𝜆^𝑖₀}. Let 𝜇<𝐵, and let Re(𝛾)>−𝜇: 𝐴=𝐴𝐵+𝛾≤𝐴𝐵−𝜇+𝜇+𝛾𝛾≤𝐴𝐵−𝜇+𝜇+Re2𝐵−𝜇𝛾.𝜇+Re (97) Fix 𝜎>0. For 𝐵 sufficiently large there exist some 𝜇 and 𝜇<𝜇 that 𝐴2𝐵−𝜇𝜇−𝜇<𝜎. (98) For Re(𝛾)>−𝜇 we have 𝐴⩽𝐴𝐵+𝛾2𝐵−𝜇𝛾⩽𝐴𝜇+Re2𝐵−𝜇𝜇−𝜇<𝜎. (99) In consequence for 𝛼^𝑖 and 𝜆^𝑖₀ sufficiently large, we can find that 𝜇>0 such that the modulus of the right-hand side of (68) is less than 1 for every 𝛾, for which Re(𝛾)>−𝜇.

Acknowledgment

The first author acknowledges the support from Bialystok University of Technology (Grant no. S/WI/2/2011).