Dynamics of an HTLV-I Model Incorporating Viral Replication Delay and CTL Immunity

Yinji Huang

College of Mathematics and Statistics, Northwest Normal University, Lanzhou, China.

DOI: 10.4236/jamp.2025.138156                

Abstract

Human T-cell leukemia/lymphoma virus type I (HTLV-I) is a pathogenic retrovirus, and cytotoxic T lymphocytes (CTLs) or specific CD8⁺ cytotoxic T lymphocytes are considered crucial factors in the human immune system against viral infections. Therefore, studying the impact of CTL immune response on HTLV-I infection is essential. This paper investigates a class of HTLV-I infection models that consider latent infected CD4⁺ T cells and viral replication delay (the time from initial infection of healthy CD4⁺ T cells to becoming latent CD4⁺ T infected cells). By analyzing the model, the basic reproduction number for immunological inactivation $R_0$ and the basic reproduction number for immunological activation $R_1$ are defined. Using $R_0$ and $R_1$ as thresholds, the local stability of the infection-free equilibrium $E_0$ , the immunological inactivation equilibrium point $E_1$ , and the immunological activation equilibrium point $E^{\text{*}}$ is established by analyzing the distribution of roots of the corresponding characteristic equations when $\tau =0$ . Additionally, the local stability for $\tau \ne 0$ and the existence of Hopf bifurcations with delay as a bifurcation parameter are discussed.

Keywords

HTLV-I Model, Viral Replication Delay, Stability, Hopf Bifurcation

1. Introduction

Recently, the study of population dynamics of infectious diseases has attracted widespread attention. Several mathematical models have been proposed to describe the humoral immune response during the infection process within the body, such as Human Immunodeficiency Virus (HIV), Hepatitis B Virus (HBV), and Human T-cell Leukemia/Lymphotropic Virus Type I (HTLV-I) [1]-[5]. These intrahost models can capture some basic characteristics of the immune system and generate various immune responses, many of which have been observed in experiments and clinics. The findings from intrahost modeling can also be used to guide the development of efficient antiviral drug therapies [6].

Human T-cell leukemia/lymphoma virus type I (HTLV-I) is a pathogenic retrovirus that can cause slowly developing neurological diseases, known as HTLV-I-associated myelopathy (HAM), which is a chronic inflammatory disease of the central nervous system (CNS) and is also called tropical spastic paraparesis (TSP) [7] [8]. HTLV-I infection can also lead to an aggressive blood cancer called adult T-cell leukemia/lymphoma (ATL) [9] [10]. Globally, about 20 to 40 million people are infected with HTLV-I, mainly distributed in the Caribbean, southern Japan, Central and South America, the Middle East, and equatorial Africa [11]. Most HTLV-I-infected individuals remain asymptomatic throughout their lives and are carriers of the virus (ACs). Approximately 0.25% - 3.8% of infected individuals develop HAM/TSP, and another 2% - 3% develop ATL [12].

HTLV-I infection is achieved through cellular contact between healthy CD4⁺ T cells and infected CD4⁺ T cells, whereas viral particles are passed from cell to cell via viral synapses [13]. Cytotoxic T-lymphocytes or specific CD8⁺T cytotoxic T-lymphocytes (CTLs) are thought to be important in the human immune system to fight against viral infections by inhibiting viral reproduction and killing potentially infected cells. HTLV-I infection triggers a strong cytotoxic T lymphocyte (CTL) response in individuals . CTL recognize and kill CD4⁺ T cells actively infected with HTLV-I, further reducing proviral load [14] [15]. However, experiments have shown that the cytotoxicity of CTL ultimately leads to demyelination of the HAM/TSP CNS and progression of HAM/TSP disease . HTLV-I infection is also vertically transmitted via mitosis of healthy CD4⁺ T cells carrying HTLV-I provirus . Vertical transmission permits viral transmission in the absence of the HTLV-I genome and describes a low mutation rate of the HTLV-I genome . It follows that the impact of the CTL immune response on HTLV-I infection is much more complex. Therefore, consideration of CTL immune responses in HTLV-I infection models is important for studying the development and treatment of ATL or HAM/TSP.

Due to the importance of biological significance, it is urgent to study the role of CTL immune response in HTLV-I infection models. A number of scholars have investigated the kinetic behavior of HTLV-I infection by building three-dimensional mathematical models with healthy cells, infected cells, and CTL immune responses [17]-[22]. Some scholars have also introduced viral replication delays or CTL immune response delays in 3D models to portray the kinetic behavior of HTLV-I infection models [6] [23]-[26]. However, for studying the mathematical model of HTLV-I infection, it is crucial to have latently infected cells. In 2011, M.Y. Li and A.G. Lim classified the CD4⁺ T-cell population into three different classes [27], denoting the healthy, latently infected (Tax-), and active at time $t$ , by $x\left(t\right)$ , $u\left(t\right)$ , and $y\left(t\right)$ , respectively. The number of infected (Tax+) CD4⁺ helper T cells at time $t$ was used to construct a mathematical model exploring the dynamic interaction of Tax expression in transcriptional latency and viral activation

$\{\begin{array}{l}{x}^{\prime }=\lambda -\beta xy-{\mu }_{1}x,\\ u^{\prime }=\sigma \beta xy+ϵry\left(1-\frac{x+u}{k}\right)-\left(\tau +{\mu }_2\right)u,\\ y^{\prime }=\tau u-{\mu }_{3}y,\end{array}$ (1.1)

the authors analyzed the stability of the model equilibrium and the existence of Backward Bifurcation by calculating the fundamental regeneration number $R_0$ . In 2014, A.G. Lim and P.K. Maini extended the model (1.1) [28] to incorporate two key features of HTLV-I in vivo infections: 1) viral latency, 2) HTLV-I-specific CTL responses. The following model was developed

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=\lambda -\beta xy-{\mu }_{1}x,\\ \frac{\text{d}u}{\text{d}t}=\beta xy+\gamma y-\left(\tau +{\mu }_2\right)u,\\ \frac{\text{d}y}{\text{d}t}=\tau u-\alpha yz-{\mu }_{3}y,\\ \frac{\text{d}z}{\text{d}t}=\upsilon y-{\mu }_{4}z,\end{array}$ (1.2)

where $x$ , $u$ , $y$ , $z$ denote the densities of healthy CD4⁺ T cells, latently infected CD4⁺ T cells, actively infected CD4⁺ T cells, and HTLV-I specific CD8⁺ CTLs, respectively, at the time $t$ .

2016 S.M. Li, Y.C. Zhou constructed a mathematical model [29] which considered spontaneous expression of the HTLV-I antigen Tax, CTL immune response, intercellular propagation, and mitotic propagation, defined two parameters ($R_0$ ) and ($R_1$ ) to study the dynamics of the model, and discussed the Backward Bifurcation. 2021 S. Khajanchi, S. Bera, T.K. Roy extended the model proposed by and [27] to study a four-dimensional model containing (healthy cells, latently infected cells, actively infected cells, HTLV-I specific CTLs) , compared to the model (1.2), the authors considered the activated CD4⁺ T cell clearance term due to HTLV-I infection and used $\upsilon yz$ to describe CTL proliferation. The authors analyzed the local and global stability of the equilibrium point by calculating $R_0,R_1$ , while consistent persistence was described using geometric methods and a global asymptotic stability analysis was performed. In the same year, C.W. Song, R. Xu, inspired by the literature and [28], considered a model of HTLV-I infection with delayed CTL immune response [31], replacing $\upsilon y$ with $\upsilon y\left(t-\tau \right)z\left(t-\tau \right)$ , which denotes that the amount of CTL produced at time $t$ depends on the concentration and activity of CTL at time $t-\tau$ the concentration of CTL at time $t$ and the concentration of actively infected CD4⁺ T cells. The local stability of the equilibrium point is analyzed by calculating $R_0,R_1$ and the existence of Hopf branches is discussed, followed by analysis of the global stability of the equilibrium point of the model by constructing a Lyapunov function. 2022 S. Bera, S. Khajanchi, T. K. Roy in the literature modeled the delay in CTL immunity by adding the CTL immunity delay $\upsilon y\left(t-\tau \right)z\left(t-\tau \right){\text{e}}^{-m\tau }$ [32], to analyze the local stability of the equilibrium point by calculating $R_0,R_1$ and the global stability of the equilibrium point of the model by constructing the Lyapunov function, followed by a discussion of the existence of the Hopf branch as well as the orientation of the Hopf branch and period of the solution. In the same year, R. Xu and Y. Yang considered the saturated occurrence rate of CTL proliferation, the time delay of immune response, and the immune impairment term in the model presented in reference , discussing both the local and global dynamics of the model [33]. In 2022 A.M. Elaiw, A.S. Shflota, A.D. Hobinya [34] inspired by (Katri and Ruan) [35] proposed a HTLV-I kinetic model with two distributed time delays based on the literature ; this was followed in 2023 by the formulation and development of a generic HTLV-I mathematical model [36].

Inspired by the above, in this paper, we consider the HTLV-I model with a time delay between initial infection of uninfected CD4⁺ T cells and becoming latent CD4⁺ T-infected cells and CTL immune response, we first give the model without considering the time delay as follows

