*Article* **Global Stability Analysis of a Five-Dimensional Unemployment Model with Distributed Delay**

**Eva Kaslik 1,2,\*,†, Mihaela Neam¸tu 3,\*,† and Loredana Flavia Vesa 1,2,\*,†**


**Abstract:** The present paper proposes a five-dimensional mathematical model for studying the labor market, focusing on unemployment, migration, fixed term contractors, full time employment and the number of available vacancies. The distributed time delay is considered in the rate of change of available vacancies that depends on the past regular employment levels. The non-dimensional mathematical model is introduced and the existence of the equilibrium points is analyzed. The positivity and boundedness of solutions are provided and global asymptotic stability findings are presented both for the employment free equilibrium and the positive equilibrium. The numerical simulations support the theoretical results.

**Keywords:** unemployment model; global stability; Lyapunov function; distributed delay

**1. Introduction**

Over the time, one of the many challenges that a country faces is economical downturn fuelled by unemployment. The causes of appearance of this phenomena are countless and are based on country specificity. Problems ranging from the increasing population to slow economic growth are some of the leading factors linked to unemployment spanning out of control. One of the negative impact of unemployment pertains to the social side that puts the affected population at high risk, such as psychological and mental health problems, just to name a few [1]. Any responsible government cannot ignore the signals coming from the labour market and should take appropriate and immediate actions to improve the overall situation.

Therefore, the need to manage and handle the tendency of unwanted spread of unemployment brings about the necessity of studying the behaviour of complex mathematical models in conjunction with numerical simulations. In order to comprehend the dynamics of unemployment, using some concepts from [2], the following variables have been considered in a previous mathematical model investigated in [3,4]: number of unemployed persons, employed persons and new vacancies; time delay has been incorporated into the rate of change for creation of new vacancies. Moreover, by resorting to the skill development programs, Misra et al. [5] demonstrated the link between the betterment of workers' capabilities and the reduction of unemployment. The effect of training programs has also been studied in [6]. Furthermore, Harding and Neam¸tu [7] considered the migration as a contributing factor when defining policies pertaining to unemployment, including distributed time delay. The optimal control analysis has been completed in [8,9].

The motivation of the present paper is centered around the existing mathematical models, enabling the development of new ways for studying unemployment, based on

**Citation:** Kaslik, E.; Neam¸tu M.; Vesa, L.F. Global Stability Analysis of a Five-Dimensional Unemployment Model with Distributed Delay. *Mathematics* **2021**, *9*, 3037. https:// doi.org/10.3390/math9233037

Academic Editor: Snezhana Hristova

Received: 11 November 2021 Accepted: 23 November 2021 Published: 26 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the past history of the state variables. The potential impact of our findings is related to upcoming biological and economical mathematical models connected to population dynamics. The most important benefit could be related to the strategic policy makers of our society.

The novelty of this paper is the investigation of the interaction among the number of unemployed persons, immigrants, temporary employed persons, regularly employed persons and the number of available vacancies, in the set up of delay differential equations. The distributed time delay has been incorporated to reflect the dependence of the rate of change of available vacancies on the past regular employment levels. Basically, this type of delay is highlighted for acquiring the realistic sides of the economic process [10–16]. It is important to emphasize that there are many mathematical models including distributed time delays, originating from population biology, epidemiology and economics [17–33].

The paper is structured as follows: the mathematical model and its non-dimensional version are introduced in Section 2; the existence of the equilibrium points of the model is discussed in Section 3 and the positivity and boundedness of solutions are proved in Section 4; global asymptotic stability results are proved for the employment free equilibrium and the positive equilibrium in Sections 5 and 6, respectively; Section 7 includes numerical simulations to highlight the theoretical findings, followed by conclusions which are formulated in Section 8.

#### **2. The Mathematical Model and Its Non-Dimensional Version**

We consider the following five variables which describe the mathematical model for to control the unemployment: the number of unemployed persons *U*(*t*), the number of immigrants *M*(*t*), the number of temporary employed persons *T*(*t*), the number of regularly employed persons *R*(*t*), respectively and the number of available vacancies *V*(*t*) at time *t*.

The system of differential equations is:

$$\begin{cases} \dot{L}(t) = a\_1 - a\_2 \mathcal{U}(t)V(t) + a\_3 \mathcal{R}(t) - a\_4 \mathcal{U}(t) + a\_5 \mathcal{T}(t) - b\_1 \mathcal{U}(t) \\ \dot{M}(t) = m\_1 - m\_2 \mathcal{M}(t)V(t) - b\_2 \mathcal{M}(t) \\ \dot{T}(t) = a\_4 \mathcal{U}(t) - a\_5 \mathcal{T}(t) - c\_1 \mathcal{T}(t)V(t) - b\_3 \mathcal{T}(t) \\ \mathcal{R}(t) = a\_2 \mathcal{U}(t)V(t) + m\_2 \mathcal{M}(t)V(t) + c\_1 \mathcal{T}(t)V(t) - a\_3 \mathcal{R}(t) - b\_4 \mathcal{R}(t) \\ \dot{V}(t) = b\_4 \int\_0^\infty h(s)\mathcal{R}(t-s)ds - b\_5 V(t) \end{cases} \tag{1}$$

where *a*1, *a*2, *a*3, *a*<sup>4</sup> and *a*<sup>5</sup> are positive constants: *a*1the constant growth rate of unemployed persons entering the labor market, *a*<sup>2</sup> the rate of hiring, *a*<sup>3</sup> the rate of firing, *a*<sup>4</sup> the rate of move to unemployment, *a*<sup>5</sup> the rate of move to temporary employment; *b*<sup>1</sup> the rate of migration of unemployed persons; *b*<sup>2</sup> is the rate of return or death, *b*<sup>3</sup> the rate of temporary employment, *b*<sup>4</sup> the rate of retirement, migration or death of employed persons; *b*<sup>5</sup> the rate of available vacancies; *c*<sup>1</sup> is the rate of good job finding; *m*<sup>1</sup> represent the exogenous increase in migration, *m*<sup>2</sup> is the migrant employment rate.

The bounded, piecewise continuous function *h* : [0, ∞) → [0, ∞) is the delay kernel with average time delay *τ* for available vacancies based on past employment levels. Hence, *h* is a probability density function satisfying the following properties:

$$\int\_0^\infty h(s)ds = 1, \quad \tau = \int\_0^\infty sh(s)ds < \infty. \tag{2}$$

With the changes of variables:

$$\begin{cases} \begin{aligned} \mathbf{x}\_{1}(t) &= \frac{a\_{2}b\_{4}}{a\_{5}^{2}} U \left(\frac{t}{a\_{5}}\right) \\ \mathbf{x}\_{2}(t) &= \frac{a\_{2}b\_{4}}{a\_{5}^{2}} M \left(\frac{t}{a\_{5}}\right) \\ \mathbf{x}\_{3}(t) &= \frac{a\_{2}b\_{4}}{a\_{5}^{2}} T \left(\frac{t}{a\_{5}}\right) \\ \mathbf{x}\_{4}(t) &= \frac{a\_{2}b\_{4}}{a\_{5}^{2}} R \left(\frac{t}{a\_{5}}\right) \\ \mathbf{x}\_{5}(t) &= \frac{a\_{2}}{a\_{5}} V \left(\frac{t}{a\_{5}}\right) \end{aligned} \tag{3}$$

system (1) becomes a non-dimensional system, with fewer dimensionless parameters and dimensionless state variables:

$$\begin{cases}
\dot{\mathbf{x}}\_{1}(t) = \gamma\_{1} - \mathbf{x}\_{1}(t)\mathbf{x}\_{5}(t) + a\_{3}\mathbf{x}\_{4}(t) - a\_{4}\mathbf{x}\_{1}(t) + \mathbf{x}\_{3}(t) - \beta\_{1}\mathbf{x}\_{1}(t), \\
\dot{\mathbf{x}}\_{2}(t) = \gamma\_{2} - a\_{2}\mathbf{x}\_{2}(t)\mathbf{x}\_{5}(t) - \beta\_{2}\mathbf{x}\_{2}(t), \\
\dot{\mathbf{x}}\_{3}(t) = a\_{4}\mathbf{x}\_{1}(t) - \mathbf{x}\_{3}(t) - a\_{1}\mathbf{x}\_{3}(t)\mathbf{x}\_{5}(t) - \beta\_{3}\mathbf{x}\_{3}(t), \\
\dot{\mathbf{x}}\_{4}(t) = \mathbf{x}\_{1}(t)\mathbf{x}\_{5}(t) + a\_{2}\mathbf{x}\_{2}(t)\mathbf{x}\_{5}(t) + a\_{1}\mathbf{x}\_{3}(t)\mathbf{x}\_{5}(t) - a\_{3}\mathbf{x}\_{4}(t) - \beta\_{4}\mathbf{x}\_{4}(t), \\
\dot{\mathbf{x}}\_{5}(t) = \int\_{0}^{\infty} k(\mathbf{s})\mathbf{x}\_{4}(t-s)ds - \beta\_{5}\mathbf{x}\_{5}(t),
\end{cases} (4)$$

