Qualitative Analysis of a Model of Renewable Resources and Population with Distributed Delays

: This work generalizes Matsumoto et al.’s dynamic model of population and renewable resources by substituting a distributed delay for the time delay. It is proved that the equilibrium point may lose or gain local stability, allowing for the observation of alternating stability/instability areas if some conditions hold.


Introduction
Following the findings of archeological studies, many Pacific islands follow similar evolutionary patterns in terms of natural resources and population dynamics. The research by Brander and Taylor [1] (BT henceforth) examines the archeological and anthropological findings from Easter Island as a case in point from an economic standpoint. They present a general equilibrium model of renewable resources and population dynamics to explain how and why Easter Island grew and fell throughout the 1400 years between the 4th century and the middle of the 18th century. Their findings suggest that an economic model linking resources and population dynamics may be able to explain not only the sources of past historical evolution discovered on these small islands, but also the possibility of sustainable growth for our global economy in an era in which a rapidly increasing population and a rapidly degrading environment are becoming serious problems. The analysis included inside the BT model has been expanded in several ways. Coats and Dalton [2] investigated the influence of market institutions and various property-rights frameworks on the distribution of wealth. According to Reuveny and Decker [3], the long-run dynamics of Easter Island were altered by technical advancement and population management reform in the past. However, there has not been much disclosed about the "history" of Easter Island. It is thought that a small group of Polynesians came on the island about the year 400, that deforestation began around the year 1000, that the majority of the sculptures were carved between the years 1000 and 1400, and so on. BT and other academics are attempting to replicate a dynamic pattern of natural resources and population based on this "common knowledge" in their study. However, "wisdom" is still just one of the numerous plausible possibilities, and it has not yet been proven beyond reasonable doubt. For instance, new reconsiderations of archeological material on the island, according to Intoh [4], indicate that the arrival time of the Polynesians on the island cannot be determined with certainty. It is merely a guess that a small number of people arrived between 410 and 1270 AD. Compared to traditional knowledge, this discovery is a jarring departure. In other words, we may have a different history and an evolutionary pattern of Easter Island than what is now known, although the historical data is consistent. It is thus essential to develop a model for small islands capable of generating a variety of dynamic patterns to cope with the confusing qualities of archeological data. Agricultural output may be critical to economic activity in a preindustrial economic civilization. It is well-known that a vital aspect of this kind of agricultural production is the substantial time lag between when producers decide to sow seeds in the fields and when they harvest the crops. Thus, it is logical to ask how delayed output has influenced the evolutionary pattern of a small island economy. Following this consideration, Matsumoto et al. [5] incorporate a delay in manufacturing into the BT model and discover the destabilizing effect caused by the delay in production on the evolution of a small island economy; namely, there is a critical value of the delay for which a loss of stability occurs. In the current literature, time delays are treated as fixed or continuously distributed (distributed delay henceforth). The former relates to economic situations when there is a specified time gap for the actors involved. This latter is suited for economic scenarios where actors' delays vary. A fundamental issue is the unknown nature of time delays. Distributed delays, on the other hand, are the weighted average of all prior data from time zero to the present. Thus, distributed delays better describe time-delayed economic systems (see MacDonald [6]). Additionally, Caperon [7] shows there is some experimental evidence that they are more precise than those with instantaneous time delays. Cushing [8] pioneered and popularized distributed delays in mathematical biology, while Invernizzi and Medio [9] brought distributed delays to mathematical economics. In this paper, the authors reconsider Matsumoto et al.'s model [5] replacing time delay with distributed delay and investigate the dynamic effects of such lags on the adjustment process in population and stock of natural resources. Analytically, our system of functional delay differential equations has a single equilibrium point. The qualitative study demonstrates that this critical point loses or gains local stability when the lag duration increases or decreases. If certain circumstances are met, a series of intervals in which zones of stability and instability alternate may occur. This indicates that it may have significant challenges stabilizing the economic system. Another significant possibility that emerges from the qualitative investigation is the occurrence of Hopf bifurcation-induced limit cycles. The paper is organised as follows. Section 2 considers the distributed delay version of the time-delayed dynamic model of Matsumoto et al. [5]. Section 3 discusses the dynamics of weak delays in which the growth rate's greatest weighted response is to present population density, whereas previous densities have an exponentially declining effect. Section 4 deals with strong delays in the sense that the largest effect on growth rate response at any time t is related to population density at the prior time t-T. Sections 5 and 6 analyze the cases of a weak delay together with a fixed delay given by the average length of the continuous delay. Finally, Section 7 concludes.