$\{\begin{array}{l}\frac{\text{d}U}{\text{d}t}=\lambda -\overline{v}U-\beta UA,\\ \frac{\text{d}L}{\text{d}t}=\beta UA+{\epsilon }^{\text{*}}A-\left(\pi +\mu \right)L,\\ \frac{\text{d}A}{\text{d}t}=\pi L-{\delta }^{\text{*}}A-\alpha AC,\\ \frac{\text{d}C}{\text{d}t}=\rho AC-\gamma C.\end{array}$ (1.3)

where $U$ , $L$ , $A$ , $C$ denote the densities of healthy CD4⁺ T-cells, latently infected CD4⁺ T-cells, actively infected CD4⁺ T-cells, and HTLV-I specific CD8⁺ CTLs at the time of $t$ , $\lambda$ denotes the rate of healthy CD4⁺ T-cells production, respectively; $\overline{v}$ , $\mu$ , ${\delta }^{\text{*}}$ , $\gamma$ denote the natural mortality rate of healthy CD4⁺ T cells, latently infected CD4⁺ T cells, actively infected CD4⁺ T cells, and HTLV-I-specific CD8⁺ CTL, respectively; $\beta$ denotes the coefficient of infectivity; ${\epsilon }^{\text{*}}$ mitotic rate of infected CD4⁺ T cells; $\pi$ denotes the rate of conversion of latently infected to infected cells; $\alpha$ denotes the rate at which the CTL immune response removes infected cells; and $\rho$ denotes the strength of the CTL immune response or CTL proliferation.

The rest of the paper is as follows, in Section II, the model is analyzed by discretizing the model and then introducing the viral replication delay, analyzing the positivity and boundedness of the model, analyzing the existence of the equilibrium point by calculating $R_0,R_1$ , and in Section III, discussing the stability of the equilibrium point and verifying the existence of the Hopf branch for $\tau =0$ and $\tau \ne 0$ , respectively, Section IV performs numerical simulations of the results in the paper using Matlab, and finally summarizes the conclusions of the paper.

2. Model Analysis and Equilibrium Point

2.1. Simplify the Model and Introduce Time Delay

In this section we first invariably steel the model (1.3) such that $u=aU$ , $v=bL$ , $w=\eta A$ , $z=dC$ , $\tau =et$ , it can be obtained

$\begin{array}{l}\frac{\text{d}u}{\text{d}\tau }=\frac{\text{d}u}{\text{d}t}\frac{\text{d}t}{\text{d}\tau }=\frac{a}{e}\frac{\text{d}U}{\text{d}t}=\frac{a}{e}\left(\lambda -\frac{\overline{v}}{a}u-\frac{\beta }{a\eta }uw\right)=\frac{a\lambda }{e}-\frac{\overline{v}}{e}u-\frac{\beta }{e\eta }uw,\\ \frac{\text{d}v}{\text{d}\tau }=\frac{\text{d}v}{\text{d}t}\frac{\text{d}t}{\text{d}\tau }=\frac{b}{e}\frac{\text{d}L}{\text{d}t}=\frac{b}{e}\left(\frac{\beta }{a\eta }uw+\frac{{\epsilon }^{\text{*}}}{\eta }w-\frac{\pi +\mu }{b}v\right)=\frac{b\beta }{ea\eta }uw+\frac{b{\epsilon }^{\text{*}}}{e\eta }w-\frac{\pi +\mu }{e}v,\\ \frac{\text{d}w}{\text{d}\tau }=\frac{\text{d}w}{\text{d}t}\frac{\text{d}t}{\text{d}\tau }=\frac{\eta }{e}\frac{\text{d}A}{\text{d}t}=\frac{\eta }{e}\left(\frac{\pi }{b}v-\frac{{\delta }^{\text{*}}}{\eta }w-\frac{\alpha }{d\eta }wz\right)=\frac{\eta \pi }{eb}v-\frac{{\delta }^{\text{*}}}{e}w-\frac{\alpha }{de}wz,\\ \frac{\text{d}z}{\text{d}\tau }=\frac{\text{d}z}{\text{d}t}\frac{\text{d}t}{\text{d}\tau }=\frac{d}{e}\frac{\text{d}C}{\text{d}t}=\frac{d}{e}\left(\frac{\rho }{d\eta }wz-\frac{\gamma }{d}z\right)=\frac{\rho }{e\eta }wz-\frac{\gamma }{e}z.\end{array}$

Parameter normalization:

$\frac{\pi +\mu }{e}=1,\frac{a\lambda }{e}=1,\frac{\beta }{e\eta }=1,\frac{b{\epsilon }^{\text{*}}}{e\eta }=1,\frac{\alpha }{de}=1,$

and replace $\tau$ by $t$ to obtain the following model

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=1-hu-uw,\\ \frac{\text{d}v}{\text{d}t}=fuw+w-v,\\ \frac{\text{d}w}{\text{d}t}=lv-mw-wz,\\ \frac{\text{d}z}{\text{d}t}=nwz-kz,\end{array}$ (2.1)

where

$\begin{array}{l}a=\frac{\pi +\mu }{\lambda },b=\frac{\beta }{{\epsilon }^{\text{*}}},\eta =\frac{\beta }{\pi +\mu },d=\frac{\alpha }{\pi +\mu },e=\pi +\mu ,h=\frac{\overline{v}}{\pi +\mu },\\ f=\frac{\beta \lambda }{\overline{v}\left(\pi +\mu \right)},l=\frac{\overline{v}\pi }{{\left(\pi +\mu \right)}^2},m=\frac{{\delta }^{\text{*}}}{\pi +\mu },n=\frac{\rho }{\beta },k=\frac{\gamma }{\pi +\mu }.\end{array}$

Inspired by (Katri and Ruan) [35] and considering the effect of time delay on viral replication, we introduce a viral replication delay in model (2.1), which describes the time between the initial infection of a healthy CD4⁺ T cell and the time when it becomes a latent CD4⁺ T-infected cell, obtaining the following model

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=1-hu-uw,\\ \frac{\text{d}v}{\text{d}t}=fu\left(t-\tau \right)w\left(t-\tau \right)+w-v,\\ \frac{\text{d}w}{\text{d}t}=lv-mw-wz,\\ \frac{\text{d}z}{\text{d}t}=nwz-kz,\end{array}$ (2.2)

where $u,v,w,z$ denote the densities of healthy CD4⁺ T cells, latently infected CD4⁺ T cells, actively infected CD4⁺ T cells, and HTLV-I specific CD8⁺ CTLs, respectively.

The latency period $\tau$ primarily reflects the time delay required for healthy CD4⁺ T cells ($u$ ) to be infected and converted into latent infected cells ($v$ ). Specifically, it refers to the time lag in the infection process: When the virus (HTLV-I) infects healthy cells ($u$ ), it does not immediately convert them into latent infected cells ($v$ ). The virus requires time to complete the steps of entering the cell $\to$ reverse transcription $\to$ integration into the host genome $\to$ gene expression, a process that may take several hours to several days. Therefore, the number of latent infected cells ($v$ ) at time $t$ depends on the interaction between healthy cells ($u$ ) and infected cells ($w$ ) at time $t-\tau$ , i.e., it is determined by $fu\left(t-\tau \right)w\left(t-\tau \right)$ .

2.2. The Existence of Model Equilibrium Points and Basic Reproduction Numbers

In this section, we define $R_0$ as the basic reproduction number for viral infection. The basic reproduction number $R_0$ indicates the average number of secondary infections produced by an infected individual over the entire course of their infection in a population of completely susceptible individuals. Additionally, to study the long-term dynamics of HTLV, we need to determine $R_1$ , the basic reproduction number for CTL responses, which represents the number of CTL immune stimuli and determines when CTL immunity can be activated. The value of $R_1$ describes the risk of developing HAM/TSP in a population where most individuals are lifelong asymptomatic carriers (AC) [30].

$R_0$ : Describes the virus’s ability to spread in a completely susceptible population (without immune response). $R_1$ : Describes the virus’s net replication capacity in the presence of an immune response (CTL activity). $R_0>1$ indicates that the virus can sustain transmission within the host, leading to an increase in the number of infected cells $\left(v,w\right)$ and the establishment of chronic infection, corresponding to the asymptomatic carrier phase of HTLV-I infection (where the virus is present but does not cause disease); $R_0<1$ indicates that the virus cannot effectively replicate, and the infection is spontaneously cleared (e.g., by innate immune suppression). $R_1>1$ indicates that even in the presence of CTL responses, the virus can evade immune control, with infected cells persisting (e.g., immune escape), potentially corresponding to the active disease phase (e.g., HAM/TSP or ATL); $R_1<1$ indicates that CTLs effectively suppress the virus, with the system tending toward immune control homeostasis (low viral load, e.g., long-term asymptomatic).

Model (2.2) always has an uninfected equilibrium point $E_0=\left(\frac{1}{h},0,0,0\right)$ . In the following, we compute the basic regeneration number $R_0$ .

Firstly, we isolate the infection subsystem, and decompose the right term of the infection subsystem into $ℱ-\mathcal{V}$ , the new infection term ($ℱ$ ): describes the generation of a new infected cell, and the transfer term ($\mathcal{V}$ ): describes the removal or transfer of an infected cell

$\left(\begin{array}{l}\frac{\text{d}u}{\text{d}t}\hfill \\ \frac{\text{d}w}{\text{d}t}\hfill \end{array}\right)=\underset{ℱ}{\underbrace{\left(\begin{array}{c}fu\left(t-\tau \right)w\left(t-\tau \right)+w\\ 0\end{array}\right)}}-\underset{\mathcal{V}}{\underbrace{\left(\begin{array}{c}v\\ -lv+mw\end{array}\right)}}.$

Solve for the Jacobi matrix at $E_0$ . Assume that the time lag term can be approximated near the steady state as $w\left(t-\tau \right)\approx w\left(t\right)$ .

$F={\frac{\partial ℱ}{\partial \left(v,w\right)}|}_{E_0}=\left(\begin{array}{cc}0& f{u}_0+1\\ 0& 0\end{array}\right),V={\frac{\partial \mathcal{V}}{\partial \left(v,w\right)}|}_{E_0}=\left(\begin{array}{cc}1& 0\\ -l& m\end{array}\right).$

Calculate the next generation matrix $F{V}^{-1}$

$F{V}^{-1}=\left(\begin{array}{cc}0& f{u}_0+1\\ 0& 0\end{array}\right)\cdot \frac{1}{m}\left(\begin{array}{cc}m& 0\\ l& 1\end{array}\right)=\frac{1}{m}\left(\begin{array}{cc}l\left(f{u}_0+1\right)& f{u}_0+1\\ 0& 0\end{array}\right).$

The basic regeneration number $R_0$ is the spectral radius $\rho \left(F{V}^{-1}\right)$ , i.e., the largest eigenvalue of $F{V}^{-1}$ .

$R_0=\rho \left(F{V}^{-1}\right)=\frac{l\left(f{u}_0+1\right)}{m},$

substitute $u_0=1/h$ into the equation to obtain:

$R_0=\frac{l\left(f+h\right)}{mh}.$

In addition to the infection-free equilibrium, if $R_0>1$ and $m>l$ , there exists an immune inactivation equilibrium $E_1=\left(u_1,v_1,w_1,0\right)$ for model (2.2), where

$\begin{array}{l}{u}_1=\frac{m-l}{fl},v_1=\frac{m\left(lf+lh-mh\right)}{l\left(m-l\right)}=\frac{m^{2}h\left(R_0-1\right)}{l\left(m-l\right)},\\ w_1=\frac{lf+lh-mh}{m-l}=\frac{mh\left(R_0-1\right)}{m-l}.\end{array}$

Additionally, through calculation, we obtained the immune activation reproduction number

$R_1=\frac{l\left(nf+nh+k\right)}{m\left(nh+k\right)},$

if $R_1>1$ , in addition to $E_0$ , $E_1$ , there exists an immune activation equilibrium $E^{\text{*}}\left(u^{\text{*}},v^{\text{*}},w^{\text{*}},z^{\text{*}}\right)$ for model (2.2), where

$\begin{array}{l}{u}^{\text{*}}=\frac{n}{nh+k},v^{\text{*}}=\frac{k}{n}\frac{nf+nh+k}{nh+k}=\frac{mk}{nl}{R}_1,\\ w^{\text{*}}=\frac{k}{n},z^{\text{*}}=\frac{nfl+nh\left(l-m\right)+k\left(l-m\right)}{nh+k}=m\left(R_1-1\right).\end{array}$

2.3. Positive and Boundedness Analysis

Model (2.2) must be studied under the following initial conditions. Spatially define $\phi =\left({\phi }_1,{\phi }_2,{\phi }_3,{\phi }_4\right)$ , the

$\begin{array}{l}{C}_{+}=\{\phi \in C\left(\left[-\tau ,0\right],{ℝ}_{+}^4\right):u\left(\theta \right)={\phi }_1\left(\theta \right),v\left(\theta \right)={\phi }_2\left(\theta \right),\\ w\left(\theta \right)={\phi }_3\left(\theta \right),z\left(\theta \right)={\phi }_4\left(\theta \right)\},\end{array}$ (2.3)

where${\phi }_i\left(\theta \right)\ge 0$ , $\theta \in \left[-\tau ,0\right]$ , ${\phi }_i>0$ , $i=1,2,3,4$ , $\phi =\left({\phi }_1,{\phi }_2,{\phi }_3,{\phi }_4\right)\in C\left(\left[-\tau ,0\right],{ℝ}_{+}^4\right)$ , $C$ denotes the Banach space of continuous real-valued functions from $\left[-\tau ,0\right]\left(\tau >0\right)$ to ${ℝ}_{+}^4$ with sub-norm

$\Vert \phi \Vert =\text{max}\left\{\underset{-\tau \le \theta \le 0}{\text{sup}}\left|{\phi }_1\left(\theta \right)\right|,\underset{-\tau \le \theta \le 0}{\text{sup}}\left|{\phi }_2\left(\theta \right)\right|,\underset{-\tau \le \theta \le 0}{\text{sup}}\left|{\phi }_3\left(\theta \right)\right|,\underset{-\tau \le \theta \le 0}{\text{sup}}\left|{\phi }_4\left(\theta \right)\right|\right\},\phi \in C_{+},$

${ℝ}_{+}^4=\left\{\left(u,v,w,z\right):u\ge 0,u\ge 0,w\ge 0,z\ge 0\right\}.$

Theorem 2.1 Any solution of model (2.2) under initial conditions (2.3) is defined on $\left[0,+\infty \right)$ and remains positive for all $t>0$ .

Proof: Let $\left(u\left(t\right),v\left(t\right),w\left(t\right),z\left(t\right)\right)$ be an arbitrary solution of model (2.2) under the initial condition (2.3), firstly, according to the first equation of model (2.2), we have

${\frac{\text{d}u}{\text{d}t}|}_{u=0}=1>0.$

Its solution is

$u\left(t\right)=u\left(0\right){\text{e}}^{-{\displaystyle {\int }_0^t\left(h+w\left(s\right)\right)\text{d}s}}+{\displaystyle {\int }_0^t{\text{e}}^{-{\displaystyle {\int }_s^t\left(h+w\left(\xi \right)\right)\text{d}\xi }}\text{d}s}>0.$

For $u\left(t\right)>0$ , the equation $\frac{\text{d}v}{\text{d}t}=fu\left(t-\tau \right)w\left(t-\tau \right)+w-v$ has derivative $fu\left(t-\tau \right)w\left(t-\tau \right)+w\ge 0$ for $v=0$ , and the initial value of $v\left(0\right)>0$ , so $v\left(t\right)>0$ . The solution is

$v\left(t\right)=v\left(0\right){\text{e}}^{-t}+{\displaystyle {\int }_0^t\left(fu\left(s-\tau \right)w\left(s-\tau \right)+w\left(s\right)\right){\text{e}}^{-\left(t-s\right)}\text{d}s}>0.$

For $w\left(t\right)>0$ , the equation $\frac{\text{d}w}{\text{d}t}=lv-mw-wz$ has derivative $lv\ge 0$ at $w=0$ and initial value $w\left(0\right)>0$ , so $w\left(t\right)>0$ . The solution is

$w\left(t\right)=w\left(0\right){\text{e}}^{-{\displaystyle {\int }_0^t\left(m+z\left(s\right)\right)\text{d}s}}+{\displaystyle {\int }_0^{t}lv\left(s\right){\text{e}}^{-{\displaystyle {\int }_s^t\left(m+z\left(\xi \right)\right)\text{d}\xi }}\text{d}s}>0.$

For $z\left(t\right)>0$ , the equation $\frac{\text{d}z}{\text{d}t}=z\left(nw-k\right)$ is solved by

$z\left(t\right)=z\left(0\right){\text{e}}^{{\displaystyle {\int }_0^t\left(nw\left(s\right)-k\right)\text{d}s}}>0,$

Because $z\left(0\right)>0$ and the exponential function is positive.

Theorem 2.2 There exists a constant $M>0$ such that all solutions satisfy

$\underset{t\to \infty }{\lim \sup }u\left(t\right)\le M,\underset{t\to \infty }{\lim \sup }v\left(t\right)\le M,\underset{t\to \infty }{\lim \sup }w\left(t\right)\le M,\underset{t\to \infty }{\lim \sup }z\left(t\right)\le M.$

Proof Since $\frac{\text{d}u}{\text{d}t}=1-hu-uw\le 1-hu$ , solve the differential inequality:

$u\left(t\right)\le \frac{1}{h}+\left(u\left(0\right)-\frac{1}{h}\right){\text{e}}^{-ht}\Rightarrow \underset{t\to \infty }{\lim \sup }u\left(t\right)\le \frac{1}{h}.$

Construct the composite function $N\left(t\right)$ , Definition:

$N\left(t\right)=u\left(t\right)+v\left(t\right)+w\left(t\right)+\frac{f}{n}z\left(t+\tau \right),$

differentiation

$\begin{array}{c}\frac{\text{d}N}{\text{d}t}=\underset{\le 1}{\underbrace{1-hu}}+\underset{\left(\text{i}\right)}{\underbrace{fu\left(t-\tau \right)w\left(t-\tau \right)-uw}}+\underset{\left(\text{ii}\right)}{\underbrace{w-v+lv-mw}}\\ +\underset{\left(\text{iii}\right)}{\underbrace{fw\left(t+\tau \right)z\left(t+\tau \right)-\frac{fk}{n}z\left(t+\tau \right)}}.\end{array}$

Since $u\left(t\right)\le \frac{1}{h}$ and $w\ge 0$ , the term 1) satisfies