where the delay kernel is *k*(*s*) = <sup>1</sup> *a*5 *h <sup>s</sup> a*5 and the coefficients are expressed as

$$\begin{aligned} \gamma\_1 &= \frac{a\_1 a\_2 b\_4}{a\_5^3}, & \gamma\_2 &= \frac{a\_2 b\_4 m\_1}{a\_5^3} \\ \alpha\_1 &= \frac{c\_1}{a\_2}, & a\_2 &= \frac{m\_2}{a\_2}, & a\_3 &= \frac{a\_3}{a\_5}, & a\_4 &= \frac{a\_4}{a\_5} \\ \beta\_1 &= \frac{b\_1}{a\_5}, & \beta\_2 &= \frac{b\_2}{a\_5}, & \beta\_3 &= \frac{b\_3}{a\_5}, & \beta\_4 &= \frac{b\_4}{a\_5}, & \beta\_5 &= \frac{b\_5}{a\_5} \end{aligned}$$

Initial conditions for system (4) are considered of the form

$$\mathfrak{x}\_{i}(\theta) = \mathfrak{q}\_{i}(\theta) \quad , \,\forall \,\theta \in (-\infty, 0], \,\forall \, i = 1, 5, \,\forall i$$

where *<sup>ϕ</sup><sup>i</sup>* belong to the Banach space *<sup>C</sup>*0,*μ*(R−, R) (where *<sup>μ</sup>* > 0) of continuous real valued functions defined on (−∞, 0] such that lim *<sup>t</sup>*→−<sup>∞</sup> *<sup>e</sup>μtϕ*(*t*) = 0, considered with respect to the norm:

$$||q||\_{\infty,\mu} = \sup\_{t \in (-\infty,0]} \varepsilon^{\mu t} |q(t)|.$$

The existence and uniqueness of solutions, as well as the continuous dependence of solutions on initial conditions, in the framework of the distributed delay system (4), are guaranteed by theoretical results presented in [34].

#### **3. Equilibrium Points of the Model**

The equilibrium points are constant solutions of system (4) and hence they satisfy the following algebraic equations:

$$\begin{cases} \gamma\_1 - \mathbf{x}\_1 \mathbf{x}\_5 + \mathbf{a}\_3 \mathbf{x}\_4 - \mathbf{a}\_4 \mathbf{x}\_1 + \mathbf{x}\_3 - \beta\_1 \mathbf{x}\_1 = \mathbf{0} \\ \gamma\_2 - \mathbf{a}\_2 \mathbf{x}\_2 \mathbf{x}\_5 - \beta\_2 \mathbf{x}\_2 = \mathbf{0} \\ \mathbf{a}\_4 \mathbf{x}\_1 - \mathbf{x}\_3 - \mathbf{a}\_1 \mathbf{x}\_3 \mathbf{x}\_5 - \beta\_3 \mathbf{x}\_3 = \mathbf{0} \\ \mathbf{x}\_1 \mathbf{x}\_5 + \mathbf{a}\_2 \mathbf{x}\_2 \mathbf{x}\_5 + \mathbf{a}\_1 \mathbf{x}\_3 \mathbf{x}\_5 - \mathbf{a}\_3 \mathbf{x}\_4 - \beta\_4 \mathbf{x}\_4 = \mathbf{0} \\ \mathbf{x}\_4 - \beta\_5 \mathbf{x}\_5 = \mathbf{0} \end{cases} \tag{5}$$

On one hand, we observe that system (5) has at least one solution, which corresponds to the case when *x*<sup>5</sup> = 0. Hence, we obtain the equilibrium point *S*<sup>0</sup> given by:

$$S^0 := (\delta\_1 (1 + \beta\_3), \delta\_2, \delta\_1 \alpha\_{4\prime} 0, 0)$$

where

$$
\delta\_1 = \frac{\gamma\_1}{\beta\_1 + \beta\_1 \beta\_3 + \alpha\_4 \beta\_3} \quad \text{and} \quad \delta\_2 = \frac{\gamma\_2}{\beta\_2}.
$$

This equilibrium point corresponds to the state of no regular employment and no available vacancies.

We now introduce the basic reproduction number *R*0, which has the role of a threshold parameter that prognosticates whether the unemployment, immigration and temporary employment problems will increase or decrease. Using the next generation matrix method, we deduce:

$$R\_0 = \frac{\delta\_1 (1 + \beta\_3 + \alpha\_1 \alpha\_4) + \delta\_2 \alpha\_2}{(\alpha\_3 + \beta\_4)\beta\_5}$$

On the other hand, system (5) has at least one solution with the last component *x*<sup>5</sup> > 0, if and only if *x*<sup>5</sup> is the solution of the following cubic equation:

$$E\_3 \mathbf{x}\_5^3 + E\_2 \mathbf{x}\_5^2 + E\_1 \mathbf{x}\_5 + E\_0 = \mathbf{0},\tag{6}$$

where

$$\begin{aligned} E\_3 &= \alpha\_1 \beta\_4 \beta\_5\\ E\_2 &= \alpha\_1 (\nu\_2 \beta\_4 \beta\_5 - \gamma\_1 - \gamma\_2) + \beta\_4 \beta\_5 (1 + \beta\_3 + a\_1 a\_4) + a\_1 \beta\_1 \mu\\ E\_1 &= (\nu\_2 \beta\_4 \beta\_5 - \gamma\_1 - \gamma\_2)(1 + \beta\_3 + a\_1 a\_4) - \nu\_2 \gamma\_1 a\_1 + \mu \beta\_1 \left(1 + \beta\_3 + a\_1 a\_4 \frac{\nu\_1}{\beta\_1}\right) + (\mu \nu\_2 - \gamma\_2) \beta\_1 a\_1\\ E\_0 &= \frac{\mu \gamma\_1 \gamma\_2}{a\_2 \delta\_1 \delta\_2} (1 - R\_0),\end{aligned}$$

where *R*0, *δ*<sup>1</sup> and *δ*<sup>2</sup> are given above and

$$\nu\_1 = \frac{\beta\_3}{\alpha\_1}, \qquad \nu\_2 = \frac{\beta\_2}{\alpha\_2}, \qquad \mu = (\alpha\_3 + \beta\_4)\beta\_5.$$

In this case, system (4) has at least one equilibrium point *S*<sup>+</sup> of the form

$$S^{+} := \left( -\frac{Q(\mathbf{x\_{5}}) + \nu\_{1}d(\mathbf{x\_{5}})}{\left(\beta\_{1} - \nu\_{1}\right)(\mathbf{x\_{5}} + \nu\_{2})}, \frac{\gamma\_{2}}{a\_{2}(\nu\_{2} + \mathbf{x\_{5}})}, \frac{Q(\mathbf{x\_{5}}) + \beta\_{1}d(\mathbf{x\_{5}})}{a\_{1}(\beta\_{1} - \nu\_{1})(\nu\_{2} + \mathbf{x\_{5}})}, \beta\_{5}\mathbf{x\_{5}}, \mathbf{x\_{5}} \right).$$

where

$$\begin{aligned} Q(\mathbf{x}\mathbf{s}) &= (\beta\_4 \beta\_5 \mathbf{x}\mathbf{s} - \gamma\_1 - \gamma\_2)(\mathbf{x}\mathbf{s} + \nu\_2) + \gamma\_2 \nu\_{2,0}, \\ d(\mathbf{x}\_5) &= \mu(\mathbf{x}\_5 + \nu\_2) - \gamma\_2. \end{aligned}$$

In fact, we distinguish two cases:

**Case 1:** If *R*<sup>0</sup> > 1 then *E*<sup>0</sup> < 0 and system (4) has at least one positive equilibrum point. Moreover, if either *E*<sup>1</sup> < 0, or *E*<sup>1</sup> > 0 and *E*<sup>2</sup> > 0, Descartes' rule of signs guarantees the existence of a unique positive equilibrium point *S*+.

**Case 2:** If *R*<sup>0</sup> < 1 then *E*<sup>0</sup> > 0, we may have either two or zero positive equilibrium points. For instance, if *E*<sup>1</sup> > 0 and *E*<sup>2</sup> > 0, system (4) does not have positive equilibrium points.

#### **4. Positivity and Boundedness of Solutions**

**Theorem 1.** *The open positive (nonnegative) orthant* R<sup>5</sup> <sup>+</sup> = (0, ∞)<sup>5</sup> *is invariant to the flow of system* (4)*. Furthermore, denoting β<sup>m</sup>* = min(*β*1, *β*2, *β*3, *β*4)*, the set*

$$\Omega = \left\{ (\mathbf{x}\_1, \mathbf{x}\_2, \mathbf{x}\_3, \mathbf{x}\_4, \mathbf{x}\_5) : 0 \le \mathbf{x}\_1 + \mathbf{x}\_2 + \mathbf{x}\_3 + \mathbf{x}\_4 \le \frac{\gamma\_1 + \gamma\_2}{\beta\_m}, \ 0 \le \mathbf{x}\_5 \le \frac{\gamma\_1 + \gamma\_2}{\beta\_m \beta\_5} \right\}$$

