Electronic Journal of Qualitative Theory of Differential Equations

In this work we study a Nicholson-type periodic system with variable delay, density-dependent mortality and linear harvesting rate. Using the topological degree and Lyapunov stability theories, we obtain sufficient conditions that allow us to demonstrate the existence of periodic solutions for the Nicholson-type system and, under suitable conditions, the uniqueness and local exponential stability of the periodic solution is established. We illustrate our results with an example and numerical simulations.


Introduction
In recent years, the question of the existence of periodic solutions for Nicholson-type systems with periodic coefficients has received the attention of many researchers. This class of systems of differential equations with delays was introduced as a coupled patch population model for marine protected areas and B-cell chronic lymphocytic leukemia [7]. However, it has been pointed out that the new models applied to the fishery must consider nonlinear density-dependent mortality rates [6]. Consequently, research on Nicholson-type equations and systems with density-dependent mortality has developed rapidly. But despite that, few studies have considered periodic Nicholson models with density-dependent mortality and harvesting. The goal of this article is to investigate the existence and stability of positive periodic solutions for a m-dimensional Nicholson-type system with periodic coefficients, nonlinear mortality rates, and linear harvesting.

The Nicholson models
In [16] Gurney, Blythe and Nisbet proposed a model to describe the behavior of a population of flies that had been studied in the 1950s by Nicholson [27]. The model corresponds to the following delayed differential equatioṅ where x is the density of the adult population, m is the per capita mortality rate, b the maximum birth rate, τ is the time to maturity and γ indicates where the unimodal function reaches its maximum. Equation (1.1) is known as the Nicholson model. In [7] Berezansky, Idels and Troib studied the dynamics of metapopulation models with migration between two patches. Within the models studied, the authors considered a model of a marine population, with an age structure that inhabits two areas, one protected and the other for extraction. From this model, they obtained the system of differential equations with delay:ẋ (1.2) where x i corresponds to the densities of adult populations, m i are the per capita mortality rates, d i are the diffusion rates between patches, b i are the maximum birth rates, γ i indicates where the unimodal functions reaches its maximum, τ is the time to maturity, and h is the harvesting rate. Due to the presence of a nonlinear birth rate that considers delay, models similar to (1.2) are known as Nicholson-type systems. The model (1.2) has been extended to the non-autonomous case to consider variations due to the passage of time, such as the seasons of the year, which has led to the study of periodic and almost periodic solutions, see [14,15,22,28,29,35].
Since the model (1.2) allows predicting the dynamics of an adult population, it is relevant to include some types of harvesting in them so that they can be applied in models of fishery or agricultural livestock production. Different authors have considered Nicholson-type equations and systems with linear harvesting [13,24,38] and nonlinear harvesting [1,4,5] among others.
Berezansky, Braverman, and Idels in [6] mention that for marine populations at low densities it is appropriate a linear model of density-dependent mortality and that new fishery models must consider nonlinear density-dependent mortality rates. Afterward, research on Nicholson-type equations and systems with density-dependent mortality has been developing rapidly, see [3,8,9,19,23,25,30,33]. However, the study of periodic Nicholson models with density-dependent nonlinear mortality and harvesting terms have not yet been sufficiently explored and this work aims to contribute in this direction.

Novelty of this work
We consider a Nicholson-type system with nonlinear density-dependent mortality, linear harvesting terms, and several concentrated delays of the form where r(x) = x exp(−x), and δ ij , c ij , b ij , τ ij , h i : R → (0, +∞) , i = 1, . . . , m, j = 1, . . . , n, are bounded, continuous and ω-periodic functions. Note that the above system includes the case where each patch considers a different Rickertype function, namely r i (y i ) = y i e −γ −1 i y i . In fact, in this case the system (1.3) is obtained by making the change of variable y i = γ i x i .
Our objective is to apply topological degree and Lyapunov stability theory to the system (1.3) to determine the conditions that guarantee the existence and exponential stability of periodic solutions of the system.

Outline
Section 2 deals with fundamental preliminary aspects of this work, particularly the theory of differential equations with delay and a theorem of continuation of the topological degree; In addition, a result of the existence of solutions and a priori estimates are obtained. Section 3 establishes the main results of this work: Theorem 3.1 provides sufficient conditions for the existence of positive periodic solutions, while Theorems 3.3 and 3.5 prove the local asymptotic and exponential stability, respectively. Section 4 focuses on an example and its numerical simulations. Section 5 is dedicated to the conclusions and discussion of the results, particularly the possible extension of the present study to one involving nonlinear harvesting terms previously considered in population models, see [18,34].

