Discretization and Analysis of HIV-1 and HTLV-I Coinfection Model with Latent Reservoirs

Abstract

This article formulates and analyzes a discrete-time Human immunodeficiency virus type 1 (HIV-1) and human T-lymphotropic virus type I (HTLV-I) coinfection model with latent reservoirs. We consider that the HTLV-I infect the CD$4^{+}$T cells, while HIV-1 has two classes of target cells—CD$4^{+}$T cells and macrophages. The discrete-time model is obtained by discretizing the original continuous-time by the non-standard finite difference (NSFD) approach. We establish that NSFD maintains the positivity and boundedness of the model’s solutions. We derived four threshold parameters that determine the existence and stability of the four equilibria of the model. The Lyapunov method is used to examine the global stability of all equilibria. The analytical findings are supported via numerical simulation. The impact of latent reservoirs on the HIV-1 and HTLV-I co-dynamics is discussed. We show that incorporating the latent reservoirs into the HIV-1 and HTLV-I coinfection model will reduce the basic HIV-1 single-infection and HTLV-I single-infection reproductive numbers. We establish that neglecting the latent reservoirs will lead to overestimation of the required HIV-1 antiviral drugs. Moreover, we show that lengthening of the latent phase can suppress the progression of viral coinfection. This may draw the attention of scientists and pharmaceutical companies to create new treatments that prolong the latency period.

1. Introduction

Over the past decades, many scientists and researchers from the disciplines of mathematics, biology, and medicine have been interested in modeling the dynamics of viral infection within the host. Mathematical modeling of viral infection has a long history of helping to provide insight that is difficult to obtain through pure experiments. Examples of viral single-infection that have been modeled and studied are as follows: (i) chronic viral infections such as human immunodeficiency virus type 1 (HIV-1) [1], human T-lymphotropic virus type I (HTLV-I) [2], hepatitis B virus (HBV) [3] and hepatitis C virus (HCV) [4], (ii) respiratory viral infections such as influenza A virus (IAV) [5] and severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) [6,7], (iii) vector-borne viral infections such as dengue virus [8], chikungunya virus [9] and Zika virus [10]. The in-host viral coinfections were also modeled in recent years such as Zika/dengue [11], HIV-1/HTLV-I [12,13], IAV/SARS-CoV-2 [14,15], SARS-CoV-2/HIV-1 [16], SARS-CoV-2/HTLV-I [17], HIV-1/HCV [18] and HIV-1/HBV [19]. HIV-1 and HTLV-I are two dangerous retroviruses that attack the central component of the immune system, CD$4^{+}$T cells and can cause chronic diseases. HIV-1 causes acquired immunodeficiency syndrome (AIDS), while HTLV-I causes HTLV-I-associated myelopathy/tropical spastic paraparesis (HAM/TSP) and adult T-cell leukemia (ATL) diseases. Cytotoxic T lymphocytes (CTLs) and B cells play important roles in the immune response against viral infection. CTLs attack and kill the viral-infected cells, while B cells produce antibodies to neutralize viruses.

Since formulating the basic model of HIV-1 single-infection in [1], many modifications and developments have been made to make the model more accurate in describing the dynamics of the virus within the host. Some biological factors have been taken into account such as: (i) time delay [20,21], (ii) drug therapies [20,22], (iii) CTL immunity [1,23], (iv) antibody immunity [24], (v) reaction–diffusion [25], (vi) stochastic effects [26], and (vii) two target cells (CD4${}^{+}$T cells and macrophages) [22,27,28,29].