$fu\left(t-\tau \right)w\left(t-\tau \right)-uw\le f\cdot \frac{1}{h}w\left(t-\tau \right),$

the term 2) satisfies

$\left(1-m\right)w+\left(l-1\right)v\le -m^{\prime }\left(v+w\right),m^{\prime }=\text{min}\left\{m,1-l\right\},$

the term 3) satisfies

$fw\left(t+\tau \right)z\left(t+\tau \right)-\frac{fk}{n}z\left(t+\tau \right)\le -\frac{fk}{n}z\left(t+\tau \right),$

therefore, it can be concluded

$\frac{\text{d}N}{\text{d}t}\le 1-m^{″}N\left(t\right),m^{″}=\text{min}\left\{h,m^{\prime },k\right\}.$

Solve the inequality

$N\left(t\right)\le \frac{1}{m^{″}}+\left(N\left(0\right)-\frac{1}{m^{″}}\right){\text{e}}^{-m^{″}t}\underset{t\to \infty }{\lim \sup }N\left(t\right)\le \frac{1}{m^{″}},$

from the boundedness of $N\left(t\right)$ , we directly obtain:

$u\left(t\right),v\left(t\right),w\left(t\right),z\left(t\right)\le \frac{1}{m^{″}}.$

Therefore, the solution is globally bounded by:

$\text{Ω}=\left\{\left(u,v,w,z\right)\in {ℝ}_{+}^4:u\left(t\right)\le \frac{1}{h},u\left(t\right)+v\left(t\right)+w\left(t\right)+\frac{f}{n}z\left(t+\tau \right)\le \frac{1}{m^{″}}\right\}.$

3. Stability of Equilibrium Points and Existence of Hopf Bifurcations

To study the local asymptotic stability of the model (2.2) at each equilibrium point, we compute the Jacobi matrix, which is expressed at $E\left(u,v,w,z\right)$ as

$J=\left(\begin{array}{cccc}-\left(h+w\right)& 0& -u& 0\\ fw{\text{e}}^{-\lambda \tau }& -1& 1+fu{\text{e}}^{-\lambda \tau }& 0\\ 0& l& -m-z& -w\\ 0& 0& nz& nw-k\end{array}\right).$

Theorem 3.1 When $\tau =0$ , if the basic regeneration number $R_0<1$ , the disease-free equilibrium point $E_0$ of model (2.2) is locally asymptotically stable; if $R_0>1$ , $E_0$ is unstable. When $\tau >0$ , $E_0$ is locally asymptotically stable if $R_0<1$ .

Proof The characteristic equation of model (2.2) at the equilibrium point $E_0$ when $\tau =0$ is

${h\left(\lambda ,0\right)|}_{E_0}=\left(\lambda +k\right)\left(\lambda +h\right)\left[{\lambda }^2+\left(1+m\right)\lambda +m\left(1-R_0\right)\right]=0,$ (3.1)

clearly (3.1) has negative real roots ${\lambda }_1=-k,{\lambda }_2=-h$ , and the other roots of (3.1) are determined by the following equation

$f\left(\lambda ,0\right)={\lambda }^2+\left(1+m\right)\lambda +m\left(1-R_0\right)=0.$

If $R_0<1$ , it is easy to show that all roots of $f\left(\lambda \right)$ have negative real parts. Therefore, the equilibrium point $E_0$ is locally asymptotically stable.