Delay differential equations
Time delays occur naturally in many population dynamical models and their presence is due, among others, to factors like sexual maturity or gestation. Mathematical models with timedelays has a significant role in population dynamics, we refer the reader to [12,26,32,36]. Delayed differential equations may exhibit more complex dynamics than ODE's because of the presence of delay may induce a Hopf bifurcation, periodic and oscillatory solutions or chaos, see [17,21,36].
We introduce some definitions and notation for delay differential equations. For τ ≥ 0, we consider C = C([−τ, 0], R m ) the Banach space with the norm ϕ τ = sup −τ≤θ≤0 ϕ(θ) , where · is the maximum norm in R m . Any vector v ∈ R m is identified in C with the constant function v(θ) = v for θ ∈ [−τ, 0]. A general system of functional differential equations take the formẋ where f : R × C ⊃ D → R m and x t corresponds to the translation of a function x(t) on the interval [t − τ, t] to the interval [−τ, 0], more precisely x t ∈ C is given by we say x(t; 0, φ) is a solution of system (2.1) with initial value φ at 0 if there is an A > 0 such that x(t; 0, φ) is a solution of equation (2.1) on [−τ, A) and x 0 (t; 0, φ) = φ. In addition, for a given continuous and bounded function f ∈ C(R, R) we will denote by f + and f − respectively, the supremum and infimum of f over R. Now, for system (1.3) we consider Since nonnegative solutions are significant for population models, the following subsets of C are often introduced : Proof. We will denote by F i (t, x(t), x(t − τ i1 (t)), . . . , x(t − τ ij (t))) the right hand side of system (1.3) and x(t) = (x 1 (t), . . . , x m (t)) T , then (1.3) can be written as, . Now, applying Theorems 3.1 and 3.2 of [36], it follows that the system (1.3) has a unique solution defined over a maximal interval, for each initial condition φ ∈ C + . In order to show that x(t; 0, φ) takes nonnegative values, we fix i ∈ {1, . . . , m} and t in the maximal interval, in addition we assume that entries of the function F are nonnegative vectors while x ∈ R m + is such that x i = 0, then Consequently, each nonnegative initial condition φ has a corresponding solution x(t; 0, φ) that takes nonnegative values for t in the maximal interval. Now we will prove that the solutions of (1.3), corresponding to nonnegative initial conditions, are defined for all t ≥ 0. Otherwise, they would be defined over an interval Whence, integrating the above estimation we obtain This estimates ensure that A = +∞, because if A < +∞ then |x(t)| → ∞ as t → A, contradicting the estimates.

Topological degree and periodic functions
We begin this subsection by recalling some definitions and notations that will be used in this work. The closure and the boundary of a subset A of a topological space will be denoted respectively by A and ∂A. Let the Banach space of the continuous vector functions ω periodic with the norm It is useful consider the usual notation for the natural embedding R m → C ω given by y → y, where y(t) = y for t ∈ R. Given a continuous function and ω periodic f ∈ C(R, R) notice that f + and f − coincide, respectively, with the maximum and the minimum value of f over the interval [0, ω].
The existence of periodic solutions of the system (1.3) will be proved as a consequence of a general continuation theorem, see [2, Theorem 6.3], in our case we consider: has no solutions on ∂Ω for λ ∈ (0, 1).
Then there exist at least one solution of (1.3) in Ω.
To study conditions ii) and iii) is useful introduce additional notation, let the i-th opposite faces. Condition iii) of the lemma 2.2 will be obtained by the construction of an affine isomorphism homotopic to g combined with the homotopy invariance property of the Brouwer degree.