The Model
Matsumoto et al. [5] introduce a discrete-time delay into a dynamic economic model for a small island based on the Easter Island study by BT [1]. The model is an economy with two items and three economic agents (two producers and one consumer). The agricultural good is the harvest of renewable resources, while the manufactured good is something else. The natural resource stock S(t) and population L(t) are supplied. Producers set labor and supply needs to optimize profitability. A manufacturing producer offers the manufactured item created by labor alone. BT presume that per capita use of the resource good increases fertility and/or reduces death in accordance with Malthusian population dynamics. Including a delay τ ≥ 0 in production to account for the reality that agricultural items require time to mature from sowing to harvesting, they obtain the following twodimensional system of delay differential equations: where α is a positive constant indicating the harvesting efficiency; β ∈ (0, 1) indicates the agricultural good's preference, φ is a positive constant. K is the maximum feasible size of the resource stock, while r denotes the natural resource's intrinsic growth rate, and both are positive constants. The difference between the underlying birth and death rates determines population change. The net rate, represented as c, is expected to be negative. When modeling the delay distribution across a large population, using a discrete delay might be considered a crude approximation at times. Of course, if the delay is continuously distributed by a continuous distribution function with a mean delay equal to the discrete delay and a positive variance to account for the delay difference across people, it may seem to be much more realistic. As a result, we suggest a generalizing model (1) with distributed delays. The dynamics of the economy are then governed by the following integro-differential equations system where the delay kernel g is given by a gamma distribution function, i.e., with (j, l) = (1, m) or (j, l) = (2, n), m, n positive integers, and T ≥ 0 is a parameter associated with the mean time delay of the distribution. The weighting function's form is determined by the parameter l. The two cases l = 0 and l = 1 in (3) are the weak delay kernel and the strong delay kernel, respectively. Notice that the distribution function approaches the Dirac distribution as T → 0, while it converges to a fixed delay T as the shape parameters m and n approach infinity. It is clear that the delayed system has the same equilibrium point as the basic system (1). The coordinates of an interior equilibrium are obtained by solving and so they are To analyze the stability of system (2), we should first convert its coordinates to create a new system centered on the equilibrium point (S e , L e ), and then linearize the resulting system at the origin to obtain its characteristic equation. In order to render the eigenvalue analysis analytically tractable, we shall focus on some special cases and apply the linear chain trick technique (see MacDonald [6]), which allows an equation with gamma distributed kernels to be replaced by an equivalent system of ordinary differential equations.