If $R_0>1$ , then we have $f\left(0\right)<0$ and $f\left(\lambda \right)\to +\infty \left(\lambda \to +\infty \right)$ . Therefore, $f\left(\lambda \right)$ has at least one positive real root. Accordingly, $E_0$ is unstable.

The characteristic equation of model (2.2) at the equilibrium point $E_0$ when $\tau >0$ is

${h\left(\lambda ,\tau \right)|}_{E_0}=\left(\lambda +k\right)\left(\lambda +h\right)\left[{\lambda }^2+\left(1+m\right)\lambda +m-l-\frac{lf}{h}{\text{e}}^{-\lambda \tau }\right]=0,$ (3.2)

clearly (3.2) has negative real roots ${\lambda }_1=-k,{\lambda }_2=-h$ , and the other roots of (3.2) are determined by the following equation

$f\left(\lambda ,\tau \right)={\lambda }^2+\left(1+m\right)\lambda +m-l-\frac{lf}{h}{\text{e}}^{-\lambda \tau }=0.$ (3.3)

Assuming that $\text{i}\omega$ is a solution to Equation (3.3), bringing in (3.3) yields

$-{\omega }^2+\left(1+m\right)\text{i}\omega +m-l-\frac{lf}{h}\left(\cos \omega \tau -\text{i}\sin \omega \tau \right)=0,$

separating the real and imaginary parts of the above equation yields

$\{\begin{array}{l}-{\omega }^2+m-l=\frac{lf}{h}\cos \omega \tau ,\\ \left(1+m\right)\omega =-\frac{lf}{h}\sin \omega \tau ,\end{array}$

this is equivalent to

${\omega }^4+\left[{\left(m+1\right)}^2+2\left(l-m\right)\right]{\omega }^2+{\left(m-l\right)}^2-{\left(\frac{lf}{h}\right)}^2=0,$ (3.4)

therefore

${\omega }^2=\frac{-\left[{\left(m+1\right)}^2+2\left(l-m\right)\right]\pm \sqrt{{\left[{\left(m+1\right)}^2+2\left(l-m\right)\right]}^2-4\left[{\left(m-l\right)}^2-{\left(\frac{lf}{h}\right)}^2\right]}}{2},$

when $\left(m-l\right)>\frac{lf}{h}$ , i.e., when $R_0<1$ , ${\omega }^2<0$ , it means that Equation (3.4) has no positive roots, i.e., Equation (3.3) has no pure imaginary roots. According to [[37], Corollary 2.4], all roots of the characteristic Equation (3.3) have negative real parts, at which point $E_0$ is locally asymptotically stable for the time lag $\tau >0$ .

Theorem 3.2 When $\tau =0$ , if $R_1<1<R_0$ , then model (2.2) is locally asymptotically stable at the immune inactivation equilibrium point $E_1=\left(u_1,v_1,w_1,0\right)$ ; and if $R_1>1$ , $E_1$ is unstable.

Proof When $\tau =0$ , the characteristic equation for model (2.2) at $E_1$ is

${h\left(\lambda ,0\right)|}_{E_1}=\left[\lambda -\left(n{w}_1-k\right)\right]\left({\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3\right)=0,$ (3.5)

where

$\begin{array}{l}{A}_1=h+1+m+w_1>0,\\ A_2=h+w_1+m\left(h+w_1\right)>0,\\ A_3=\left(m-l\right)w_1>0.\end{array}$

When $R_1<1$ , $w_1<\frac{k}{n}$ , we know that Equation (3.5) has a negative real root $\lambda =n{w}_1-k<0$ , and for

${\lambda }^3+A_1{\lambda }^2+A_2\lambda +A_3=0,$ (3.6)

the calculation shows that ${\text{Δ}}_2=A_1{A}_2-A_3>0$ . By the Routh-Hurwitz criterion, we know the sufficient condition for stability, i.e. the root of (3.6) is negative or has a negative real part if $A_1>0,A_3>0$ and $A_1{A}_2-A_3>0$ . Thus when $R_1<1<R_0$ , model (2.2) is locally asymptotically stable at the immune inactivation equilibrium point $E_1$ .

Correspondingly, when $R_1>1$ , $w_1>\frac{k}{n}$ , we know that Equation (3.5) has a positive real root $\lambda =n{w}_1-k>0$ , and hence $E_1$ is unstable.

In the following, we discuss the stability of model (2.2) at $E_1$ when $\tau >0$ .

When $\tau >0$ , the characteristic equation of model (2.2) at $E_1$ is

${h\left(\lambda ,\tau \right)|}_{E_1}=\left[\lambda -\left(n{w}_1-k\right)\right]\left[\left({\lambda }^3+a_{11}{\lambda }^2+a_{12}\lambda +a_{13}\right)+\left(b_{11}\lambda +b_{12}\right){\text{e}}^{-\lambda \tau }\right]=0,$ (3.7)

where

$\begin{array}{l}{a}_{11}=h+1+m+w_1>0,\\ a_{12}=h+w_1+m\left(h+w_1\right)+m-l>0,\\ a_{13}=\left(m-l\right)\left(h+w_1\right)>0,\\ b_{11}=-lf{u}_1<0,\\ b_{12}=-lhf{u}_1<0.\end{array}$

When $R_1>1$ , then $w_1>\frac{k}{n}$ , we know that Equation (3.5) has a positive real root $\lambda =n{w}_1-k>0$ , and hence $E_1$ is unstable. Correspondingly, when $R_1<1$ , then $w_1<\frac{k}{n}$ , we know that Equation (3.5) has a negative real root $\lambda =n{w}_1-k<0$ , and for

$\left({\lambda }^3+a_{11}{\lambda }^2+a_{12}\lambda +a_{13}\right)+\left(b_{11}\lambda +b_{12}\right){\text{e}}^{-\lambda \tau }=0,$ (3.8)

suppose $\text{i}\omega$ is a root of Equation (3.8) if and only if $\omega$ satisfies

$-\text{i}{\omega }^3-a_{11}{\omega }^2+a_{12}\text{i}\omega +a_{13}+\left(b_{11}\text{i}\omega +b_{12}\right)\left(\cos \omega \tau -\text{i}\sin \omega \tau \right)=0,$

separating the real and imaginary parts gives

$\{\begin{array}{l}{\omega }^3-a_{12}\omega =b_{11}\omega \cos \omega \tau -b_{12}\sin \omega \tau ,\\ a_{11}{\omega }^2-a_{13}=b_{12}\cos \omega \tau +b_{11}\omega \sin \omega \tau ,\end{array}$ (3.9)

the squares of the two Equations of (3.9) are added together to give

${\omega }^6+\left(a_{11}-2a_{12}\right){\omega }^4+\left(a_{12}^2-2a_{11}{a}_{13}-b_{11}^2\right){\omega }^2+a_{13}^2-b_{12}^2=0.$ (3.10)

Let $x={\omega }^2$ , (3.10) be reduced to

$x^3+g_1{x}^2+g_{2}x+g_3=0,$ (3.11)

where

$g_1=a_{11}-2a_{12},g_2=a_{12}^2-2a_{11}{a}_{13}-b_{11}^2,g_3=a_{13}^2-b_{12}^2.$

Find below the condition for (3.11) to have at least one positive root, let

$g\left(x\right)=x^3+g_1{x}^2+g_{2}x+g_3,$

Case 1: $g_3<0$ . Since ${\lim }_{x\to +\infty }g\left(x\right)=+\infty$ , $g\left(0\right)=g_3<0$ , then $g\left(x\right)=0$ has at least one positive root.

Case 2: Since

$g^{\prime }\left(x\right)=3x^2+2g_{1}x+g_2,$

and $\text{Δ}=4g_1^2-12g_2$ , so if $\text{Δ}\le 0$ , then $g\left(x\right)=0$ does not have a positive real root; if $\text{Δ}>0$ , then (3.11) has two zeros

$x_1=\frac{-2g_1+\sqrt{\text{Δ}}}{6},x_2=\frac{-2g_1-\sqrt{\text{Δ}}}{6}.$

Since $g^{″}\left(x\right)=6x+2g_1$ , then $g^{″}\left(x_1\right)>0$ , $g^{″}\left(x_2\right)<0$ . That is, $x_1$ is a local minima of $g\left(x\right)$ , and $x_2$ is a local maxima of $g\left(x\right)$ .

Lemma 3.1 [37] Let $g_3=a_{13}^2-b_{12}^2\ge 0$ and $\text{Δ}=4g_1^2-12g_2>0$ , then there exists a positive root of the equation $g\left(x\right)=0$ if and only if $x_1>0$ and $g\left(x_1\right)\le 0$ .

Proof Sufficiency: When $g_3\ge 0$ and $\text{Δ}>0$ , if $x_1>0$ and $g\left(x_1\right)\le 0$ , then $x=x_1$ is a root of $g\left(x\right)=0$ . If $g\left(x\right)<0$ , since $\underset{x\to +\infty }{\text{lim}}g\left(x\right)=+\infty$ , it follows that there exists $x_3>x_1$ such that $g\left(x_3\right)>0$ , and $g\left(x\right)=0$ has at least one positive root in $\left(x_1,x_3\right)$ .

Necessity: If $g_3\ge 0$ and $\text{Δ}>0$ , there exists a positive root of $g\left(x\right)=0$ , and by the monotonicity of $g\left(x\right)$ , $g\left(x\right)$ is an increasing function on $\left(-\infty ,x_2\right)$ and $\left(x_1,+\infty \right)$ and a decreasing function on $\left(x_2,x_1\right)$ . Assuming that Assuming that $x_1\le 0$ , then $g\left(x_1\right)\le g\left(0\right)=g_3$ , by the monotonicity of $g\left(x\right)$ , $g\left(x\right)>g\left(0\right)$ holds for any $x>0$ . At this point, $g\left(x\right)$ has no root in $\left(0,+\infty \right)$ , which contradicts the assumption. Therefore $x_1>0$ . Assuming that $g\left(x_1\right)>0$ , it follows from the monotonicity of $g\left(x\right)$ that $g\left(x_2\right)>g\left(x_1\right)$ and $g\left(0\right)=g_3\ge 0$ , the equation does not have any positive roots on $\left(0,+\infty \right)$ , which is a contradiction in terms, hence $g\left(x_1\right)\le 0$ .

Lemma 3.2 ([37] Lemma 2.1) For Equation (3.11), the following conclusions are drawn.

1) If $g_3<0$ , Equation (3.11) has at least one positive root.

2) If $g_3\ge 0$ and $\text{Δ}\le 0$ , Equation (3.11) has no positive roots.

3) If $g_3\ge 0$ and $\text{Δ}>0$ , Equation (3.11) has positive roots if and only if $x_1>0$ and $g\left(x_1\right)\le 0$ .