Latent HIV-1 reservoirs are considered a significant obstacle to HIV-1 elimination [30]. Latent HIV-1-infected cells can be activated after a period of time and then become viral production cells when drug therapies are stopped. These cells contains the HIV-1 virions; however, they do not produce them until they are stimulated. This fact causes an escapement of latent HIV-1-infected cells from the immune response. HIV-1 dynamics models with latent infected cells were developed in several works (see e.g., [23,26,30,31,32].

Stilianakis and Seydel [2] formulated a mathematical model for in-host HTLV-I dynamics. After that, several HTLV-I dynamics models were developed. HTLV-I infection models with CTL immune response were addressed in [33,34,35,36]. HTLV-I infection models have been incorporated with intracellular delay in [37], and with immune response delay in [35,37]. Reaction-diffusion HTLV-I infection models were investigated in [34].

HIV-1 and HTLV-1 share the methods of transmission between people through sexual relationships, infected sharp objects, blood transfusions and organ transplantation. Therefore, some nonlinear continuous-time models were recently formulated that describe the dynamics of HIV-1 and HTLV-1 coinfection in-host [12,13]. In [12,13], it was assumed that HIV-1 has one type of target cells, CD4${}^{+}$T cells. In fact, HIV-1 can also infect macrophages [31]. In [38], an HIV-1 and HTLV-I coinfection model was formulated by considering two types of target cells for HIV-1, CD4${}^{+}$T cells and macrophages:

$$\left\{\begin{array}{c}\frac{dx}{dt}=\stackrel{{\mathrm{CD}4}^{+}\mathrm{T}\mathrm{cells}\mathrm{production}}{\overbrace{{\zeta }_1}}-\stackrel{\mathrm{death}}{\overbrace{{\varrho }_{1}x}}-\stackrel{\mathrm{HIV}-1\mathrm{infectious}\mathrm{transmission}}{\overbrace{{\rho }_{1}xv}}-\stackrel{\mathrm{HTLV}-\mathrm{I}\mathrm{infectious}\mathrm{transmission}}{\overbrace{{\rho }_{3}xu}},\hfill \\ \frac{dy}{dt}=\stackrel{\mathrm{HIV}-1\mathrm{infectious}\mathrm{transmission}}{\overbrace{{\rho }_{1}xv}}-\stackrel{\mathrm{death}}{\overbrace{{\alpha }_{1}y}},\hfill \\ \frac{dw}{dt}=\stackrel{\mathrm{macrophages}\mathrm{production}}{\overbrace{{\zeta }_2}}-\stackrel{\mathrm{death}}{\overbrace{{\varrho }_{2}w}}-\stackrel{\mathrm{HIV}-1\mathrm{infectious}\mathrm{transmission}}{\overbrace{{\rho }_{2}wv}},\hfill \\ \frac{dz}{dt}=\stackrel{\mathrm{HIV}-1\mathrm{infectious}\mathrm{transmission}}{\overbrace{{\rho }_{2}wv}}-\stackrel{\mathrm{death}}{\overbrace{{\alpha }_{2}z}},\hfill \\ \frac{dv}{dt}=\stackrel{\mathrm{generation}\mathrm{of}\mathrm{HIV}-1}{\overbrace{{\beta }_1{\alpha }_{1}y+{\beta }_2{\alpha }_{2}z}}-\stackrel{\mathrm{death}}{\overbrace{\theta v}},\hfill \\ \frac{du}{dt}=\stackrel{\mathrm{HTLV}-\mathrm{I}\mathrm{infectious}\mathrm{transmission}}{\overbrace{{\rho }_{3}xu}}-\stackrel{\mathrm{death}}{\overbrace{\delta u}},\hfill \end{array}\right.$$

(1)

where x, y, w, z, v and u denote the concentrations of uninfected CD4${}^{+}$T cells, HIV-1-infected CD4${}^{+}$T cells, uninfected macrophages, HIV-1-infected macrophages, HIV-1 particles and HTLV-I-infected CD4${}^{+}$T cells, respectively. Model (1) was discretized by the NSFD method in [39]. Global stability of the discretized model was established using the Lyapunov technique.

In model (1), latent HIV-1 and HTLV-I reservoirs were not included. Therefore, model (1) has been extended in [40] by including three additional populations, latent HIV-1-infected CD4${}^{+}$T cells (g), latent HIV-1-infected macrophages (s) and latent HTLV-I-infected CD4${}^{+}$T cells (q) as:

$$\left\{\begin{array}{c}\frac{dx}{dt}={\zeta }_1-{\varrho }_{1}x-{\rho }_{1}xv-{\rho }_{3}xu,\hfill \\ \frac{dg}{dt}={\rho }_{1}xv-({\pi }_1+{\mu }_1)g,\hfill \\ \frac{dy}{dt}={\pi }_{1}g-{\alpha }_{1}y,\hfill \\ \frac{dw}{dt}={\zeta }_2-{\varrho }_{2}w-{\rho }_{2}wv,\hfill \\ \frac{ds}{dt}={\rho }_{2}wv-({\pi }_2+{\mu }_2)s,\hfill \\ \frac{dz}{dt}={\pi }_{2}s-{\alpha }_{2}z,\hfill \\ \frac{dv}{dt}={\beta }_1{\alpha }_{1}y+{\beta }_2{\alpha }_{2}z-\theta v,\hfill \\ \frac{dq}{dt}={\rho }_{3}xu-({\pi }_3+{\mu }_3)q,\hfill \\ \frac{du}{dt}={\pi }_{3}q-\delta u,\hfill \end{array}\right.$$

(2)

where the latent HIV-1-infected CD4${}^{+}$T cells, latent HIV-1-infected macrophages and latent HTLV-I-infected CD4${}^{+}$T cells are activated by rates ${\pi }_{1}g$, ${\pi }_{2}s$ and ${\pi }_{3}q$, respectively, while they die at rates ${\mu }_{1}g$, ${\mu }_{2}s$ and ${\mu }_{3}q$, respectively.

We note that model (2) is nonlinear continuous-time, and its exact analytical solution is unknown; therefore, discretization is unavoidable. Further, blood measurements from infected patients can only be available at discrete-time instants. As a result, an adequate discretization approach has to be chosen such that the basic and global properties of the original model is maintained. Mickens [41] introduced a non-standard finite difference (NSFD) scheme for solving different types of differential equations. NSFD was successfully utilized in discretizing several within-host virus dynamics models [42,43,44,45]. HIV-1 continuous-time models were discretized via the NSFD approach in [46,47,48,49]. A stability analysis of a discrete HIV-1 dynamics model with the Beddington–DeAngelis incidence and cure rate was studied in [47]. In [49], the global stability of discrete HIV-1 dynamics models with three classes of HIV-1-infected cells was studied. In [48], the HIV-1 dynamics model given by PDEs was discretized via the NSFD method. The Lyapunov method was used to prove the global stability of equilibria.

We mention that all mathematical models for HTLV-I single-infection and HIV-1/HTLV-I coinfection presented in the literature are given as continuous-time systems. The only exception is our recent paper [39], where a discrete HIV-1 and HTLV-I coinfection model is considered. In [39], the presence of latent reservoirs has not been modeled. The aim of the present article is to use the NSFD method to discretize an HIV-1 and HTLV-I coinfection model with latent reservoirs. We first establish the positivity and ultimate boundedness of the discrete-time model’s solutions, then calculate all equilibria and deduce their existence conditions. We examine the global stability of the four equilibria using the Lyapunov approach. We present some numerical simulations to clarify the theoretical results. We discuss the impact of latent reservoirs on the HIV-1 and HTLV-I co-dynamics.

2. The Discrete-Time Model

Classical numerical methods, such as Euler, Runge–Kutta and others, when used in solving nonlinear differential equations, suffer from numerical instability and bias when large step sizes are used in the numerical simulation [50]. In this situation, these numerical methods may provide non-physical solutions and can produce ‘false’ or ‘spurious’ fixed points, which are not fixed points of the original continuous-time model [51,52]. The NSFD method preserves the essential qualitative features of the original continuous-time model such as equilibria, positivity, boundedness and global behaviors of solutions independently of the selected step-size.

Applying the NSFD approach on system (2), we obtain

$$\begin{array}{cc}\hfill \frac{x_{n+1}-x_n}{{\rm Y}\left(h\right)}& ={\zeta }_1-{\varrho }_1{x}_{n+1}-{\rho }_1{x}_{n+1}{v}_n-{\rho }_3{x}_{n+1}{u}_n,\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \frac{g_{n+1}-g_n}{{\rm Y}\left(h\right)}& ={\rho }_1{x}_{n+1}{v}_n-({\pi }_1+{\mu }_1)g_{n+1},\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \frac{y_{n+1}-y_n}{{\rm Y}\left(h\right)}& ={\pi }_1{g}_{n+1}-{\alpha }_1{y}_{n+1},\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \frac{w_{n+1}-w_n}{{\rm Y}\left(h\right)}& ={\zeta }_2-{\varrho }_2{w}_{n+1}-{\rho }_2{w}_{n+1}{v}_n,\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \frac{s_{n+1}-s_n}{{\rm Y}\left(h\right)}& ={\rho }_2{w}_{n+1}{v}_n-({\pi }_2+{\mu }_2)s_{n+1},\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \frac{z_{n+1}-z_n}{{\rm Y}\left(h\right)}& ={\pi }_2{s}_{n+1}-{\alpha }_2{z}_{n+1},\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \frac{v_{n+1}-v_n}{{\rm Y}\left(h\right)}& ={\beta }_1{\alpha }_1{y}_{n+1}+{\beta }_2{\alpha }_2{z}_{n+1}-\theta v_{n+1},\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \frac{q_{n+1}-q_n}{{\rm Y}\left(h\right)}& ={\rho }_3{x}_{n+1}{u}_n-({\pi }_3+{\mu }_3)q_{n+1},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \frac{u_{n+1}-u_n}{{\rm Y}\left(h\right)}& ={\pi }_3{q}_{n+1}-\delta u_{n+1},\hfill \end{array}$$

(11)

where $h>0$ is the time step, and $(x_{n,},g_n,y_n,w_n,s_n,z_n,v_n,q_n,u_n)$ are the approximation of the solution $(x\left(t_n\right),g\left(t_n\right),y\left(t_n\right),w\left(t_n\right),s\left(t_n\right),z\left(t_n\right),v\left(t_n\right),q\left(t_n\right),u\left(t_n\right))$ of the system (2) at the discrete time point $t_n=nh,$ $n\in N=\{0,1,2,\dots \}$. The denominator function ${\rm Y}\left(h\right)$ is selected such that ${\rm Y}\left(h\right)=h+O\left(h^2\right)$. We consider the following form of ${\rm Y}\left(h\right)$

$${\rm Y}\left(h\right)=\frac{1-e^{-{\varrho }_{1}h}}{{\varrho }_1}.$$

(12)

The initial conditions of system (3)–(11) are

$$\begin{array}{cc}\hfill (x_0,g_0,y_0,w_0,s_0,z_0,v_0,q_0,u_0)& \in {\mathbb{R}}_{+}^9=\{(x,g,y,w,s,z,v,q,u)\mid x>0,g>0,y>0,w>0,s>0,\hfill \\ \hfill z& >0,v>0,q>0,u>0\}.\hfill \end{array}$$

(13)

3. Preliminaries

Let $\sigma =min\{{\varrho }_1,{\alpha }_1,\delta ,{\varrho }_2,{\alpha }_2,{\mu }_1,{\mu }_2,{\mu }_3\},$ and ${\zeta }_{12}={\zeta }_1+{\zeta }_2$ and define the regions

$$\begin{array}{cc}\hfill {\Gamma }_1& =\left\{(x,g,y,w,s,z,v,q,u)\in {\mathbb{R}}_{+}^9:x\le \frac{{\zeta }_1}{{\varrho }_1},w\le \frac{{\zeta }_2}{{\varrho }_2}\right.\hfill \\ \hfill & \left.0<x+g+y+w+s+z+q+u\le \frac{{\zeta }_{12}}{\sigma },v\le \frac{({\beta }_1{\alpha }_1+{\beta }_2{\alpha }_2){\zeta }_{12}}{\theta \sigma }\right\},\hfill \\ \hfill {\Gamma }_0& =\left\{(x,0,0,w,0,0,0,0,0)\in {\mathbb{R}}_{+}^9:x\ge 0,w\ge 0\right\}.\hfill \end{array}$$

Lemma 1.

Any solution ($x,g,y,w,s,z,v,q,u$) of model (3)–(11) with initial conditions (13) is positive and ultimately bounded.

Proof. 

Equations (3)–(11) imply that

$$\begin{array}{cc}\hfill x_{n+1}& =\frac{{\rm Y}\left(h\right){\zeta }_1+x_n}{1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_1{v}_n+{\rho }_3{u}_n)},\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill g_{n+1}& =\frac{{\rm Y}\left(h\right){\rho }_1{x}_{n+1}{v}_n+g_n}{1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1)},\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill y_{n+1}& =\frac{{\rm Y}\left(h\right){\pi }_1{g}_{n+1}+y_n}{1+{\rm Y}\left(h\right){\alpha }_1},\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill w_{n+1}& =\frac{{\rm Y}\left(h\right){\zeta }_2+w_n}{1+{\rm Y}\left(h\right)({\varrho }_2+{\rho }_2{v}_n)},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill s_{n+1}& =\frac{{\rm Y}\left(h\right){\rho }_2{w}_{n+1}{v}_n+s_n}{1+{\rm Y}\left(h\right)({\pi }_2+{\mu }_2)},\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill z_{n+1}& =\frac{{\rm Y}\left(h\right){\pi }_2{s}_{n+1}+z_n}{1+{\rm Y}\left(h\right){\alpha }_2},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill v_{n+1}& =\frac{{\rm Y}\left(h\right)({\beta }_1{\alpha }_1{y}_{n+1}+{\beta }_2{\alpha }_2{z}_{n+1})+v_n}{1+{\rm Y}\left(h\right)\theta },\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill q_{n+1}& =\frac{{\rm Y}\left(h\right){\rho }_3{x}_{n+1}{u}_n+q_n}{1+{\rm Y}\left(h\right)({\pi }_3+{\mu }_3)},\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill u_{n+1}& =\frac{{\rm Y}\left(h\right){\pi }_3{q}_{n+1}+u_n}{1+{\rm Y}\left(h\right)\delta },\hfill \end{array}$$

(22)

Since all parameters of model (2) are positive and the initial values are also positive, then by induction we obtain $x_n>0,g_n>0,y_n>0,w_n>0,s_n>0,z_n>0,v_n>0,q_n>0,$ and $u_n>0$ for all $n\in N$.

From Equations (3) and (6), we have

$$\begin{array}{cc}\hfill x_{n+1}-x_n& ={\rm Y}\left(h\right)\left[{\zeta }_1-{\varrho }_1{x}_{n+1}-{\rho }_1{x}_{n+1}{v}_n-{\rho }_3{x}_{n+1}{u}_n\right]\hfill \\ \hfill & \le {\rm Y}\left(h\right)\left[{\zeta }_1-{\varrho }_1{x}_{n+1}\right],\hfill \\ \hfill w_{n+1}-w_n& ={\rm Y}\left(h\right)\left[{\zeta }_2-{\varrho }_2{w}_{n+1}-{\rho }_2{w}_{n+1}{v}_n\right],\hfill \\ \hfill & \le {\rm Y}\left(h\right)\left[{\zeta }_2-{\varrho }_2{w}_{n+1}\right].\hfill \end{array}$$

It follows that

$$x_{n+1}\le \frac{x_n}{1+{\varrho }_1{\rm Y}\left(h\right)}+\frac{{\zeta }_1{\rm Y}\left(h\right)}{1+{\varrho }_1{\rm Y}\left(h\right)}\mathrm{and}{w}_{n+1}\le \frac{w_n}{1+{\varrho }_2{\rm Y}\left(h\right)}+\frac{{\zeta }_2{\rm Y}\left(h\right)}{1+{\varrho }_2{\rm Y}\left(h\right)}.$$

Using Lemma 2.2 in [53], we obtain

$$\begin{array}{cc}\hfill x_n& \le {\left(\frac{1}{1+{\rm Y}\left(h\right){\varrho }_1}\right)}^n{x}_0+\frac{{\zeta }_1}{{\varrho }_1}\left[1-{\left(\frac{1}{1+{\rm Y}\left(h\right){\varrho }_1}\right)}^n\right],\hfill \\ \hfill w_n& \le {\left(\frac{1}{1+{\rm Y}\left(h\right){\varrho }_2}\right)}^n{w}_0+\frac{{\zeta }_2}{{\varrho }_2}\left[1-{\left(\frac{1}{1+{\rm Y}\left(h\right){\varrho }_2}\right)}^n\right].\hfill \end{array}$$

Consequently, $\underset{n\to \infty }{lim\; sup}{x}_n\le \frac{{\zeta }_1}{{\varrho }_1}$ and $\underset{n\to \infty }{lim\; sup}{w}_n\le \frac{{\zeta }_2}{{\varrho }_2}$. Define a sequence $K_n$ as:

$$K_n=x_n+g_n+y_n+w_n+s_n+z_n+q_n+u_n.$$

Hence

$$\begin{array}{cc}\hfill K_{n+1}-K_n& =(x_{n+1}-x_n)+(g_{n+1}-g_n)+(y_{n+1}-y_n)+(w_{n+1}-w_n)+(s_{n+1}-s_n)+(z_{n+1}-z_n)\hfill \\ \hfill & +(q_{n+1}-q_n)+(u_{n+1}-u_n)\hfill \\ \hfill & ={\rm Y}\left(h\right)\left[{\zeta }_1-{\varrho }_1{x}_{n+1}-{\mu }_1{g}_{n+1}-{\alpha }_1{y}_{n+1}+{\zeta }_2-{\varrho }_2{w}_{n+1}-{\mu }_2{s}_{n+1}-{\alpha }_2{z}_{n+1}-{\mu }_3{q}_{n+1}\right.\hfill \\ \hfill & \left.-\delta u_{n+1}\right]\hfill \\ \hfill & \le {\rm Y}\left(h\right){\zeta }_{12}-{\rm Y}\left(h\right)\sigma [x_{n+1}+g_{n+1}+y_{n+1}+w_{n+1}+s_{n+1}+z_{n+1}+u_{n+1}]\hfill \\ \hfill & ={\rm Y}\left(h\right){\zeta }_{12}-{\rm Y}\left(h\right)\sigma K_{n+1}.\hfill \end{array}$$

Hence

$$K_{n+1}\le \frac{K_n}{1+{\rm Y}\left(h\right)\sigma }+\frac{{\rm Y}\left(h\right){\zeta }_{12}}{1+{\rm Y}\left(h\right)\sigma }.$$

Lemma 2.2 in [53] gives

$$K_n\le {\left(\frac{1}{1+{\rm Y}\left(h\right)\sigma }\right)}^n{K}_0+\frac{{\zeta }_{12}}{\sigma }\left[1-{\left(\frac{1}{1+{\rm Y}\left(h\right)\sigma }\right)}^n\right].$$

Then, $\underset{n\to \infty }{lim\; sup}{K}_n\le \frac{{\zeta }_{12}}{\sigma }$. From Equation (9), we obtain

$$\begin{array}{cc}\hfill v_{n+1}-v_n& ={\rm Y}\left(h\right)\left[{\beta }_1{\alpha }_1{y}_{n+1}+{\beta }_2{\alpha }_2{z}_{n+1}-\theta v_{n+1}\right]\hfill \\ \hfill & \le {\rm Y}\left(h\right)\left[{\beta }_1{\alpha }_1\frac{{\zeta }_{12}}{\sigma }+{\beta }_2{\alpha }_2\frac{{\zeta }_{12}}{\sigma }-\theta v_{n+1}\right].\hfill \end{array}$$

Hence

$$v_{n+1}\le \frac{v_n}{1+{\rm Y}\left(h\right)\theta }+\frac{{\rm Y}\left(h\right)({\beta }_1{\alpha }_1+{\beta }_2{\alpha }_2){\zeta }_{12}}{\left(1+{\rm Y}\left(h\right)\theta \right)\sigma }.$$

By induction, we obtain

$$v_n\le {\left(\frac{1}{1+{\rm Y}\left(h\right)\theta }\right)}^n{v}_0+\frac{({\beta }_1{\alpha }_1+{\beta }_2{\alpha }_2){\zeta }_{12}}{\theta \sigma }\left[1-{\left(\frac{1}{1+{\rm Y}\left(h\right)\theta }\right)}^n\right].$$

Consequently, $\underset{n\to \infty }{lim\; sup}{v}_n\le \frac{({\beta }_1{\alpha }_1+{\beta }_2{\alpha }_2){\zeta }_{12}}{\theta \sigma }.$ Therefore, $\left(x_{n,},g_n,y_n,w_n,s_n,z_n,v_n,q_n,u_n\right)$ converge to $\Gamma$ as $n\to \infty .$ □

4. Equilibria

Here, we calculate the model’s equilibria and deduce their existence conditions.

Lemma 2.

Model (3)–(11) has four equilibria that are determined by four threshold parameters $R_j>0,$ $j=0,1,2,3$:

(1) Infection-free equilibrium $E{Q}_0=(x^0,0,0,w^0,0,0,0,0,0)$, which always exists.

(2) Chronic HIV-1 single-infection equilibrium $E{Q}_1=(\widehat{x},\widehat{g},\widehat{y},\widehat{w},\widehat{s},\widehat{z},\widehat{v},0,0)$ exists when $R_0=R_{01}+R_{02}>1.$

(3) Chronic HTLV-I single-infection equilibrium $E{Q}_2=(\tilde{x},0,0,\tilde{w},0,0,0,\tilde{q},\tilde{u})$ exists when $R_1>1$.

(4) Chronic HIV-1/HTLV-I coinfection infection equilibrium $E{Q}_3=(\overline{x},\overline{g},\overline{y},\overline{w},\overline{s},\overline{z},\overline{v},\overline{q},\overline{u})$ exists when $\frac{R_1}{R_{01}}>1,$ $R_2>1$ and $R_3>1.$

Proof. 

Any equilibrium $EQ=(x,g,y,w,s,z,v,q,u)$ satisfies

$$\begin{array}{cc}\hfill 0& ={\zeta }_1-{\varrho }_{1}x-{\rho }_{1}xv-{\rho }_{3}xu,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill 0& ={\rho }_{1}xv-({\pi }_1+{\mu }_1)g\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill 0& ={\pi }_{1}g-{\alpha }_{1}y\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill 0& ={\zeta }_2-{\varrho }_{2}w-{\rho }_{2}wv,\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill 0& ={\rho }_{2}wv-({\pi }_2+{\mu }_2)s\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill 0& ={\pi }_{2}s-{\alpha }_{2}z,\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill 0& ={\beta }_1{\alpha }_{1}y+{\beta }_2{\alpha }_{2}z-\theta v,\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill 0& ={\rho }_{3}xu-({\pi }_3+{\mu }_3)q\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill 0& ={\pi }_{3}q-\delta u.\hfill \end{array}$$

(31)

From Equations (30) and (31), we obtain two options $u=0$ and $x=\frac{({\pi }_3+{\mu }_3)\delta }{{\rho }_3{\pi }_3}.$ First, we consider $u=0,$ then $q=0.$

From Equations (23) and (26), we obtain

$$x=\frac{{\zeta }_1}{{\varrho }_1+{\rho }_{1}v},w=\frac{{\zeta }_2}{{\varrho }_2+{\rho }_{2}v}.$$

(32)

and from Equations (25) and (28), we obtain

$$g=\frac{{\alpha }_1}{{\pi }_1}y\mathrm{and}s=\frac{{\alpha }_2}{{\pi }_2}z.$$

(33)

Now substituting in Equations (24) and (27), we obtain

$$y=\frac{{\rho }_{1}xv{\pi }_1}{{\alpha }_1({\pi }_1+{\mu }_1)}\mathrm{and}z=\frac{{\rho }_{2}wv{\pi }_2}{{\alpha }_2({\pi }_2+{\mu }_2)}.$$

(34)

Now substituting in Equation (29), we obtain

$$\left(\frac{{\beta }_1{\rho }_{1}x{\pi }_1}{{\pi }_1+{\mu }_1}+\frac{{\beta }_2{\rho }_{2}w{\pi }_2}{{\pi }_2+{\mu }_2}-\theta \right)v=0.$$

(35)

There are two solutions for Equation (35), $v=0$ and $\left(\frac{{\beta }_1{\rho }_1{\pi }_1}{{\pi }_1+{\mu }_1}x+\frac{{\beta }_2{\rho }_2{\pi }_2}{{\pi }_2+{\mu }_2}w-\theta \right)=0$. When $v=0$, we obtain $y=z=g=s=0$, $x=\frac{{\zeta }_1}{{\varrho }_1}$ and $w=\frac{{\zeta }_2}{{\varrho }_2}$, which gives the infection-free equilibrium $E{Q}_0=(x^0,0,0,w^0,0,0,0,0,0),$ where

$$x^0=\frac{{\zeta }_1}{{\varrho }_1}\mathrm{and}{w}^0=\frac{{\zeta }_2}{{\varrho }_2}.$$

When $v\ne 0$ and $\frac{{\beta }_1{\rho }_1{\pi }_1}{{\pi }_1+{\mu }_1}x+\frac{{\beta }_2{\rho }_2{\pi }_2}{{\pi }_2+{\mu }_2}w-\theta =0,$ then from Equation (32), we obtain

$$\frac{{\beta }_1{\rho }_1{\pi }_1{\zeta }_1}{({\pi }_1+{\mu }_1)({\varrho }_1+{\rho }_{1}v)}+\frac{{\beta }_2{\rho }_2{\pi }_2{\zeta }_2}{({\pi }_2+{\mu }_2)({\varrho }_2+{\rho }_{2}v)}-\theta =0.$$

We define a function H as:

$$H\left(v\right)=\frac{{\beta }_1{\rho }_1{\pi }_1{\zeta }_1}{({\pi }_1+{\mu }_1)({\varrho }_1+{\rho }_{1}v)}+\frac{{\beta }_2{\rho }_2{\pi }_2{\zeta }_2}{({\pi }_2+{\mu }_2)({\varrho }_2+{\rho }_{2}v)}-\theta =0.$$

Then

$$H\left(0\right)=\frac{{\beta }_1{\rho }_1{\pi }_1{\zeta }_1}{{\varrho }_1({\pi }_1+{\mu }_1)}+\frac{{\beta }_2{\rho }_2{\pi }_2{\zeta }_2}{{\varrho }_2({\pi }_2+{\mu }_2)}-\theta =\theta \left(\frac{{\beta }_1{\rho }_1{\pi }_1{\zeta }_1}{\theta {\varrho }_1({\pi }_1+{\mu }_1)}+\frac{{\beta }_2{\rho }_2{\pi }_2{\zeta }_2}{\theta {\varrho }_2({\pi }_2+{\mu }_2)}-1\right)=\theta (R_0-1),$$

where

$$R_0=R_{01}+R_{02},R_{01}=\frac{{\beta }_1{\rho }_1{\pi }_1{\zeta }_1}{\theta {\varrho }_1({\pi }_1+{\mu }_1)}\mathrm{and}{R}_{02}=\frac{{\beta }_2{\rho }_2{\pi }_2{\zeta }_2}{\theta {\varrho }_2({\pi }_2+{\mu }_2)}.$$

(36)

Thus, $H\left(0\right)>0,$ when $R_0>1$. The parameter $R_0$ represents the basic HIV-1 single-infection reproductive number.

$$\underset{v\to \infty }{lim}H\left(v\right)=-\theta .$$

Further,

$$H^{\prime }\left(v\right)=-\left(\frac{{\beta }_1{\zeta }_1{\pi }_1{\rho }_1^2}{({\pi }_1+{\mu }_1){({\varrho }_1+{\rho }_{1}v)}^2}+\frac{{\beta }_2{\zeta }_2{\pi }_2{\rho }_2^2}{({\pi }_2+{\mu }_2){({\varrho }_2+{\rho }_{2}v)}^2}\right)<0.$$

Hence, H is a strictly decreasing function of v, and thus, there exists a unique $\widehat{v}\in (0,\infty )$ such that $H\left(\widehat{v}\right)=0$. It follows that

$$\widehat{x}=\frac{{\zeta }_1}{{\varrho }_1+{\rho }_1\widehat{v}}>0\mathrm{and}\widehat{w}=\frac{{\zeta }_2}{{\varrho }_2+{\rho }_2\widehat{v}}>0.$$

Then, Equations (33) and (34) give

$$\widehat{y}=\frac{{\pi }_1{\rho }_1\widehat{x}\widehat{v}}{{\alpha }_1({\pi }_1+{\mu }_1)}>0,\widehat{z}=\frac{{\pi }_2{\rho }_2\widehat{w}\widehat{v}}{{\alpha }_2({\pi }_2+{\mu }_2)}>0,\widehat{g}=\frac{{\alpha }_1}{{\pi }_1}\widehat{y}>0\mathrm{and}\widehat{s}=\frac{{\alpha }_2}{{\pi }_2}\widehat{z}>0$$

Here, $\widehat{v}$ satisfies the following quadratic equation:

$$A{\widehat{v}}^2+B\widehat{v}+C=0,$$

(37)

with

$$A=\theta {\rho }_1{\rho }_2({\pi }_1+{\mu }_1)({\pi }_2+{\mu }_2),$$

$$B=\theta ({\pi }_1+{\mu }_1)({\pi }_2+{\mu }_2)({\varrho }_1{\rho }_2+{\varrho }_2{\rho }_1)-{\rho }_1{\rho }_2({\beta }_1{\zeta }_1{\pi }_1({\pi }_2+{\mu }_2)+{\beta }_2{\zeta }_2{\pi }_2({\pi }_1+{\mu }_1)),$$

$$\begin{array}{cc}\hfill C& =\theta {\varrho }_1{\varrho }_2({\pi }_1+{\mu }_1)({\pi }_2+{\mu }_2)-{\beta }_1{\varrho }_2{\rho }_1{\zeta }_1{\pi }_1({\pi }_2+{\mu }_2)-{\beta }_2{\varrho }_1{\rho }_2{\zeta }_2{\pi }_2({\pi }_1+{\mu }_1)\hfill \\ \hfill & =\theta {\varrho }_1{\varrho }_2({\pi }_1+{\mu }_1)({\pi }_2+{\mu }_2)\left(1-\frac{{\beta }_1{\rho }_1{\zeta }_1{\pi }_1}{{\varrho }_1\theta ({\pi }_1+{\mu }_1)}-\frac{{\beta }_2{\rho }_2{\zeta }_2{\pi }_2}{{\varrho }_2\theta ({\pi }_2+{\mu }_2)}\right)\hfill \\ \hfill & =-\theta {\varrho }_1{\varrho }_2({\pi }_1+{\mu }_1)({\pi }_2+{\mu }_2)(R_0-1).\hfill \end{array}$$

Obviously, $C<0$ when $R_0>1$. Equation (37) has a positive root as:

$$\widehat{v}=\frac{-B+\sqrt{B^2-4AC}}{2A}>0.$$

Hence, the chronic HIV-1 single-infection equilibrium $E{Q}_1=(\widehat{x},\widehat{g},\widehat{y},\widehat{w},\widehat{s},\widehat{z},\widehat{v},0,0)$ exists when $R_0>1$.

Now consider $\tilde{x}=\frac{({\pi }_3+{\mu }_3)\delta }{{\rho }_3{\pi }_3}$ and $u\ne 0$. Solving Equations (26)–(30), we obtain two equilibria: The chronic HTLV-I single-infection equilibrium $E{Q}_2=(\tilde{x},0,0,\tilde{w},0,0,0,\tilde{q},\tilde{u})$, where

$$\tilde{x}=\frac{({\pi }_3+{\mu }_3)\delta }{{\rho }_3{\pi }_3},\tilde{w}=\frac{{\zeta }_2}{{\varrho }_2}=w^0,\tilde{u}=\frac{{\varrho }_1}{{\rho }_3}(R_1-1),\tilde{q}=\frac{\delta {\varrho }_1}{{\pi }_3{\rho }_3}(R_1-1),$$

where

$$R_1=\frac{{\rho }_3{\zeta }_1{\pi }_3}{{\varrho }_1\delta ({\pi }_3+{\mu }_3)},$$

(38)

Parameter $R_1$ is the basic HTLV-I single-infection reproductive number. Consequently, $E{Q}_2$ exists when $R_1>1$. The other equilibrium is the chronic HIV-1/HTLV-I coinfection equilibrium $E{Q}_3=(\overline{x},\overline{g},\overline{y},\overline{w},\overline{s},\overline{z},\overline{v},\overline{q},\overline{u}),$ where

$$\begin{array}{cc}\hfill \overline{x}& =\frac{({\pi }_3+{\mu }_3)\delta }{{\rho }_3{\pi }_3}=\tilde{x},\overline{g}=\frac{{\varrho }_2{\rho }_1\delta ({\pi }_3+{\mu }_3)}{{\rho }_2{\rho }_3{\pi }_3({\pi }_1+{\mu }_1)}(R_2-1),\overline{y}=\frac{{\pi }_1{\varrho }_2{\rho }_1\delta ({\pi }_3+{\mu }_3)}{{\alpha }_1{\rho }_2{\rho }_3{\pi }_3({\pi }_1+{\mu }_1)}(R_2-1),\hfill \\ \hfill \overline{w}& =\frac{{\pi }_1{\beta }_1{\rho }_1\delta ({\pi }_2+{\mu }_2)({\pi }_3+{\mu }_3)}{{\beta }_2{\rho }_2{\rho }_3{\pi }_2{\pi }_3({\pi }_1+{\mu }_1)}\left(\frac{R_1}{R_{01}}-1\right),\hfill \\ \hfill \overline{s}& =\frac{{\varrho }_2{\pi }_1{\beta }_1{\rho }_1\delta ({\pi }_3+{\mu }_3)}{{\beta }_2{\rho }_2{\rho }_3{\pi }_2{\pi }_3({\pi }_1+{\mu }_1)}\left(\frac{R_1}{R_{01}}-1\right)(R_2-1)=\frac{{\varrho }_2\theta }{{\beta }_2{\rho }_2{\pi }_2}\left(\frac{R_{01}+R_{02}{R}_1}{R_1}-1\right),\hfill \\ \hfill \overline{z}& =\frac{{\varrho }_2{\beta }_1{\rho }_1{\pi }_1\delta ({\pi }_3+{\mu }_3)}{{\alpha }_2{\beta }_2{\rho }_2{\rho }_3{\pi }_3({\pi }_1+{\mu }_1)}\left(\frac{R_1}{R_{01}}-1\right)(R_2-1)=\frac{{\varrho }_2\theta }{{\beta }_2{\rho }_2{\alpha }_2}\left(\frac{R_{01}+R_{02}{R}_1}{R_1}-1\right),\hfill \\ \hfill \overline{v}& =\frac{{\varrho }_2}{{\rho }_2}(R_2-1),\overline{q}=\frac{\delta {\rho }_1{\varrho }_2}{{\rho }_2{\rho }_3{\pi }_3}\left[(R_2-1)(R_3-1)\right],\overline{u}=\frac{{\rho }_1{\varrho }_2}{{\rho }_2{\rho }_3}\left[(R_2-1)(R_3-1)\right].\hfill \end{array}$$

and

$$R_2=\frac{{\zeta }_2{\beta }_2{\rho }_2{\rho }_3{\pi }_2{\pi }_3({\pi }_1+{\mu }_1)}{{\varrho }_2{\beta }_1{\rho }_1\delta {\pi }_1({\pi }_2+{\mu }_2)({\pi }_3+{\mu }_3)(\frac{R_1}{R_{01}}-1)},R_3=\frac{{\varrho }_1{\rho }_2}{{\varrho }_2{\rho }_1}\left(\frac{R_1-1}{R_2-1}\right).$$

We can see that $E{Q}_3$ exists when $\frac{R_1}{R_{01}}>1,$ $R_2>1$ and $R_3>1$. □

5. Global Stability

In this section, we demonstrate the global asymptotic stability of all equilibria by establishing appropriate Lyapunov functions. Define a function $G\left(x\right)\ge 0$ as $G\left(x\right)=x-1-lnx$. We have

$$lnx\le x-1.$$

(39)

Theorem 1.

If $R_0\le 1$ and $R_1\le 1$, then $E{Q}_0=(x^0,0,0,w^0,0,0,0,0,0)$ is globally asymptotically stable (GAS) in ${\Gamma }_1$.

Proof. 

Define a discrete Lyapunov function ${\Lambda }_n(x_n,g_n,y_n,w_n,s_n,z_n,v_n,q_n,u_n)$ as

$$\begin{array}{cc}\hfill {\Lambda }_n& =\frac{1}{{\rm Y}\left(h\right)}\left[x^{0}G\left(\frac{x_n}{x^0}\right)+g_n+\frac{{\pi }_1+{\mu }_1}{{\pi }_1}{y}_n+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{w}^{0}G\left(\frac{w_n}{w^0}\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{s}_n\right.\hfill \\ \hfill & \left.+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{z}_n+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(1+{\rm Y}\left(h\right)\theta )v_n+q_n+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}(1+{\rm Y}\left(h\right)\delta )u_n\right].\hfill \end{array}$$

Clearly, ${\Lambda }_n>0$ for all $x_n>0,g_n>0,y_n>0,w_n>0,s_n>0,z_n>0,v_n>0,q_n>0,u_n>0$. In addition, ${\Lambda }_n(x^0,0,0,w^0,0,0,0,0,0)=0$. Evaluating the difference $\Delta {\Lambda }_n={\Lambda }_{n+1}-{\Lambda }_n$ as:

$$\begin{array}{cc}\hfill \Delta {\Lambda }_n& ={\Lambda }_{n+1}-{\Lambda }_n=\frac{1}{{\rm Y}\left(h\right)}\left[x^{0}G\left(\frac{x_{n+1}}{x^0}\right)+g_{n+1}+\frac{{\pi }_1+{\mu }_1}{{\pi }_1}{y}_{n+1}+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{w}^{0}G\left(\frac{w_{n+1}}{w^0}\right)\right.\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{s}_{n+1}+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{z}_{n+1}+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(1+{\rm Y}\left(h\right)\theta )v_{n+1}+q_{n+1}\hfill \\ \hfill & +\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}(1+{\rm Y}\left(h\right)\delta )u_{n+1}-x^{0}G\left(\frac{x_n}{x^0}\right)-g_n-\frac{{\pi }_1+{\mu }_1}{{\pi }_1}{y}_n-\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{w}^{0}G\left(\frac{w_n}{w^0}\right)\hfill \\ \hfill & -\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{s}_n-\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{z}_n-\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(1+{\rm Y}\left(h\right)\theta )v_n-q_n\hfill \\ \hfill & \left.-\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}(1+{\rm Y}\left(h\right)\delta )u_n\right]\hfill \\ \hfill & =\frac{1}{{\rm Y}\left(h\right)}\left[x^0\left(\frac{x_{n+1}-x_n}{x^0}+ln\left(\frac{x_n}{x_{n+1}}\right)\right)+\left(g_{n+1}-g_n\right)+\frac{{\pi }_1+{\mu }_1}{{\pi }_1}\left(y_{n+1}-y_n\right)\right.\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{w}^0\left(\frac{w_{n+1}-w_n}{w^0}+ln\left(\frac{w_n}{w_{n+1}}\right)\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(s_{n+1}-s_n\right)\hfill \\ \hfill & +\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(z_{n+1}-z_n\right)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(1+{\rm Y}\left(h\right)\theta )\left(v_{n+1}-v_n\right)+\left(q_{n+1}-q_n\right)\hfill \\ \hfill & \left.+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}(1+{\rm Y}\left(h\right)\delta )\left(u_{n+1}-u_n\right)\right].\hfill \end{array}$$