Case m = n = 1
System (2) with weak kernels, i.e., Introducing the new variable then, according to the law of solving the derivative for an integral with parameterized variables, system (4) can be rewritten as the following three-dimensional system of ODEs: The equilibrium (S e , L e ) of system (2) is transformed into the equilibrium (S e , x e , L e ) of (5), where x e = S e . Thus, the stability study of equilibrium (S e , L e ) of (2) is equivalent to the stability study of equilibrium (S e , x e , L e ) of (5). Linearizing the system at the equilibrium point, one finds that the associated characteristic equation is where and To study the stability of the equilibrium, we need to investigate the distribution of roots in the complex plane of the characteristic Equation (6). The stability will be determined by the real parts of the roots of Equation (6). If all roots of Equation (6) are located in the left-half complex plane, then the equilibrium is locally asymptotically stable. If Equation (6) has a root with positive real part, the equilibrium is unstable. According to the Routh-Hurwitz stability criterion, the necessary and sufficient conditions of asymptotic stability are a 1 > 0, a 3 > 0 and a 1 a 2 > a 3 . From (7) and (8), these conditions reduce to a 1 a 2 > a 3 . Hence, we must have Notice that the inequality (9) is verified for all T if M ≤ 0 and for T < 1/M if M > 0. Using the Hopf bifurcation theorem, we check the possibility of the emergence of a limit cycle at T = T * , where At this value of T, one has a * 1 a * 2 = a * 3 , where a * j = a j (T * ) (j = 1, 2, 3). As a result, Equation (6) becomes (λ + a * 1 )(λ 2 + a * 2 ) = 0, yielding the existence of a pair of purely imaginary roots λ 1,2 = ±iω * , with ω * = a * 2 , and a real root λ 3 = −a * 1 . Differentiating Equation (6) with respect to T, we obtain where Recalling that ω 2 * = a * 2 , from (11) we derive Then, Consequently, one pair of complex roots of (6) crosses through the imaginary axis transversally at T = T * from the left half-plane to the right half. Finally, notice from (11) that λ = iω * is a simple root of (6). Otherwise, (11) would imply a * 1 (iω * ) 2 + a * 2 (iω * ) + a * 3 = 0, and so the absurd ω * = 0. By virtue of the previous analysis, we have the following conclusion: Theorem 1. Let T * be defined as in (10).
If λ = iω * is a repeated root of (14), then we get from (15) that a 1 λ 3 + a 2 λ 2 + a 3 λ + a 4 evaluated at T = T * must be zero. Thus, separating the real and imaginary parts, we find a * 1 ω 2 * = a * 3 , a * 2 ω 2 * = a * 4 , where a * j = a j (T * ) (j = 1, 2, 3, 4). However, λ = iω * is a root of (14), and so a * 1 ω 2 * = a * 3 . In other words, one obtains a * 3 a * 1 = a * 1 a * 3 , and so the leading to where . A negative sign of ψ (T * ), and so a positive sign of (16), implies that one pair of complex roots of (14) crosses through the imaginary axis transversally at T = T * from the left half plane to the right half plane. On the other hand, a positive sign of ψ (T * ) yields that the roots can cross the imaginary axis only from right to left as T increases. Summarizing the above analysis, we have the following results.

Theorem 2. Let T *
1 and T * 2 be defined as in the previous Lemma. (1) Let αβφK + 2c = 0 and 16c + r < 0. The equilibrium point (S e , x e , y e , L e ) of (13) is locally asymptotically stable for T < T * 1 , unstable for T > T * 1 , and it bifurcates to chaos at T = T * 1 .

Case m = 1 and n → ∞ System (2) becomes
Letting the dynamical system (17) assumes the form , The characteristic equation of the linearized system around the interior equilibrium point (S e , x e , L e ), where x e = S e is i.e., where We now turn our attention to the issue of stability switching. The stability analysis concerns with whether all roots are located in the left half of the complex plane. When the delay changes, the study of stability switching is concerned with whether the roots cross the imaginary axis. Obviously, Since e −iωT = cos ωT − i sin ωT, separating the real and imaginary parts, we become Squaring and adding these two equations yields where Let z = ω 2 , then Equation (22) rewrites as follows: Since f (0) = s < 0 and f (+∞) = +∞, we conclude that Equation (23) has at least one positive real root.
Proof. Since s < 0, in virtue of Descartes rule, the three enunciated cases follow immediately.
To determine the values of T at which Equation (19) has purely imaginary roots λ = ±iω 0 , we derive from (20) and (21) sin Therefore, from (24) we obtain where j = 0, 1, 2, . . . According to the Hopf bifurcation theorem, to verify stability switching, we need to determine the sign of the derivative of Re(λ) at T = T