The following determines the positive and negative of $g_1,g_2,g_3$ .

$\begin{array}{c}{g}_1=a_{11}-2a_{12}\\ =h+1+m+w_1-2\left(h+w_1+m\left(h+w_1\right)+m-l\right),\end{array}$

$\begin{array}{c}{g}_2=a_{12}^2-2a_{11}{a}_{13}-b_{11}^2=\left(a_{12}+b_{11}\right)\left(a_{12}-b_{11}\right)-2a_{11}{a}_{13}\\ =\left[h+w_1+m\left(h+w_1\right)+m-l-lf{u}_1\right]\left[h+w_1+m\left(h+w_1\right)+m-l+lf{u}_1\right]\\ -2\left(h+1+m+w_1\right)\left(\left(m-l\right)\left(h+w_1\right)\right)\\ =\left[h+w_1+m\left(h+w_1\right)\right]\left[h+w_1+m\left(h+w_1\right)+2\left(m-l\right)\right]\\ -2\left(h+1+m+w_1\right)\left(\left(m-l\right)\left(h+w_1\right)\right),\end{array}$

$\begin{array}{c}{g}_3=a_{13}^2-b_{12}^2=\left(a_{13}+b_{12}\right)\left(a_{13}-b_{12}\right)\\ =\left[\left(m-l\right)\left(h+w_1\right)-lhf{u}_1\right]\left[\left(m-l\right)\left(h+w_1\right)+lhf{u}_1\right]\\ =\left[\left(m-l\right)\left(h+w_1\right)-h\left(m-l\right)\right]\left[\left(m-l\right)\left(h+w_1\right)+h\left(m-l\right)\right]>0,\end{array}$

we know that $g_3>0$ , and $g_1,g_2$ is not good for positive or negative judgment.

Without loss of generality, assuming that there are three positive roots (3.11), defined as $x_1^{\text{*}},x_2^{\text{*}},x_3^{\text{*}}$ , Equation (3.12) has three positive roots ${\omega }_1=\sqrt{x_1^{\text{*}}}$ , ${\omega }_2=\sqrt{x_2^{\text{*}}}$ , ${\omega }_3=\sqrt{x_3^{\text{*}}}$ .

From (3.9)

$\cos \omega \tau =\frac{b_{11}\omega \left({\omega }^3-a_{12}\omega \right)+b_{12}\left(a_{11}{\omega }^2-a_{13}\right)}{b_{11}^2{\omega }^2+b_{12}^2},$

let

${\tau }_k^{\left(j\right)}=\frac{1}{{\omega }_k}\text{arccos}\left[\frac{b_{11}{\omega }_k\left({\omega }_k^3-a_{12}{\omega }_k\right)+b_{12}\left(a_{11}{\omega }_k^2-a_{13}\right)}{b_{11}^2{\omega }_k^2+b_{12}^2}+2j\pi \right],$

where $k=1,2,3$ ; $j=0,1,2,\cdots$ . Therefore, $\pm \text{i}{\omega }_k$ is a pair of pure imaginary roots of Equation (3.8) at $\tau ={\tau }_k^{\left(j\right)}$ .

Define ${\tau }_0={\tau }_k^{\left(0\right)}=\underset{k=1,2,3}{\text{min}}{\tau }_k^{\left(0\right)}$ , ${\omega }_0={\omega }_{k_0}$ , i.e.

${\tau }_0=\text{min}\left\{\frac{1}{{\omega }_0}\arccos \frac{b_{11}{\omega }_0\left({\omega }_0^3-a_{12}{\omega }_0\right)+b_{12}\left(a_{11}{\omega }_0^2-a_{13}\right)}{b_{11}^2{\omega }_0^2+b_{12}^2}\right\}$ (3.14)

Suppose $\lambda \left(\tau \right)=\alpha \left(\tau \right)+\text{i}\omega \left(\tau \right)$ is the eigenvalue of the system (2.2) at $E_1$ , i.e., there exists a root $\lambda \left(\tau \right)$ of Equation (3.8), and satisfies $\alpha \left({\tau }_0\right)=0$ , $\omega \left({\tau }_0\right)={\omega }_0$ .

The derivation of the characteristic Equation (3.8) with respect to $\tau$ gives

$\frac{\text{d}\lambda }{\text{d}\tau }=\frac{\lambda \left(b_{11}\lambda +b_{12}\right){\text{e}}^{-\lambda \tau }}{\left(3{\lambda }^2+2a_{11}\lambda +a_{12}\right)+b_{11}{\text{e}}^{-\lambda \tau }-\tau {\text{e}}^{-\lambda \tau }\left(b_{11}\lambda +b_{12}\right)},$

from the characteristic Equation (3.8), it can be solved that

${\text{e}}^{-\lambda \tau }=-\frac{{\lambda }^3+a_{11}{\lambda }^2+a_{12}\lambda +a_{13}}{b_{11}\lambda +b_{12}},$

bringing this into $\frac{\text{d}\lambda }{\text{d}\tau }$ , yields

${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=-\frac{3{\lambda }^2+2a_{11}\lambda +a_{12}}{\lambda \left({\lambda }^3+a_{11}{\lambda }^2+a_{12}\lambda +a_{13}\right)}+\frac{b_{11}}{\lambda \left(b_{11}\lambda +b_{12}\right)}-\frac{\tau }{\lambda },$

the two equations of Equation (3.9) are squared and then added together to satisfy the following equation at ${\omega }_0$ .

${\left({\omega }_0^3-a_{12}{\omega }_0\right)}^2+{\left(a_{11}{\omega }_0^2-a_{13}\right)}^2=b_{11}^2{\omega }_0^2+b_{12}^2,$

solving yields the real part of ${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}$ at $\tau ={\tau }_0$ asSolving yields the real part of ${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}$ at $\tau ={\tau }_0$ as

$\begin{array}{c}{{\left(\frac{\text{d}\left(\mathrm{Re}\lambda \right)}{\text{d}\tau }\right)}^{-1}|}_{\tau ={\tau }_0}={\text{Re}{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}|}_{\tau ={\tau }_0}\\ =\text{Re}\left({-\frac{3{\lambda }^2+2a_{11}\lambda +a_{12}}{\lambda \left({\lambda }^3+a_{11}{\lambda }^2+a_{12}\lambda +a_{13}\right)}|}_{\lambda ={\omega }_0}\right)\\ +\text{Re}\left({-\frac{b_{11}}{\lambda \left(b_{11}\lambda +b_{12}\right)}|}_{\lambda ={\omega }_0}\right)\\ =\frac{3{\omega }_0^4+2\left(a_{11}-2a_{12}\right){\omega }_0^2+a_{12}^2-2a_{11}{a}_{13}-b_{11}^2}{b_{11}^2{\omega }_0^2+b_{12}^2}\\ =\frac{g^{\prime }\left({\omega }_0^2\right)}{b_{11}^2{\omega }_0^2+b_{12}^2},\end{array}$

since $b_{11}^2{\omega }_0^2+b_{12}^2>0$ , therefore

${\text{sign}\left(\frac{\text{d}\left(\mathrm{Re}\lambda \right)}{\text{d}\tau }\right)|}_{\tau ={\tau }_0}={\text{sign}{\left(\frac{\text{d}\left(\mathrm{Re}\lambda \right)}{\text{d}\tau }\right)}^{-1}|}_{\tau ={\tau }_0}=\text{sign}\left(g^{\prime }\left({\omega }_0^2\right)\right).$

Thus, if $g^{\prime }\left({\omega }_0^2\right)\ne 0$ , then ${\frac{\text{d}\left(\mathrm{Re}\lambda \right)}{\text{d}\tau }|}_{\tau ={\tau }_0}\ne 0$ , Hopf’s bifurcation point at $\tau ={\tau }_0$ at which the transversality condition is satisfied. Based on the above discussion and literature [37] [38], we arrive at the following result.

Theorem 3.3 Let ${\tau }_0,{\omega }_0$ be given by (3.14), and the following conclusions about the stability of model (2.2) at $E_1$ are drawn

1) When $R_1<1<R_0$ and $g^{\prime }\left({\omega }_0^2\right)\ne 0$ holds, then there are

a) When $g_3\ge 0$ and $\text{Δ}\le 0$ , then for any $\tau \ge 0$ , the equilibrium point $E_1$ is locally asymptotically stable.

b) Assuming that $g_3\ge 0$ , $\text{Δ}>0$ , $x_1>0$ and $g\left(x_1\right)\le 0$ holds, the equilibrium point $E_1$ is locally asymptotically stable when $\tau \in \left[0,{\tau }_0\right)$ ; the model (2.2) generates a Hopf bifurcation at $E_1$ when $\tau ={\tau }_0$ ; $E_1$ is unstable when $\tau >{\tau }_0$ .

2) When $R_1>1$ , then the equilibrium point $E_1$ is unstable for any $\tau \ge 0$ .

In the following we discuss the stability of model (2.2) at the immune activation equilibrium $E^{\text{*}}\left(u^{\text{*}},v^{\text{*}},w^{\text{*}},z^{\text{*}}\right)$ , which, from Section 2, exists for $E^{\text{*}}$ when $R_1>1$ .

The Jacobi matrix of model (2.2) at $E^{\text{*}}$ is

${J|}_{E^{\text{*}}}=\left(\begin{array}{cccc}-\left(h+w^{\text{*}}\right)& 0& -u^{\text{*}}& 0\\ f{w}^{\text{*}}{\text{e}}^{-\lambda \tau }& -1& 1+f{u}^{\text{*}}{\text{e}}^{-\lambda \tau }& 0\\ 0& l& -m-z^{\text{*}}& -w^{\text{*}}\\ 0& 0& n{z}^{\text{*}}& n{w}^{\text{*}}-k\end{array}\right),$

the corresponding characteristic equation is

$\begin{array}{c}{h\left(\lambda ,\tau \right)|}_{E^{\text{*}}}=\left|{\lambda I-J|}_{E^{\text{*}}}\right|\\ ={\lambda }^4+c_1{\lambda }^3+c_2{\lambda }^2+c_3\lambda +c_4+\left(d_1{\lambda }^2+d_2\lambda \right){\text{e}}^{-\lambda \tau }\\ =0,\end{array}$ (3.15)

where