Using inequality (39), we obtain

$$\begin{array}{cc}\hfill \Delta {\Lambda }_n& \le \frac{1}{{\rm Y}\left(h\right)}\left(\left(x_{n+1}-x_n+x^0\left(\frac{x_n}{x_{n+1}}-1\right)\right)+\left(g_{n+1}-g_n\right)+\frac{{\pi }_1+{\mu }_1}{{\pi }_1}\left(y_{n+1}-y_n\right)\right.\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(w_{n+1}-w_n+w^0\left(\frac{w_n}{w_{n+1}}-1\right)\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(s_{n+1}-s_n\right)\hfill \\ \hfill & +\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(z_{n+1}-z_n\right)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(1+{\rm Y}\left(h\right)\theta )\left(v_{n+1}-v_n\right)+\left(q_{n+1}-q_n\right)\hfill \\ \hfill & \left.+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}(1+{\rm Y}\left(h\right)\delta )\left(u_{n+1}-u_n\right)\right)\hfill \\ \hfill & =\frac{1}{{\rm Y}\left(h\right)}\left(\left(1-\frac{x^0}{x_{n+1}}\right)\left(x_{n+1}-x_n\right)+\left(g_{n+1}-g_n\right)+\frac{{\pi }_1+{\mu }_1}{{\pi }_1}\left(y_{n+1}-y_n\right)\right.\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{w^0}{w_{n+1}}\right)\left(w_{n+1}-w_n\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(s_{n+1}-s_n\right)\hfill \\ \hfill & +\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(z_{n+1}-z_n\right)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(1+{\rm Y}\left(h\right)\theta )\left(v_{n+1}-v_n\right)+\left(q_{n+1}-q_n\right)\hfill \\ \hfill & \left.+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}(1+{\rm Y}\left(h\right)\delta )\left(u_{n+1}-u_n\right)\right).\hfill \end{array}$$

From Equations (3)–(11), we have

$$\begin{array}{cc}\hfill \Delta {\Lambda }_n& \le \left(1-\frac{x^0}{x_{n+1}}\right)\left({\zeta }_1-{\varrho }_1{x}_{n+1}-{\rho }_1{x}_{n+1}{v}_n-{\rho }_3{x}_{n+1}{u}_n\right)+({\rho }_1{x}_{n+1}{v}_n-({\pi }_1+{\mu }_1)g_{n+1})\hfill \\ \hfill & +\frac{{\pi }_1+{\mu }_1}{{\pi }_1}\left({\pi }_1{g}_{n+1}-{\alpha }_1{y}_{n+1}\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{w^0}{w_{n+1}}\right)\left({\zeta }_2-{\varrho }_2{w}_{n+1}-{\rho }_2{w}_{n+1}{v}_n\right)\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left({\rho }_2{w}_{n+1}{v}_n-({\pi }_2+{\mu }_2)s_{n+1}\right)+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left({\pi }_2{s}_{n+1}-{\alpha }_2{z}_{n+1}\right)\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left({\beta }_1{\alpha }_1{y}_{n+1}+{\beta }_2{\alpha }_2{z}_{n+1}-\theta v_{n+1}\right)+\frac{\theta }{{\beta }_1}\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}\left(v_{n+1}-v_n\right)\hfill \\ \hfill & +({\rho }_3{x}_{n+1}{u}_n-({\pi }_3+{\mu }_3)q_{n+1})+\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\left(u_{n+1}-u_n\right)+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\left({\pi }_3{q}_{n+1}-\delta u_{n+1}\right).\hfill \end{array}$$

Collecting terms yields

$$\begin{array}{cc}\hfill \Delta {\Lambda }_n& \le \left(1-\frac{x^0}{x_{n+1}}\right)\left({\zeta }_1-{\varrho }_1{x}_{n+1}\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{w^0}{w_{n+1}}\right)({\zeta }_2-{\varrho }_2{w}_{n+1})\hfill \\ \hfill & +\left(\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{\rho }_2{w}^0+{\rho }_1{x}^0-\frac{\theta }{{\beta }_1}\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}\right)v_n+\left({\rho }_3{x}^0-\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\right)u_n.\hfill \end{array}$$

We have ${\zeta }_1={\varrho }_1{x}^0,$ ${\zeta }_2={\varrho }_2{w}^0$, then we obtain

$$\begin{array}{cc}\hfill \Delta {\Lambda }_n& \le \left(1-\frac{x^0}{x_{n+1}}\right)({\varrho }_1{x}^0-{\varrho }_1{x}_{n+1})+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{w^0}{w_{n+1}}\right)({\varrho }_2{w}^0-{\varrho }_2{w}_{n+1})\hfill \\ \hfill & +\left(\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{\rho }_2{w}^0+{\rho }_1{x}^0-\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\right)v_n+\left({\rho }_3{x}^0-\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\right)u_n\hfill \\ \hfill & =-\frac{{\varrho }_1{\left(x_{n+1}-x^0\right)}^2}{x_{n+1}}-\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\frac{{\varrho }_2{(w_{n+1}-w^0)}^2}{w_{n+1}}\hfill \\ \hfill & +\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(\frac{{\rho }_1{\beta }_1{\zeta }_1{\pi }_1}{\theta {\varrho }_1\left({\pi }_1+{\mu }_1\right)}+\frac{{\rho }_2{\beta }_2{\zeta }_2{\pi }_2}{\theta {\varrho }_2({\pi }_2+{\mu }_2)}-1\right)v_n\hfill \\ \hfill & +\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\left(\frac{{\rho }_3{\zeta }_1{\pi }_3}{\delta {\varrho }_1({\pi }_3+{\mu }_3)}-1\right)u_n.\hfill \end{array}$$

From Equations (36) and (38), we can write

$$\begin{array}{cc}\hfill \Delta {\Lambda }_n& \le -{\varrho }_1\frac{{\left(x_{n+1}-x^0\right)}^2}{x_{n+1}}-\frac{{\beta }_2{\pi }_2{\varrho }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\frac{{(w_{n+1}-w^0)}^2}{w_{n+1}}+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(R_0-1)v_n\hfill \\ \hfill & +\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}(R_1-1)u_n.\hfill \end{array}$$

Since $R_0\le 1$ and $R_1\le 1,$ then ${\Lambda }_n$ is monotonically decreasing. Clearly ${\Lambda }_n\ge 0,$ and hence, there is a limit ${lim}_{n\to \infty }{\Lambda }_n\ge 0$ and thus ${lim}_{n\to \infty }\Delta {\Lambda }_n=0,$ which gives ${lim}_{n\to \infty }{x}_n=x^0,$ ${lim}_{n\to \infty }{w}_n=w^0,$ ${lim}_{n\to \infty }(R_0-1)v_n=0$ and ${lim}_{n\to \infty }(R_1-1)u_n=0.$ We consider four cases:

(i) $R_0=1$ and $R_1=1,$ and then from Equation (6),

$$0={\zeta }_2-{\varrho }_2{w}^0-{\rho }_2{w}^0\underset{n\to \infty }{lim}{v}_n\Rightarrow \underset{n\to \infty }{lim}{v}_n=0.$$

(40)

In addition, from Equations (3), (4), (7) and (9), we obtain

$$\begin{array}{cc}\hfill 0& ={\zeta }_1-{\varrho }_1{x}^0-{\rho }_1{x}^0\underset{n\to \infty }{lim}{v}_n-{\rho }_3{x}^0\underset{n\to \infty }{lim}{u}_n\Rightarrow \underset{n\to \infty }{lim}{u}_n=0,\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill 0& ={\rho }_1{x}^0\underset{n\to \infty }{lim}{v}_n-({\pi }_1+{\mu }_1)\underset{n\to \infty }{lim}{g}_{n+1}\Rightarrow \underset{n\to \infty }{lim}{g}_n=0,\hfill \end{array}$$

(42)

$$\begin{array}{cc}\hfill 0& ={\rho }_2{w}^0\underset{n\to \infty }{lim}{v}_n-({\pi }_2+{\mu }_2)\underset{n\to \infty }{lim}{s}_{n+1}\Rightarrow \underset{n\to \infty }{lim}{s}_n=0,\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill 0& ={\beta }_1{\alpha }_1\underset{n\to \infty }{lim}{y}_{n+1}+{\beta }_2{\alpha }_2\underset{n\to \infty }{lim}{z}_{n+1}-\theta \underset{n\to \infty }{lim}{v}_{n+1}\Rightarrow \underset{n\to \infty }{lim}{y}_n=\underset{n\to \infty }{lim}{z}_n=0.\hfill \end{array}$$

(44)

Therefore, from Equation (11), we obtain

$$0={\pi }_3\underset{n\to \infty }{lim}{q}_{n+1}-\delta \underset{n\to \infty }{lim}{u}_{n+1}\Rightarrow \underset{n\to \infty }{lim}{q}_n=0.$$

(45)

(ii) $R_0=1$, $R_1<1$ and ${lim}_{n\to \infty }{u}_n=0.$ Equations (40) and (42)–(45) yield ${lim}_{n\to \infty }{v}_n=0$, ${lim}_{n\to \infty }{g}_n=0,$ ${lim}_{n\to \infty }{y}_n=0$, ${lim}_{n\to \infty }{s}_n=0$, ${lim}_{n\to \infty }{z}_n=0$ and ${lim}_{n\to \infty }{q}_n=0$.

(iii) $R_0<1$, $R_1=1$ and ${lim}_{n\to \infty }{v}_n=0.$ Equations (41)–(45), give ${lim}_{n\to \infty }{u}_n=0,$ ${lim}_{n\to \infty }{g}_n=0,$ ${lim}_{n\to \infty }{y}_n=0$, ${lim}_{n\to \infty }{s}_n=0,$ ${lim}_{n\to \infty }{z}_n=0$ and ${lim}_{n\to \infty }{q}_n=0$.

(iv) $R_0<1$, $R_1<1,$ ${lim}_{n\to \infty }{v}_n=0$ and ${lim}_{n\to \infty }{u}_n=0$. From Equations (42)–(45), we obtain ${lim}_{n\to \infty }{g}_n=0,$ ${lim}_{n\to \infty }{y}_n=0,$ ${lim}_{n\to \infty }{s}_n=0,$ ${lim}_{n\to \infty }{z}_n=0$ and ${lim}_{n\to \infty }{q}_n=0$.

Consequently, if $R_0\le 1$ and $R_1\le 1$, then ${lim}_{n\to \infty }(x_n,g_n,y_n,w_n,s_n,z_n,v_n,q_n,u_n)=(x^0,0,0,w^0,0,0,0,0,0)$. This gives that $E{Q}_0$ is GAS. □

Theorem 2.

If $R_0>1$ and $R_1\le 1$ then $E{Q}_1=(\widehat{x},\widehat{g},\widehat{y},\widehat{w},\widehat{s},\widehat{z},\widehat{v},0,0)$ is GAS in ${\Gamma }_1\setminus {\Gamma }_0$.