A priori bounds
To prove the existence of a periodic solution of (1.3) by using the theory of topological degree we need to find some a priori bounds for any ω-periodic solution of the system (2.3). Next, we will state some propositions related to upper and lower a priori bounds that will be useful when proving the existence of positive periodic solutions of (1.3). To obtain the existence of upper bounds for the solutions of the system (2.3) we consider the following assumption: (H1) The coefficients of the system satisfy: Proposition 2.3. If (H1) holds, then every non-negative ω-periodic solution of (2.3) is bounded above for any λ ∈ (0, 1).
Now, combining the monotonicity of the map u → δu c+u , the assumptions over the functions b ij (·), δ ij (·), c ij (·), h i (·) and, the fact that r(u) ≤ 1 e for u ∈ R + we obtain Next, adding and subtracting the terms The above inequality implies On the other hand, (H1) and the continuity of the coefficients imply that there is ζ > 0 such that Now, for R 0 taking the minimum in (2.4), by using the estimations (2.5) and (2.6) we obtain the contradiction Consequently there is a positive number R 0 such that x i (t) < R 0 , for t ∈ R and i = 1, 2, . . . , m. (2.7) To study the a priori lower bounds for the solutions of the system (2.3) we will proceed in a similar way to the proof of the proposition 2.3, but this time the key hypothesis is: (H2) For i = 1, 2, . . . , m we have: Proof. Consider ε = min{x − 1 , x − 2 , . . . , x − m } and, without loss of generality, we suppose that Since (H1) holds, proposition 2.3 implies that the periodic solutions of (2.3) are bounded from above by R 0 . We assume that R 0 ≥ 1 and consider ρ 0 as the unique value in (0, 1] such that r(ρ 0 ) = r(R 0 ). We may suppose that ε ≤ ρ 0 since otherwise, we have trivially a lower bounds for the solutions of (2.3), from ρ 0 < x i (t), for t ∈ R. Now, since ε ≤ ρ 0 , it follows By adding and subtracting the terms Since ε > 0, the above inequality is equivalent to On the other hand, (H2) and the continuity of the coefficients imply that there is ζ > 0 such that Note that there exists 0 < ε 1 such that Therefore, for ε > 0 arbitrarily small values we obtain Consequently there is a positive number ε 0 such that ε 0 < x i (t) < R 0 , for t ∈ R and i = 1, 2, . . . , m.

Results
In this section, we address the problem of the existence and local stability of positive periodic solution for (1.3). We prove the existence of at least one periodic solution of the system (1.3) under assumptions (H1) and (H2) by using the degree topological theory. Proof. The proof of this result is supported by lemma 2.2. Since (H1) and (H2) hold, we apply propositions 2.3 and 2.4 to obtain lower and upper bounds for the periodic solutions of (2.3) for all λ ∈ (0, 1). Next define the set Ω ⊂ C ω as where the positive constants R 0 and ε 0 are, respectively, the upper and lower bounds given by propositions 2.3 and 2.4, we note that Ω ∩ R m = (ε 0 , R 0 ) m . As a consequence of these propositions, it follows that the system (2.3) has no solution in ∂Ω for any λ ∈ (0, 1). We will prove that there are positive constants ε and R such that g(x) = 0 for x ∈ ∂I, where I = [ε, R] m . We recall that, for i = 1, 2, . . . , m and x = (x i ) ∈ R m , we have From the definition of g i (x), considering the notation 1 = (1, 1, . . . , 1), it follows that for z ∈ I − i we obtain Analogously to the estimates made in the proof of proposition 2.4, we deduce that From (H2), it follows that there exists some 0 < ε 1 such that Therefore, there exists a positive number ε 1 such that if ε ≤ ε 1 we have On the other hand, if z ∈ I + i then Since r(R) ≤ 1 e for R ∈ R + and analogously to the estimates made in the proof of proposition 2.3, for z ∈ I + i we obtain From (H1), it follows that there exists some R > R 0 such that Hence there is R 1 > 0 such that if R ≥ R 1 , then We have proved that if ε < ε 1 and R > R 1 , then g(x) = 0 for x ∈ ∂I, where I = [ε, R] m . We claim that g is homotopic to an affine isomorphism. In fact we consider A : where b ∈ R m and the diagonal matrix M ∈ M m×m are completely defined by the systems of linear equation It follows immediately that which is a homotopy between A and g. Since sign g(I + i ) = sign A(I + i ) and sign g( σ) does not vanish on ∂I for any σ ∈ [0, 1], and we conclude that g is homotopic to the affine isomorphism A. The homotopy invariance property of Brouwer degree implies that and by the definition of Brouwer degree it follows that Finally we apply Lemma 2.2 to conclude that the system (1.3) has at least one solution x(t) ∈ Ω.