$\begin{array}{l}{c}_1=h+w^{\text{*}}+m+z^{\text{*}}+1,\\ c_2=h+w^{\text{*}}+m+z^{\text{*}}+\left(h+w^{\text{*}}\right)\left(m+z^{\text{*}}\right)+n{z}^{\text{*}}{w}^{\text{*}}-l,\\ c_3=\left(h+w^{\text{*}}\right)\left(m+z^{\text{*}}\right)+n{z}^{\text{*}}{w}^{\text{*}}+n{z}^{\text{*}}{w}^{\text{*}}\left(h+w^{\text{*}}\right)-l\left(h+w^{\text{*}}\right),\\ c_4=n{z}^{\text{*}}{w}^{\text{*}}\left(h+w^{\text{*}}\right),\\ d_1=-lf{u}^{\text{*}},\\ d_2=-lfh{u}^{\text{*}}.\end{array}$

When $\tau =0$ , Equation (3.15) becomes

$\begin{array}{l}\lambda \left(\lambda +h+w^{\text{*}}\right)\left(\lambda +1\right)\left(\lambda +m+z^{\text{*}}\right)+lf{u}^{\text{*}}{w}^{\text{*}}\lambda +\left(\lambda +h+w^{\text{*}}\right)\left(\lambda +1\right)k{z}^{\text{*}}\\ =\lambda \left(\lambda +h+w^{\text{*}}\right)l\left(f{u}^{\text{*}}+1\right)=\lambda \left(\lambda +h+w^{\text{*}}\right)\left(m+z^{\text{*}}\right).\end{array}$ (3.16)

Now, we require that all roots of (3.16) have negative real parts. Otherwise, (3.16) has at least one root ${\lambda }_1=x_1+\text{i}{y}_1$ , where $x_1\ge 0$ . Notice that

$\begin{array}{l}\left|{\lambda }_1\left({\lambda }_1+h+w^{\text{*}}\right)\left({\lambda }_1+1\right)\left({\lambda }_1+m+z^{\text{*}}\right)+lf{u}^{\text{*}}{w}^{\text{*}}{\lambda }_1+\left({\lambda }_1+h+w^{\text{*}}\right)\left({\lambda }_1+1\right)k{z}^{\text{*}}\right|\\ >\left|{\lambda }_1\left({\lambda }_1+h+w^{\text{*}}\right)\left(m+z^{\text{*}}\right)\right|,\end{array}$

which contradicts (3.16), Thus, all roots of (3.16) have negative real parts. Hence, when $\tau =0$ , the equilibrium point $E^{\text{*}}$ is locally asymptotically stable.

Consider the case where $\tau >0$ . Assuming that $\text{i}\omega$ is a root of Equation (3.15), substituting it into Equation (3.15) yields

${\omega }^4-c_1\text{i}{\omega }^3-c_2{\omega }^2+c_3\text{i}\omega +c_4+\left(-d_1{\omega }^2+d_2\text{i}\omega \right)\left(\text{cos}\omega \tau -\text{i}\sin \omega \tau \right)=0,$

separate the real and imaginary parts of the pair of superscripts

$\{\begin{array}{l}{\omega }^4-c_2{\omega }^2+c_4=d_1{\omega }^2\cos \omega \tau -d_2\omega \sin \omega \tau ,\\ c_1{\omega }^3-c_3\omega =d_2\omega \cos \omega \tau +d_1{\omega }^2\sin \omega \tau .\end{array}$ (3.17)

Squaring both ends of (3.17) and adding them together yields

${\left({\omega }^4-c_2{\omega }^2+c_4\right)}^2+{\left(c_1{\omega }^3-c_3\omega \right)}^2={\left(d_1{\omega }^2\right)}^2+{\left(d_2\omega \right)}^2,$ (3.18)

i.e.

${\omega }^8+s_1{\omega }^6+s_2{\omega }^4+s_3{\omega }^2+s_4=0,$ (3.19)

where

$\begin{array}{l}{s}_1=c_1^2-2c_2,s_2=2c_4+c_2^2-2c_1{c}_3-d_1^2,\\ s_3=c_3^2-2c_2{c}_4-d_2^2,s_4=c_4^2.\end{array}$

Let $r={\omega }^2$ , (3.19) become

$s\left(r\right)\triangleq r^4+s_1{r}^3+s_2{r}^2+s_{3}r+s_4=0,$ (3.20)

then

$s^{\prime }\left(r\right)=4r^3+3s_1{r}^2+2s_{2}r+s_3.$ (3.21)

Define

$P_0=\frac{8s_2-3s_1^2}{16},Q_0=\frac{s_1^3-4s_1{s}_2+8s_3}{32},D_0=\frac{Q_0^2}{4}+\frac{P_0^3}{27},\sigma =\frac{-1+\sqrt{3}\text{i}}{2}.$

From Kadam’s formula, the largest real root of Equation (3.21) is found in the following cases

when $D_0>0$ ,

$r_1^{\text{*}}=-\frac{s_1}{4}+\sqrt[3]{-\frac{Q_0}{2}+\sqrt{D_0}}+\sqrt[3]{-\frac{Q_0}{2}-\sqrt{D_0}};$

when $D_0=0$ ,

$r_2^{\text{*}}=\text{max}\left\{-\frac{s_1}{4}-2\sqrt[3]{-\frac{Q_0}{2}},-\frac{s_1}{4}+\sqrt[3]{-\frac{Q_0}{2}}\right\};$

when $D_0<0$ ,

$r_3^{*}=\max \left\{-\frac{s_1}{4}+2\mathrm{Re}\left\{\xi \right\},-\frac{s_1}{4}+2\mathrm{Re}\left\{\xi \sigma \right\},-\frac{s_1}{4}+2\mathrm{Re}\left\{\xi \overline{\sigma }\right\}\right\},$

where $\xi =\sqrt[3]{-\frac{Q_0}{2}+\sqrt{D_0}}$ .

Since $s_4>0$ , therefore, similar to the discussion in literature [39], the following lemma about equation (3.20) can be obtained.

Lemma 3.3

1) If $s_4\ge 0$ , then Equation (3.20) has no positive real roots when one of the following conditions hold

a) $D_0>0$ , $r_1^{\text{*}}<0$ ;

b) $D_0=0$ , $r_2^{\text{*}}<0$ ;

c) $D_0<0$ , $r_3^{\text{*}}<0$ .

2) If $s_4\ge 0$ , then there exists at least one positive real root of Equation (3.20) when one of the following conditions hold

a) $D_0>0$ , $r_1^{\text{*}}>0$ , $s\left(r_1^{\text{*}}\right)<0$ ;

b) $D_0=0$ , $r_2^{\text{*}}>0$ , $s\left(r_2^{\text{*}}\right)<0$ ;

c) $D_0<0$ , $r_3^{\text{*}}>0$ , $s\left(r_3^{\text{*}}\right)<0$ .

Without loss of generality, we let Equation (3.20) have four positive roots defined as $r_1,r_2,r_3,r_4$ , then Equation (3.19) have four positive real roots ${\omega }_1=\sqrt{r_1}$ ,${\omega }_2=\sqrt{r_2}$ , ${\omega }_3=\sqrt{r_3}$ , ${\omega }_4=\sqrt{r_4}$ . From (3.17) we get

$\cos \omega \tau =\frac{d_1\left({\omega }^4-c_2{\omega }^2+c_4\right)+d_2\left(c_1{\omega }^2-c_3\right)}{d_1^2{\omega }^2+d_2^2},$

solving for $\tau$ the expression is

${\tau }_k^{\left(j\right)}=\frac{1}{{\omega }_k}\text{arccos}\left[\frac{d_1\left({\omega }^4-c_2{\omega }^2+c_4\right)+d_2\left(c_1{\omega }^2-c_3\right)}{d_1^2{\omega }^2+d_2^2}+2j\pi \right],$

where $k=1,2,3,4$ ; $j=0,1,2,\cdots$ . Thus, $\pm \text{i}{\omega }_k$ are a pair of pure imaginary roots of Equation (3.15) at $\tau ={\tau }_k^{\left(j\right)}$ .

Define ${\tau }_0={\tau }_k^{\left(0\right)}=\underset{k=1,2,3,4}{\text{min}}{\tau }_k^{\left(0\right)}$ , ${\omega }_0={\omega }_{k_0}$ , i.e.

${\tau }_0={\tau }_k^{\left(0\right)}=\text{min}\left\{\frac{1}{{\omega }_k}\arccos \frac{d_1\left({\omega }^4-c_2{\omega }^2+c_4\right)+d_2\left(c_1{\omega }^2-c_3\right)}{d_1^2{\omega }^2+d_2^2}\right\}.$ (3.22)

Suppose $\lambda \left(\tau \right)=\alpha \left(\tau \right)+\text{i}\omega \left(\tau \right)$ is the eigenvalue of the model (2.2) at $E^{\text{*}}$ , i.e., there exists a root $\lambda \left(\tau \right)$ of Equation (3.15), and satisfies $\alpha \left({\tau }_0\right)=0$ ,$\omega \left({\tau }_0\right)={\omega }_0$ .

The derivation of the characteristic Equation (3.15) with respect to $\tau$ yields

$\frac{\text{d}\lambda }{\text{d}\tau }=\frac{\lambda \left(d_1{\lambda }^2+d_2\lambda \right){\text{e}}^{-\lambda \tau }}{\left(4{\lambda }^3+3c_1{\lambda }^2+2c_2\lambda +c_3\right)+\left(2d_1\lambda +d_2\right){\text{e}}^{-\lambda \tau }-\tau {\text{e}}^{-\lambda \tau }\left(d_1{\lambda }^2+d_2\lambda \right)},$

from the characteristic Equation (3.15), we can solve for

${\text{e}}^{-\lambda \tau }=-\frac{{\lambda }^4+c_1{\lambda }^3+c_2{\lambda }^2+c_3\lambda +c_4}{d_1{\lambda }^2+d_2\lambda },$

bringing this into $\frac{\text{d}\lambda }{\text{d}\tau }$ , we get