Proof. 

Define

$$\begin{array}{cc}\hfill {\Theta }_n& =\frac{1}{{\rm Y}\left(h\right)}\left[\widehat{x}G\left(\frac{x_n}{\widehat{x}}\right)+\widehat{g}G\left(\frac{g_n}{\widehat{g}}\right)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}\widehat{y}G\left(\frac{y_n}{\widehat{y}}\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\widehat{w}G\left(\frac{w_n}{\widehat{w}}\right)\right.\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\widehat{s}G\left(\frac{s_n}{\widehat{s}}\right)+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\widehat{z}G\left(\frac{z_n}{\widehat{z}}\right)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\widehat{v}(1+{\rm Y}\left(h\right)\theta )G\left(\frac{v_n}{\widehat{v}}\right)+q_n\hfill \\ \hfill & \left.+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}(1+{\rm Y}\left(h\right)\delta )u_n\right].\hfill \end{array}$$

Clearly ${\Theta }_n>0$ for all $x_n>0,g_n>0,y_n>0,w_n>0,s_n>0,z_n>0,v_n>0,q_n>0,u_n>0$. In addition ${\Theta }_n(\widehat{x},\widehat{g},\widehat{y},\widehat{w},\widehat{s},\widehat{z},\widehat{v},0,0)=0$. Computing the difference $\Delta {\Theta }_n={\Theta }_{n+1}-{\Theta }_n$ as:

$$\begin{array}{cc}\hfill \Delta {\Theta }_n& =\frac{1}{{\rm Y}(h}\left[\widehat{x}\left(\frac{x_{n+1}}{\widehat{x}}-\frac{x_n}{\widehat{x}}+ln\left(\frac{x_n}{x_{n+1}}\right)\right)+\widehat{g}\left(\frac{g_{n+1}}{\widehat{g}}-\frac{g_n}{\widehat{g}}+ln\left(\frac{g_n}{g_{n+1}}\right)\right)\right.\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}\widehat{y}\left(\frac{y_{n+1}}{\widehat{y}}-\frac{y_n}{\widehat{y}}+ln\left(\frac{y_n}{y_{n+1}}\right)\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\widehat{w}\left(\frac{w_{n+1}}{\widehat{w}}-\frac{w_n}{\widehat{w}}+ln\left(\frac{w_n}{w_{n+1}}\right)\right)\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\widehat{s}\left(\frac{s_{n+1}}{\widehat{s}}-\frac{s_n}{\widehat{s}}+ln\left(\frac{s_n}{s_{n+1}}\right)\right)+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\widehat{z}\left(\frac{z_{n+1}}{\widehat{z}}-\frac{z_n}{\widehat{z}}+ln\left(\frac{z_n}{z_{n+1}}\right)\right)\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(1+{\rm Y}\left(h\right)\theta )\widehat{v}\left(\frac{v_{n+1}}{\widehat{v}}-\frac{v_n}{\widehat{v}}+ln\left(\frac{v_n}{v_{n+1}}\right)\right)+\left(q_{n+1}-q_n\right)\hfill \\ \hfill & \left.+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}(1+{\rm Y}\left(h\delta \right)\left(u_{n+1}-u_n\right)\right].\hfill \end{array}$$

Using inequality (39), we obtain

$$\begin{array}{cc}\hfill \Delta {\Theta }_n& \le \frac{1}{{\rm Y}\left(h\right)}\left[x_{n+1}-x_n+\widehat{x}\left(\frac{x_n}{x_{n+1}}-1\right)+g_{n+1}-g_n+\widehat{g}\left(\frac{g_n}{g_{n+1}}-1\right)\right.\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}\left(y_{n+1}-y_n+\widehat{y}\left(\frac{y_n}{y_{n+1}}-1\right)\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(w_{n+1}-w_n+\widehat{w}\left(\frac{w_n}{w_{n+1}}-1\right)\right)\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(s_{n+1}-s_n+\widehat{s}\left(\frac{s_n}{s_{n+1}}-1\right)\right)+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(z_{n+1}-z_n+\widehat{z}\left(\frac{z_n}{z_{n+1}}-1\right)\right)\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(v_{n+1}-v_n+\widehat{v}\left(\frac{v_n}{v_{n+1}}-1\right)\right)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\rm Y}\left(h\right)\theta \left(v_{n+1}-v_n+\widehat{v}ln\left(\frac{v_n}{v_{n+1}}\right)\right)\hfill \\ \hfill & \left.+\left(q_{n+1}-q_n\right)+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}(1+{\rm Y}\left(h\right)\delta )\left(u_{n+1}-u_n\right)\right].\hfill \end{array}$$

(46)

Inequality (46) can be written as:

$$\begin{array}{cc}\hfill \Delta {\Theta }_n& \le \frac{1}{{\rm Y}\left(h\right)}\left[\left(1-\frac{\widehat{x}}{x_{n+1}}\right)\left(x_{n+1}-x_n\right)+\left(1-\frac{\widehat{g}}{g_{n+1}}\right)\left(g_{n+1}-g_n\right)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}\left(1-\frac{\widehat{y}}{y_{n+1}}\right)\left(y_{n+1}-y_n\right)\right.\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\widehat{w}}{w_{n+1}}\right)\left(w_{n+1}-w_n\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\widehat{s}}{s_{n+1}}\right)\left(s_{n+1}-s_n\right)\hfill \\ \hfill & +\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(1-\frac{\widehat{z}}{z_{n+1}}\right)\left(z_{n+1}-z_n\right)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(1-\frac{\widehat{v}}{v_{n+1}}\right)\left(v_{n+1}-v_n\right)+\left(q_{n+1}-q_n\right)\hfill \\ \hfill & \left.+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}(1+{\rm Y}\left(h\right)\delta )\left(u_{n+1}-u_n\right)\right]+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{v}_{n+1}-\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{v}_n+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\widehat{v}ln\left(\frac{v_n}{v_{n+1}}\right).\hfill \end{array}$$

From Equations (3)–(11), we have

$$\begin{array}{cc}\hfill \Delta {\Theta }_n& \le \left(1-\frac{\widehat{x}}{x_{n+1}}\right)({\zeta }_1-{\varrho }_1{x}_{n+1}-{\rho }_1{x}_{n+1}{v}_n-{\rho }_3{x}_{n+1}{u}_n)+\left(1-\frac{\widehat{g}}{g_{n+1}}\right)({\rho }_1{x}_{n+1}{v}_n\hfill \\ \hfill & -({\pi }_1+{\mu }_1)g_{n+1})+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}\left(1-\frac{\widehat{y}}{y_{n+1}}\right)\left({\pi }_1{g}_{n+1}-{\alpha }_1{y}_{n+1}\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\widehat{w}}{w_{n+1}}\right)\hfill \\ \hfill & \times \left({\zeta }_2-{\varrho }_2{w}_{n+1}-{\rho }_2{w}_{n+1}{v}_n\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\widehat{s}}{s_{n+1}}\right)\left({\rho }_2{w}_{n+1}{v}_n-({\pi }_2+{\mu }_2)s_{n+1}\right)\hfill \\ \hfill & +\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(1-\frac{\widehat{z}}{z_{n+1}}\right)({\pi }_2{s}_{n+1}-{\alpha }_2{z}_{n+1})+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(1-\frac{\widehat{v}}{v_{n+1}}\right)\hfill \\ \hfill & \times \left({\beta }_1{\alpha }_1{y}_{n+1}+{\beta }_2{\alpha }_2{z}_{n+1}-\theta v_{n+1}\right)+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{v}_{n+1}-\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{v}_n+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\widehat{v}ln\left(\frac{v_n}{v_{n+1}}\right)\hfill \\ \hfill & +{\rho }_3{x}_{n+1}{u}_n-({\pi }_3+{\mu }_3)q_{n+1}+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\left({\pi }_3{q}_{n+1}-\delta u_{n+1}\right)+\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}(u_{n+1}-u_n).\hfill \end{array}$$

Collecting terms, we obtain

$$\begin{array}{cc}\hfill \Delta {\Theta }_n& \le \left(1-\frac{\widehat{x}}{x_{n+1}}\right)({\zeta }_1-{\varrho }_1{x}_{n+1})+\left({\rho }_1\widehat{x}+\frac{{\beta }_2{\rho }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\widehat{w}-\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\right)v_n+{\rho }_3\widehat{x}\widehat{u}\frac{u_n}{\widehat{u}}\hfill \\ \hfill & -{\rho }_1\widehat{x}\widehat{v}\frac{x_{n+1}{v}_n\widehat{g}}{\widehat{x}\widehat{v}{g}_{n+1}}+({\pi }_1+{\mu }_1)\widehat{g}-({\pi }_1+{\mu }_1)\widehat{g}\frac{g_{n+1}\widehat{y}}{\widehat{g}{y}_{n+1}}+\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\widehat{w}}{w_{n+1}}\right)({\zeta }_2-{\varrho }_2{w}_{n+1})-\frac{{\beta }_2{\rho }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\widehat{w}\widehat{v}\frac{\widehat{s}{w}_{n+1}{v}_n}{s_{n+1}\widehat{w}\widehat{v}}+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\widehat{s}\hfill \\ \hfill & -\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\widehat{s}\frac{\widehat{z}{s}_{n+1}}{z_{n+1}\widehat{s}}+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}-\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}{\alpha }_1\widehat{y}\frac{\widehat{v}{y}_{n+1}}{v_{n+1}\widehat{y}}-\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}\frac{\widehat{v}{z}_{n+1}}{v_{n+1}\widehat{z}}\hfill \\ \hfill & +\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\widehat{v}+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\widehat{v}ln\left(\frac{v_n}{v_{n+1}}\right)-\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}{u}_n.\hfill \end{array}$$

Utilizing the following conditions for $E{Q}_1$:

$$\begin{array}{cc}\hfill {\zeta }_1& ={\varrho }_1\widehat{x}+{\rho }_1\widehat{x}\widehat{v},{\rho }_1\widehat{x}\widehat{v}=({\pi }_1+{\mu }_1)\widehat{g},\hfill \\ \hfill {\pi }_1\widehat{g}& ={\alpha }_1\widehat{y},{\zeta }_2={\varrho }_2\widehat{w}+{\rho }_2\widehat{w}\widehat{v},\hfill \\ \hfill {\rho }_2\widehat{w}\widehat{v}& =({\pi }_2+{\mu }_2)\widehat{s},{\pi }_2\widehat{s}={\alpha }_2\widehat{z},\hfill \\ \hfill \theta \widehat{v}& ={\beta }_1{\alpha }_1\widehat{y}+{\beta }_2{\alpha }_2\widehat{z},\hfill \end{array}$$

we obtain

$${\zeta }_1={\varrho }_1\widehat{x}+\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y},{\zeta }_2={\varrho }_2\widehat{w}+\frac{({\pi }_2+{\mu }_2)}{{\pi }_2}{\alpha }_2\widehat{z}.$$

Then

$$\left({\rho }_1\widehat{x}+\frac{{\beta }_2{\rho }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\widehat{w}-\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\right)v_n=0,$$

and

$$\begin{array}{cc}\hfill \Delta {\Theta }_n& \le \left(1-\frac{\widehat{x}}{x_{n+1}}\right)({\varrho }_1\widehat{x}-{\varrho }_1{x}_{n+1})+\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}-\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}\frac{\widehat{x}}{x_{n+1}}+{\rho }_3\widehat{x}\widehat{u}\frac{u_n}{\widehat{u}}\hfill \\ \hfill & -\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}\frac{x_{n+1}{v}_n\widehat{g}}{\widehat{x}\widehat{v}{g}_{n+1}}+\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}-\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}\frac{g_{n+1}\widehat{y}}{\widehat{g}{y}_{n+1}}+\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\widehat{w}}{w_{n+1}}\right)({\varrho }_2\widehat{w}-{\varrho }_2{w}_{n+1})+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}-\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}\frac{\widehat{w}}{w_{n+1}}\hfill \\ \hfill & -\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}\frac{\widehat{s}{w}_{n+1}{v}_n}{s_{n+1}\widehat{w}\widehat{v}}+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}-\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}\frac{\widehat{z}{s}_{n+1}}{z_{n+1}\widehat{s}}+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}\hfill \\ \hfill & -\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}\frac{\widehat{v}{y}_{n+1}}{v_{n+1}\widehat{y}}-\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}\frac{\widehat{v}{z}_{n+1}}{v_{n+1}\widehat{z}}+\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}\hfill \\ \hfill & +\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}ln\left(\frac{v_n}{v_{n+1}}\right)+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}ln\left(\frac{v_n}{v_{n+1}}\right)-\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}{u}_n.\hfill \end{array}$$

(47)

Inequality (47) takes the form

$$\begin{array}{cc}\hfill \Delta {\Theta }_n& \le -{\varrho }_1\frac{{(x_{n+1}-\widehat{x})}^2}{x_{n+1}}-\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\frac{{\varrho }_2{(w_{n+1}-\widehat{w})}^2}{w_{n+1}}+\left({\rho }_3\widehat{x}-\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\right)u_n\hfill \\ \hfill & +\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}\left[4-\frac{\widehat{x}}{x_{n+1}}-\frac{x_{n+1}{v}_n\widehat{g}}{\widehat{x}\widehat{v}{g}_{n+1}}-\frac{g_{n+1}\widehat{y}}{\widehat{g}{y}_{n+1}}-\frac{\widehat{v}{y}_{n+1}}{v_{n+1}\widehat{y}}+ln\left(\frac{v_n}{v_{n+1}}\right)\right]\hfill \\ \hfill & +\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}\left[4-\frac{\widehat{w}}{w_{n+1}}-\frac{\widehat{s}{w}_{n+1}{v}_n}{s_{n+1}\widehat{w}\widehat{v}}-\frac{\widehat{z}{s}_{n+1}}{z_{n+1}\widehat{s}}-\frac{\widehat{v}{z}_{n+1}}{v_{n+1}\widehat{z}}+ln\left(\frac{v_n}{v_{n+1}}\right)\right].\hfill \end{array}$$

Using the following equalities:

$$\begin{array}{cc}\hfill ln\left(\frac{v_n}{v_{n+1}}\right)& =ln\left(\frac{\widehat{x}}{x_{n+1}}\right)+ln\left(\frac{x_{n+1}{v}_n\widehat{g}}{\widehat{x}\widehat{v}{g}_{n+1}}\right)+ln\left(\frac{g_{n+1}\widehat{y}}{\widehat{g}{y}_{n+1}}\right)+ln\left(\frac{\widehat{v}{y}_{n+1}}{v_{n+1}\widehat{y}}\right),\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill ln\left(\frac{v_n}{v_{n+1}}\right)& =ln\left(\frac{\widehat{w}}{w_{n+1}}\right)+ln\left(\frac{\widehat{s}{w}_{n+1}{v}_n}{s_{n+1}\widehat{w}\widehat{v}}\right)+ln\left(\frac{\widehat{z}{s}_{n+1}}{z_{n+1}\widehat{s}}\right)+ln\left(\frac{\widehat{v}{z}_{n+1}}{v_{n+1}\widehat{z}}\right).\hfill \end{array}$$

(49)

We obtain

$$\begin{array}{cc}\hfill \Delta {\Theta }_n& \le -{\varrho }_1\frac{{(x_{n+1}-\widehat{x})}^2}{x_{n+1}}-\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\frac{{\varrho }_2{(w_{n+1}-\widehat{w})}^2}{w_{n+1}}+\left(\frac{{\rho }_3{\zeta }_1}{{\varrho }_1+{\rho }_1\widehat{v}}-\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\right)u_n\hfill \\ \hfill & -\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}\left[G\left(\frac{\widehat{x}}{x_{n+1}}\right)+G\left(\frac{x_{n+1}{v}_n\widehat{g}}{\widehat{x}\widehat{v}{g}_{n+1}}\right)+G\left(\frac{g_{n+1}\widehat{y}}{\widehat{g}{y}_{n+1}}\right)+G\left(\frac{\widehat{v}{y}_{n+1}}{v_{n+1}\widehat{y}}\right)\right]\hfill \\ \hfill & -\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}\left[G\left(\frac{\widehat{w}}{w_{n+1}}\right)+G\left(\frac{\widehat{s}{w}_{n+1}{v}_n}{s_{n+1}\widehat{w}\widehat{v}}\right)+G\left(\frac{\widehat{z}{s}_{n+1}}{z_{n+1}\widehat{s}}\right)+G\left(\frac{\widehat{v}{z}_{n+1}}{v_{n+1}\widehat{z}}\right)\right].\hfill \end{array}$$

Since $R_0>1$, then $\widehat{v}>0$ and $\frac{{\rho }_3{\zeta }_1}{{\varrho }_1+{\rho }_1\widehat{v}}<\frac{{\rho }_3{\zeta }_1}{{\varrho }_1}$. Therefore, we obtain

$$\begin{array}{cc}\hfill \Delta {\Theta }_n& \le -{\varrho }_1\frac{{(x_{n+1}-\widehat{x})}^2}{x_{n+1}}-\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\frac{{\varrho }_2{(w_{n+1}-\widehat{w})}^2}{w_{n+1}}+\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\left(R_1-1\right)u_n\hfill \\ \hfill & -\frac{({\pi }_1+{\mu }_1)}{{\pi }_1}{\alpha }_1\widehat{y}\left[G\left(\frac{\widehat{x}}{x_{n+1}}\right)+G\left(\frac{x_{n+1}{v}_n\widehat{g}}{\widehat{x}\widehat{v}{g}_{n+1}}\right)+G\left(\frac{g_{n+1}\widehat{y}}{\widehat{g}{y}_{n+1}}\right)+G\left(\frac{\widehat{v}{y}_{n+1}}{v_{n+1}\widehat{y}}\right)\right]\hfill \\ \hfill & -\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\widehat{z}\left[G\left(\frac{\widehat{w}}{w_{n+1}}\right)+G\left(\frac{\widehat{s}{w}_{n+1}{v}_n}{s_{n+1}\widehat{w}\widehat{v}}\right)+G\left(\frac{\widehat{z}{s}_{n+1}}{z_{n+1}\widehat{s}}\right)+G\left(\frac{\widehat{v}{z}_{n+1}}{v_{n+1}\widehat{z}}\right)\right].\hfill \end{array}$$