*is a region of attraction for the system* (1)*, attracting all the solutions with initial conditions in the open positive orthant* R<sup>5</sup> +*.*

**Proof.** As a first step, we will show that the open positive (nonnegative) orthant of R<sup>5</sup> is invariant to the flow of system (4). Considering initial conditions *ϕ<sup>i</sup>* : (−∞, 0] → (0, ∞), for *i* = 1, 5, based on the continuity of the solutions, there exists *T* > 0 such that *xi*(*t*) > 0, for any *t* ∈ (0, *T*) and any *i* = 1, 5.

Therefore, the fifth equation of system (4) provides

$$
\dot{\mathfrak{x}}\_5(t) \ge -\beta\_5 \mathfrak{x}\_5(t), \quad \forall \ t \in (0, T)
$$

and an integration over the interval (0, *T*) leads to:

$$\varkappa\_5(T) \ge \varrho\_5(0)e^{-b\_5T} > 0.1$$

In a similar manner, it can be show that *xi*(*T*) > 0, for any *i* = 1, 4. Moreover, from the first four equations of system (4) we observe that

$$\frac{d}{dt}[\mathbf{x}\_1(t) + \mathbf{x}\_2(t) + \mathbf{x}\_3(t) + \mathbf{x}\_4(t)] = \gamma\_1 + \gamma\_2 - \beta\_1 \mathbf{x}\_1(t) - \beta\_2 \mathbf{x}\_2(t) - \beta\_3 \mathbf{x}\_3(t) - \beta\_4 \mathbf{x}\_4(t),$$

and taking into account that *β<sup>m</sup>* = min(*β*1, *β*2, *β*3, *β*4) we have

$$\frac{d}{dt}[\mathbf{x}\_1(t) + \mathbf{x}\_2(t) + \mathbf{x}\_3(t) + \mathbf{x}\_4(t)] \le \gamma\_1 + \gamma\_2 - \beta\_m[\mathbf{x}\_1(t) + \mathbf{x}\_2(t) + \mathbf{x}\_3(t) + \mathbf{x}\_4(t)], \quad \forall \ t > 0.$$

Therefore, using basic differential inequality techniques, we obtain:

$$\mathbf{x}\_1(t) + \mathbf{x}\_2(t) + \mathbf{x}\_3(t) + \mathbf{x}\_4(t) \le e^{-\beta\_m t} \left( \mathbf{x}\_1(0) + \mathbf{x}\_2(0) + \mathbf{x}\_3(0) + \mathbf{x}\_4(0) - \frac{\gamma\_1 + \gamma\_2}{\beta\_m} \right) + \frac{\gamma\_1 + \gamma\_2}{\beta\_m}$$

for any *t* ≥ 0.