${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=-\frac{4{\lambda }^3+3c_1{\lambda }^2+2c_2\lambda +c_3}{\lambda \left({\lambda }^4+c_1{\lambda }^3+c_2{\lambda }^2+c_3\lambda +c_4\right)}+\frac{2d_1\lambda +d_2}{\lambda \left(d_1{\lambda }^2+d_2\lambda \right)}-\frac{\tau }{\lambda },$

combining with (3.18), we solve for the real part of ${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}$ at $\tau ={\tau }_0$ as

$\begin{array}{l}{{\left(\frac{\text{d}\left(\mathrm{Re}\lambda \right)}{\text{d}\tau }\right)}^{-1}|}_{\tau ={\tau }_0}={\text{Re}{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}|}_{\tau ={\tau }_0}\\ =\text{Re}\left({-\frac{4{\lambda }^3+3c_1{\lambda }^2+2c_2\lambda +c_3}{\lambda \left({\lambda }^4+c_1{\lambda }^3+c_2{\lambda }^2+c_3\lambda +c_4\right)}|}_{\lambda ={\omega }_0}\right)\\ +\text{Re}\left({-\frac{2d_1\lambda +d_2}{\lambda \left(d_1{\lambda }^2+d_2\lambda \right)}-\frac{\tau }{\lambda }|}_{\lambda ={\omega }_0}\right)\\ =\frac{4{\omega }_0^6+3\left(c_1^2-2c_2\right){\omega }_0^4+2\left(2c_4+c_2^2-2c_1{c}_3-d_1^2\right){\omega }_0^2+\left(c_3^2-2c_2{c}_4-d_2^2\right)}{d_1^2{\omega }_0^4+d_2^2{\omega }_0^2}\\ =\frac{s^{\prime }\left({\omega }_0^2\right)}{d_1^2{\omega }_0^4+d_2^2{\omega }_0^2},\end{array}$

since $d_1^2{\omega }_0^4+d_2^2{\omega }_0^2>0$ , therefore

${\text{sign}\left(\frac{\text{d}\left(\mathrm{Re}\lambda \right)}{\text{d}\tau }\right)|}_{\tau ={\tau }_0}={\text{sign}{\left(\frac{\text{d}\left(\mathrm{Re}\lambda \right)}{\text{d}\tau }\right)}^{-1}|}_{\tau ={\tau }_0}=\text{sign}\left(s^{\prime }\left({\omega }_0^2\right)\right).$

Thus, if $s^{\prime }\left({\omega }_0^2\right)\ne 0$ , then ${\frac{\text{d}\left(\mathrm{Re}\lambda \right)}{\text{d}\tau }|}_{\tau ={\tau }_0}\ne 0$ , the Hopf bifurcation at $\tau ={\tau }_0$ the transversality condition is satisfied. Based on the above discussion, we arrive at the following result.

Theorem 3.4 When $R_1>1$ holds, ${\tau }_0,{\omega }_0$ defined by (3.22), then there are:

1) If Equation (3.20) has no positive real roots, $E^{\text{*}}$ is locally asymptotically stable for all $\tau \ge 0$ .

2) If there exists at least one positive real root of Equation (3.20) with $s^{\prime }\left({\omega }_0^2\right)\ne 0$ , then $E^{\text{*}}$ is locally asymptotically stable when $0\le \tau <{\tau }_0$ ; $E^{\text{*}}$ is unstable when $\tau >{\tau }_0$ ; when $\tau ={\tau }_0$ , the model (2.2) generates a Hopf bifurcation at $E^{\text{*}}$ .

4. Numerical Simulation

In this section, we select several sets of parameters to numerically simulate the above theoretical results.

1) Two sets of parameters are selected to verify the stability of $E_0$ ,

a) selection

$h=0.14,f=0.05,l=0.04,m=0.36,n=1.2,k=1,$

obtain $R_0=0.1508<1$ , $E_0$ is locally asymptotically stable for all $\tau \ge 0$ .

b) selection

$h=0.14,f=0.15,l=0.08,m=0.036,n=1.2,k=0.8,$

Obtain $R_0=4.6032>1$ , $E_0$ is unstable for all $\tau =0$ . As shown in Figure 1 and Figure 2.

Figure 1. $R_0=0.1508<1$ , $E_0$ is locally asymptotically stable for all $\tau \ge 0$ , choose $\tau =1$ .

Figure 2. $R_0=4.6032>1$ , $E_0$ is unstable for all $\tau =0$ .

2) Selection of parameters a)

$h=0.27,f=0.15,l=0.34,m=0.36,n=0.03,k=0.18,$

obtain

$\begin{array}{l}{R}_0=1.4691>1>R_1=0.9670,\\ E_1=\left(0.3924,2.4141,2.2800,0\right),\end{array}$

Equation (3.7) has a negative real root ${\lambda }_1=-0.1116$ .

$\begin{array}{l}{g}_3=0.0016>0,\\ \text{Δ}=4g_1^2-12g_2=-3.9447<0,\end{array}$

therefore, Equation (3.11) has positive roots, according to theorem 3.3, When $g_3\ge 0$ and $\text{Δ}\le 0$ , then for any $\tau \ge 0$ , the equilibrium point $E_1$ is locally asymptotically stable. As shown in Figure 3.

Figure 3. $R_0=1.4691>1>R_1=0.9670$ , $E_1=\left(0.3924,2.4141,2.2800,0\right)$ , $E_1$ is locally asymptotically stable for all $\tau \ge 0$ , choose $\tau =10$ .

b)

$h=0.37,f=0.9,l=0.7,m=0.8,n=0.3,k=0.18,$

obtain

$\begin{array}{l}{R}_0=3.0034>1,\\ R_1=1.6869>1,\\ E_1=\left(0.1587,6.7771,5.9300,0\right),\end{array}$

Equation (3.7) has a positive real root ${\lambda }_1=1.5990>0$ . Then for any $\tau \ge 0$ , the equilibrium point $E_1$ is unstable. As shown in Figure 4.

Figure 4. $R_0=3.0034>1$ , $R_1=1.6869>1$ , $E_1=\left(0.1587,6.7771,5.9300,0\right)$ , ${\lambda }_1=1.5990>0$ . Then for any $\tau \ge 0$ , the equilibrium point $E_1$ is unstable, choose $\tau =50$ .

Figure 5. $R_1=2.4298>1$ , $E^{\text{*}}$ is locally asymptotically stable for all $\tau \ge 0$ , choose $\tau =2$ .

3) Selection of parameters

$h=0.34,f=0.05,l=0.84,m=0.36,n=0.92,k=0.8,$

obtain $R_1=2.4298>1$ , and Equation (3.20) has no positive real roots, then $E^{\text{*}}$ is locally asymptotically stable for all $\tau \ge 0$ . As shown in Figure 5.

5. Conclusions

In this paper, we study a class of HTLV-I infection models with delayed viral replication and CTL immune response. Through rigorous mathematical analysis, we establish the threshold dynamics of the model, which can be determined by the immune inactivation reproduction ratio $R_0$ and the immune activation reproduction ratio $R_1$ . Firstly, the model is dimensionless, the positive invariance and boundedness of the model (2.2) are proved, the existence of the equilibrium point of the model is determined by $R_0$ and $R_1$ , and secondly, the stability of the equilibrium point of the model as well as the existence of Hopf bifurcation are discussed. It can be seen that the time lag does not affect the stability of the disease-free equilibrium point $E_0$ , but changes the stability of the immune inactivation equilibrium point $E_1$ and the immune activation equilibrium point $E^{\text{*}}$ . If $R_0<1$ , $E_0$ is locally asymptotically stable. If $R_0>1$ , the immune inactivation equilibrium point $E_1$ exists, and when $R_1<1<R_0$ , certain conditions are satisfied, and the model (2.2) undergoes a Hopf bifurcation at $E_1$ , and when $R_1>1$ , $E_1$ is unstable for any $\tau \ge 0$ . When $R_1>1$ , immune activation equilibrium $E^{\text{*}}$ exists, if $\tau =0$ immune activation equilibrium $E^{\text{*}}$ is locally asymptotically stabilized, if $\tau >0$ , take $\tau$ as the number of bifurcations, when the condition $\tau ={\tau }_0$ is satisfied, the system undergoes Hopf in $E^{\text{*}}$ . When $\tau <{\tau }_0$ , the equilibrium point $E^{\text{*}}$ is locally asymptotically stable, and when $\tau >{\tau }_0$ , the equilibrium point $E^{\text{*}}$ is unstable.

$R_0$ is an indicator of the intrinsic transmissibility of a virus, determining whether an infection can establish itself. $R_1$ is an indicator of the immune-virus equilibrium, determining whether an infection can be controlled. In the HTLV-I infection model, when the time delay ($\tau$ ) is selected as the bifurcation parameter, the stability changes, instability, and the emergence of Hopf bifurcations at the equilibrium points of immune inactivation (e.g., no CTL response) and immune activation (presence of CTL response) have the following biological significance:

1) Immune inactivation equilibrium point (no CTL response, $z=0$ )

If the delay $\tau$ is small, the model equilibrium point tends to be stable, indicating the absence of CTL response, leading to uncontrolled viral proliferation (simulating asymptomatic carriers or the latent period); immune tolerance: the body fails to effectively activate CTL, resulting in the persistent presence of infected cells $\left(v,w\right)$ . Instability (increasing $\tau$ ): Viral dynamic imbalance: When the time delay $\tau$ exceeds the critical value, stability disruption may lead to: a) Sudden increase in viral load: Explosive growth of actively infected cells ($w$ ) (e.g., entering the pre-leukemia stage); b) Activation of latent infection: Delayed $\tau$ reflects the efficiency of latent cells ($v$ ) converting to active infection ($w$ ); increased $\tau$ may trigger viral reactivation. Hopf bifurcation leads to periodic oscillations: Without CTL ($z=0$ ), oscillations caused by time delay reflect intermittent control by endogenous immunity (e.g., non-specific immunity).

2) Immune activation equilibrium point (CTL response present, $z>0$ )

Stability is indicated by: a) Immune control: If $\tau$ is small, CTL ($z$ ) effectively suppresses infected cells $\left(w,v\right)$ , and the system tends toward a stable state with low viral load (simulating long-term asymptomatic control); b) Balanced immune response: CTL proliferation ($n$ ) and decay ($k$ ) reach a dynamic equilibrium, avoiding excessive immune damage. Instability (increased $\tau$ ) indicates immune escape. An increased time delay $\tau$ may disrupt stability, leading to: a) Delayed CTL response: $\tau$ reflects the delay in CTL activation; an excessively large $\tau$ allows infected cells ($w$ ) to proliferate before CTL becomes effective; b) Disease progression: the system tends toward a high viral load state (e.g., ATL progression). Hopf bifurcation indicates periodic immune oscillations: a tug-of-war between the immune system and the virus. The periodic solution resulting from Hopf bifurcation corresponds to the reciprocal increase and decrease of CTL ($z$ ) and infected cells $\left(w,v\right)$ , simulating recurrent chronic inflammation (such as neuroinflammation in HAM/TSP).