Since $R_0>1$ and if $R_1\le 1$, then ${\Theta }_n$ is monotonically decreasing. We have ${\Theta }_n\ge 0,$ and then there is a limit ${lim}_{n\to \infty }{\Theta }_n\ge 0$ and hence ${lim}_{n\to \infty }\Delta {\Theta }_n=0,$ which implies ${lim}_{n\to \infty }{x}_n=\widehat{x}$, ${lim}_{n\to \infty }{g}_n=\widehat{g},$ ${lim}_{n\to \infty }{y}_n=\widehat{y}$, ${lim}_{n\to \infty }{w}_n=\widehat{w}$, ${lim}_{n\to \infty }{s}_n=\widehat{s}$, ${lim}_{n\to \infty }{z}_n=\widehat{z}$, ${lim}_{n\to \infty }{v}_n=\widehat{v}$, and ${lim}_{n\to \infty }\left(R_1-1\right)u_n=0$. Now, we address two cases:

(i) $R_1=1$. From Equation (3), we obtain

$$0={\zeta }_1-{\varrho }_1\widehat{x}-{\rho }_1\widehat{x}\widehat{v}-{\rho }_3\widehat{x}\underset{n\to \infty }{lim}{u}_n\Rightarrow \underset{n\to \infty }{lim}{u}_n=0.$$

From Equation (11), we obtain

$$0={\pi }_3\underset{n\to \infty }{lim}{q}_{n+1}-\delta \underset{n\to \infty }{lim}{u}_{n+1}\Rightarrow \underset{n\to \infty }{lim}{q}_n=0.$$

(50)

(ii) $R_1<1$ and then ${lim}_{n\to \infty }{u}_n=0$. From Equation (50), we obtain ${lim}_{n\to \infty }{q}_n=0.$ Hence, $E{Q}_1$ is GAS. □

Theorem 3.

If $R_1>1$ and $R_{02}+\frac{R_{01}}{R_1}\le 1$, then $E{Q}_2=(\tilde{x},0,0,\tilde{w},0,0,0,\tilde{q},\tilde{u})$ is GAS in ${\Gamma }_1\setminus {\Gamma }_0$.

Proof. 

Consider a discrete Lyapunov function ${\Phi }_n$ as:

$$\begin{array}{cc}\hfill {\Phi }_n& =\frac{1}{{\rm Y}\left(h\right)}\left[\tilde{x}G\left(\frac{x_n}{\tilde{x}}\right)+g_n+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}{y}_n+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\tilde{w}G\left(\frac{w_n}{\tilde{w}}\right)\right.\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{s}_n+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{z}_n+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(1+{\rm Y}\left(h\right)\theta )v_n\hfill \\ \hfill & \left.+\tilde{q}G\left(\frac{q_n}{\tilde{q}}\right)+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}(1+{\rm Y}\left(h\right)\delta )\tilde{u}G\left(\frac{u_n}{\tilde{u}}\right)\right].\hfill \end{array}$$

Computing the difference $\Delta {\Phi }_n={\Phi }_{n+1}-{\Phi }_n$ as:

$$\begin{array}{cc}\hfill \Delta {\Phi }_n& =\frac{1}{{\rm Y}\left(h\right)}\left[\tilde{x}\left(\frac{x_{n+1}}{\tilde{x}}-\frac{x_n}{\tilde{x}}+ln\left(\frac{x_n}{x_{n+1}}\right)\right)+(g_{n+1}-g_n)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}(y_{n+1}-y_n)\right.\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\tilde{w}\left(\frac{w_{n+1}}{\tilde{w}}-\frac{w_n}{\tilde{w}}+ln\left(\frac{w_n}{w_{n+1}}\right)\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}(s_{n+1}-s_n)\hfill \\ \hfill & +\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(z_{n+1}-z_n)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(v_{n+1}-v_n)+\tilde{q}\left(\frac{q_{n+1}}{\tilde{q}}-\frac{q_n}{\tilde{q}}+ln\left(\frac{q_n}{q_{n+1}}\right)\right)\hfill \\ \hfill & \left.+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\tilde{u}\left(\frac{u_{n+1}}{\tilde{u}}-\frac{u_n}{\tilde{u}}+ln\left(\frac{u_n}{u_{n+1}}\right)\right)\right]+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta \tilde{u}\left(G\left(\frac{u_{n+1}}{\tilde{u}}\right)-G\left(\frac{u_n}{\tilde{u}}\right)\right)\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\theta v_{n+1}-\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\theta v_n.\hfill \end{array}$$

Using inequality (39), we obtain

$$\begin{array}{cc}\hfill \Delta {\Phi }_n& \le \frac{1}{{\rm Y}\left(h\right)}\left[\tilde{x}\left(\frac{x_{n+1}-x_n}{\tilde{x}}+\frac{x_n}{x_{n+1}}-1\right)+(g_{n+1}-g_n)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}(y_{n+1}-y_n)\right.\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\tilde{w}\left(\frac{w_{n+1}-w_n}{\tilde{w}}+\frac{w_n}{w_{n+1}}-1\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}(s_{n+1}-s_n)\hfill \\ \hfill & +\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(z_{n+1}-z_n)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(v_{n+1}-v_n)+\tilde{q}\left(\frac{q_{n+1}-q_n}{\tilde{q}}+\frac{q_n}{q_{n+1}}-1\right)\hfill \\ \hfill & \left.+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\tilde{u}\left(\frac{u_{n+1}-u_n}{\tilde{u}}+\frac{u_n}{u_{n+1}}-1\right)\right]+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta \tilde{u}\left(G\left(\frac{u_{n+1}}{\tilde{u}}\right)-G\left(\frac{u_n}{\tilde{u}}\right)\right)\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\theta v_{n+1}-\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\theta v_n.\hfill \end{array}$$

(51)

We write inequality (51) as:

$$\begin{array}{cc}\hfill \Delta {\Phi }_n& \le \frac{1}{{\rm Y}\left(h\right)}\left[\left(1-\frac{\tilde{x}}{x_{n+1}}\right)(x_{n+1}-x_n)+(g_{n+1}-g_n)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}(y_{n+1}-y_n)\right.\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\tilde{w}}{w_{n+1}}\right)\left(w_{n+1}-w_n\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}(s_{n+1}-s_n)\hfill \\ \hfill & +\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(z_{n+1}-z_n)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(v_{n+1}-v_n)+\left(1-\frac{\tilde{q}}{q_{n+1}}\right)\left(q_{n+1}-q_n\right)\hfill \\ \hfill & \left.+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\left(1-\frac{\tilde{u}}{u_{n+1}}\right)\left(u_{n+1}-u_n\right)\right]+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta \tilde{u}\left(G\left(\frac{u_{n+1}}{\tilde{u}}\right)-G\left(\frac{u_n}{\tilde{u}}\right)\right)\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\theta v_{n+1}-\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\theta v_n\hfill \end{array}$$

From Equations (3)–(11), we have

$$\begin{array}{cc}\hfill \Delta {\Phi }_n& \le \left(1-\frac{\tilde{x}}{x_{n+1}}\right)({\zeta }_1-{\varrho }_1{x}_{n+1}-{\rho }_1{x}_{n+1}{v}_n-{\rho }_3{x}_{n+1}{u}_n)+({\rho }_1{x}_{n+1}{v}_n-({\pi }_1+{\mu }_1)g_{n+1})\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}({\pi }_1{g}_{n+1}-{\alpha }_1{y}_{n+1})+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\tilde{w}}{w_{n+1}}\right)\left({\zeta }_2-{\varrho }_2{w}_{n+1}-{\rho }_2{w}_{n+1}{v}_n\right.\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}({\rho }_2{w}_{n+1}{v}_n-({\pi }_2+{\mu }_2)s_{n+1})+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left({\pi }_2{s}_{n+1}-{\alpha }_2{z}_{n+1})\right.\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}({\beta }_1{\alpha }_1{y}_{n+1}+{\beta }_2{\alpha }_2{z}_{n+1}-\theta v_{n+1})+\left(1-\frac{\tilde{q}}{q_{n+1}}\right)\left({\rho }_3{x}_{n+1}{u}_n\right.\hfill \\ \hfill & \left.-({\pi }_3+{\mu }_3)q_{n+1}\right)+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\left(1-\frac{\tilde{u}}{u_{n+1}}\right)\left({\pi }_3{q}_{n+1}-\delta u_{n+1}\right)\hfill \\ \hfill & +\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta \tilde{u}\left(\frac{u_{n+1}}{\tilde{u}}-\frac{u_n}{\tilde{u}}+ln\left(\frac{u_n}{u_{n+1}}\right)\right)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\theta v_{n+1}-\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\theta v_n.\hfill \end{array}$$

By collecting the terms, we obtain

$$\begin{array}{cc}\hfill \Delta {\Phi }_n& \le \left(1-\frac{\tilde{x}}{x_{n+1}}\right)({\zeta }_1-{\varrho }_1{x}_{n+1})+{\rho }_3\tilde{x}{u}_n+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\tilde{w}}{w_{n+1}}\right)\left({\zeta }_2-{\varrho }_2{w}_{n+1}\right)\hfill \\ \hfill & -{\rho }_3{x}_{n+1}{u}_n\frac{\tilde{q}}{q_{n+1}}+({\pi }_3+{\mu }_3)\tilde{q}-({\pi }_3+{\mu }_3)q_{n+1}\frac{\tilde{u}}{u_{n+1}}+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta \tilde{u}-\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta u_n\hfill \\ \hfill & +\left({\rho }_1\tilde{x}+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{\rho }_2\tilde{w}-\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\theta \right)v_n+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta \tilde{u}ln\left(\frac{u_n}{u_{n+1}}\right).\hfill \end{array}$$

Using the equilibria conditions of $E{Q}_2$

$${\zeta }_1={\varrho }_1\tilde{x}+{\rho }_3\tilde{x}\tilde{u},{\zeta }_2={\varrho }_2\tilde{w},{\rho }_3\tilde{x}\tilde{u}=\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta \tilde{u}.$$

We obtain

$$\begin{array}{cc}\hfill \Delta {\Phi }_n& \le -\frac{{\varrho }_1}{x_{n+1}}{(x_{n+1}-\tilde{x})}^2-\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\frac{{\varrho }_2}{w_{n+1}}{(w_{n+1}-\tilde{w})}^2\hfill \\ \hfill & +\left({\rho }_1\tilde{x}+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{\rho }_2\tilde{w}-\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\right)v_n\hfill \\ \hfill & +\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta \tilde{u}\left(3-\frac{\tilde{x}}{x_{n+1}}-\frac{\tilde{q}{x}_{n+1}{u}_n}{\tilde{x}\tilde{u}{q}_{n+1}}-\frac{\tilde{u}{q}_{n+1}}{\tilde{q}{u}_{n+1}}+ln\left(\frac{u_n}{u_{n+1}}\right)\right).\hfill \end{array}$$

(52)

Since we have

$$\begin{array}{cc}\hfill & {\rho }_1\tilde{x}+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{\rho }_2\tilde{w}-\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\hfill \\ \hfill & =\frac{{\rho }_1({\pi }_3+{\mu }_3)\delta }{{\rho }_3{\pi }_3}+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right){\rho }_2{\zeta }_2}{{\beta }_1{\pi }_1{\varrho }_2({\pi }_2+{\mu }_2)}-\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\hfill \\ \hfill & =\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(\frac{{\rho }_1{\beta }_1{\pi }_1\delta ({\pi }_3+{\mu }_3)}{{\rho }_3{\pi }_3\theta \left({\pi }_1+{\mu }_1\right)}+\frac{{\beta }_2{\pi }_2{\rho }_2{\zeta }_2}{{\varrho }_2\theta ({\pi }_2+{\mu }_2)}-1\right)\hfill \\ \hfill & =\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(\frac{R_{01}}{R_1}+R_{02}-1\right),\hfill \end{array}$$

and using the following equality:

$$ln\left(\frac{u_n}{u_{n+1}}\right)=ln\left(\frac{\tilde{x}}{x_{n+1}}\right)+ln\left(\frac{\tilde{q}{x}_{n+1}{u}_n}{\tilde{x}\tilde{u}{q}_{n+1}}\right)+ln\left(\frac{q_{n+1}\tilde{u}}{\tilde{q}{u}_{n+1}}\right).$$

(53)

Then Equation (52) becomes

$$\begin{array}{cc}\hfill \Delta {\Phi }_n& \le -{\varrho }_1\frac{{(x_{n+1}-\tilde{x})}^2}{x_{n+1}}-\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\frac{{\varrho }_2{(w_{n+1}-\tilde{w})}^2}{w_{n+1}}+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(R_{02}+\frac{R_{01}}{R_1}-1\right)v_n\hfill \\ \hfill & -\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta \tilde{u}\left(G\left(\frac{\tilde{x}}{x_{n+1}}\right)+G\left(\frac{\tilde{q}{x}_{n+1}{u}_n}{\tilde{x}\tilde{u}{q}_{n+1}}\right)+G\left(\frac{q_{n+1}\tilde{u}}{\tilde{q}{u}_{n+1}}\right)\right).\hfill \end{array}$$

Since, $R_{02}+\frac{R_{01}}{R_1}\le 1,$ then $\Delta {\Phi }_n\le 0$, for all $n\ge 0.$ Hence, the sequence ${\Phi }_n$ is monotonically decreasing. Since ${\Phi }_n\ge 0,$ then ${lim}_{n\to \infty }{\Phi }_n\ge 0$ and thus ${lim}_{n\to \infty }\Delta {\Phi }_n=0.$ Thus, ${lim}_{n\to \infty }{x}_n=\tilde{x},$ ${lim}_{n\to \infty }{w}_n=\tilde{w},$ ${lim}_{n\to \infty }{q}_n=\tilde{q}$, ${lim}_{n\to \infty }{u}_n=\tilde{u}$ and ${lim}_{n\to \infty }\left(R_{02}+\frac{R_{01}}{R_1}-1\right)v_n=0.$ We have two cases:

(i) $R_{02}+\frac{R_{01}}{R_1}=1,$ and from Equation (3)

$$0={\zeta }_1-{\varrho }_1\tilde{x}-{\rho }_1\tilde{x}\underset{n\to \infty }{lim}{v}_n-{\rho }_3\tilde{x}\tilde{u}⟹\underset{n\to \infty }{lim}{v}_n=0.$$

(54)

Moreover, from Equations (9), (5) and (8),

$$\begin{array}{cc}\hfill 0& ={\beta }_1{\alpha }_1\underset{n\to \infty }{lim}{y}_{n+1}+{\beta }_2{\alpha }_2\underset{n\to \infty }{lim}{z}_{n+1}=0⟹\underset{n\to \infty }{lim}{y}_n=0\mathrm{and}\underset{n\to \infty }{lim}{z}_n=0,\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill 0& ={\pi }_1\underset{n\to \infty }{lim}{g}_{n+1}⟹\underset{n\to \infty }{lim}{g}_n=0,\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill 0& ={\pi }_2\underset{n\to \infty }{lim}{s}_{n+1}⟹\underset{n\to \infty }{lim}{s}_n=0.\hfill \end{array}$$

(57)

(ii) $R_{02}+\frac{R_{01}}{R_1}<1$ and ${lim}_{n\to \infty }{v}_n=0$. Equations (55)–(57) imply that ${lim}_{n\to \infty }{y}_n=0,$ ${lim}_{n\to \infty }{z}_n=0,$ ${lim}_{n\to \infty }{g}_n=0$ and ${lim}_{n\to \infty }{s}_n=0$. This proves that $E{Q}_2$ is GAS. □

Theorem 4.

If $\frac{R_1}{R_{01}}>1$, $R_2>1$ and $R_3>1,$ then $E{Q}_3=(\overline{x},\overline{g},\overline{y},\overline{w},\overline{s},\overline{z},\overline{v},\overline{q},\overline{u})$ is GAS in the interior of ${\Gamma }_1$.

Proof. 

Consider

$$\begin{array}{cc}\hfill {\Psi }_n& =\frac{1}{{\rm Y}\left(h\right)}\left[\overline{x}G\left(\frac{x_n}{\overline{x}}\right)+\overline{g}G\left(\frac{g_n}{\overline{g}}\right)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}\overline{y}G\left(\frac{y_n}{\overline{y}}\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\overline{w}G\left(\frac{w_n}{\overline{w}}\right)\right.\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\overline{s}G\left(\frac{s_n}{\overline{s}}\right)+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{z}G\left(\frac{z_n}{\overline{z}}\right)+\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}(1+{\rm Y}\left(h\right)\theta )\overline{v}G\left(\frac{v_n}{\overline{v}}\right)\hfill \\ \hfill & \left.+\tilde{q}G\left(\frac{q_n}{\tilde{q}}\right)+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}(1+{\rm Y}\left(h\right)\delta )\tilde{u}G\left(\frac{u_n}{\tilde{u}}\right)\right].\hfill \end{array}$$

Computing the difference $\Delta {\Psi }_n={\Psi }_{n+1}-{\Psi }_n$ as:

$$\begin{array}{cc}\hfill \Delta {\Psi }_n& =\frac{1}{{\rm Y}\left(h\right)}\left[\overline{x}\left(\frac{x_{n+1}}{\overline{x}}-\frac{x_n}{\overline{x}}+ln\left(\frac{x_n}{x_{n+1}}\right)\right)+\overline{g}\left(\frac{g_{n+1}}{\overline{g}}-\frac{g_n}{\overline{g}}+ln\left(\frac{g_n}{g_{n+1}}\right)\right)\right.\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}\overline{y}\left(\frac{y_{n+1}}{\overline{y}}-\frac{y_n}{\overline{y}}+ln\left(\frac{y_n}{y_{n+1}}\right)\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\overline{w}\left(\frac{w_{n+1}}{\overline{w}}-\frac{w_n}{\overline{w}}+ln\left(\frac{w_n}{w_{n+1}}\right)\right)\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\overline{s}\left(\frac{s_{n+1}}{\overline{s}}-\frac{s_n}{\overline{s}}+ln\left(\frac{s_n}{s_{n+1}}\right)\right)+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{z}\left(\frac{z_{n+1}}{\overline{z}}-\frac{z_n}{\overline{z}}+ln\left(\frac{z_n}{z_{n+1}}\right)\right)\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{v}\left(\frac{v_{n+1}}{\overline{v}}-\frac{v_n}{\overline{v}}+ln\left(\frac{v_n}{v_{n+1}}\right)\right)+\overline{q}\left(\frac{q_{n+1}}{\overline{q}}-\frac{q_n}{\overline{q}}+ln\left(\frac{q_n}{q_{n+1}}\right)\right)\hfill \\ \hfill & \left.+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\overline{u}\left(\frac{u_{n+1}}{\overline{u}}-\frac{u_n}{\overline{u}}+ln\left(\frac{u_n}{u_{n+1}}\right)\right)\right]+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{v}\left(G\left(\frac{v_{n+1}}{\overline{v}}\right)-G\left(\frac{v_n}{\overline{v}}\right)\right)\hfill \\ \hfill & +\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\overline{u}\left(G\left(\frac{u_{n+1}}{\overline{u}}\right)-G\left(\frac{u_n}{\overline{u}}\right)\right).\hfill \end{array}$$

Using inequality (39), we obtain

$$\begin{array}{cc}\hfill \Delta {\Psi }_n& \le \frac{1}{{\rm Y}\left(h\right)}\left[\overline{x}\left(\frac{x_{n+1}-x_n}{\overline{x}}+\frac{x_n}{x_{n+1}}-1\right)+\overline{g}\left(\frac{g_{n+1}-g_n}{\overline{g}}+\frac{g_n}{g_{n+1}}-1\right)\right.\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}\overline{y}\left(\frac{y_{n+1}-y_n}{\overline{y}}+\frac{y_n}{y_{n+1}}-1\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\overline{w}\left(\frac{w_{n+1}-w_n}{\overline{w}}+\frac{w_n}{w_{n+1}}-1\right)\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\overline{s}\left(\frac{s_{n+1}-s_n}{\overline{s}}+\frac{s_n}{s_{n+1}}-1\right)+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{z}\left(\frac{z_{n+1}-z_n}{\overline{z}}+\frac{z_n}{z_{n+1}}-1\right)\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{v}\left(\frac{v_{n+1}-v_n}{\overline{v}}+\frac{v_n}{v_{n+1}}-1\right)+\overline{q}\left(\frac{q_{n+1}-q_n}{\overline{q}}+\frac{q_n}{q_{n+1}}-1\right)\hfill \\ \hfill & \left.+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\overline{u}\left(\frac{u_{n+1}-u_n}{\overline{u}}+\frac{u_n}{u_{n+1}}-1\right)\right]+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{v}\left(G\left(\frac{v_{n+1}}{\overline{v}}\right)-G\left(\frac{v_n}{\overline{v}}\right)\right)\hfill \\ \hfill & +\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\overline{u}\left(G\left(\frac{u_{n+1}}{\overline{u}}\right)-G\left(\frac{u_n}{\overline{u}}\right)\right).\hfill \end{array}$$

(58)

We write inequality (58) as:

$$\begin{array}{cc}\hfill \Delta {\Psi }_n& \le \frac{1}{{\rm Y}\left(h\right)}\left[\left(1-\frac{\overline{x}}{x_{n+1}}\right)\left(x_{n+1}-x_n\right)+\left(1-\frac{\overline{g}}{g_{n+1}}\right)\left(g_{n+1}-g_n\right)\right.\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}\left(1-\frac{\overline{y}}{y_{n+1}}\right)\left(y_{n+1}-y_n\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\overline{w}}{w_{n+1}}\right)\left(w_{n+1}-w_n\right)\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\overline{s}}{s_{n+1}}\right)\left(s_{n+1}-s_n\right)+\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(1-\frac{\overline{z}}{z_{n+1}}\right)\left(z_{n+1}-z_n\right)\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(1-\frac{\overline{v}}{v_{n+1}}\right)\left(v_{n+1}-v_n\right)+\left(1-\frac{\overline{q}}{q_{n+1}}\right)\left(q_{n+1}-q_n\right)\hfill \\ \hfill & \left.+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\left(1-\frac{\overline{u}}{u_{n+1}}\right)\left(u_{n+1}-u_n\right)\right]+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{v}\left(\frac{v_{n+1}}{\overline{v}}-\frac{v_n}{\overline{v}}+ln\left(\frac{v_n}{v_{n+1}}\right)\right)\hfill \\ \hfill & +\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\overline{u}\left(\frac{u_{n+1}}{\overline{u}}-\frac{u_n}{\overline{u}}+ln\left(\frac{u_n}{u_{n+1}}\right)\right).\hfill \end{array}$$

From Equations (3)–(11), we have

$$\begin{array}{cc}\hfill \Delta {\Psi }_n& \le \left(1-\frac{\overline{x}}{x_{n+1}}\right)({\zeta }_1-{\varrho }_1{x}_{n+1}-{\rho }_1{x}_{n+1}{v}_n-{\rho }_3{x}_{n+1}{u}_n)+\left(1-\frac{\overline{g}}{g_{n+1}}\right)\hfill \\ \hfill & \times ({\rho }_1{x}_{n+1}{v}_n-({\pi }_1+{\mu }_1)g_{n+1})+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}\left(1-\frac{\overline{y}}{y_{n+1}}\right)({\pi }_1{g}_{n+1}-{\alpha }_1{y}_{n+1})\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\overline{w}}{w_{n+1}}\right)\left({\zeta }_2-{\varrho }_2{w}_{n+1}-{\rho }_2{w}_{n+1}{v}_n\right)\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\overline{s}}{s_{n+1}}\right)({\rho }_2{w}_{n+1}{v}_n-({\pi }_2+{\mu }_2)s_{n+1})\hfill \\ \hfill & +\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(1-\frac{\overline{z}}{z_{n+1}}\right)({\pi }_2{s}_{n+1}-{\alpha }_2{z}_{n+1})\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\left(1-\frac{\overline{v}}{v_{n+1}}\right)({\beta }_1{\alpha }_1{y}_{n+1}+{\beta }_2{\alpha }_2{z}_{n+1}-\theta v_{n+1})\hfill \\ \hfill & +\left(1-\frac{\overline{q}}{q_{n+1}}\right)({\rho }_3{x}_{n+1}{u}_n-({\pi }_3+{\mu }_3)q_{n+1})+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\left(1-\frac{\overline{u}}{u_{n+1}}\right)\left({\pi }_3{q}_{n+1}-\delta u_{n+1}\right)\hfill \\ \hfill & +\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{v}_{n+1}-\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{v}_n+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{v}ln\left(\frac{v_n}{v_{n+1}}\right)+\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}{u}_{n+1}\hfill \\ \hfill & -\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}{u}_n+\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\overline{u}ln\left(\frac{u_n}{u_{n+1}}\right).\hfill \end{array}$$

After collecting the terms, we obtain

$$\begin{array}{cc}\hfill \Delta {\Psi }_n& \le \left(1-\frac{\overline{x}}{x_{n+1}}\right)({\zeta }_1-{\varrho }_1{x}_{n+1})+{\rho }_1\overline{x}{v}_n+{\rho }_3\overline{x}{u}_n-{\rho }_1{x}_{n+1}{v}_n\frac{\overline{g}}{g_{n+1}}+({\pi }_1+{\mu }_1)\overline{g}\hfill \\ \hfill & -\left({\pi }_1+{\mu }_1\right)g_{n+1}\frac{\overline{y}}{y_{n+1}}+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}{\alpha }_1\overline{y}+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\overline{w}}{w_{n+1}}\right)\left({\zeta }_2-{\varrho }_2{w}_{n+1}\right)\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{\rho }_2\overline{w}{v}_n-\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{\rho }_2{w}_{n+1}{v}_n\frac{\overline{s}}{s_{n+1}}+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{s}\hfill \\ \hfill & -\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{s}_{n+1}\frac{\overline{z}}{z_{n+1}}+\frac{{\beta }_2{\alpha }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{z}-\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}{\alpha }_1{y}_{n+1}\frac{\overline{v}}{v_{n+1}}\hfill \\ \hfill & -\frac{{\beta }_2{\alpha }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\frac{z_{n+1}\overline{v}}{v_{n+1}}+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{v}-{\rho }_3{x}_{n+1}{u}_n\frac{\overline{q}}{q_{n+1}}\hfill \\ \hfill & +({\pi }_3+{\mu }_3)\overline{q}-({\pi }_3+{\mu }_3)\frac{q_{n+1}\overline{u}}{u_{n+1}}+\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\overline{u}-\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{v}_n\hfill \\ \hfill & +\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{v}ln\left(\frac{v_n}{v_{n+1}}\right)-\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}{u}_n+\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\overline{u}ln\left(\frac{u_n}{u_{n+1}}\right).\hfill \end{array}$$

Using the equilibrium conditions for $E{Q}_3$

$$\begin{array}{cc}\hfill {\rho }_1\overline{x}\overline{v}& =\left({\pi }_1+{\mu }_1\right)\overline{g},{\rho }_3\overline{x}\overline{u}=({\pi }_3+{\mu }_3)\overline{q},\hfill \\ \hfill {\rho }_2\overline{w}\overline{v}& =({\pi }_2+{\mu }_2)\overline{s},\theta \overline{v}={\beta }_1{\alpha }_1\overline{y}+{\beta }_2{\alpha }_2\overline{z},\hfill \\ \hfill {\zeta }_1& ={\varrho }_1\overline{x}+{\rho }_1\overline{x}\overline{v}+{\rho }_3\overline{x}\overline{u},{\zeta }_2={\varrho }_2\overline{w}+{\rho }_2\overline{w}\overline{v},\hfill \\ \hfill \overline{q}& =\frac{\delta \overline{u}}{{\pi }_3},\overline{g}=\frac{{\alpha }_1\overline{y}}{{\pi }_1},\overline{s}=\frac{{\alpha }_2\overline{z}}{{\pi }_2},\hfill \end{array}$$

we obtain

$${\zeta }_1={\varrho }_1\overline{x}+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}{\alpha }_1\overline{y}+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta \overline{u},{\zeta }_2={\varrho }_2\overline{w}+\frac{({\pi }_2+{\mu }_2)}{{\pi }_2}{\alpha }_2\overline{z}.$$

Then

$$\begin{array}{cc}\hfill \Delta {\Psi }_n& \le \left(1-\frac{\overline{x}}{x_{n+1}}\right)\left({\varrho }_1\overline{x}+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}{\alpha }_1\overline{y}+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta \overline{u}-{\varrho }_1{x}_{n+1}\right)+{\rho }_1\overline{x}{v}_n\hfill \\ \hfill & +{\rho }_3\overline{x}{u}_n-{\rho }_1\overline{x}\overline{v}\frac{x_{n+1}{v}_n\overline{g}}{\overline{x}\overline{v}{g}_{n+1}}+({\pi }_1+{\mu }_1)\overline{g}-\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}{\alpha }_1\overline{y}\frac{g_{n+1}\overline{y}}{\overline{g}{y}_{n+1}}+\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}{\alpha }_1\overline{y}\hfill \\ \hfill & +\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\left(1-\frac{\overline{w}}{w_{n+1}}\right)\left({\varrho }_2\overline{w}+\frac{({\pi }_2+{\mu }_2)}{{\pi }_2}{\alpha }_2\overline{z}-{\varrho }_2{w}_{n+1}\right)+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{\rho }_2\overline{w}{v}_n\hfill \\ \hfill & -\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}{\rho }_2\overline{w}\overline{v}\frac{w_{n+1}{v}_n\overline{s}}{\overline{w}\overline{v}{s}_{n+1}}+\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{s}-\frac{{\beta }_2{\pi }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{s}\frac{\overline{z}{s}_{n+1}}{\overline{s}{z}_{n+1}}\hfill \\ \hfill & +\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\overline{z}-\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}{\alpha }_1\overline{y}\frac{\overline{v}{y}_{n+1}}{\overline{y}{v}_{n+1}}-\frac{{\beta }_2{\alpha }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{z}\frac{z_{n+1}\overline{v}}{\overline{z}{v}_{n+1}}\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}{\alpha }_1\overline{y}+\frac{{\beta }_2{\alpha }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{z}-{\rho }_3\overline{x}\overline{u}\frac{x_{n+1}{u}_n\overline{q}}{\overline{x}\overline{u}{q}_{n+1}}+({\pi }_3+{\mu }_3)\overline{q}\hfill \\ \hfill & -({\pi }_3+{\mu }_3)\overline{q}\frac{q_{n+1}\overline{u}}{\overline{q}{u}_{n+1}}+\frac{({\pi }_3+{\mu }_3)}{{\pi }_3}\delta \overline{u}-\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{v}_n+\frac{\theta \left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}\overline{v}ln\left(\frac{v_n}{v_{n+1}}\right)\hfill \\ \hfill & -\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}{u}_n+\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\overline{u}ln\left(\frac{u_n}{u_{n+1}}\right),\hfill \end{array}$$

and we obtain

$$\begin{array}{cc}\hfill \Delta {\Psi }_n& \le -{\varrho }_1\frac{{(x_{n+1}-\overline{x})}^2}{x_{n+1}}-\frac{{\beta }_2{\pi }_2{\varrho }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\frac{{(w_{n+1}-\overline{w})}^2}{w_{n+1}}\hfill \\ \hfill & +\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}{\alpha }_1\overline{y}\left[4-\frac{\overline{x}}{x_{n+1}}-\frac{x_{n+1}{v}_n\overline{g}}{\overline{x}\overline{v}{g}_{n+1}}-\frac{g_{n+1}\overline{y}}{\overline{g}{y}_{n+1}}+\frac{y_{n+1}\overline{v}}{\overline{y}{v}_{n+1}}+ln\left(\frac{v_n}{v_{n+1}}\right)\right]\hfill \\ \hfill & +\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\overline{z}\left[4-\frac{\overline{w}}{w_{n+1}}-\frac{w_{n+1}{v}_n\overline{s}}{\overline{w}\overline{v}{s}_{n+1}}-\frac{\overline{z}{s}_{n+1}}{\overline{s}{z}_{n+1}}+\frac{z_{n+1}\overline{v}}{\overline{z}{v}_{n+1}}+ln\left(\frac{v_n}{v_{n+1}}\right)\right]\hfill \\ \hfill & +\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\overline{u}\left[3-\frac{\overline{x}}{x_{n+1}}-\frac{x_{n+1}{u}_n\overline{q}}{\overline{x}\overline{u}{q}_{n+1}}-\frac{q_{n+1}\overline{u}}{\overline{q}{u}_{n+1}}+ln\left(\frac{u_n}{u_{n+1}}\right)\right].\hfill \end{array}$$

Using equalities similar to Equations (48), (49) and (53), we obtain

$$\begin{array}{cc}\hfill \Delta {\Psi }_n& \le -{\varrho }_1\frac{{(x_{n+1}-\overline{x})}^2}{x_{n+1}}-\frac{{\beta }_2{\pi }_2{\varrho }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1({\pi }_2+{\mu }_2)}\frac{{(w_{n+1}-\overline{w})}^2}{w_{n+1}}\hfill \\ \hfill & -\frac{\left({\pi }_1+{\mu }_1\right)}{{\pi }_1}{\alpha }_1\overline{y}\left[G\left(\frac{\overline{x}}{x_{n+1}}\right)+G\left(\frac{x_{n+1}{v}_n\overline{g}}{\overline{x}\overline{v}{g}_{n+1}}\right)+G\left(\frac{g_{n+1}\overline{y}}{\overline{g}{y}_{n+1}}\right)+G\left(\frac{y_{n+1}\overline{v}}{\overline{y}{v}_{n+1}}\right)\right]\hfill \\ \hfill & -\frac{{\beta }_2\left({\pi }_1+{\mu }_1\right)}{{\beta }_1{\pi }_1}{\alpha }_2\overline{z}\left[G\left(\frac{\overline{w}}{w_{n+1}}\right)+G\left(\frac{w_{n+1}{v}_n\overline{s}}{\overline{w}\overline{v}{s}_{n+1}}\right)+G\left(\frac{\overline{z}{s}_{n+1}}{\overline{s}{z}_{n+1}}\right)+G\left(\frac{z_{n+1}\overline{v}}{\overline{z}{v}_{n+1}}\right)\right]\hfill \\ \hfill & -\frac{\delta ({\pi }_3+{\mu }_3)}{{\pi }_3}\overline{u}\left[G\left(\frac{\overline{x}}{x_{n+1}}\right)+G\left(\frac{x_{n+1}{u}_n\overline{q}}{\overline{x}\overline{u}{q}_{n+1}}\right)+G\left(\frac{q_{n+1}\overline{u}}{\overline{q}{u}_{n+1}}\right)\right].\hfill \end{array}$$

We note that $\Delta {\Psi }_n\le 0.$ Hence, the sequence ${\Psi }_n$ is monotonically decreasing. Since ${\Psi }_n\ge 0,$ then ${lim}_{n\to \infty }{\Psi }_n\ge 0$ and thus, ${lim}_{n\to \infty }\Delta {\Psi }_n=0.$ Thus, ${lim}_{n\to \infty }(x_n,g_n,y_n,w_n,s_n,z_n,v_n,q_n,u_n)$ $=(\overline{x},\overline{g},\overline{y},\overline{w},\overline{s},\overline{z},\overline{v},\overline{q},\overline{u}).$ Hence, $E{Q}_3$ is GAS. □

6. Numerical Simulations

In this section, we execute numerical simulations for the discrete-time model (3)–(11). Moreover, we discussed the impact of latent reservoirs on the HIV-1 and HTLV-I co-dynamics. We use the values of the parameters given in

Table 1. Model parameters.