On one hand, if *<sup>x</sup>*1(0) + *<sup>x</sup>*2(0) + *<sup>x</sup>*3(0) + *<sup>x</sup>*4(0) <sup>≤</sup> *<sup>γ</sup>*1+*γ*<sup>2</sup> *<sup>β</sup><sup>m</sup>* it follows that (*x*1(0) + *x*2(0) + *<sup>x</sup>*3(0) + *<sup>x</sup>*4(0) <sup>≤</sup> *<sup>γ</sup>*1+*γ*<sup>2</sup> *<sup>β</sup><sup>m</sup>* , for any *<sup>t</sup>* <sup>&</sup>gt; 0. Otherwise, if *<sup>x</sup>*1(0) + *<sup>x</sup>*2(0) + *<sup>x</sup>*3(0) + *<sup>x</sup>*4(0) <sup>&</sup>gt; *<sup>γ</sup>*1+*γ*<sup>2</sup> *<sup>β</sup><sup>m</sup>* , we have

,

$$\limsup\_{t \to \infty} \left[ \mathbf{x}\_1(t) + \mathbf{x}\_2(t) + \mathbf{x}\_3(t) + \mathbf{x}\_4(t) \right] \le \frac{\gamma\_1 + \gamma\_2}{\beta\_m}.$$

Solving the last equation of (4) for *x*5(*t*) leads to

$$\mathbf{x}\_{\mathfrak{F}}(t) = e^{-\beta\_{\mathfrak{F}}t} \mathbf{x}(0) + e^{-\beta\_{\mathfrak{F}}t} \int\_{0}^{t} e^{\beta\_{\mathfrak{F}}u} \left( \int\_{0}^{\infty} k(s) \mathbf{x}\_{\mathfrak{k}}(u-s) ds \right) du.$$

Using the generalized l'Hospital rule (see Lemma 1.1 in [35]) in conjunction with Lemma 1 from [4], we deduce:

$$\begin{split} \limsup\_{t \to \infty} \mathbf{x}\_5(t) &\leq \frac{1}{\beta\_5} \limsup\_{t \to \infty} \int\_0^\infty k(s) \mathbf{x}\_4(t-s) ds \\ &= \frac{1}{\beta\_5} \limsup\_{t \to \infty} \left( \int\_0^t k(s) \mathbf{x}\_4(t-s) ds + \int\_t^\infty k(s) \boldsymbol{q}\_4(t-s) ds \right) \\ &\leq \frac{1}{\beta\_5} \left( \limsup\_{t \to \infty} \mathbf{x}\_4(t) + \| \boldsymbol{q}\_4 \|\_{\infty, \mu} \limsup\_{t \to \infty} \int\_t^\infty k(s) e^{-\mu(t-s)} ds \right) \\ &= \frac{1}{\beta\_5} \left( \limsup\_{t \to \infty} \mathbf{x}\_4(t) + \| \boldsymbol{q}\_4 \|\_{\infty, \mu} \limsup\_{t \to \infty} e^{-\mu t} \int\_t^\infty k(s) e^{\mu s} ds \right) \\ &\leq \frac{\gamma\_1 + \gamma\_2}{\beta\_5 \beta\_{\text{m}}}, \end{split}$$

which concludes the proof.

#### **5. Global Stability Analysis for** *S***<sup>0</sup>**

**Theorem 2.** *The equilibrium point S*<sup>0</sup> *is globally asymptotically stable in the open positive orthant* R<sup>5</sup> <sup>+</sup>*, regardless of the delay kernel k*(*s*)*, if the following inequality holds:*

$$R\_0 < \frac{\beta\_4}{\alpha\_3 + \beta\_4}.\tag{7}$$

**Proof.** Let us consider an arbitrary solution *xi*(*t*), *i* = 1, 5, of system (4), with initial conditions *ϕ<sup>i</sup>* : (−∞, 0] → (0, ∞). From Theorem 1 it follows that *xi*(*t*) > 0, for any *t* > 0 and *i* = 1, 5. Hence, we consider the functions

$$L\_k(t) = x\_k(t) - x\_k^0 - x\_k^0 \ln \frac{x\_k(t)}{x\_k^0} \text{ , } k = \overline{1,3} \text{ ,}$$

and we further denote:

$$E\_k(t) = 1 - \frac{\mathbf{x}\_k(t)}{\mathbf{x}\_k^0} + \ln \frac{\mathbf{x}\_k(t)}{\mathbf{x}\_k^0} < 0 \quad , \quad E\_k(t) = 1 - \frac{\mathbf{x}\_k^0}{\mathbf{x}\_k(t)} + \ln \frac{\mathbf{x}\_k^0}{\mathbf{x}\_k(t)} < 0 \; , \; k = \overline{1,3}.$$

Therefore, we have:

*L* <sup>1</sup>(*t*) = *x*˙1 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 1 *x*1 = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 1 *x*1 [*γ*<sup>1</sup> − *x*1*x*<sup>5</sup> + *α*3*x*<sup>4</sup> + *x*<sup>3</sup> − (*α*<sup>4</sup> + *β*1)*x*1] = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 1 *x*1 [*x*0 1*x*0 <sup>5</sup> <sup>−</sup> *<sup>α</sup>*3*x*<sup>0</sup> <sup>4</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> <sup>3</sup> + (*α*<sup>4</sup> + *<sup>β</sup>*1)*x*<sup>0</sup> <sup>1</sup> − *x*1*x*<sup>5</sup> + *α*3*x*<sup>4</sup> + *x*<sup>3</sup> − (*α*<sup>4</sup> + *β*1)*x*1] = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 1 *x*1 (*α*<sup>4</sup> <sup>+</sup> *<sup>β</sup>*1)*x*<sup>0</sup> 1 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>1</sup> *x*0 1 <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 3 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*0 3 − *x*1*x*<sup>5</sup> + *α*3*x*<sup>4</sup> = (*α*<sup>4</sup> + *β*1)*x*<sup>0</sup> 1 <sup>2</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 1 *x*1 <sup>−</sup> *<sup>x</sup>*<sup>1</sup> *x*0 1 <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 3 <sup>1</sup> <sup>+</sup> *<sup>x</sup>*<sup>0</sup> 1 *x*1 · *x*3 *x*0 3 <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*0 3 − *x*0 1 *x*1 − *x*1*x*<sup>5</sup> + *α*3*x*<sup>4</sup> − *α*3*x*<sup>4</sup> *x*0 1 *x*1 + *x*5*x*<sup>0</sup> 1 <sup>≤</sup> (*α*<sup>4</sup> <sup>+</sup> *<sup>β</sup>*1)*x*<sup>0</sup> 1 <sup>2</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 1 *x*1 <sup>−</sup> *<sup>x</sup>*<sup>1</sup> *x*0 1 <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 3 <sup>2</sup> <sup>+</sup> ln *<sup>x</sup>*<sup>0</sup> 1 *x*1 <sup>+</sup> ln *<sup>x</sup>*<sup>3</sup> *x*0 3 <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*0 3 − *x*0 1 *x*1 <sup>−</sup> *<sup>x</sup>*1*x*<sup>5</sup> <sup>+</sup> *<sup>α</sup>*3*x*<sup>4</sup> <sup>+</sup> *<sup>x</sup>*5*x*<sup>0</sup> 1 = (*α*<sup>4</sup> + *β*1)*x*<sup>0</sup> <sup>1</sup>(*E*<sup>1</sup> + *<sup>E</sup>*˜ <sup>1</sup>) <sup>−</sup> *<sup>x</sup>*<sup>0</sup> <sup>3</sup>(*E*˜ <sup>1</sup> <sup>+</sup> *<sup>E</sup>*3) <sup>−</sup> *<sup>x</sup>*1*x*<sup>5</sup> <sup>+</sup> *<sup>x</sup>*<sup>0</sup> <sup>1</sup>*x*<sup>5</sup> + *α*3*x*<sup>4</sup> − *α*3*x*<sup>4</sup> *x*0 1 *x*1 .

Furthermore, we evaluate

$$\begin{split} L\_{2}^{\prime}(t) &= \dot{\mathbf{x}}\_{2} \left( 1 - \frac{\mathbf{x}\_{2}^{0}}{\mathbf{x}\_{2}} \right) = \left( 1 - \frac{\mathbf{x}\_{2}^{0}}{\mathbf{x}\_{2}} \right) (\gamma\_{2} - \alpha\_{2}\mathbf{x}\_{2}\mathbf{x}\_{5} - \beta\_{2}\mathbf{x}\_{2}) \\ &= \left( 1 - \frac{\mathbf{x}\_{2}^{0}}{\mathbf{x}\_{2}} \right) (\alpha\_{2}\mathbf{x}\_{2}^{0}\mathbf{x}\_{5}^{0} + \beta\_{2}\mathbf{x}\_{2}^{0} - \alpha\_{2}\mathbf{x}\_{2}\mathbf{x}\_{5} - \beta\_{2}\mathbf{x}\_{2}) \\ &= \beta\_{2}\mathbf{x}\_{2}^{0} \left( 2 - \frac{\mathbf{x}\_{2}^{0}}{\mathbf{x}\_{2}} - \frac{\mathbf{x}\_{2}}{\mathbf{x}\_{2}^{0}} \right) - \alpha\_{2}\mathbf{x}\_{2}\mathbf{x}\_{5} + \alpha\_{2}\mathbf{x}\_{2}^{0}\mathbf{x}\_{5} \\ &= \beta\_{2}\mathbf{x}\_{2}^{0} (E\_{2} + \vec{E}\_{2}) - \alpha\_{2}\mathbf{x}\_{2}\mathbf{x}\_{5} + \alpha\_{2}\mathbf{x}\_{2}^{0}\mathbf{x}\_{5} \end{split}$$

and

*L* <sup>3</sup>(*t*) = *x*˙3 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 3 *x*3 = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 3 *x*3 (*α*4*x*<sup>1</sup> − *x*<sup>3</sup> − *α*1*x*3*x*<sup>5</sup> − *β*3*x*3) = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 3 *x*3 *α*4(*x*<sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> <sup>1</sup>) <sup>−</sup> (<sup>1</sup> <sup>+</sup> *<sup>β</sup>*3)(*x*<sup>3</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> <sup>3</sup>) <sup>−</sup> *<sup>α</sup>*1(*x*3*x*<sup>5</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 3*x*0 5) = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 3 *x*3 *α*4*x*<sup>0</sup> 1( *x*1 *x*0 1 <sup>−</sup> <sup>1</sup>) <sup>−</sup> *<sup>x</sup>*<sup>0</sup> <sup>3</sup>(<sup>1</sup> <sup>+</sup> *<sup>β</sup>*3)( *<sup>x</sup>*<sup>3</sup> *x*0 3 − 1) − *α*1*x*3*x*<sup>5</sup> <sup>=</sup> <sup>−</sup>*α*4*x*<sup>0</sup> 1 <sup>1</sup> <sup>+</sup> *<sup>x</sup>*<sup>0</sup> 3 *x*3 · *x*1 *x*0 1 − *x*0 3 *x*3 <sup>−</sup> *<sup>x</sup>*<sup>1</sup> *x*0 1 + *x*<sup>0</sup> <sup>3</sup>(<sup>1</sup> <sup>+</sup> *<sup>β</sup>*3)(<sup>2</sup> <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*0 3 − *x*0 3 *x*3 ) <sup>−</sup> *<sup>α</sup>*1*x*3*x*<sup>5</sup> <sup>+</sup> *<sup>α</sup>*1*x*<sup>0</sup> <sup>3</sup>*x*<sup>5</sup> ≤ −*α*4*x*<sup>0</sup> 1 <sup>2</sup> <sup>+</sup> ln *<sup>x</sup>*<sup>0</sup> 3 *x*3 <sup>+</sup> ln *<sup>x</sup>*<sup>1</sup> *x*0 1 − *x*0 3 *x*3 <sup>−</sup> *<sup>x</sup>*<sup>1</sup> *x*0 1 + *x*<sup>0</sup> <sup>3</sup>(<sup>1</sup> <sup>+</sup> *<sup>β</sup>*3)(<sup>2</sup> <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*0 3 − *x*0 3 *x*3 ) <sup>−</sup> *<sup>α</sup>*1*x*3*x*<sup>5</sup> <sup>+</sup> *<sup>α</sup>*1*x*<sup>0</sup> <sup>3</sup>*x*<sup>5</sup> <sup>=</sup> <sup>−</sup>*α*4*x*<sup>0</sup> <sup>1</sup>(*E*<sup>1</sup> + *<sup>E</sup>*˜ <sup>3</sup>)+(1 + *β*3)*x*<sup>0</sup> <sup>3</sup>(*E*<sup>3</sup> + *<sup>E</sup>*˜ <sup>3</sup>) <sup>−</sup> *<sup>α</sup>*1*x*3*x*<sup>5</sup> <sup>+</sup> *<sup>α</sup>*1*x*<sup>0</sup> <sup>3</sup>*x*5.

Moreover, considering

$$L\_{4,5}(t) = x\_4 + c x\_5 + c \int\_0^\infty \left( k(s) \int\_{t-s}^t x\_4(r) dr \right) ds$$

we have:

$$\begin{split} L'\_{4,5}(t) &= \dot{\mathbf{x}}\_4 + c\dot{\mathbf{x}}\_5 + c \int\_0^\infty \dot{\mathbf{k}}(\mathbf{s}) [\mathbf{x}\_4(t) - \mathbf{x}\_4(t-s)] ds \\ &= \mathbf{x}\_1 \mathbf{x}\_5 + \mathbf{a}\_2 \mathbf{x}\_2 \mathbf{x}\_5 + \mathbf{a}\_1 \mathbf{x}\_3 \mathbf{x}\_5 - (\mathbf{a}\_3 + \beta\_4 - c) \mathbf{x}\_4 - c\beta\_5 \mathbf{x}\_5. \end{split}$$

Further, considering:

$$L(t) = L\_1(t) + L\_2(t) + L\_3(t) + L\_{4,5}(t)$$

we compute

*L* (*t*) <sup>≤</sup> (*α*<sup>4</sup> <sup>+</sup> *<sup>β</sup>*1)*x*<sup>0</sup> <sup>1</sup>(*E*<sup>1</sup> + *<sup>E</sup>*˜ <sup>1</sup>) <sup>−</sup> *<sup>x</sup>*<sup>0</sup> <sup>3</sup>(*E*˜ <sup>1</sup> <sup>+</sup> *<sup>E</sup>*3) <sup>−</sup> *<sup>x</sup>*1*x*<sup>5</sup> <sup>+</sup> *<sup>x</sup>*<sup>0</sup> <sup>1</sup>*x*<sup>5</sup> + *α*3*x*<sup>4</sup> + *β*2*x*<sup>0</sup> <sup>2</sup>(*E*<sup>2</sup> + *<sup>E</sup>*˜ <sup>2</sup>) <sup>−</sup> *<sup>α</sup>*2*x*2*x*<sup>5</sup> <sup>+</sup> *<sup>α</sup>*2*x*<sup>0</sup> <sup>2</sup>*x*<sup>5</sup> <sup>−</sup> *<sup>α</sup>*4*x*<sup>0</sup> <sup>1</sup>(*E*<sup>1</sup> + *<sup>E</sup>*˜ <sup>3</sup>)+(1 + *β*3)*x*<sup>0</sup> <sup>3</sup>(*E*<sup>3</sup> + *<sup>E</sup>*˜ <sup>3</sup>) <sup>−</sup> *<sup>α</sup>*1*x*3*x*<sup>5</sup> <sup>+</sup> *<sup>α</sup>*1*x*<sup>0</sup> <sup>3</sup>*x*<sup>5</sup> + *x*1*x*<sup>5</sup> + *α*2*x*2*x*<sup>5</sup> + *α*1*x*3*x*<sup>5</sup> − (*α*<sup>3</sup> + *β*<sup>4</sup> − *c*)*x*<sup>4</sup> − *cβ*5*x*<sup>5</sup> = *β*1*x*<sup>0</sup> <sup>1</sup>*E*<sup>1</sup> + (*α*<sup>4</sup> + *β*1)*x*<sup>0</sup> <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> 3 *E*˜ <sup>1</sup> + *β*2*x*<sup>0</sup> <sup>2</sup>(*E*<sup>2</sup> + *<sup>E</sup>*˜ <sup>2</sup>) + *β*3*x*<sup>0</sup> <sup>3</sup>*E*<sup>3</sup> + (1 + *β*3)*x*<sup>0</sup> <sup>3</sup> <sup>−</sup> *<sup>α</sup>*4*x*<sup>0</sup> 1 *E*˜ 3 <sup>−</sup> (*β*<sup>4</sup> <sup>−</sup> *<sup>c</sup>*)*x*<sup>4</sup> <sup>−</sup> (*cβ*<sup>5</sup> <sup>−</sup> *<sup>x</sup>*<sup>0</sup> <sup>1</sup> <sup>−</sup> *<sup>α</sup>*2*x*<sup>0</sup> <sup>2</sup> <sup>−</sup> *<sup>α</sup>*1*x*<sup>0</sup> <sup>3</sup>)*x*<sup>5</sup> = *β*1*x*<sup>0</sup> <sup>1</sup>*E*<sup>1</sup> + *<sup>β</sup>*2*x*<sup>0</sup> <sup>2</sup>(*E*<sup>2</sup> + *<sup>E</sup>*˜ <sup>2</sup>) + *β*3*x*<sup>0</sup> <sup>3</sup>*E*<sup>3</sup> − (*β*<sup>4</sup> − *c*)*x*<sup>4</sup> − *β*5[*c* − *R*0(*α*<sup>3</sup> + *β*4)]*x*5.

Therefore, as inequality *<sup>R</sup>*<sup>0</sup> < *<sup>β</sup>*<sup>4</sup> *<sup>α</sup>*3+*β*<sup>4</sup> holds, there exists *<sup>c</sup>* > 0 such that *<sup>L</sup>* (*t*) < 0. By means of LaSalle's invariance principle [12,36], we deduce that the equilibrium equilibrium point *S*<sup>0</sup> of system (4) is globally asymptotically stable in R<sup>5</sup> +.

#### **6. Global Stability Analysis for** *S***<sup>+</sup>**

**Theorem 3.** *The positive equilibrium point S*<sup>+</sup> *is globally asymptotically stable in the open positive orthant* R<sup>5</sup> <sup>+</sup>*, regardless of the delay kernel k*(*s*)*, if a*<sup>3</sup> = 0*.*

**Proof.** Let us consider an arbitrary solution *xi*(*t*), *i* = 1, 5, of system (4), with initial conditions *ϕ<sup>i</sup>* : (−∞, 0] → (0, ∞). From Theorem 1 it follows that *xi*(*t*) > 0, for any *t* > 0 and *i* = 1, 5.

Taking into account that *x* − 1 − ln *x* ≥ 0, for any *x* > 0, we construct the Lyapunov function as follows:

$$L(t) = L\_1(t) + L\_2(t) + L\_3(t) + L\_4(t) + cL\_5(t) + cL\_{4,d}(t),\tag{8}$$

where

$$L\_k(t) = x\_k(t) - x\_k^+ - x\_k^+ \ln \frac{x\_k(t)}{x\_k^+}\ ,\; k = \overline{1,5}\ . \tag{9}$$

and

$$L\_{4,d}(t) = \int\_0^\infty \left( k(s) \int\_{t-s}^t \left( \mathbf{x}\_4(r) - \mathbf{x}\_4^+ - \mathbf{x}\_4^+ \ln \frac{\mathbf{x}\_4(r)}{\mathbf{x}\_4^+} \right) dr \right) ds$$

and we further denote:

$$E\_k(t) = 1 - \frac{\mathbf{x}\_k(t)}{\mathbf{x}\_k^+} + \ln \frac{\mathbf{x}\_k(t)}{\mathbf{x}\_k^+} \le 0 \quad , \quad \tilde{E}\_k(t) = 1 - \frac{\mathbf{x}\_k^+}{\mathbf{x}\_k(t)} + \ln \frac{\mathbf{x}\_k^+}{\mathbf{x}\_k(t)} \le 0 \quad , \quad k = \overline{1, 5}.$$

With the notations defined above, taking into account that *S*<sup>+</sup> is an equilibrium point of system (4) we have:

*L* <sup>1</sup>(*t*) = *x*˙1 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 1 *x*1 = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 1 *x*1 [*γ*<sup>1</sup> − *x*1*x*<sup>5</sup> + *α*3*x*<sup>4</sup> + *x*<sup>3</sup> − (*α*<sup>4</sup> + *β*1)*x*1] = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 1 *x*1 [*x*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> <sup>−</sup> *<sup>α</sup>*3*x*<sup>+</sup> <sup>4</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>3</sup> + (*α*<sup>4</sup> <sup>+</sup> *<sup>β</sup>*1)*x*<sup>+</sup> <sup>1</sup> − *x*1*x*<sup>5</sup> + *α*3*x*<sup>4</sup> + *x*<sup>3</sup> − (*α*<sup>4</sup> + *β*1)*x*1] = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 1 *x*1 *x*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>1</sup> *x*+ 1 *x*5 *x*+ 5 <sup>−</sup> *<sup>α</sup>*3*x*<sup>+</sup> 4 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>4</sup> *x*+ 4 <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 3 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*+ 3 + (*α*<sup>4</sup> + *β*1)*x*<sup>+</sup> 1 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>1</sup> *x*+ 1 = *x*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>1</sup> *x*+ 1 *x*5 *x*+ 5 − *x*+ 1 *x*1 <sup>+</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>−</sup> *<sup>α</sup>*3*x*<sup>+</sup> 4 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>4</sup> *x*+ 4 − *x*+ 1 *x*1 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 1 *x*1 *x*4 *x*+ 4 <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 3 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*+ 3 − *x*+ 1 *x*1 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 1 *x*1 *x*3 *x*+ 3 + (*α*<sup>4</sup> + *β*1)*x*<sup>+</sup> 1 <sup>2</sup> <sup>−</sup> *<sup>x</sup>*<sup>1</sup> *x*+ 1 − *x*+ 1 *x*1 <sup>≤</sup> *<sup>x</sup>*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>1</sup> *x*+ 1 *x*5 *x*+ 5 − *x*+ 1 *x*1 <sup>+</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>−</sup> *<sup>α</sup>*3*x*<sup>+</sup> 4 <sup>2</sup> <sup>+</sup> ln *<sup>x</sup>*<sup>+</sup> 1 *x*1 <sup>+</sup> ln *<sup>x</sup>*<sup>4</sup> *x*+ 4 <sup>−</sup> *<sup>x</sup>*<sup>4</sup> *x*+ 4 − *x*+ 1 *x*1 − *x*<sup>3</sup> <sup>2</sup> <sup>+</sup> ln *<sup>x</sup>*<sup>+</sup> 1 *x*1 <sup>+</sup> ln *<sup>x</sup>*<sup>3</sup> *x*+ 3 − *x*+ 1 *x*1 <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*+ 3 + (*α*<sup>4</sup> + *β*1)*x*<sup>+</sup> 1 *E*<sup>1</sup> + *E*˜ 1 = *x*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>1</sup> *x*+ 1 *x*5 *x*+ 5 − *x*+ 1 *x*1 <sup>+</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>−</sup> *<sup>α</sup>*3*x*<sup>+</sup> 4 *E*˜ <sup>1</sup> + *E*<sup>4</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 3 *E*˜ <sup>1</sup> + *E*<sup>3</sup> + (*α*<sup>4</sup> + *β*1)*x*<sup>+</sup> 1 *E*<sup>1</sup> + *E*˜ 1 .

Furthermore, we get:

*L* <sup>2</sup>(*t*) = *x*˙2 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 2 *x*2 = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 2 *x*2 (*γ*<sup>2</sup> − *α*2*x*2*x*<sup>5</sup> − *β*2*x*2) = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 2 *x*2 (*α*2*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> <sup>+</sup> *<sup>β</sup>*2*x*<sup>+</sup> <sup>2</sup> − *α*2*x*2*x*<sup>5</sup> − *β*2*x*2) = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 2 *x*2 *α*2*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>2</sup> *x*+ 2 *x*5 *x*+ 5 + *β*2*x*<sup>+</sup> 2 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>2</sup> *x*+ 2 = *α*2*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>2</sup> *x*+ 2 *x*5 *x*+ 5 <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 2 *x*2 <sup>+</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 + *β*2*x*<sup>+</sup> 2 <sup>2</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 2 *x*2 <sup>−</sup> *<sup>x</sup>*<sup>2</sup> *x*+ 2 = *α*2*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>2</sup> *x*+ 2 *x*5 *x*+ 5 <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 2 *x*2 <sup>+</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 + *β*2*x*<sup>+</sup> 2 *E*<sup>2</sup> + *E*˜ 2 .

Moreover, we obtain:

*L* <sup>3</sup>(*t*) = *x*˙3 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 3 *x*3 = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 3 *x*3 (*α*4*x*<sup>1</sup> − *x*<sup>3</sup> − *α*1*x*3*x*<sup>5</sup> − *β*3*x*3) = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 3 *x*3 *α*4(*x*<sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>1</sup> ) <sup>−</sup> (<sup>1</sup> <sup>+</sup> *<sup>β</sup>*3)(*x*<sup>3</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>3</sup> ) <sup>−</sup> *<sup>α</sup>*1(*x*3*x*<sup>5</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 ) = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 3 *x*3 *α*4*x*<sup>+</sup> 1 *x*1 *x*+ 1 − 1 <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>3</sup> (1 + *β*3) *x*3 *x*+ 3 − 1 <sup>−</sup> *<sup>α</sup>*1*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 *x*<sup>3</sup> *x*+ 3 *x*5 *x*+ 5 − 1 = *α*4*x*<sup>+</sup> 1 *x*1 *x*+ 1 <sup>−</sup> <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 3 *x*3 · *x*1 *x*+ 1 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 3 *x*3 + *x*<sup>+</sup> <sup>3</sup> (1 + *β*3) <sup>2</sup> <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*+ 3 − *x*+ 3 *x*3 + *α*1*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 3 *x*3 <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*+ 3 *x*5 *x*+ 5 <sup>+</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>≤</sup> *<sup>α</sup>*4*x*<sup>+</sup> 1 <sup>−</sup><sup>2</sup> <sup>−</sup> ln *<sup>x</sup>*<sup>+</sup> 3 *x*3 <sup>−</sup> ln *<sup>x</sup>*<sup>1</sup> *x*+ 1 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 3 *x*3 <sup>+</sup> *<sup>x</sup>*<sup>1</sup> *x*+ 1 + *x*<sup>+</sup> <sup>3</sup> (1 + *β*3) <sup>2</sup> <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*+ 3 − *x*+ 3 *x*3 + *α*1*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 3 *x*3 <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*+ 3 *x*5 *x*+ 5 <sup>+</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 = *α*4*x*<sup>+</sup> <sup>1</sup> (−*E*<sup>1</sup> <sup>−</sup> *<sup>E</sup>*˜ <sup>3</sup>)+(1 + *β*3)*x*<sup>+</sup> <sup>3</sup> (*E*<sup>3</sup> + *<sup>E</sup>*˜ <sup>3</sup>) + *α*1*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 3 *x*3 <sup>−</sup> *<sup>x</sup>*<sup>3</sup> *x*+ 3 *x*5 *x*+ 5 <sup>+</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 ,

and finally:

*L* <sup>4</sup>(*t*) = *x*˙4 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 4 *x*4 = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 4 *x*4 (*x*1*x*<sup>5</sup> + *α*2*x*2*x*<sup>5</sup> + *α*1*x*3*x*<sup>5</sup> − *α*3*x*<sup>4</sup> − *β*4*x*4) = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 4 *x*4 [*x*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> 5 *x*1*x*<sup>5</sup> *x*+ <sup>1</sup> *<sup>x</sup>*<sup>+</sup> 5 − 1 + *α*2*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> 5 *x*2*x*<sup>5</sup> *x*+ <sup>2</sup> *<sup>x</sup>*<sup>+</sup> 5 − 1 + *α*1*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 *x*3*x*<sup>5</sup> *x*+ <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 − 1 <sup>−</sup> (*α*<sup>3</sup> <sup>+</sup> *<sup>β</sup>*4)*x*<sup>+</sup> 4 *x*4 *x*+ 4 − 1 ] = *x*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> 5 *x*1*x*<sup>5</sup> *x*+ <sup>1</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>−</sup> <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>4</sup> *x*1*x*<sup>5</sup> *x*4*x*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 4 *x*4 + *α*2*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> 5 *x*2*x*<sup>5</sup> *x*+ <sup>2</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>−</sup> <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>4</sup> *x*2*x*<sup>5</sup> *x*4*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 4 *x*4 + *α*1*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 *x*3*x*<sup>5</sup> *x*+ <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>−</sup> <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>4</sup> *x*3*x*<sup>5</sup> *x*4*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 4 *x*4 + (*α*<sup>3</sup> + *β*4)*x*<sup>+</sup> 4 <sup>2</sup> <sup>−</sup> *<sup>x</sup>*<sup>4</sup> *x*+ 4 − *x*+ 4 *x*4 <sup>≤</sup> *<sup>x</sup>*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> 5 *x*1*x*<sup>5</sup> *x*+ <sup>1</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>−</sup> <sup>2</sup> <sup>−</sup> ln *<sup>x</sup>*<sup>+</sup> 4 *x*4 <sup>−</sup> ln *<sup>x</sup>*<sup>1</sup> *x*+ 1 <sup>−</sup> ln *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 4 *x*4 + *α*2*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> 5 *x*2*x*<sup>5</sup> *x*+ <sup>2</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>−</sup> <sup>2</sup> <sup>−</sup> ln *<sup>x</sup>*<sup>+</sup> 4 *x*4 <sup>−</sup> ln *<sup>x</sup>*<sup>2</sup> *x*+ 2 <sup>−</sup> ln *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 4 *x*4 + *α*1*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 *x*3*x*<sup>5</sup> *x*+ <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>−</sup> <sup>2</sup> <sup>−</sup> ln *<sup>x</sup>*<sup>+</sup> 4 *x*4 <sup>−</sup> ln *<sup>x</sup>*<sup>3</sup> *x*+ 3 <sup>−</sup> ln *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 4 *x*4 + (*α*<sup>3</sup> + *β*4)*x*<sup>+</sup> <sup>4</sup> (*E*<sup>4</sup> <sup>+</sup> *<sup>E</sup>*˜ 4).

Furthermore, taking into account the algebraic relations satisfied by the coordinates of the equilibrium point *S*+, we obtain:

*L* <sup>1</sup>(*t*) + *L* <sup>2</sup>(*t*) + *L* <sup>3</sup>(*t*) + *L* <sup>4</sup>(*t*) ≤ <sup>≤</sup> *<sup>x</sup>*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 1 *x*1 <sup>+</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>−</sup> <sup>2</sup> <sup>+</sup> ln *<sup>x</sup>*<sup>+</sup> 1 *x*1 <sup>−</sup> ln *<sup>x</sup>*<sup>+</sup> 4 *x*4 <sup>−</sup> ln *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 4 *x*4 + *α*2*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 2 *x*2 <sup>+</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>−</sup> <sup>2</sup> <sup>−</sup> ln *<sup>x</sup>*<sup>+</sup> 4 *x*4 <sup>+</sup> ln *<sup>x</sup>*<sup>+</sup> 2 *x*2 <sup>−</sup> ln *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 4 *x*4 + *α*1*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 3 *x*3 <sup>+</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>−</sup> <sup>2</sup> <sup>−</sup> ln *<sup>x</sup>*<sup>+</sup> 4 *x*4 <sup>+</sup> ln *<sup>x</sup>*<sup>+</sup> 3 *x*3 <sup>−</sup> ln *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>+</sup> *<sup>x</sup>*<sup>+</sup> 4 *x*4 <sup>−</sup> *<sup>α</sup>*3*x*<sup>+</sup> <sup>4</sup> (*E*˜ <sup>1</sup> <sup>+</sup> *<sup>E</sup>*4) <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>3</sup> (*E*˜ <sup>1</sup> + *E*3)+(*α*<sup>4</sup> + *β*1)*x*<sup>+</sup> <sup>1</sup> (*E*<sup>1</sup> <sup>+</sup> *<sup>E</sup>*˜ 1) + *β*2*x*<sup>+</sup> <sup>2</sup> (*E*˜ <sup>2</sup> <sup>+</sup> *<sup>E</sup>*2) <sup>−</sup> *<sup>α</sup>*4*x*<sup>+</sup> <sup>1</sup> (*E*<sup>1</sup> <sup>+</sup> *<sup>E</sup>*˜ <sup>3</sup>)+(1 + *β*3)*x*<sup>+</sup> <sup>3</sup> (*E*<sup>3</sup> + *<sup>E</sup>*˜ 3) + (*α*<sup>3</sup> + *β*4)*x*<sup>+</sup> <sup>4</sup> (*E*<sup>4</sup> <sup>+</sup> *<sup>E</sup>*˜ 4) = *x*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> (*E*˜ <sup>1</sup> <sup>−</sup> *<sup>E</sup>*<sup>5</sup> <sup>−</sup> *<sup>E</sup>*˜ <sup>4</sup>) + *α*2*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> (*E*˜ <sup>2</sup> <sup>−</sup> *<sup>E</sup>*<sup>5</sup> <sup>−</sup> *<sup>E</sup>*˜ 4) + *α*1*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> (*E*˜ <sup>3</sup> <sup>−</sup> *<sup>E</sup>*<sup>5</sup> <sup>−</sup> *<sup>E</sup>*˜ <sup>4</sup>) <sup>−</sup> *<sup>α</sup>*3*x*<sup>+</sup> <sup>4</sup> (*E*˜ <sup>1</sup> <sup>+</sup> *<sup>E</sup>*4) <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>3</sup> (*E*˜ <sup>1</sup> + *E*3) + (*α*<sup>4</sup> + *β*1)*x*<sup>+</sup> <sup>1</sup> (*E*<sup>1</sup> <sup>+</sup> *<sup>E</sup>*˜ <sup>1</sup>) + *β*2*x*<sup>+</sup> <sup>2</sup> (*E*˜ <sup>2</sup> <sup>+</sup> *<sup>E</sup>*2) <sup>−</sup> *<sup>α</sup>*4*x*<sup>+</sup> <sup>1</sup> (*E*<sup>1</sup> <sup>+</sup> *<sup>E</sup>*˜ 3) + (1 + *β*3)*x*<sup>+</sup> <sup>3</sup> (*E*<sup>3</sup> + *<sup>E</sup>*˜ <sup>3</sup>)+(*α*<sup>3</sup> + *β*4)*x*<sup>+</sup> <sup>4</sup> (*E*<sup>4</sup> <sup>+</sup> *<sup>E</sup>*˜ 4) = *E*˜ 1 *x*+ <sup>1</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> <sup>−</sup> *<sup>α</sup>*3*x*<sup>+</sup> <sup>4</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>3</sup> + (*α*<sup>4</sup> <sup>+</sup> *<sup>β</sup>*1)*x*<sup>+</sup> 1 + *β*1*x*<sup>+</sup> <sup>1</sup> *E*<sup>1</sup> + (*α*2*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> <sup>+</sup> *<sup>β</sup>*2*x*<sup>+</sup> <sup>2</sup> )*E*˜ <sup>2</sup> + *β*2*x*<sup>+</sup> <sup>2</sup> *E*<sup>2</sup> + *E*˜ 3[*α*1*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> <sup>−</sup> *<sup>α</sup>*4*x*<sup>+</sup> <sup>1</sup> + (<sup>1</sup> <sup>+</sup> *<sup>β</sup>*3)*x*<sup>+</sup> <sup>3</sup> ] + *<sup>β</sup>*3*x*<sup>+</sup> <sup>3</sup> *E*<sup>3</sup> + *E*˜ 4[−*x*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> <sup>−</sup> *<sup>α</sup>*2*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> <sup>−</sup> *<sup>α</sup>*1*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> + (*α*<sup>3</sup> <sup>+</sup> *<sup>β</sup>*4)*x*<sup>+</sup> <sup>4</sup> ] + *<sup>β</sup>*4*x*<sup>+</sup> <sup>4</sup> *E*<sup>4</sup> <sup>+</sup> *<sup>E</sup>*5[−*x*<sup>+</sup> <sup>1</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> <sup>−</sup> *<sup>α</sup>*2*x*<sup>+</sup> <sup>2</sup> *<sup>x</sup>*<sup>+</sup> <sup>5</sup> <sup>−</sup> *<sup>α</sup>*1*x*<sup>+</sup> <sup>3</sup> *<sup>x</sup>*<sup>+</sup> 5 ] = *γ*1*E*˜ <sup>1</sup> + *γ*2*E*˜ <sup>2</sup> + *β*1*x*<sup>+</sup> <sup>1</sup> *<sup>E</sup>*<sup>1</sup> <sup>+</sup> *<sup>β</sup>*2*x*<sup>+</sup> <sup>2</sup> *<sup>E</sup>*<sup>2</sup> <sup>+</sup> *<sup>β</sup>*3*x*<sup>+</sup> <sup>3</sup> *<sup>E</sup>*<sup>3</sup> <sup>+</sup> *<sup>β</sup>*4*x*<sup>+</sup> <sup>4</sup> *E*<sup>4</sup> <sup>−</sup> (*α*<sup>3</sup> <sup>+</sup> *<sup>β</sup>*4)*x*<sup>+</sup> <sup>4</sup> *E*5.

Denoting *<sup>D</sup>*[*x*4] = ∞ 0 *k*(*s*)*x*4(*t* − *s*)*ds* and employing the inequality *x* ≥ ln *x* + 1, for any *x* > 0, we have:

*L* <sup>5</sup>(*t*) = *x*˙5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 5 *x*5 = <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 5 *x*5 (*D*[*x*4] − *β*5*x*5) = *D*[*x*4] − *β*5*x*<sup>5</sup> − *D*[*x*4] *x*+ 5 *x*5 + *β*5*x*<sup>+</sup> 5 = (*D*[*x*4] <sup>−</sup> *<sup>x</sup>*4) + *<sup>x</sup>*<sup>+</sup> 4 *x*4 *x*+ 4 + *β*5*x*<sup>+</sup> 5 <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 4 *D*[*x*4]*x*<sup>+</sup> 5 *x*+ <sup>4</sup> *x*<sup>5</sup> <sup>≤</sup> (*D*[*x*4] <sup>−</sup> *<sup>x</sup>*4) + *<sup>x</sup>*<sup>+</sup> 4 *x*4 *x*+ 4 <sup>+</sup> <sup>1</sup> <sup>−</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>−</sup> *<sup>x</sup>*<sup>+</sup> 4 ln *<sup>D</sup>*[*x*4] *x*4 <sup>+</sup> ln *<sup>x</sup>*<sup>4</sup> *x*+ 4 <sup>+</sup> ln *<sup>x</sup>*<sup>+</sup> 5 *x*5 + 1 = (*D*[*x*4] <sup>−</sup> *<sup>x</sup>*4) + *<sup>x</sup>*<sup>+</sup> 4 *x*4 *x*+ 4 <sup>−</sup> ln *<sup>x</sup>*<sup>4</sup> *x*+ 4 <sup>−</sup> *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>+</sup> ln *<sup>x</sup>*<sup>5</sup> *x*+ 5 <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>4</sup> ln *<sup>D</sup>*[*x*4] *x*4 = (*D*[*x*4] <sup>−</sup> *<sup>x</sup>*4) + *<sup>x</sup>*<sup>+</sup> <sup>4</sup> (*E*<sup>5</sup> <sup>−</sup> *<sup>E</sup>*4) <sup>−</sup> *<sup>x</sup>*<sup>+</sup> <sup>4</sup> ln *<sup>D</sup>*[*x*4] *x*4 .

Moreover:

$$\begin{aligned} L\_{4,d}'(t) &= \int\_0^\infty k(s) \left( \mathbf{x}\_4(t) - \mathbf{x}\_4^+ - \mathbf{x}\_4^+ \ln \frac{\mathbf{x}\_4(t)}{\mathbf{x}\_4^+} - \mathbf{x}\_4(t-s) + \mathbf{x}\_4^+ + \mathbf{x}\_4^+ \ln \frac{\mathbf{x}\_4(t-s)}{\mathbf{x}\_4^+} \right) ds \\ &= (\mathbf{x}\_4 - D[\mathbf{x}\_4]) + \mathbf{x}\_4^+ \int\_0^\infty k(s) \ln \frac{\mathbf{x}\_4(t-s)}{\mathbf{x}\_4(t)} ds \end{aligned}$$

Therefore, combining the previous relations, we have:

$$\begin{split} L'(t) &= L\_1'(t) + L\_2'(t) + L\_3'(t) + L\_4'(t) + c \cdot L\_5'(t) + c \cdot L\_{4,d}'(t) \\ &\leq \gamma\_1 \tilde{E}\_1 + \gamma\_2 \tilde{E}\_2 + \beta\_1 \mathbf{x}\_1^+ \mathbf{E}\_1 + \beta\_2 \mathbf{x}\_2^+ \mathbf{E}\_2 + \beta\_3 \mathbf{x}\_3^+ \mathbf{E}\_3 + \beta\_4 \mathbf{x}\_4^+ \mathbf{E}\_4 - (\alpha\_3 + \beta\_4) \mathbf{x}\_4^+ \mathbf{E}\_5 \\ &+ c \mathbf{x}\_4^+ (E\_5 - \mathbf{E}\_4) - c \mathbf{x}\_4^+ \left( \ln \frac{D[\mathbf{x}\_4]}{\mathbf{x}\_4} - \int\_0^\infty k(s) \ln \frac{\mathbf{x}\_4(t-s)}{\mathbf{x}\_4(t)} ds \right) \\ &= \gamma\_1 \tilde{E}\_1 + \gamma\_2 \tilde{E}\_2 + \beta\_1 \mathbf{x}\_1^+ \mathbf{E}\_1 + \beta\_2 \mathbf{x}\_2^+ \mathbf{E}\_2 + \beta\_3 \mathbf{x}\_3^+ \mathbf{E}\_3 + (\beta\_4 - c) \mathbf{x}\_4^+ \mathbf{E}\_4 + (c - a\mathbf{x}\_3 - \beta\_4) \mathbf{x}\_4^+ \mathbf{E}\_5 \\ &- c \mathbf{x}\_4^+ \left( \ln \int\_0^\infty k(s) \mathbf{x}\_4(t-s) ds - \int\_0^\infty k(s) \ln \mathbf{x}\_4(t-s) ds \right). \end{split}$$

For *α*<sup>3</sup> = 0, choosing *c* = *β*4, it follows that:

$$\begin{split} L'(t) &\leq \gamma\_1 \ddot{E}\_1 + \gamma\_2 \ddot{E}\_2 + \beta\_1 \mathbf{x}\_1^+ E\_1 + \beta\_2 \mathbf{x}\_2^+ E\_2 + \beta\_3 \mathbf{x}\_3^+ E\_3 \\ &\quad - \beta\_4 \mathbf{x}\_4^+ \left( \ln \int\_0^\infty k(s) \mathbf{x}\_4(t-s) ds - \int\_0^\infty k(s) \ln \mathbf{x}\_4(t-s) ds \right). \end{split}$$

Jensen's inequality for probability density functions applied to the concave logarithmic function [37] provides that

$$\ln \int\_0^\infty k(s)x\_4(t-s)ds \ge \int\_0^\infty k(s)\ln x\_4(t-s)ds,$$

and hence, we obtain that *L* (*t*) ≤ 0, for any *t* ≥ 0. Hence, by means of LaSalle's invariance principle [12,36], we deduce that the equilibrium equilibrium point *S*<sup>+</sup> of system (4) is globally asymptotically stable in R<sup>5</sup> +.

#### **7. Numerical Simulations**

*7.1. Scenario 1: Global Asymptotic Stability of S*<sup>0</sup>

In this case, we first consider the following parameter values: *a*<sup>1</sup> = 687.5, *a*<sup>2</sup> = 0.0000152588, *a*<sup>3</sup> = 1, *a*<sup>4</sup> = 0.9375, *a*<sup>5</sup> = 0.125, *b*<sup>1</sup> = 0.5, *b*<sup>2</sup> = 0.51, *b*<sup>3</sup> = 0.03125, *b*<sup>4</sup> = 0.5, *b*<sup>5</sup> = 0.5, *c*<sup>1</sup> = 0.0000152588, *m*<sup>1</sup> = 25.5, *m*<sup>2</sup> = 0.0078125. For these parameter values, system (1) has a unique equilibrium point:

$$S^0 = (\mathcal{U}\_{0\prime} \mathcal{M}\_{0\prime} \mathcal{T}\_{0\prime} \mathcal{R}\_{0\prime} \mathcal{V}\_0) = (1000, 50, 6000, 0, 0).$$

We remark that the basic reproduction number *R*<sup>0</sup> satisfies inequality (7), and hence, based on Theorem 2, the equilibrium point *S*<sup>0</sup> is globally asymptotically stable, for any delay kernel considered in system (1).

In fact, in this scenario, the very low hiring rates for unemployed persons, immigrants and temporarily employed persons and the low value of the basic reproduction number, justify the existence of only one globally asymptotically stable equilibrium point, *S*0. For the numerical simulations shown in Figure 1, the following initial condition has been considered: *U*(0) = 280, *M*(0) = 50, *T*(0) = 2000, *R*(0) = 6400, *V*(0) = 40, which represents a plausible situation in a stable economic environment. However, due to the low hiring rates mentioned above, the trajectories of system (1) quickly converge to the "crisis" equilibrium *S*0, corresponding to no regular employment and no available vacancies.

**Figure 1.** Evolution of the state variables *U*(*t*), *M*(*t*), *T*(*t*), *R*(*t*) and *V*(*t*) in Scenario 1, with fixed initial conditions: *U*(0) = 280, *M*(0) = 50, *T*(0) = 2000, *R*(0) = 6400, *V*(0) = 40 and a discrete time delay *<sup>τ</sup>* <sup>∈</sup> [0, 100]; the trajectories converge to the equilibrium point *<sup>S</sup>*0=(1000, 50, 6000, 0, 0).

*7.2. Scenario 2: Global Asymptotic Stability of S*<sup>+</sup>

For the numerical simulations we considered: *a*<sup>1</sup> = 300, *a*<sup>2</sup> = 0.0303571, *a*<sup>3</sup> = 0.0582031, *a*<sup>4</sup> = 0.5, *a*<sup>5</sup> = 0.035, *b*<sup>1</sup> = 0.9375, *b*<sup>2</sup> = 0.5, *b*<sup>3</sup> = 0.02625, *b*<sup>4</sup> = 0.003125, *b*<sup>5</sup> = 0.5, *c*<sup>1</sup> = 0.00021875, *m*<sup>1</sup> = 60, *m*<sup>2</sup> = 0.0175. For these values of the parameters, we obtain a unique positive equilibrium point

$$S^{+} = (\mathcal{U}\_{+}, \mathcal{M}\_{+}, T\_{+}, \mathcal{R}\_{+}, V\_{+}) = (280, 50, 2000, 6400, 40).$$

In fact, in this scenario, we observe that the hiring rate *a*<sup>3</sup> is very small, and hence, global asymptotic stability can be expected for the equilibrium *S*+, in line with the theoretical results presented in Theorem 3. Indeed, in the numerical simulations shown in Figure 2, convergence to the positive equilibrium *S*<sup>+</sup> is observed, starting from a "crisis" initial condition: *U*(0) = 260, *M*(0) = 120, *T*(0) = 2100, *R*(0) = 1, *V*(0) = 1, with almost no regular employment or available vacancies. We emphasize that in this scenario, the convergence is very slow. Moreover, a larger time delay is associated with an even slower convergence to the positive equilibrium.

**=0 =20 =40 =60 =80 =100**

**Figure 2.** Evolution of the state variables *U*(*t*), *M*(*t*), *T*(*t*), *R*(*t*) and *V*(*t*) in Scenario 2, with initial conditions: *U*(0) = 260, *M*(0) = 120, *T*(0) = 2100, *R*(0) = 1, *V*(0) = 1 and a discrete time delay *<sup>τ</sup>* <sup>∈</sup> [0, 100]; the trajectories converge to *<sup>S</sup>*<sup>+</sup> <sup>=</sup> (280, 50, 2000, 6400, 40).

#### **8. Conclusions**

The present paper introduced a five-dimensional mathematical model that facilitates the understanding of new ways for studying the labour market while observing the levels of unemployment, migration, fixed term contractors, full time employment and the number of available vacancies. The distributed time delay has been incorporated to reflect the dependence of the rate of change of available vacancies on the past regular employment levels.

The positivity and boundedness of solutions have been provided. Using the basic reproduction number, the existence of equilibrium points has been discussed. As a consequence, the employment free equilibrium and positive equilibrium came into focus, for which we have undertaken a global stability analysis, regardless of the delay kernel included in the mathematical model.

Through the numerical simulations, the theoretical findings have been emphasized as well, further revealing that in an economic crisis scenario (Scenario 1), a fast convergence to the "crisis" equilibrium (with no regular employment or vacancies) is observed. However, in Scenario 2, if the initial conditions correspond to a crisis situation and the parameters of the system are adjusted, the convergence to the positive equilibrium is very slow, which reflects the slow recovery of the job market.

The findings of the paper could be used as input for future developments linked to population dynamics, and the decision making entities of the society could also use these findings in various strategies.

The proposed mathematical model is open ended that allows the introduction of future variables depicting the effects of the unemployment, such as poverty and unsocial behaviours that can lead to crime, as in [38].

**Author Contributions:** Conceptualization, E.K., M.N. and L.F.V.; methodology, E.K., M.N. and L.F.V.; software, E.K., M.N. and L.F.V.; validation, E.K., M.N. and L.F.V.; formal analysis, E.K., M.N. and L.F.V.; investigation, E.K., M.N. and L.F.V.; resources, E.K., M.N. and L.F.V.; data curation, E.K., M.N. and L.F.V.; writing—original draft preparation, E.K., M.N. and L.F.V.; writing—review and editing, E.K., M.N. and L.F.V.; visualization, E.K., M.N. and L.F.V.; supervision, E.K., M.N. and L.F.V.; project administration, E.K., M.N. and L.F.V.; funding acquisition, E.K., M.N. and L.F.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