Remark 3.2.
Several types of delayed harvesting terms have been considered for the Nicholson scalar equation. If we modify the harvesting terms h i (t)x i (t) in our model to delayed terms similar to those used in the work of Qiyuan Zhou in [38], then we obtain the system Next, we will address the asymptotic and exponential stability of the system (1.3). As is common in the literature on Nicholson-type models, our results are obtained by constructing appropriate Lyapunov functions. We define the region of stability of the solutions of our system as the set (3.6) To achieve our stability results, we assume the following: (H3) The delays involve in the model (1.3) are continuously differentiable and satisfy: (H4) For i = 1, 2, . . . , m we have Now we state and prove our first stability theorem. Proof. Let x(t) = (x i (t)) and y(t) = (y i (t)) two solutions in B of system (1.3). We consider the functions: Calculating the upper right Dini derivative of V i (t) along the solutions of (1.3), since 0 ≤ x i (t), y i (t) ≤ K i and |r (x)| ≤ 1 for x ∈ [0, +∞), then proceeding similarly to theorem 2 in [31] we have Notice that assumption (H3) implies that hence we obtain the following estimate , and by a straightforward computation of the corresponding sums it follows Hypothesis (H4) ensure the existence of a positive constant µ such that then we get Therefore, all solution of the system (1.3) in B converge to an ω-periodic solution, hence there is a unique periodic solution of (1.3) in B.

Remark 3.4.
Note that in the proof of theorem (3.3), we use arguments similar to those presented in the proof of theorem (4.5) of [37]. Both results are supported by considering the derivative of Dini and the definition of an adequate Lyapunov functional, in addition to the uniform continuity of the integrands of (3.7) of our proof, equivalent to the integrand given in (4.13) of the proof used in [37]. These are key aspects in the literature on stability in Nicholson-type models, see for instance [13] and references therein.
In order to state and prove our second stability theorem we define, for i = 1, . . . , m, the continuous functions G i : R → R given by Notice that hypothesis (H4) ensures that G i (0) > 0 for each i = 1, . . . , m, furthermore, the continuity of G i guarantees the existence of positive constants r i , such that and we define λ 0 := min 1≤i≤m {r i }, so G i (λ 0 ) > 0 for i = 1, . . . , m. Proof. We consider x(t) = (x i (t)) and y(t) = (y i (t)) two arbitrary solutions in B of system (1.3) and we define the functions: Calculating the upper right Dini derivative of W i (t) along the solutions of model (1.3) we have Replacing x i and y i given in the system, applying triangular inequality, considering (H3), 0 ≤ x i (t), y i (t) ≤ K i , |r (x)| ≤ 1 for x ∈ [0, +∞) and grouping we obtain Extending the sum for i = 1 to m and grouping terms we obtain that the Lyapunov functional We fix λ = λ 0 = min 1≤i≤m {r i }, since (3.8) and (3.9) hold we deduce that It follows that W(t) is decreasing for all t > 0 along the solutions of system (1.3), consequently we have and the exponential convergence it is obtained for solutions of (1.3) in B.

Examples
In this section we show an example of the asymptotic stability of the solution and include numerical simulations performed in R software using the library PBSddesolve, see for instance [11]. In this example x i is the density of biomass in patch i, s(t) = sin(2πt/365), c(t) = cos(2πt/365), and i ∈ {1, 2, 3}.

Conclusion and further work
A Nicholson-type system with nonlinear density-dependent mortality and linear harvesting has been studied in this paper. Based on the theory of topological degree, has been obtained sufficient conditions for the existence of a positive periodic solution of the model. In addition, by using the Lyapunov-Krasovskii functional method, the uniqueness, stability, and exponential stability of the Nicholson-type system were addressed. Numerical simulations were performed based on an example to illustrate the results obtained.
Among the projections of this work, we will focus on the possible extension of the present study to one involving nonlinear harvesting terms. We recall that in the works [1,4,5] advances in this direction have been developed. However, from the point of view of applications, it seems more realistic to consider the harvesting terms, proposed by Clark and Mangel in [10], of the form where q is the catch coefficient, E is the external effort dedicated to the harvest, c and are constants. Population models with terms of this type have been studied in [18,34]. Thus, a new version of the system (1.3) naturally arises, this time with these nonlinear harvesting terms as a new research goal. We anticipate that the main aspects to take into account when applying the methods presented in this work to these nonlinear terms is to search for alternative hypotheses to (H2) and (H4), which can be deduced after a careful reading of this work.