Parameter	Value	Source
${\zeta }_1$	10 cells mm${}^{-3}$ day${}^{-1}$	[27,35,54]
${\zeta }_2$	$0.03198$ cells mm${}^{-3}$ day${}^{-1}$	[28,29]
${\alpha }_1$	$0.5$ day${}^{-1}$	[20,22]
${\alpha }_2$	$0.1$ day${}^{-1}$	[13,55]
${\gamma }_1$	$0.01$ day${}^{-1}$	[27,56]
${\gamma }_2$	$0.01$ day${}^{-1}$	[28,29]
${\beta }_1$	6 viruses cells${}^{-1}$	[55]
${\beta }_2$	6 viruses cells${}^{-1}$	[55]
$\theta$	2 day${}^{-1}$	[13,57]
$\delta$	$0.2$ day${}^{-1}$	[13,34]
h	$0.1$	[58]
${\rho }_1$	(varied) viruses${}^{-1}$ mm${}^3$ day${}^{-1}$	Assumed
${\rho }_2$	(varied) viruses${}^{-1}$ mm${}^3$ day${}^{-1}$	Assumed
${\rho }_3$	(varied) cells${}^{-1}$ mm${}^3$ day${}^{-1}$	Assumed
${\mu }_1$	$0.02$ day${}^{-1}$	[54,59]
${\mu }_2$	$0.01$ day${}^{-1}$	[55]
${\mu }_3$	$0.01$ day${}^{-1}$	[35]
${\pi }_1$	$0.2$ day${}^{-1}$	[59]
${\pi }_2$	$0.01$ day${}^{-1}$	Assumed
${\pi }_3$	$0.03$ day${}^{-1}$	[60]

. Some of these values are taken from previous works for HIV-1 and HTLV-I single-infections. The values of other parameters are assumed just to execute the numerical simulations. The reason for that is the unavailability of real data from HIV-1 and HTLV-I coinfection patients. However, when the real data are available, then the values of the parameters of the coinfection model can be estimated.

6.1. Stability of the Equilibria

To support the global stability results given Theorems 1–4, we show that the solutions of the discrete-time system starting from any point (any disease stage) in the feasible region will tend to one of the four equilibria. Let us consider three different initial values as:

$$\begin{array}{cc}\hfill \mathrm{IV}1:& (x_0,g_0,y_0,w_0,s_0,z_0,v_0,q_0,u_0)=(850,25,5.5,2,1.5,0.1,20,55,35),\\ \hfill \mathrm{IV}2:& (x_0,g_0,y_0,w_0,s_0,z_0,v_0,q_0,u_0)=(650,20,3.5,1.5,1,0.15,15,45,25),\\ \hfill \mathrm{IV}3:& (x_0,g_0,y_0,w_0,s_0,z_0,v_0,q_0,u_0)=(350,15,2,1,0.5,0.2,0.4,35,15).\end{array}$$

We choose ${\rho }_1$, ${\rho }_2$ and ${\rho }_3$ as:

Case (I) ${\rho }_1=0.0002,$ ${\rho }_2=0.001$ and ${\rho }_3=0.0001$. This gives $R_0=0.550252\le 1$ and $R_1=0.375\le 1$. Figure 1 illustrates that the concentrations of uninfected CD4${}^{+}$T cells and uninfected macrophages increase and tend to the healthy values $x^0=1000$ and $w^0=3.1980$, while the concentrations of other populations decrease and converge to zero. Therefore, $E{Q}_0$ is GAS, and this agrees with the result of Theorem 1. In this case, both HIV-1 and HTLV-I are cleared from the human body, regardless of the starting states.

Case (II) ${\rho }_1=0.0007,$ ${\rho }_2=0.001$ and ${\rho }_3=0.0001$. These values give $R_0=1.91389>1$ and $R_1=0.375\le 1$. From Figure 2, we see that the solutions of the discrete-time model tend to the equilibrium $E{Q}_1=(522.71,21.695,8.678,1.388,0.905,0.091,13.044,0,0)$. As a result, $E{Q}_1$ exists, and based on Theorem 2, it is GAS. This result shows that the HIV-1 single-infection can be reached for all initial states.

Case (III) ${\rho }_1=0.0003,$ ${\rho }_2=0.0001$ and ${\rho }_3=0.00045$, and then $R_1=1.6875>1$ and $R_{02}+(R_{01}/R_1)=0.489645\le 1$. Figure 3 demonstrates that the solutions of the discrete-time model reach the equilibrium $E{Q}_2=(592.593,0,0,3.198,0,0,0,101.852,15.278)$ for all the initial states. According to Lemma 2 and Theorem 3, $E{Q}_2$ exists and it is GAS. This result shows that the HTLV-I single-infection can be reached for all starting states.

Case (IV) ${\rho }_1=0.00054,$ ${\rho }_2=0.03$ and ${\rho }_3=0.0004$, and thus, $R_1/R_{01}=1.01852>1,$ $R_2=7.91505>1,$ and $R_3=4.017>1$. Figure 4 illustrates that the solutions of the discrete-time model starting with initial values IV1-IV3 converge to the equilibrium $E{Q}_3=(666.667,3.772,1.509,0.404,1.397,0.1396,2.305,62.588,9.388)$. Based on Lemma 2 and Theorem 4, $E{Q}_3$ exists and it is GAS. This result shows that the HIV-1 and HTLV-I coinfection can be reached for all starting states.

For more confirmation, we examine the local stability of the equilibria of the discrete-time model in Cases (I)–(IV). The Jacobian matrix $J=J(x,g,y,w,s,z,v,q,u)$ of model (9)–(25) is calculated as:

$$J=\left(\begin{array}{ccccccccc}{J}_{11}& 0& 0& 0& 0& 0& J_{17}& 0& J_{19}\\ J_{21}& J_{22}& 0& 0& 0& 0& J_{27}& 0& J_{29}\\ J_{31}& J_{32}& J_{33}& 0& 0& 0& J_{37}& 0& J_{39}\\ 0& 0& 0& J_{44}& 0& 0& J_{47}& 0& 0\\ 0& 0& 0& J_{54}& J_{55}& 0& J_{57}& 0& 0\\ 0& 0& 0& J_{64}& J_{65}& J_{66}& J_{67}& 0& 0\\ J_{71}& J_{72}& J_{73}& J_{74}& J_{75}& J_{76}& J_{77}& 0& J_{79}\\ J_{81}& 0& 0& 0& 0& 0& J_{87}& J_{88}& J_{89}\\ J_{91}& 0& 0& 0& 0& 0& J_{97}& J_{98}& J_{99}\end{array}\right),$$

where

$$\begin{array}{cc}\hfill J_{11}& =\frac{1}{1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u)},\hfill \\ \hfill J_{17}& =-\frac{{\rho }_1{\rm Y}\left(h\right)(x+{\zeta }_1{\rm Y}\left(h\right))}{{(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))}^2},\hfill \\ \hfill J_{19}& =-\frac{{\rm Y}\left(h\right)(x+{\zeta }_1{\rm Y}\left(h\right)){\rho }_3}{{(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))}^2},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{21}& =\frac{v{\rm Y}\left(h\right){\rho }_1}{(1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1))(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))},\hfill \\ \hfill J_{22}& =\frac{1}{1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1)},\hfill \\ \hfill J_{27}& =\frac{{\rho }_1{\rm Y}\left(h\right)(x+{\zeta }_1{\rm Y}\left(h\right))(1+{\rm Y}\left(h\right){\varrho }_1+{\rm Y}\left(h\right)u{\rho }_3)}{(1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1)){(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))}^2},\hfill \\ \hfill J_{29}& =-\frac{v{{\rm Y}}^2\left(h\right)(x+{\zeta }_1{\rm Y}\left(h\right)){\rho }_1{\rho }_3}{(1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1)){(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))}^2},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{31}& =\frac{v{\pi }_1{{\rm Y}}^2\left(h\right){\rho }_1}{(1+{\alpha }_1{\rm Y}\left(h\right))(1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1))(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))},\hfill \\ \hfill J_{32}& =\frac{{\pi }_1{\rm Y}\left(h\right)}{(1+{\alpha }_1{\rm Y}\left(h\right))(1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1))},\hfill \\ \hfill J_{33}& =\frac{1}{1+{\alpha }_1{\rm Y}\left(h\right)},\hfill \\ \hfill J_{37}& =\frac{{\pi }_1{{\rm Y}}^2\left(h\right)(x+{\zeta }_1{\rm Y}\left(h\right)){\rho }_1(1+{\rm Y}\left(h\right){\varrho }_1+{\rm Y}\left(h\right){\rho }_{3}u)}{(1+{\alpha }_1{\rm Y}\left(h\right))(1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1)){(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))}^2},\hfill \\ \hfill J_{39}& =-\frac{v{\pi }_1{{\rm Y}}^3\left(h\right)(x+{\zeta }_1{\rm Y}\left(h\right)){\rho }_1{\rho }_3}{(1+{\alpha }_1{\rm Y}\left(h\right))(1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1)){(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))}^2},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{44}& =\frac{1}{1+{\rm Y}\left(h\right){\varrho }_2+{\rm Y}\left(h\right){\rho }_{2}v},\hfill \\ \hfill J_{47}& =-\frac{{\rm Y}\left(h\right)(w+{\rm Y}\left(h\right){\zeta }_2){\rho }_2}{{(1+{\rm Y}\left(h\right){\varrho }_2+{\rm Y}\left(h\right){\rho }_{2}v)}^2},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{54}& =\frac{{\rm Y}\left(h\right){\rho }_{2}v}{(1+{\rm Y}\left(h\right)({\pi }_2+{\mu }_2))(1+{\rm Y}\left(h\right){\varrho }_2+{\rm Y}\left(h\right){\rho }_{2}v)},\hfill \\ \hfill J_{55}& =\frac{1}{1+{\rm Y}\left(h\right)({\pi }_2+{\mu }_2)},\hfill \\ \hfill J_{57}& =\frac{{\rm Y}\left(h\right)(1+{\rm Y}\left(h\right){\varrho }_2)(w+{\zeta }_2{\rm Y}\left(h\right)){\rho }_2}{(1+{\rm Y}\left(h\right)({\pi }_2+{\mu }_2)){(1+{\rm Y}\left(h\right){\varrho }_2+{\rm Y}\left(h\right){\rho }_{2}v)}^2},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{64}& =\frac{{{\rm Y}}^2\left(h\right)v{\pi }_2{\rho }_2}{(1+{\rm Y}\left(h\right){\alpha }_2)(1+{\rm Y}\left(h\right)({\pi }_2+{\mu }_2))(1+{\rm Y}\left(h\right){\varrho }_2+{\rm Y}\left(h\right){\rho }_{2}v)},\hfill \\ \hfill J_{65}& =\frac{{\rm Y}\left(h\right){\pi }_2}{(1+{\rm Y}\left(h\right){\alpha }_2)(1+{\rm Y}\left(h\right)({\pi }_2+{\mu }_2))},\hfill \\ \hfill J_{66}& =\frac{1}{(1+{\rm Y}\left(h\right){\alpha }_2)},\hfill \\ \hfill J_{67}& =\frac{{{\rm Y}}^2\left(h\right){\pi }_2(w+{\rm Y}\left(h\right){\zeta }_2)({\rho }_2+{\rm Y}\left(h\right){\varrho }_2{\rho }_2)}{(1+{\rm Y}\left(h\right){\alpha }_2)(1+{\rm Y}\left(h\right)({\pi }_2+{\mu }_2)){(1+{\rm Y}\left(h\right){\varrho }_2+{\rm Y}\left(h\right){\rho }_{2}v)}^2},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{71}& =\frac{{{\rm Y}}^3\left(h\right)v{\alpha }_1{\beta }_1{\pi }_1{\rho }_1}{(1+{\rm Y}\left(h\right){\alpha }_1)(1+{\rm Y}\left(h\right)\theta )(1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1))(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))},\hfill \\ \hfill J_{72}& =\frac{{{\rm Y}}^2\left(h\right){\alpha }_1{\beta }_1{\pi }_1}{(1+{\rm Y}\left(h\right){\alpha }_1)(1+{\rm Y}\left(h\right)\theta )(1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1))},\hfill \\ \hfill J_{73}& =\frac{{\rm Y}\left(h\right){\alpha }_1{\beta }_1}{(1+{\rm Y}\left(h\right){\alpha }_1)(1+{\rm Y}\left(h\right)\theta )},\hfill \\ \hfill J_{74}& =\frac{{{\rm Y}}^3\left(h\right)v{\alpha }_2{\beta }_2{\pi }_2{\rho }_2}{(1+{\rm Y}\left(h\right){\alpha }_2)(1+{\rm Y}\left(h\right)\theta )(1+{\rm Y}\left(h\right)({\pi }_2+{\mu }_2))(1+{\rm Y}\left(h\right){\varrho }_2+{\rm Y}\left(h\right){\rho }_{2}v)},\hfill \\ \hfill J_{75}& =\frac{{{\rm Y}}^2\left(h\right){\alpha }_2{\beta }_2{\pi }_2}{(1+{\rm Y}\left(h\right){\alpha }_2)(1+{\rm Y}\left(h\right)\theta )(1+{\rm Y}\left(h\right)({\pi }_2+{\mu }_2))},\hfill \\ \hfill J_{76}& =\frac{{\rm Y}\left(h\right){\alpha }_2{\beta }_2}{(1+{\rm Y}\left(h\right){\alpha }_2)(1+{\rm Y}\left(h\right)\theta )},\hfill \\ \hfill J_{77}& =\frac{1}{1+{\rm Y}\left(h\right)\theta }+\frac{{{\rm Y}}^3\left(h\right)}{1+{\rm Y}\left(h\right)\theta }\left[\frac{{\alpha }_2{\beta }_2{\pi }_2{\rho }_2(w+{\rm Y}\left(h\right){\zeta }_2)(1+{\rm Y}\left(h\right){\varrho }_2)}{(1+{\rm Y}\left(h\right){\alpha }_2)(1+{\rm Y}\left(h\right)({\pi }_2+{\mu }_2)){(1+{\rm Y}\left(h\right){\varrho }_2+{\rm Y}\left(h\right){\rho }_{2}v)}^2}\right.\hfill \\ \hfill & \left.+\frac{{\alpha }_1{\beta }_1{\pi }_1(x+{\rm Y}\left(h\right){\zeta }_1){\rho }_1(1+{\rm Y}\left(h\right){\varrho }_1+{\rm Y}\left(h\right){\rho }_{3}u)}{(1+{\rm Y}\left(h\right){\alpha }_1)(1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1)){(1+{\rm Y}\left(h\right)({\varrho }_1+v{\rho }_1+{\rm Y}\left(h\right){\rho }_{2}u))}^2}\right],\hfill \\ \hfill J_{79}& =-\frac{{{\rm Y}}^4\left(h\right)v{\alpha }_1{\beta }_1{\pi }_1(x+{\rm Y}\left(h\right){\zeta }_1){\rho }_1{\rho }_3}{(1+{\rm Y}\left(h\right){\alpha }_1)(1+{\rm Y}\left(h\right)\theta )(1+{\rm Y}\left(h\right)({\pi }_1+{\mu }_1)){(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))}^2},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{81}& =\frac{{\rm Y}\left(h\right)u{\rho }_3}{(1+{\rm Y}\left(h\right)({\pi }_3+{\mu }_3))(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))},\hfill \\ \hfill J_{87}& =-\frac{{{\rm Y}}^2\left(h\right)(x+{\rm Y}\left(h\right){\zeta }_1){\rho }_1{\rho }_3}{(1+{\rm Y}\left(h\right)({\pi }_3+{\mu }_3)){(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))}^2},\hfill \\ \hfill J_{88}& =\frac{1}{1+{\rm Y}\left(h\right)({\pi }_3+{\mu }_3)},\hfill \\ \hfill J_{89}& =\frac{{\rho }_3{\rm Y}\left(h\right)(x+{\rm Y}\left(h\right){\zeta }_1)(1+{\rm Y}\left(h\right){\varrho }_1+{\rm Y}\left(h\right){\rho }_{1}v)}{(1+{\rm Y}\left(h\right)({\pi }_3+{\mu }_3)){(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))}^2},\hfill \end{array}$$

$$\begin{array}{cc}\hfill J_{91}& =\frac{{{\rm Y}}^2\left(h\right)u{\pi }_3{\rho }_3}{(1+{\rm Y}\left(h\right)\delta )(1+{\rm Y}\left(h\right)({\pi }_3+{\mu }_3))(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))},\hfill \\ \hfill J_{97}& =-\frac{{{\rm Y}}^3\left(h\right)u{\pi }_3(x+{\zeta }_1{\rm Y}\left(h\right)){\rho }_1{\rho }_3}{(1+{\rm Y}\left(h\right)\delta )(1+{\rm Y}\left(h\right)({\pi }_3+{\mu }_3)){(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))}^2},\hfill \\ \hfill J_{98}& =\frac{{\rm Y}\left(h\right){\pi }_3}{(1+{\rm Y}\left(h\right)\delta )(1+{\rm Y}\left(h\right)({\pi }_3+{\mu }_3))},\hfill \\ \hfill J_{99}& =\frac{1}{1+{\rm Y}\left(h\right)\delta }\left[1+\frac{{{\rm Y}}^2\left(h\right){\pi }_3(x+{\zeta }_1{\rm Y}\left(h\right))(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v)){\rho }_3}{(1+{\rm Y}\left(h\right)({\pi }_2+{\mu }_2)){(1+{\rm Y}\left(h\right)({\varrho }_1+{\rho }_{1}v+{\rho }_{3}u))}^2}\right],\hfill \end{array}$$

Then, we compute the eigenvalues ${\lambda }_j,$ $j=1,2,\dots ,9$ of the matrix $J,$ at each equilibrium. An equilibrium point of the discrete-time model is locally asymptotically stable (LAS) when $\left|{\lambda }_j\right|<1$, for all $j=1,2,\dots ,9$. We compute the eigenvalues of all nonnegative equilibria using the values of ${\rho }_1$, ${\rho }_2$ and ${\rho }_3$ given in Cases (I)–(IV).

Table 2. Local stability of equilibria.