3) Time delay $\tau$ is a key indicator of immune response efficiency: a smaller $\tau$ facilitates rapid immune control of the virus, while a larger $\tau$ may lead to immune failure or periodic pathological damage. Hopf bifurcation oscillations: These reflect the dynamic fluctuations in viral load and immune response observed clinically, providing a theoretical basis for explaining the diversity of HTLV-I infection (e.g., latency period, HAM/TSP, ATL). Therapeutic targets: Reducing latent cells (lowering $\tau$ ) or enhancing CTL efficacy (e.g., checkpoint inhibitors) through drug therapy may improve prognosis [32].

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Wasserstein-Robbins, F. (2010) A Mathematical Model of HIV Infection: Simulating T4, T8, Macrophages, Antibody, and Virus via Specific Anti-HIV Response in the Presence of Adaptation and Tropism. Bulletin of Mathematical Biology, 72, 1208-1253.[CrossRef] [PubMed]
[2]	Wang, K., Fan, A. and Torres, A. (2010) Global Properties of an Improved Hepatitis B Virus Model. Nonlinear Analysis: Real World Applications, 11, 3131-3138.[CrossRef]
[3]	Elaiw, A.M. (2010) Global Properties of a Class of HIV Models. Nonlinear Analysis: Real World Applications, 11, 2253-2263.[CrossRef]
[4]	Vieira, I.T., Cheng, R.C.H., Harper, P.R. and de Senna, V. (2009) Small World Network Models of the Dynamics of HIV Infection. Annals of Operations Research, 178, 173-200.[CrossRef]
[5]	Yu, Y., Nieto, J.J., Torres, A. and Wang, K. (2009) A Viral Infection Model with a Nonlinear Infection Rate. Boundary Value Problems, 2009, 1-19.[CrossRef]
[6]	Li, M.Y. and Shu, H. (2012) Global Dynamics of a Mathematical Model for HTLV-I Infection of CD4+ T Cells with Delayed CTL Response. Nonlinear Analysis: Real World Applications, 13, 1080-1092.[CrossRef]
[7]	Osame, M., Usuku, K. and Izumo, S. (1986) HTLV-I Associated Myelopathy, a New Clinical Entity. The Lancet, 327, 1031-1032.[CrossRef] [PubMed]
[8]	Osame, M., Janssen, R., Kubota, H., Nishitani, H., Igata, A., Nagataki, S., et al. (1990) Nationwide Survey of HTLV-I-Associated Myelopathy in Japan: Association with Blood Transfusion. Annals of Neurology, 28, 50-56.[CrossRef] [PubMed]
[9]	Gallo, R.C. (2005) History of the Discoveries of the First Human Retroviruses: HTLV-1 and HTLV-2. Oncogene, 24, 5926-5930.[CrossRef] [PubMed]
[10]	Kubota, R. (2000) Retrovirus: Human T-Cell Lymphotropic Virus Type I-Associated Diseases and Immune Dysfunction. Effects of Microbes on the Immune System, 349-371.
[11]	Coffin, J.M., Hughes, S.H. and Varmus, H.E. (1997) Retroviruses. Cold Spring Harbor Laboratory Press.
[12]	Bangham, C.R.M. (2000) The Immune Response to HTLV-I. Current Opinion in Immunology, 12, 397-402.[CrossRef] [PubMed]
[13]	Bangham, C.R.M. (2003) The Immune Control and Cell-to-Cell Spread of Human T-Lymphotropic Virus Type 1. Journal of General Virology, 84, 3177-3189.[CrossRef] [PubMed]
[14]	Asquith, B. and Bangham, C.R.M. (2007) Quantifying HTLV-I Dynamics. Immunology & Cell Biology, 85, 280-286.[CrossRef] [PubMed]
[15]	Bangham, C.R.M. (2009) CTL Quality and the Control of Human Retroviral Infections. European Journal of Immunology, 39, 1700-1712.[CrossRef] [PubMed]
[16]	Wattel, E., Cavrois, M., Gessain, A. and Wain-Hobson, S. (1996) Clonal Expansion of Infected Cells: A Way of Life for HTLV-I. Journal of Acquired Immune Deficiency Syndromes and Human Retrovirology, 13, S92-S99.[CrossRef] [PubMed]
[17]	Jacobson, S. (2002) Immunopathogenesis of Human T Cell Lymphotropic Virus Type I–Associated Neurologic Disease. The Journal of Infectious Diseases, 186, S187-S192.[CrossRef] [PubMed]
[18]	Eshima, N., Tabata, M., Okada, T., et al. (2003) Population Dynamics of HTLV-I Infection: A Discrete-Time Mathematical Epidemic Model Approach. Mathematical Medicine and Biology, 20, 29-45.[CrossRef]
[19]	Gómez-Acevedo, H., Li, M.Y. and Jacobson, S. (2010) Multistability in a Model for CTL Response to HTLV-I Infection and Its Implications to HAM/TSP Development and Prevention. Bulletin of Mathematical Biology, 72, 681-696.[CrossRef] [PubMed]
[20]	Qi, K., Jiang, D., Hayat, T. and Alsaedi, A. (2019) The Stationary Distribution and Extinction of a Double Thresholds HTLV-I Infection Model with Nonlinear CTL Immune Response Disturbed by White Noise. International Journal of Biomathematics, 12, Article 1950058.[CrossRef]
[21]	Kuang, D.P., Yin, Q. and Li, J.L. (2021) Dynamics of Stochastic HTLV‐I Infection Model with Nonlinear CTL Immune Response. Mathematical Methods in the Applied Sciences, 44, 14059-14078.[CrossRef]
[22]	Chen, S.Y., Liu, Z., Wang, L., et al. (2024) Global Dynamics Analysis for a Nonlinear HTLV-I Model with Logistic Proliferation and CTL Response. International Journal of Biomathematics, 17, Article 2350023.[CrossRef]
[23]	Song, X.Y., Wang, S.L. and Dong, J. (2011) Stability Properties and Hopf Bifurcation of a Delayed Viral Infection Model with Lytic Immune Response. Journal of Mathematical Analysis and Applications, 373, 345-355.[CrossRef]
[24]	Shamsara, E., Afsharnezhad, Z. and Javidmanesh, E. (2018) Graphical Hopf Bifurcation of a Filippov HTLV-1 Model with Delay in Cytotoxic T Cells Response. Journal of Dynamic Systems, Measurement, and Control, 140, Article 091007.[CrossRef]
[25]	Jia, X.J. and Xu, R. (2022) Global Dynamics of a Delayed HTLV-I Infection Model with Beddington-Deangelis Incidence and Immune Impairment. Chaos, Solitons & Fractals, 155, Article 111733.[CrossRef]
[26]	Chen, S.Y., Liu, Z., Wang, L. and Zhang, X. (2024) Stability and Hopf Bifurcation Analysis of a HTLV-Ⅰ Infection Model with Time-Delay CTL Immune Response. Discrete and Continuous Dynamical Systems B, 29, 812-832.[CrossRef]
[27]	Li, M.Y. and Lim, A.G. (2011) Modelling the Role of Tax Expression in HTLV-I Persistence in Vivo. Bulletin of Mathematical Biology, 73, 3008-3029.[CrossRef] [PubMed]
[28]	Lim, A.G. and Maini, P.K. (2014) HTLV-I Infection: A Dynamic Struggle between Viral Persistence and Host Immunity. Journal of Theoretical Biology, 352, 92-108.[CrossRef] [PubMed]
[29]	Li, S.M. and Zhou, Y.C. (2016) Backward Bifurcation of an HTLV-I Model with Immune Response. Discrete and Continuous Dynamical Systems, Series B, 21, 863-881.[CrossRef]
[30]	Khajanchi, S., Bera, S. and Roy, T.K. (2021) Mathematical Analysis of the Global Dynamics of a HTLV-I Infection Model, Considering the Role of Cytotoxic T-Lymphocytes. Mathematics and Computers in Simulation, 180, 354-378.[CrossRef]
[31]	Song, C.W. and Xu, R. (2021) Mathematical Analysis of an HTLV-I Infection Model with the Mitosis of CD4+ T Cells and Delayed CTL Immune Response. Nonlinear Analysis: Modelling and Control, 26, 1-20.[CrossRef]
[32]	Bera, S., Khajanchi, S. and Roy, T.K. (2022) Dynamics of an HTLV-I Infection Model with Delayed CTLs Immune Response. Applied Mathematics and Computation, 430, Article 127206.[CrossRef]
[33]	Xu, R. and Yang, Y. (2022) Dynamics of an HTLV-I Infection Model with Delayed and Saturated CTL Immune Response and Immune Impairment. Acta Mathematica Scientia, 42, 1836-1848.
[34]	Elaiw, A.M., Shflot, A.S. and Hobiny, A.D. (2022) Stability Analysis of General Delayed HTLV-I Dynamics Model with Mitosis and CTL Immunity. Mathematical Biosciences and Engineering, 19, 12693-12729.[CrossRef] [PubMed]
[35]	Katri, P. and Ruan, S. (2004) Dynamics of Human T-Cell Lymphotropic Virus I (HTLV-I) Infection of CD4+ T-Cells. Comptes Rendus. Biologies, 327, 1009-1016.[CrossRef] [PubMed]
[36]	Elaiw, A.M., Shflot, A.S. and Hobiny, A.D. (2023) Global Stability of a General HTLV-I Infection Model with Cytotoxic T-Lymphocyte Immune Response and Mitotic Transmission. Alexandria Engineering Journal, 67, 77-91.[CrossRef]
[37]	Ruan, S.G. and Wei, J.J. (2003) On the Zeros of Transcendental Functions with Applications to Stability of Delay Differential Equations with Two Delays. Dynamics of Continuous Discrete and Impulsive Systems Series A, 10, 863-874.
[38]	Culshaw, R.V. and Ruan, S. (2000) A Delay-Differential Equation Model of HIV Infection of CD4+ T-Cells. Mathematical Biosciences, 165, 27-39.[CrossRef] [PubMed]
[39]	Yan, X. and Li, W. (2006) Stability and Bifurcation in a Simplified Four-Neuron BAM Neural Network with Multiple Delays. Discrete Dynamics in Nature and Society, 2006, Article 32529.