Proof. Substituting λ = λ(T) into (19) and by differentiating with respect to T yields
For convenience, we check the sign of (dλ/dT) −1 , that is written as Plugging λ = iω 0 and T = T It remains to prove the simplicity of the root λ = iω 0 . Assuming that this root is a repeated root of (19), then (26)  As T grows, all roots cross the imaginary axis from left to right if the sign is positive, indicating that stability is lost or instability is maintained. On the other hand, the negative sign indicates that the axis is being crossed in the opposite direction, implying that stability may be achieved. Bearing all the above discussion in mind, the following theorem holds.
Theorem 3. Let f (z) and T (j) 0 (j = 0, 1, 2, . . .) be defined as in (23) and (25), respectively. If f (z) = 0 has one or two solutions, there is only one stability switch. If f (z) = 0 has three solutions, then at least a stability switch exists. System (18) may undergo a Hopf bifurcation at (S e , x e , L e ) for those values of T at which a stability switch occurs.
Proof. Suppose that Equation (23) has a unique solution and let ω 0 be the corresponding unique positive root of Equation (22). Recalling f (0) < 0 and f (+∞) = +∞, the polynomial f (z) is an increasing function in a neighborhood of ω 0 , and so its the derivative at ω 0 is positive. Consequently, λ = ±iω 0 crosses the imaginary axis from left to right. Hence, a stability switch occurs. Suppose now that Equation (23) has two solutions, with ω 1 < ω 2 the resulting positive roots of (22). Then, the polynomial f (z) is a decreasing function in a neighborhood of ω 1 and an increasing function in a neighborhood of ω 2 . Thus, one has a crossing of imaginary axis from left to right in correspondence of ω 2 and a crossing from right to left in correspondence of ω 1 . As a result, there is just one stability switch. Finally, suppose that Equation (23) has three solutions and ω 1 < ω 2 < ω 3 are the three positive roots of (22). In this case, the polynomial f (z) is an increasing function in a neighborhood of ω 1 , ω 3 , and a decreasing function in a neighborhood of ω 2 . Crossing of imaginary axis is from left to right in correspondence of ω 1 and ω 3 , from right to left in correspondence of ω 2 . There is at least one stability switch present.

Case m → ∞ and n = 1
System (2) takes the form Setting system (28) rewrites as To examine local dynamics of the above system in a neighborhood of the equilibrium point (S e , x e , L e ), where x e = S e , we consider the linearized system and derive the following characteristic equation i.e., where Let λ = iω with ω > 0 be a solution of (30). Substituting it into (30), and separating the real and imaginary parts, we have The sum of the squares of these two equations yields the sextic equation in ω where Letting z = ω 2 , Equation (33) reduces to a cubic equation Lemma 3. Let ∆ = p 2 q 2 + 18pqs − 4q 3 − 4p 3 s − 27s 2 denote the discriminant of the cubic equation. The following scenarios may be identified.
Proof. The statement follows from Descartes rule, noticing that s > 0.

The conclusion holds.
Then, we have the following results.
Theorem 4. Let g(z) and T (j) 0 (j = 0, 1, 2, ...) be defined as in (34) and (35), respectively. If g(z) = 0 has no solutions, then no stability switches exist. If g(z) = 0 has one or two solutions, then there is at least a stability switch. For those values of T at which a stability transition occurs, system (29) may exhibit a Hopf bifurcation.
Proof. If Equation (34) has no solutions, the proof is trivial. The second part of the statement follows the same justifications for the theorem presented in the previous section.

Conclusions
This paper generalizes the BT's delayed continuous-time model for small islands presented by Matsumoto et al. [5]. Time lags modeled by way of distributed delays allow the reduction of the dynamics to a set of ordinary differential equations. The gamma distribution function has been considered to have only weak and strong kernels to simplify our analysis. Sufficient conditions of existence and stability of equilibria in four different cases are derived. In contrast with the delayed model, it is found that the positive equilibrium becomes unstable for all large delay values, or the stability of equilibrium switches back, leading to the occurrence of multiple stability switches.