Case	Equilibrium Point	$\left|{\mathit{\lambda }}_{\mathit{j}}\right|,\mathit{j}=1,2,\mathit{\dots },9$	Stability
(I)	$\begin{array}{c}E{Q}_0=\left(1000,0,0,3.20,0,0,0,0,0\right)\end{array}$	$\begin{array}{c}\left(0.999,0.999,0.998,0.998,0.993,0.990,0.979,0.934,0.837\right)\end{array}$	$\begin{array}{c}\mathrm{LAS}\hfill \end{array}$
(II)	$\begin{array}{c}E{Q}_0=\left(1000,0,0,3.20,0,0,0,0,0\right)\hfill \\ E{Q}_1=(522.72,21.69,8.68,1.39,\hfill \\ 0.91,0.09,13.04,0,0)\hfill \end{array}$	$\begin{array}{c}\left(1.011,0.999,0.999,0.998,0.998,0.990,0.979,0.902,0.852\right)\hfill \\ \\ \left(0.999,0.999,0.998,0.998,0.997,0.990,0.979,0.923,0.841\right)\hfill \end{array}$	$\begin{array}{c}\mathrm{unstable}\hfill \\ \\ \mathrm{LAS}\hfill \end{array}$
(III)	$\begin{array}{c}E{Q}_0=\left(1000,0,0,3.20,0,0,0,0,0\right)\hfill \\ E{Q}_2=(592.59,0,0,3.198,0,0,0,\hfill \\ 101.85,15.28)\hfill \end{array}$	$\begin{array}{c}\left(1.002,0.999,0.999,0.998,0.997,0.990,0.974,0.927,0.840\right)\hfill \\ \multicolumn{1}{c}{}\\ \left(0.999,0.999,0.999,0.9988,0.996,0.990,0.976,0.936,0.837\right)\hfill \end{array}$	$\begin{array}{l}\mathrm{unstable}\\ \multicolumn{1}{c}{}\\ \mathrm{LAS}\end{array}$
(IV)	$\begin{array}{c}E{Q}_0=\left(1000,0,0,3.20,0,0,0,0,0\right)\hfill \\ E{Q}_1=(675.48,14.75,5.9,0.16,\hfill \\ 1.54,0.15,8.9,0,0)\hfill \\ E{Q}_2=(666.67,0,0,3.20,0,0,0,\hfill \\ 83.33,12.50)\hfill \\ E{Q}_3=(666.67,3.77,1.51,0.4,\hfill \\ 1.4,0.14,2.31,62.59,9.39)\hfill \end{array}$	$\begin{array}{c}\left(1.006,1.002,0.999,0.999,0.998,0.990,0.975,0.912,0.846\right)\hfill \\ \\ \left(1,0.999,0.999,0.998,0.990,0.976,0.973,0.923,0.841\right)\hfill \\ \\ \left(1.001,0.999,0.999,0.999,0.996,0.990,0.976,0.924,0.841\right)\hfill \\ \\ \left(0.9999,0.999,0.999,0.998,0.992,0.990,0.976,0.924,0.841\right)\hfill \end{array}$	$\begin{array}{c}\mathrm{unstable}\\ \\ \mathrm{unstable}\\ \\ \multicolumn{1}{c}{\mathrm{unstable}}\\ \\ \multicolumn{1}{c}{\mathrm{LAS}}\end{array}$

contains the nonnegative equilibria, the absolute value of the eigenvalues and whether the equilibrium point is LAS or unstable. We note that when an equilibrium point is GAS, then it is also LAS, and all the other equilibria will be unstable.

6.2. Impact of Latent Reservoirs on the HIV-1 and HTLV-I Co-Dynamics

In this part, we investigate the impact of the presence of latent reservoirs on HIV-1 and HTLV-I co-dynamics. Currently, there is no available treatment for HTLV-I [61]. Therefore, we only consider one type of HIV-1 antiviral drug, reverse transcriptase inhibitor (RTI), which prevents the HIV-1 from infecting the macrophages and CD4${}^{+}$T cells. Let $ϵ\in \left[0,1\right]$ be the efficacy of the RTI drug. Now, we write model (1) under the effect of RTI:

$$\left\{\begin{array}{c}\frac{dx}{dt}={\zeta }_1-{\varrho }_{1}x-(1-ϵ){\rho }_{1}xv-{\rho }_{3}xu,\hfill \\ \frac{dy}{dt}=(1-ϵ){\rho }_{1}xv-{\alpha }_{1}y,\hfill \\ \frac{dw}{dt}={\zeta }_2-{\varrho }_{2}w-(1-ϵ){\rho }_{2}wv,\hfill \\ \frac{dz}{dt}=(1-ϵ){\rho }_{2}wv-{\alpha }_{2}z,\hfill \\ \frac{dv}{dt}={\beta }_1{\alpha }_{1}y+{\beta }_2{\alpha }_{2}w-\theta v,\hfill \\ \frac{du}{dt}={\rho }_{3}xu-\delta u.\hfill \end{array}\right.$$

(59)

The basic reproductive numbers of model (59) are given by

$$\begin{array}{cc}\hfill R_0^{\mathrm{Without}\mathrm{latent}}\left(ϵ\right)& =\frac{(1-ϵ){\beta }_1{\rho }_1{\zeta }_1}{\theta {\varrho }_1}+\frac{(1-ϵ){\beta }_2{\rho }_2{\zeta }_2}{\theta {\varrho }_2},\hfill \\ \hfill R_1^{\mathrm{Without}\mathrm{latent}}& =\frac{{\rho }_3{\zeta }_1}{{\varrho }_1\delta }.\hfill \end{array}$$

Similarly, the model with latent (2) under the influence of RTI drug is given as:

$$\left\{\begin{array}{c}\frac{dx}{dt}={\zeta }_1-{\varrho }_{1}x-(1-ϵ){\rho }_{1}xv-{\rho }_{3}xu,\hfill \\ \frac{dg}{dt}=(1-ϵ){\rho }_{1}xv-({\pi }_1+{\mu }_1)g,\hfill \\ \frac{dy}{dt}={\pi }_{1}g-{\alpha }_{1}y,\hfill \\ \frac{dw}{dt}={\zeta }_2-{\varrho }_{2}w-(1-ϵ){\rho }_{2}wv,\hfill \\ \frac{ds}{dt}=(1-ϵ){\rho }_{2}wv-({\pi }_2+{\mu }_2)s,\hfill \\ \frac{dz}{dt}={\pi }_{2}s-{\alpha }_{2}z,\hfill \\ \frac{dv}{dt}={\beta }_1{\alpha }_{1}y+{\beta }_2{\alpha }_{2}z-\theta v,\hfill \\ \frac{dq}{dt}={\rho }_{3}xu-({\pi }_3+{\mu }_3)q,\hfill \\ \frac{du}{dt}={\pi }_{3}q-\delta u.\hfill \end{array}\right.$$

(60)

The basic reproductive numbers of model (60) are given by

$$\begin{array}{cc}\hfill R_0^{\mathrm{With}\mathrm{latent}}\left(ϵ\right)& =\frac{(1-ϵ){\beta }_1{\rho }_1{\pi }_1{\zeta }_1}{\theta {\varrho }_1({\pi }_1+{\mu }_1)}+\frac{(1-ϵ){\beta }_2{\rho }_2{\pi }_2{\zeta }_2}{\theta {\varrho }_2({\pi }_2+{\mu }_2)},\hfill \\ \hfill R_1^{\mathrm{With}\mathrm{latent}}& =\frac{{\rho }_3{\zeta }_1{\pi }_3}{{\varrho }_1\delta ({\pi }_3+{\mu }_3)}.\hfill \end{array}$$

We suppose that $R_1^{\mathrm{Without}\mathrm{latent}}\le 1$ and $R_1^{\mathrm{With}\mathrm{latent}}\le 1$. In order to stabilize the infection-free equilibrium $E{Q}_0$ for both systems (59) and (60), we determine the minimum drug efficacies ${ϵ}_{min}^{\mathrm{Without}\mathrm{latent}}$ and ${ϵ}_{min}^{\mathrm{With}\mathrm{latent}}$ that make

$$\begin{array}{cc}\hfill R_0^{\mathrm{Without}\mathrm{latent}}\left(ϵ\right)& \le 1,\mathrm{for}\mathrm{all}{ϵ}_{min}^{\mathrm{Without}\mathrm{latent}}\le ϵ\le 1,\hfill \\ \hfill R_0^{\mathrm{With}\mathrm{latent}}\left(ϵ\right)& \le 1,\mathrm{for}\mathrm{all}{ϵ}_{min}^{\mathrm{With}\mathrm{latent}}\le ϵ\le 1,\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill {ϵ}_{min}^{\mathrm{Without}\mathrm{latent}}& =max\left\{0,1-{\left(\frac{{\beta }_1{\rho }_1{\zeta }_1}{\theta {\varrho }_1}+\frac{{\beta }_2{\rho }_2{\zeta }_2}{\theta {\varrho }_2}\right)}^{-1}\right\},\hfill \\ \hfill {ϵ}_{min}^{\mathrm{With}\mathrm{latent}}& =max\left\{0,1-{\left(\frac{{\beta }_1{\rho }_1{\pi }_1{\zeta }_1}{\theta {\varrho }_1({\pi }_1+{\mu }_1)}+\frac{{\beta }_2{\rho }_2{\pi }_2{\zeta }_2}{\theta {\varrho }_2({\pi }_2+{\mu }_2)}\right)}^{-1}\right\}.\hfill \end{array}$$

We note that $R_0^{\mathrm{With}\mathrm{latent}}\le R_0^{\mathrm{Without}\mathrm{latent}}$ and ${ϵ}_{min}^{\mathrm{With}\mathrm{latent}}\le {ϵ}_{min}^{\mathrm{Without}\mathrm{latent}}$. This means that neglecting the latent reservoirs in the HIV-1 and HTLV-I coinfection model will lead to an overestimation of the required HIV-1 antiviral drugs.

Lengthening of the Latent Phase

In this part, we study the effect of lengthening the latent phase on the HIV-1 and HTLV-I co-dynamics. The average length of the latent phases are given by $1/{\pi }_i$, $i=1,2,3$. Therefore, reducing the parameter ${\pi }_i$ will increase the transition time from latent to active and then will slow the viral replication. Let us assume that there exist two control actions $ϵ$ and $\overline{ϵ}$ with the objective of reducing the transition rate (lengthening of the latent durations). These control actions may represent the efficacies of two drugs [62,63]. Let us multiply ${\pi }_1$ by $(1-ϵ)$, ${\pi }_2$ by $(1-ϵ)$ and ${\pi }_3$ by $(1-\overline{ϵ}),$ where $ϵ,\overline{ϵ}\in [0,1]$. Then, model (2) under the effect of such control actions can be written as:

$$\left\{\begin{array}{c}\frac{dx}{dt}={\zeta }_1-{\varrho }_{1}x-{\rho }_{1}xv-{\rho }_{3}xu,\hfill \\ \frac{dg}{dt}={\rho }_{1}xv-((1-ϵ){\pi }_1+{\mu }_1)g,\hfill \\ \frac{dy}{dt}=(1-ϵ){\pi }_{1}g-{\alpha }_{1}y,\hfill \\ \frac{dw}{dt}={\zeta }_2-{\varrho }_{2}w-{\rho }_{2}wv,\hfill \\ \frac{ds}{dt}={\rho }_{2}wv-((1-ϵ){\pi }_2+{\mu }_2)s,\hfill \\ \frac{dz}{dt}=(1-ϵ){\pi }_{2}s-{\alpha }_{2}z,\hfill \\ \frac{dv}{dt}={\beta }_1{\alpha }_{1}y+{\beta }_2{\alpha }_{2}z-\theta v,\hfill \\ \frac{dq}{dt}={\rho }_{3}xu-((1-\overline{ϵ}){\pi }_3+{\mu }_3)q,\hfill \\ \frac{du}{dt}=(1-\overline{ϵ}){\pi }_{3}q-\delta u.\hfill \end{array}\right.$$

(61)

Let the model’s parameters be fixed other than $ϵ$ and $\overline{ϵ}$, then the basic reproductive numbers of system (61) $R_0$ and $R_1$ as functions of $ϵ$ and $\overline{ϵ}$, respectively, are given by

$$\begin{array}{cc}\hfill R_0\left(ϵ\right)& =\frac{(1-ϵ){\beta }_1{\rho }_1{\pi }_1{\zeta }_1}{\theta {\varrho }_1((1-ϵ){\pi }_1+{\mu }_1)}+\frac{(1-ϵ){\beta }_2{\rho }_2{\pi }_2{\zeta }_2}{\theta {\varrho }_2((1-ϵ){\pi }_2+{\mu }_2)},\hfill \\ \hfill R_1\left(\overline{ϵ}\right)& =\frac{(1-\overline{ϵ}){\rho }_3{\zeta }_1{\pi }_3}{{\varrho }_1\delta ((1-\overline{ϵ}){\pi }_3+{\mu }_3)},\hfill \end{array}$$

and we have

$$\begin{array}{cc}\hfill \frac{d{R}_0}{dϵ}& =-\frac{{\beta }_1{\mu }_1{\zeta }_1{\pi }_1{\rho }_1}{{\varrho }_1\theta {((1-ϵ){\pi }_1+{\mu }_1)}^2}-\frac{{\beta }_2{\mu }_2{\zeta }_2{\pi }_2{\rho }_2}{{\varrho }_2\theta {((1-ϵ){\pi }_2+{\mu }_2)}^2}<0,\hfill \\ \hfill \frac{d{R}_1}{d\overline{ϵ}}& =-\frac{{\mu }_3{\zeta }_1{\pi }_3{\rho }_3}{{\varrho }_1\delta {((1-\overline{ϵ}){\pi }_3+{\mu }_3)}^2}<0.\hfill \end{array}$$

Hence, $R_0\left(ϵ\right)$ and $R_1\left(\overline{ϵ}\right)$ are strictly decreasing functions of $ϵ$ and $\overline{ϵ}$, respectively. Therefore, increasing the values of $ϵ$ and $\overline{ϵ}$ will decrease the parameters $R_0$ and $R_1$. Now we determine the minimum control actions ${ϵ}_{min}$ and ${\overline{ϵ}}_{min}$ that make

$$\begin{array}{cc}\hfill R_0\left(ϵ\right)& \le 1,\mathrm{for}\mathrm{all}{ϵ}_{min}\le ϵ\le 1,\hfill \\ \hfill R_1\left(\overline{ϵ}\right)& \le 1,\mathrm{for}\mathrm{all}{\overline{ϵ}}_{min}\le \overline{ϵ}\le 1,\hfill \end{array}$$

and stabilize the system around the infection-free equilibrium $E{Q}_0$. Using the values of the parameters given in Case (IV), we obtain ${ϵ}_{min}=0.853353$ and ${\overline{ϵ}}_{min}=0.666667$. Hence, both HIV-1 and HTLV-I will die out if the control actions satisfy $0.853353\le ϵ\le 1$ and $0.666667\le \overline{ϵ}\le 1$. It means that increasing the values of $ϵ$ and $\overline{ϵ}$ will lengthen the latent phase and then clear the viruses from the body. This give some impression to develop new drug therapies for lengthening the latent phase.

7. Discussion and Conclusions

In this paper, we developed and addressed a mathematical model that describes in-host co-dynamics of HIV-1 and HTLV-I with latent reservoirs. The nonlinear continuous-time model is discretized by using the NSFD method. We first proved the positivity and boundedness of the discrete-time model, and then we calculated the model’s equilibria. We found that the model has four equilibria—infection-free equilibrium $E{Q}_0$, chronic HIV-1 single-infection equilibrium $E{Q}_1$, chronic HTLV-I single-infection equilibrium $E{Q}_2$ and chronic HIV-1/HTLV-I coinfection equilibrium $E{Q}_3$. We deduced the conditions that determine the existence and stability of equilibria. These conditions were given in terms of four threshold parameters $R_j>0$, $j=0,1,2,3$. The global stability of all equilibria of the discrete-time model was examined by constructing Lyapunov functions. We proved that

• $E{Q}_0$ always exists, further, if $R_0\le 1$ and $R_1\le 1$, then $E{Q}_0$ is GAS. This result recommends that when $R_0\le 1$ and $R_1\le 1$, both HIV-1 and HTLV-I infections will be cleared regardless of the starting points.
• $E{Q}_1$ exists when $R_0>1$, and it is GAS when $R_0>1$ and $R_1\le 1$. This result recommends obtain when $R_0>1$ and $R_1\le 1$, the HIV-1 single-infection is always established regardless of the starting points.
• $E{Q}_2$ exists when $R_1>1$, and it is GAS when $R_1>1$ and $R_{02}+\frac{R_{01}}{R_1}\le 1$. This result recommends that when $R_1>1$ and $R_{02}+\frac{R_{01}}{R_1}\le 1$, the HTLV-I single-infection is always established regardless of the starting points.
• $E{Q}_3$ exists, and it is GAS when $\frac{R_1}{R_{01}}>1$, $R_2>1$ and $R_3>1$. This result recommends that the HIV-1 and HTLV-I coinfection is always established regardless of the starting points.

We simulated the discrete-time model to confirm the theoretical results. We discussed the impact of latent reservoirs on the HIV-1 and HTLV-I co-dynamics. We established that including the latent reservoirs in the HIV-1 and HTLV-I coinfection model will reduce both $R_0$ and $R_1$. Therefore, neglecting the latent reservoirs can lead to an overestimation of the required HIV-1 antiviral drugs. Further, we found that lengthening the latent phase can suppress the viral progression. This may draw the attention of scientists and pharmaceutical companies to create new treatments that lengthen the latency period.

The HIV-1 and HTLV-I coinfection model presented in this work is given by a system of ordinary differential equations. The HIV-1 and HTLV-I co-dynamics can be described by (i) delay differential equations by incorporating the time delay [21], (ii) partial differential equations by considering the mobility of cells and viruses [64], (iii) stochastic differential equations by taking into account the random fluctuations [65,66,67,68] and (iv) fractional differential equations by considering the memory effects [69]. Great efforts are needed to study such points; therefore, we leave them for future work.