Analysis of a Delayed Multiscale AIDS/HIV-1 Model Coupling Between-Host and Within-Host Dynamics

Abstract

Taking into account the effects of the immune response and delay, and complexity on HIV-1 transmission, a multiscale AIDS/HIV-1 model is formulated in this paper. The multiscale model is described by a within-host fast time model with intracellular delay and immune delay, and a between-host slow time model with latency delay. The dynamics of the fast time model is analyzed, and includes the stability of equilibria and properties of Hopf bifurcation. Further, for the coupled slow time model without an immune response, the basic reproduction number ${\mathcal{R}}_0^h$ is defined, which determines whether the model may have zero, one, or two positive equilibria under different conditions. This implies that the slow time model demonstrates more complex dynamic behaviors, including saddle-node bifurcation, backward bifurcation, and Hopf bifurcation. For the other case, that is, the coupled slow time model with an immune response, the threshold dynamics, based on the basic reproduction number ${\tilde{\mathcal{R}}}_0^h$, is rigorously investigated. More specifically, if ${\tilde{\mathcal{R}}}_0^h<1$, the disease-free equilibrium is globally asymptotically stable; if ${\tilde{\mathcal{R}}}_0^h>1$, the model exhibits a unique endemic equilibrium that is globally asymptotically stable. With regard to the coupled slow time model with an immune response and stable periodic solution, the basic reproduction number ${\mathcal{R}}_0$ is derived, which serves as a threshold value determining whether the disease will die out or lead to periodic oscillations in its prevalence. The research results suggest that the disease is more easily controlled when hosts have an extensive immune response and the time required for new immune particles to emerge in response to antigenic stimulation is within a certain range. Finally, numerical simulations are presented to validate the main results and provide some recommendations for controlling the spread of HIV-1.

1. Introduction

Acquired immune deficiency syndrome (AIDS), caused by infection with the human immunodeficiency virus (HIV), has been one of the deadliest and most persistent epidemics in humankind since the reporting of the first AIDS cases in 1981 [1]. It is well known that the common pathways for transmitting HIV include engaging in sexual intercourse with the infected individual, mother-to-child infection, sharing of injecting equipment, and exposure to contaminated blood in healthcare settings [2]. The signs and symptoms of HIV differ based on the stage of the infection. During the early stages of HIV infection, individuals may exhibit no symptoms or flu-like symptoms such as fever, fatigue, sore throat, swollen lymph nodes, and skin rashes. As the infection progresses to the chronic phase, there may be no clinical symptoms for several years, but the virus continues actively assaulting and debilitating the immune system during this stage. This is why HIV is often referred to as a silent infection. In the final stage of infection, individuals with HIV develop AIDS, wherein the immune system is severely impaired, resulting in a range of symptoms and opportunistic infections. As of the end of 2022, there were an estimated 39.0 million [33.1–45.7 million] people living with HIV, including 1.5 million [1.2–2.1 million] children (0–14 years old). HIV remains a significant global public health concern, having claimed the lives of an estimated 40.4 million [32.9–51.3 million] individuals to date [3]. This is a reminder of the need to take appropriate preventive measures for effective HIV control.

To gain a deeper understanding of AIDS/HIV mechanisms and develop more effective control strategies, numerous scholars have utilized mathematical models to investigate this disease. For example, McCluskey et al. [4], concerned that each HIV-infected individual has the potential to progress through multiple stages of infection before reaching the advanced stage of AIDS, proposed an HIV/AIDS model with multiple stages of infection, in which the infected individual is divided into different infectious groups. In the study conducted by [5], researchers focused on the varying transmission abilities of individuals in different stages of infection, and a nonautonomous stage-structured model was developed to simulate the spread of HIV/AIDS between hosts. In addition, mathematical modeling of HIV/AIDS based on within-host dynamics are increasingly being used to obtain a comprehensive view of the HIV-1 viral infection. As in the mid-1990s, an introductory model, including both healthy and infected cells as well as free virus particles, was developed by Perelson et al. [6] to study the virus–cell interactions, which laid a foundation for HIV-1 infection modeling. The immune response is the body’s defensive response to alien or mutated autologous components, which play a vital role in eliminating or controlling the disease during viral infection [7]. As a result, work on HIV-1 infection models that focus on the relationship between immune responses and invading pathogens has been developed and progressed considerably. For example, since the role of humoral immunity is crucial and cannot be replaced within the broader scope of human immunity [8,9], some scholars have incorporated humoral immunity into a basic HIV-1 virus model [10,11]. This is based on the B cells, which are produced by antigenic stimulation and are specifically programmed to kill viruses.

For many infectious diseases, such as tuberculosis, measles, AIDS, and so on, following sufficient exposure to the disease-causing agent, the susceptible individual becomes infected but is not yet capable of transmitting the infection to others, and the disease incubates within the body for a certain period of time until the exposed individual becomes infectious, and this time between infection and becoming infective is called the latent period. Capturing this feature, the introduction of the delay, $\tau$, or exposed class, E, in modeling the dynamics of infectious diseases has received widespread consideration [12,13,14]. In particular, Tipsri et al. [15] focused on examining a $SEIR$ (susceptible, exposed, infected, recovered) model that incorporated a saturated incidence rate and discrete delay in three different cases. Their theoretical analysis revealed that introducing a time delay in the transmission term could potentially disrupt the stability of the system, and periodic solutions could be generated through a phenomenon known as Hopf bifurcation when utilizing the delay as a bifurcation parameter. In [16], Jiang et al. proposed a $SEIRS$ model with two delays (where one delay is the length of the immunity period and the other delay is the latent period of the disease) and the general nonlinear incidence rate, to discuss the the impact of delay on the disease transmission dynamics.

Additionally, delay has also been incorporated into within-host virus models to investigate HIV-1 infection more accurately [17,18,19]. For instance, to characterize the time required for the infected cells to produce new virions after viral invasion and the time needed for the adaptive immune response to emerge to control viral replication, intracellular delay and immune delay were introduced by Pawelek et al. [20]. And their results suggest that incorporating the intracellular delay into the model does not affect the stability results, but introducing the immune delay generates diverse and complex dynamics. Lin et al. [21] extended the classical HIV model by adding cell-to-cell transmission and intracellular delay, and established the threshold dynamics that were determined by the immune-inactivated reproduction number, ${\mathcal{R}}_0$, and immune-activated reproduction number, ${\mathcal{R}}_1$. In fact, when HIV invades the body, it takes time for both the infection and the immune response to occur. Thus, Yang et al. [22] incorporated two discrete delays into an HIV model, respectively, modeling intracellular delay and immune delay. Conclusions show that immune delay, ${\tau }_2$, will change the positive equilibrium from stable to unstable and cause population oscillations.

As we all know, when the virus enters the human body, it replicates and survives with the help of the host cell, destroys host cellular tissues, and than induces an innate immune response. The complex interactions between the virus and host cells directly affect the infection severity of the host. Then, the susceptible individual becomes infected after close contact with the infected individual and the disease transmission between hosts begins. That is to say, two of the important processes in transmission of viral infectious diseases are the immunological process concerning the virus–cell interactions at the individual level and the epidemiological process of disease transmission at the population level [23]. Yet, traditional infectious diseases models are often focused on one of these two processes. In fact, numerous researchers have shown that the higher the viral load in the host is, the greater is ability of the disease to spread from host to host [24,25], implying that the within-host and between-host dynamics are not independent but interrelated. Initially, researchers introduced a nested modeling approach to simulate the evolutionary dynamics of host–parasite interactions [26,27,28]. Since then, several attempts have been made to integrate the dynamics of pathogen progression within a host with those of transmission between hosts in the context of infectious diseases. For instance, Feng et al. [29,30,31] proposed a coupled model that links the cell–parasite interaction and the between-host infection via the contaminated environment to analyze the transmission of Toxoplasma Gondii, and found that backward bifurcations can be generated by coupled models. Recently, in order to understand the connection between pathogens, hosts, and the environment better, Wang et al. [32] developed a multiscale model to analyze the dynamics of COVID-19, and assumed that the person-to-person transmission rates depend on the viral load within the human body, and the within-host viral load depends on the environmental pathogen concentration, which enables two-way coupling of micro within-host and macro between-host dynamics. It should be noted that, for some environment-driven infectious diseases, scholars have effectively integrated both micro and macro models bidirectionally [29,30,31,32]. However, for HIV, it remains unclear how to reflect the impact of host-to-host transmission on viral loads within the hosts, or how to build a bridge that can be coupled from epidemiological models to viral dynamics models. Under this incentive, Xue et al. [33] proposed a bidirectional coupling HIV-1 model with standard incidence, in which the entry of viruses from infected individuals into the within-host system is based on the fact that the virus can be transmitted through direct contact with HIV-positive individuals and AIDS patients.

The purpose of this paper, motivated by the above works, is to discuss the impact of immune response and delay on the multiscale HIV-1 model coupling within-host and between-host dynamics. The subsequent sections of this paper are organized as following. In Section 2, the virus infection model with humoral immunity, intracellular delay, and immune delay in the host, and the disease transmission model with latency delay between hosts, are proposed at fast and slow time scales, respectively. In Section 3, some properties and important results of the virus infection model are given. In Section 4, the positivity and boundedness of solutions for the coupled slow time model without an immune response are discussed, and criteria on the existence and local stability of equilibria are established. In Section 5, we further study the coupled slow time model with an immune response. The basic reproduction number is calculated, the existence and stability of equilibria are studied, and the uniform persistence of the model is obtained. In Section 6, the stability of the periodic solution for the coupled slow time model with an immune response is analyzed. The basic reproduction number is defined and the global dynamics are illustrated in terms of this value. In Section 7, numerical simulations are conducted to explain the main findings further, and a brief summary is provided in the last section.

2. Model Formulation

Following the transmission process of HIV-1 at the host level, the total number of hosts at time t is divided into four classes: susceptible, $S\left(t\right)$, exposed, $E\left(t\right)$, HIV-positive, $I\left(t\right)$ (individuals without clinical symptoms), and AIDS patients, $A\left(t\right)$ (individuals that have developed clinical symptoms). Inspired by the studies [32,33], it is assumed that the transmission rates and disease-induced mortality of infectious individuals depend on the viral load, V, within the host, thereby embedding micro virus dynamics into macro epidemic dynamics. Supposing that infectivity is saturated with an increasing number of infectious individuals, $\alpha$ represents the inhibitory effect exerted by the infected individuals, which determines the level of infectivity saturation. In addition, based on the transmission characteristics of HIV, assume that susceptible individuals first undergo a latent period after infection, which is described by a delay $\tau$, before becoming infected. From these, a between-host disease transmission slow time model reads

(1)$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}=& {\mathsf{\Lambda }}_h-\frac{c_1{\beta }_1\left(V\right)I\left(t\right)+c_2{\beta }_2\left(V\right)A\left(t\right)}{1+\alpha \left(I\right(t)+A(t\left)\right)}S\left(t\right)-\mu S\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}E\left(t\right)}{\mathrm{d}t}=& \frac{c_1{\beta }_1\left(V\right)I\left(t\right)+c_2{\beta }_2\left(V\right)A\left(t\right)}{1+\alpha \left(I\right(t)+A(t\left)\right)}S\left(t\right)-\mu E\left(t\right)\hfill \\ \hfill & -\frac{c_1{\beta }_1\left(V\right)I(t-\tau )+c_2{\beta }_2\left(V\right)A(t-\tau )}{1+\alpha \left(I\right(t-\tau )+A(t-\tau \left)\right)}S(t-\tau )e^{-\mu \tau },\hfill \\ \hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=& \frac{c_1{\beta }_1\left(V\right)I(t-\tau )+c_2{\beta }_2\left(V\right)A(t-\tau )}{1+\alpha \left(I\right(t-\tau )+A(t-\tau \left)\right)}S(t-\tau )e^{-\mu \tau }-(\mu +\xi +{\delta }_1\left(V\right))I\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}=& \xi I\left(t\right)-(\mu +{\delta }_2\left(V\right))A\left(t\right).\hfill \end{array}\right.$

All the parameters in Model (1) are nonnegative constants, and their biological interpretations are provided in Table 1.

Remark 1.

The transmission rate ${\beta }_i\left(V\right)$$(i=1,2)$ is assumed to satisfy ${\beta }_i\left(0\right)=0$, ${\beta }_i\left(V\right)\ge 0$, ${\beta }_i^{\prime }\left(V\right)\ge 0$, ${\beta }_i^{″}\left(V\right)\le 0$. One of the simplest function forms for ${\beta }_i\left(V\right)$ is ${\beta }_i\left(V\right)={\overline{a}}_{i}V$; other forms for ${\beta }_i\left(V\right)$ include ${\beta }_i\left(V\right)=\frac{{\overline{a}}_{i}V}{1+{\overline{b}}_{i}V}$ and ${\beta }_i\left(V\right)={\overline{a}}_i{V}^q(0<q<1)$, where ${\overline{a}}_i$ and ${\overline{b}}_i$ are positive constants.

Next, taking into account the infection process of HIV-1 at the individual host level, let $T\left(s\right)$, $T^{*}\left(s\right)$, $V\left(s\right)$, and $B\left(s\right)$ denote the densities of susceptible cells (healthy cells or uninfected cells), infected cells, viral load, and B cells at time s in the host, respectively. For HIV, in order to construct a bridge coupling the macro level to the micro scale, investigate the effect of disease transmission on the dynamics within the hosts, motivated by the works of [33], where we suppose that the viral load in the host can increase by directly contact with infected individuals, which is characterized by saturation incidence. Humoral immunity plays an essential role in the whole human immunity [34]. In fact, after the virus enters the human body, it does not immediately stimulate the production of B cells; therefore, based on [21], the immune response with time delay ${\tau }_2$ is considered. According to the assumptions above, a within-host HIV-1 infection fast time model takes the form

(2)$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}T\left(s\right)}{\mathrm{d}s}& ={\mathsf{\Lambda }}_c-{\kappa }_{1}T\left(s\right)V\left(s\right)-{\kappa }_{2}T\left(s\right)T^{*}\left(s\right)-{\mu }_{c}T\left(s\right),\hfill \\ \hfill \frac{\mathrm{d}{T}^{*}\left(s\right)}{\mathrm{d}s}& =({\kappa }_{1}T(s-{\tau }_1)V(s-{\tau }_1)+{\kappa }_{2}T(s-{\tau }_1)T^{*}(s-{\tau }_1))e^{-{\mu }_c{\tau }_1}-({\mu }_c+{\delta }_c)T^{*}\left(s\right),\hfill \\ \hfill \frac{\mathrm{d}V\left(s\right)}{\mathrm{d}s}& =\frac{c_1{\eta }_{1}I\left(t\right)+c_2{\eta }_{2}A\left(t\right)}{1+\alpha \left(I\right(t)+A(t\left)\right)}+p{T}^{*}\left(s\right)-qB\left(s\right)V\left(s\right)-{\mu }_{v}V\left(s\right),\hfill \\ \hfill \frac{\mathrm{d}B\left(s\right)}{\mathrm{d}s}& =\sigma B(s-{\tau }_2)V(s-{\tau }_2)-\zeta B\left(s\right).\hfill \end{array}\right.$

Similarly, all the parameters of Model (2) are nonnegative constants, and their biological interpretations are given in Table 2.

Remark 2.

The flowchart of Model (1) and Model (2) are shown by the Figure 1. In addition, the multiscale model is characterized by two distinct time scales, denoted by t and s, respectively. The former represents the transmission time at the population level, while the latter represents the evolution time at the individual host level. Typically, the virus–cell dynamics within the host occur on a faster time scale than those between hosts. Therefore, s is a faster time scale than t (see [29,30,31]). Without loss of generality, we suppose that $s=ϵt$$(0<ϵ\le 1)$, where ϵ is a small positive constant.

Since the second equation of Model (1) is decoupled from the other equations, in the analysis below, it suffices to study the following between-host model

(3)$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}=& {\mathsf{\Lambda }}_h-\frac{c_1{\beta }_1\left(V\left(s\right)\right)I\left(t\right)+c_2{\beta }_2\left(V\left(s\right)\right)A\left(t\right)}{1+\alpha \left(I\right(t)+A(t\left)\right)}S\left(t\right)-\mu S\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=& \frac{c_1{\beta }_1\left(V\left(s\right)\right)I(t-\tau )+c_2{\beta }_2\left(V\left(s\right)\right)A(t-\tau )}{1+\alpha \left(I\right(t-\tau )+A(t-\tau \left)\right)}S(t-\tau )e^{-\mu \tau }-(\mu +\xi +{\delta }_1\left(V\left(s\right)\right))I\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}=& \xi I\left(t\right)-(\mu +{\delta }_2\left(V\left(s\right)\right))A\left(t\right).\hfill \end{array}\right.$

3. Analysis of the Fast Time Model

In the fast time Model (2), assume that, for HIV-positive individuals, $I\left(t\right)$ is a constant I and, for AIDS individuals, $A\left(t\right)$ is a constant A, where $I=A=0$ means there is no possibility of viral transmission through direct contact with infected individuals, and I, $A>0$ represent that the virus in the process of contacting infected individuals. And we denote $\frac{c_1{\eta }_{1}I+c_2{\eta }_{2}A}{1+\alpha (I+A)}:=\mathsf{\Theta }(I,A)$. Hence, Model (2) becomes the following isolated within-host virus infection model

(4)$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}T\left(s\right)}{\mathrm{d}s}& ={\mathsf{\Lambda }}_c-{\kappa }_{1}T\left(s\right)V\left(s\right)-{\kappa }_{2}T\left(s\right)T^{*}\left(s\right)-{\mu }_{c}T\left(s\right),\hfill \\ \hfill \frac{\mathrm{d}{T}^{*}\left(s\right)}{\mathrm{d}s}& =({\kappa }_{1}T(s-{\tau }_1)V(s-{\tau }_1)+{\kappa }_{2}T(s-{\tau }_1)T^{*}(s-{\tau }_1))e^{-{\mu }_c{\tau }_1}-({\mu }_c+{\delta }_c)T^{*}\left(s\right),\hfill \\ \hfill \frac{\mathrm{d}V\left(s\right)}{\mathrm{d}s}& =\mathsf{\Theta }(I,A)+p{T}^{*}\left(s\right)-qB\left(s\right)V\left(s\right)-{\mu }_{v}V\left(s\right),\hfill \\ \hfill \frac{\mathrm{d}B\left(s\right)}{\mathrm{d}s}& =\sigma B(s-{\tau }_2)V(s-{\tau }_2)-\zeta B\left(s\right).\hfill \end{array}\right.$

Let ${\tau }_{max}=max\{{\tau }_1,{\tau }_2\}>0$, then we define $C([-{\tau }_{max},0],{\mathbb{R}}_{+}^4)$, equipped with the suitable sup-norm, as the Banach space of all continuous functions $\varphi =({\varphi }_1,{\varphi }_2,{\varphi }_3,{\varphi }_4):[-{\tau }_{max},0]\to {\mathbb{R}}_{+}^4$. From the biological background of Model (4), the initial condition for any solution is given by

(5)$\begin{array}{cc}\hfill & T\left(\vartheta \right)={\varphi }_1\left(\vartheta \right),T^{*}\left(\vartheta \right)={\varphi }_2\left(\vartheta \right),V\left(\vartheta \right)={\varphi }_3\left(\vartheta \right),B\left(\vartheta \right)={\varphi }_4\left(\vartheta \right),\hfill \\ \hfill & {\varphi }_i\left(\vartheta \right)\ge 0,\vartheta \in [-{\tau }_{max},0],{\varphi }_1\left(0\right)>0,i=1,2,3,4,\hfill \end{array}$

where $\varphi \in C([-{\tau }_{max},0],{\mathbb{R}}_{+}^4)$. It can be proved by the fundamental theory of functional differential equations [35] that Model (4) has the unique solution $(T\left(s\right),T^{*}\left(s\right),V\left(s\right),B\left(s\right))$ that satisfies the initial condition (5). Regarding the nonnegativity and boundedness of solutions for Model (4), we can state the following result.

Theorem 1.

Under the initial condition (5), the solution $(T\left(s\right),T^{*}\left(s\right),V\left(s\right),B\left(s\right))$ of Model (4) is nonnegative and bounded for all $s\ge 0$.

Following the approach in [21], the proof of nonnegativity and boundedness of solutions is straightforward and is omitted here.

Remark 3.

Theorem 1 indicates that the biologically feasible region

$\begin{array}{c}\hfill \mathsf{\Omega }=\left\{(T,T^{*},V,B)\in {\mathbb{R}}_{+}^4:T+e^{{\mu }_c{\tau }_1}{T}^{*}\le \frac{{\mathsf{\Lambda }}_c}{{\mu }_c},V+\frac{q}{\sigma }B\le \frac{p{\mathsf{\Lambda }}_c+e^{{\mu }_c{\tau }_1}{\mu }_c\mathsf{\Theta }(I,A)}{e^{{\mu }_c{\tau }_1}{\mu }_{c}min\{{\mu }_v,\zeta \}}\right\}\end{array}$

is positively invariant for Model (4). Therefore, in what follows, we only consider the solutions with initial conditions within the region Ω.

3.1. Feasible Equilibria and Basic Reproduction Number

In this section, we are dedicated to exploring the existence of equilibria and the reproduction number for Model (4). For $I=A=0$, the immune-inactivated reproduction number, ${\mathcal{R}}_0^w$, and the immune-activated reproduction number, ${\mathcal{R}}_1^w$, are defined, respectively, by

$\begin{array}{c}\hfill {\mathcal{R}}_0^w=\frac{(p{\kappa }_1+{\mu }_v{\kappa }_2){\mathsf{\Lambda }}_c}{{\mu }_c{\mu }_v({\mu }_c+{\delta }_c)}{e}^{-{\mu }_c{\tau }_1},{\mathcal{R}}_1^w=\frac{{\mathcal{R}}_0^w}{1+\frac{{\mu }_v({\mu }_c+{\delta }_c)\zeta }{\sigma p{\mathsf{\Lambda }}_c{e}^{-{\mu }_c{\tau }_1}}{\mathcal{R}}_0^w}.\end{array}$

The existence of equilibria of Model (4) in the case of $I=A=0$ has been investigated in [21], and here we can summarize the following results.

Lemma 1.

Assume that $I=A=0$ in Model (4), then the following statements are valid.

• $\left(i\right)$

  Model (4) always has infection-free equilibrium $U_0=(T_0,0,0,0)$, where $T_0={\mathsf{\Lambda }}_c/{\mu }_c$.

• $\left(ii\right)$

  If ${\mathcal{R}}_0^w>1$, Model (4) has an immunity-inactivated infection equilibrium $U_1=(T_1,T_1^{*},V_1,0)$, where

  $\begin{array}{c}\hfill T_1=\frac{T_0}{{\mathcal{R}}_0^w},T_1^{*}=\frac{{\mu }_v{\mu }_c}{{\kappa }_{1}p+{\kappa }_2{\mu }_v}({\mathcal{R}}_0^w-1),V_1=\frac{p{\mu }_c}{{\kappa }_{1}p+{\kappa }_2{\mu }_v}({\mathcal{R}}_0^w-1).\end{array}$

• $\left(iii\right)$

  If ${\mathcal{R}}_1^w>1$, except for $U_0$ and $U_1$, Model (4) also has an immunity-activated infection equilibrium $U_2=(T_2,T_2^{*},V_2,B_2)$, where

  $\begin{array}{c}\hfill T_2=\frac{{\mathsf{\Lambda }}_c}{{\kappa }_1{V}_2+{\kappa }_2{T}_2^{*}+{\mu }_c},T_2^{*}=\frac{q\zeta }{p\sigma }{B}_2+\frac{{\mu }_v\zeta }{p\sigma },V_2=\frac{\zeta }{\sigma },\end{array}$

  and $B_2$ is the positive real root of the following equation

  $\begin{array}{c}\hfill B_2^2+\left(\frac{2{\kappa }_2\zeta {\mu }_v+{\kappa }_1\zeta p+\sigma {\mu }_c{\mu }_v}{{\kappa }_2\zeta q}-\frac{\sigma p{\mathsf{\Lambda }}_c{e}^{-{\mu }_c{\tau }_1}}{({\mu }_c+{\delta }_c)\zeta q}\right)B_2+\left(\frac{{\kappa }_2\zeta {\mu }_v^2+{\kappa }_1\zeta p{\mu }_v+\sigma {\mu }_{c}p{\mu }_v}{{\kappa }_2\zeta q^2}\right)(1-{\mathcal{R}}_1^w)=0.\end{array}$

For I, $A>0$, Model (4) has no infection-free equilibrium, and the immune-activated reproduction number for Model (4) is defined by

$\begin{array}{c}\hfill {\mathcal{R}}_w=\frac{2\sigma \mathsf{\Theta }(I,A)+\sigma p\left(-b_1+\sqrt{b_1^2-4b_2}\right)}{2{\mu }_v\zeta },\end{array}$

where

$b_1=\frac{({\mu }_c+{\delta }_c)({\kappa }_1\zeta +{\mu }_c\sigma )e^{{\mu }_c{\tau }_1}-{\mathsf{\Lambda }}_c{\kappa }_2\sigma }{{\kappa }_2({\mu }_c+{\delta }_c)\sigma e^{{\mu }_c{\tau }_1}},b_2=-\frac{{\mathsf{\Lambda }}_c{\kappa }_1\zeta }{{\kappa }_2({\mu }_c+{\delta }_c)\sigma e^{{\mu }_c{\tau }_1}}.$

The existence of equilibria for Model (4) in the scenario where $I,A>0$ can be obtained based on the reproduction number ${\mathcal{R}}_w$, as shown below.

Theorem 2.

Assume that $I,A>0$ in Model (4), the following statements are valid.

• $\left(i\right)$

  If ${\mathcal{R}}_w\le 1$, Model (4) only has an immunity-inactivated infection equilibrium ${\tilde{U}}_1({\tilde{T}}_1,{\tilde{T}}_1^{*},{\tilde{V}}_1,0)$, where ${\tilde{T}}_1$, ${\tilde{T}}_1^{*}$, and ${\tilde{V}}_1$ are given below. Furthermore,

  $\underset{A,I\to 0^{+}}{lim}{\tilde{U}}_1({\tilde{T}}_1,{\tilde{T}}_1^{*},{\tilde{V}}_1,0)=\left\{\begin{array}{cc}{U}_0(T_0,0,0,0),\hfill & \hfill {\mathcal{R}}_0^w\le 1,\\ U_1(T_1,T_1^{*},V_1,0),\hfill & \hfill {\mathcal{R}}_0^w>1.\end{array}\right.$

• $\left(ii\right)$

  If ${\mathcal{R}}_w>1$, except for ${\tilde{U}}_1$, Model (4) also has a unique immunity-activated infection equilibrium, ${\tilde{U}}_2({\tilde{T}}_2,{\tilde{T}}_2^{*},{\tilde{V}}_2,{\tilde{B}}_2)$, where ${\tilde{T}}_2$, ${\tilde{T}}_2^{*}$, ${\tilde{V}}_2$, and ${\tilde{B}}_2$ are given below. Furthermore,

  $\underset{A,I\to 0^{+}}{lim}{\tilde{U}}_2({\tilde{T}}_2,{\tilde{T}}_2^{*},{\tilde{V}}_2,{\tilde{B}}_2)=\left\{\begin{array}{cc}{U}_0(T_0,0,0,0),\hfill & \hfill {\mathcal{R}}_0^w\le 1,\\ U_1(T_1,T_1^{*},V_1,0),\hfill & \hfill {\mathcal{R}}_0^w>1,{\mathcal{R}}_1^w\le 1,\\ U_2(T_2,T_2^{*},V_2,B_2),\hfill & \hfill {\mathcal{R}}_0^w>1,{\mathcal{R}}_1^w>1.\end{array}\right.$

Proof. 

Firstly, when I, $A>0$, the equilibrium ${\tilde{U}}_1=({\tilde{T}}_1,{\tilde{T}}_1^{*},{\tilde{V}}_1,0)$ satisfies the following equations

(6)$\left\{\begin{array}{c}{\mathsf{\Lambda }}_c-{\kappa }_1{\tilde{T}}_1{\tilde{V}}_1-{\kappa }_2{\tilde{T}}_1{\tilde{T}}_1^{*}-{\mu }_c{\tilde{T}}_1=0,\hfill \\ ({\kappa }_1{\tilde{T}}_1{\tilde{V}}_1+{\kappa }_2{\tilde{T}}_1{\tilde{T}}_1^{*})e^{-{\mu }_c{\tau }_1}-({\mu }_c+{\delta }_c){\tilde{T}}_1^{*}=0,\hfill \\ \mathsf{\Theta }(I,A)+p{\tilde{T}}_1^{*}-{\mu }_v{\tilde{V}}_1=0.\hfill \end{array}\right.$

We directly obtain

(7)$\begin{array}{c}\hfill {\tilde{T}}_1^{*}=\frac{({\mathsf{\Lambda }}_c-{\mu }_c{\tilde{T}}_1)e^{-{\mu }_c{\tau }_1}}{{\mu }_c+{\delta }_c},{\tilde{V}}_1=\frac{\mathsf{\Theta }(I,A)+p{\tilde{T}}_1^{*}}{{\mu }_v},\end{array}$

and ${\tilde{T}}_1$ is a solution of the following quadratic equation

(8)$\begin{array}{c}\hfill {\tilde{T}}_1^2+a_1{\tilde{T}}_1+a_2=0,\end{array}$

where $a_1=-T_0\left(1+\frac{1}{{\mathcal{R}}_0^w}+\frac{{\kappa }_1\mathsf{\Theta }(I,A)}{{\mu }_v{\mu }_c{\mathcal{R}}_0^w}\right)<0$, $a_2=\frac{T_0^2}{{\mathcal{R}}_0^w}>0$. By calculation, Equation (8) has two real roots

$\begin{array}{c}\hfill {\tilde{T}}_{1+}=\frac{1}{2}\left(-a_1+\sqrt{a_1^2-4a_2}\right)>0,{\tilde{T}}_{1-}=\frac{1}{2}\left(-a_1-\sqrt{a_1^2-4a_2}\right)>0,\end{array}$

with $a_1^2-4a_2>0$. Further, it is easy to verify ${\tilde{T}}_{1+}>\frac{{\mathsf{\Lambda }}_c}{{\mu }_c}>{\tilde{T}}_{1-}$. Hence, Model (4) has a unique immunity-inactivated infection equilibrium, ${\tilde{U}}_1=({\tilde{T}}_1,{\tilde{T}}_1^{*},{\tilde{V}}_1,0)$, with

$\begin{array}{c}\hfill {\tilde{T}}_1=\frac{1}{2}\left(-a_1-\sqrt{a_1^2-4a_2}\right),{\tilde{V}}_1=\frac{\mathsf{\Theta }(I,A)({\mu }_c+{\delta }_c)+p({\mathsf{\Lambda }}_c-{\mu }_c{\tilde{T}}_1)e^{-{\mu }_c{\tau }_1}}{{\mu }_v({\mu }_c+{\delta }_c)},{\tilde{T}}_1^{*}=\frac{({\mathsf{\Lambda }}_c-{\mu }_c{\tilde{T}}_1)e^{-{\mu }_c{\tau }_1}}{{\mu }_c+{\delta }_c}.\end{array}$

Secondly, when I, $A>0$, we obtain

$\begin{array}{c}\hfill {\tilde{T}}_2=\frac{{\mathsf{\Lambda }}_c-({\mu }_c+{\delta }_c){\tilde{T}}_2^{*}{e}^{{\mu }_c{\tau }_1}}{{\mu }_c},{\tilde{V}}_2=\frac{\zeta }{\sigma },\end{array}$

and ${\tilde{T}}_2^{*}$ satisfies ${\tilde{T}}_1^{*2}+b_1{\tilde{T}}_1^{*}+b_2=0$. Solving this equation, we obtain two real roots

$\begin{array}{c}\hfill {\tilde{T}}_{2+}^{*}=\frac{1}{2}\left(-b_1+\sqrt{b_1^2-4b_2}\right)>0,{\tilde{T}}_{2-}^{*}=\frac{1}{2}\left(-b_1-\sqrt{b_1^2-4b_2}\right)<0,\end{array}$

with $b_1^2-4b_2>0$. For the positive equilibrium, ${\tilde{U}}_2$, it should satisfy ${\tilde{T}}_2$, ${\tilde{T}}_2^{*}$, ${\tilde{V}}_2$, ${\tilde{B}}_2>0$. Thus, if ${\mathcal{R}}_w>1$, Model (4) has a unique immunity-activated infection equilibrium, ${\tilde{U}}_2=({\tilde{T}}_2,{\tilde{T}}_2^{*},{\tilde{V}}_2,{\tilde{B}}_2)$, with

$\begin{array}{cc}\hfill {\tilde{V}}_2& =\frac{\zeta }{\sigma },{\tilde{T}}_2^{*}=\frac{1}{2}\left(-b_1+\sqrt{b_1^2-4b_2}\right),{\tilde{T}}_2=\frac{{\mathsf{\Lambda }}_c}{{\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c},\hfill \\ \hfill {\tilde{B}}_2& =\frac{2\sigma \mathsf{\Theta }(I,A)+\sigma p\left(-b_1+\sqrt{b_1^2-4b_2}\right)-2{\mu }_v\zeta }{2q\zeta }.\hfill \end{array}$

The proof is finished. □

3.2. Stability of Equilibria and Hopf Bifurcation

In this section, we with the global asymptotic stability of feasible equilibria and the existence of Hopf bifurcation of Model (4). To facilitate our discussion, we introduce the function $\rho \left(x\right)=x-1-lnx$. Note that $\rho \left(x\right)\ge 0$ for $x>0$ and $\frac{\mathrm{d}\rho \left(x\right)}{\mathrm{d}x}=1-\frac{1}{x}$. Moreover, $\rho \left(x\right)=0$ if and only if $x=1$. On the global asymptotic stability of ${\tilde{U}}_1$ and ${\tilde{U}}_2$, the following result is established.

Theorem 3.

Let I, $A>0$ in Model (4),

• $\left(i\right)$

  if ${\mathcal{R}}_w\le 1$, the immunity-inactivated equilibrium, ${\tilde{U}}_1$, is globally asymptotically stable;

• $\left(ii\right)$

  if ${\mathcal{R}}_w>1$ and ${\tau }_2=0$, the immunity-activated infection equilibrium, ${\tilde{U}}_2$, is globally asymptotically stable.

Proof. 

We focus, firstly, the stability of equilibrium ${\tilde{U}}_1$. To this end, we choose

$\begin{array}{cc}\hfill L_1\left(s\right)=& {\tilde{T}}_1\rho \left(\frac{T\left(s\right)}{{\tilde{T}}_1}\right)+e^{{\mu }_c{\tau }_1}{\tilde{T}}_1^{*}\rho \left(\frac{T^{*}\left(s\right)}{{\tilde{T}}_1^{*}}\right)+\frac{{\kappa }_1{\tilde{T}}_1{\tilde{V}}_1}{p{\tilde{T}}_1^{*}}{\tilde{V}}_1\rho \left(\frac{V\left(s\right)}{{\tilde{V}}_1}\right)+{\kappa }_1{\tilde{T}}_1{\tilde{V}}_1{\int }_{s-{\tau }_1}^s\rho \left(\frac{T\left(x_1\right)V\left(x_1\right)}{{\tilde{T}}_1{\tilde{V}}_1}\right)\mathrm{d}{x}_1\hfill \\ \hfill & +\frac{q{\kappa }_1{\tilde{T}}_1{\tilde{V}}_1}{\sigma p{\tilde{T}}_1^{*}}B\left(s\right)+{\kappa }_2{\tilde{T}}_1{\tilde{T}}_1^{*}{\int }_{s-{\tau }_1}^s\rho \left(\frac{T\left(x_1\right)T^{*}\left(x_1\right)}{{\tilde{T}}_1{\tilde{T}}_1^{*}}\right)\mathrm{d}{x}_1+\frac{q{\kappa }_1{\tilde{T}}_1{\tilde{V}}_1}{p{\tilde{T}}_1^{*}}{\int }_{s-{\tau }_2}^{s}B\left(x_1\right)V\left(x_1\right)\mathrm{d}{x}_1.\hfill \end{array}$

Calculating the time derivative of $L_1\left(s\right)$ along the positive solution of Model (4) yields

$\begin{array}{cc}\hfill \frac{\mathrm{d}{L}_1\left(s\right)}{\mathrm{d}s}=& -{\kappa }_1{\tilde{T}}_1{\tilde{V}}_1\rho \left(\frac{{\tilde{T}}_1}{T\left(s\right)}\right)-{\kappa }_1{\tilde{T}}_1{\tilde{V}}_1\rho \left(\frac{T^{*}\left(s\right){\tilde{V}}_1}{{\tilde{T}}_1^{*}V\left(s\right)}\right)-{\kappa }_1{\tilde{T}}_1{\tilde{V}}_1\rho \left(\frac{T(s-{\tau }_1)V(s-{\tau }_1){\tilde{T}}_1^{*}}{{\tilde{T}}_1{\tilde{V}}_1{T}^{*}\left(s\right)}\right)\hfill \\ \hfill & -{\kappa }_2{\tilde{T}}_1{\tilde{T}}_1^{*}\rho \left(\frac{{\tilde{T}}_1}{T\left(s\right)}\right)-{\kappa }_2{\tilde{T}}_1{\tilde{T}}_1^{*}\rho \left(\frac{T(s-{\tau }_1)T^{*}(s-{\tau }_1)}{{\tilde{T}}_1{T}^{*}\left(s\right)}\right)+\frac{\zeta q{\kappa }_1{\tilde{T}}_1{\tilde{V}}_1}{\sigma p{\tilde{T}}_1^{*}}\left({\mathcal{R}}_w-1\right)B\left(s\right)\hfill \\ \hfill & -{\mu }_c\frac{{(T\left(s\right)-{\tilde{T}}_1)}^2}{T\left(s\right)}+\frac{{\kappa }_1{\tilde{T}}_1{\tilde{V}}_1}{p{\tilde{T}}_1^{*}}\mathsf{\Theta }(I,A)\left(2-\frac{V\left(s\right)}{{\tilde{V}}_1}-\frac{{\tilde{V}}_1}{V\left(s\right)}\right).\hfill \end{array}$

Therefore, under the assumption ${\mathcal{R}}_w\le 1$, we have $\frac{\mathrm{d}{L}_1\left(s\right)}{\mathrm{d}s}\le 0$. Additionally, it can be easily verified that $T\left(s\right)={\tilde{T}}_1$, $T^{*}\left(s\right)={\tilde{T}}_1^{*}$, $V\left(s\right)={\tilde{V}}_1$, and $B\left(s\right)=0$ when $\frac{\mathrm{d}{L}_1\left(s\right)}{\mathrm{d}s}=0$. Consequently, the singleton $\left\{{\tilde{U}}_1\right\}$ is the largest compact invariant set in $\{(T\left(s\right),T^{*}\left(s\right),B\left(s\right),V\left(s\right))\in {\mathbb{R}}_{+}^4:\frac{\mathrm{d}{L}_1\left(s\right)}{\mathrm{d}s}=0\}$. Therefore, by applying LaSalle’s invariance principle, the immunity-inactivated infection equilibrium, ${\tilde{U}}_1$, is globally asymptotically stable when ${\mathcal{R}}_w\le 1$.

Now, we turn to the conclusion $\left(ii\right)$. For equilibrium ${\tilde{U}}_2$, one defines

$\begin{array}{cc}\hfill L_2\left(s\right)=& {\tilde{T}}_2\rho \left(\frac{T\left(s\right)}{{\tilde{T}}_2}\right)+e^{{\mu }_c{\tau }_1}{\tilde{T}}_2^{*}\rho \left(\frac{T^{*}\left(s\right)}{{\tilde{T}}_2^{*}}\right)+\frac{{\kappa }_1{\tilde{T}}_2{\tilde{V}}_2}{p{\tilde{T}}_2^{*}}{\tilde{V}}_2\rho \left(\frac{V\left(s\right)}{{\tilde{V}}_2}\right)+\frac{q{\kappa }_1{\tilde{T}}_2{\tilde{V}}_2}{\sigma P{\tilde{T}}_2^{*}}{\tilde{B}}_2\rho \left(\frac{B\left(s\right)}{{\tilde{B}}_2}\right)\hfill \\ \hfill & +{\kappa }_1{\tilde{T}}_2{\tilde{V}}_2{\int }_{s-{\tau }_1}^s\rho \left(\frac{T\left(x_1\right)V\left(x_1\right)}{{\tilde{T}}_2{\tilde{V}}_2}\right)\mathrm{d}{x}_1+{\kappa }_2{\tilde{T}}_2{\tilde{T}}_2^{*}{\int }_{s-{\tau }_1}^s\rho \left(\frac{T\left(x_1\right)T^{*}\left(x_1\right)}{{\tilde{T}}_2{\tilde{T}}_2^{*}}\right)\mathrm{d}{x}_1.\hfill \end{array}$

Calculating the time derivative of $L_2\left(s\right)$ along the positive solution of Model (4) yields

$\begin{array}{cc}\hfill \frac{\mathrm{d}{L}_2\left(s\right)}{\mathrm{d}s}=& -{\mu }_c\frac{{(T\left(s\right)-{\tilde{T}}_2)}^2}{T\left(s\right)}-{\kappa }_1{\tilde{T}}_2{\tilde{V}}_2\rho \left(\frac{{\tilde{T}}_2}{T\left(s\right)}\right)-{\kappa }_1{\tilde{T}}_2{\tilde{V}}_2\rho \left(\frac{T^{*}\left(s\right){\tilde{V}}_2}{{\tilde{T}}_2^{*}V\left(s\right)}\right)-{\kappa }_1{\tilde{T}}_2{\tilde{V}}_2\rho \left(\frac{T(s-{\tau }_1)V(s-{\tau }_1){\tilde{T}}_2^{*}}{{\tilde{T}}_2{\tilde{V}}_2{T}^{*}\left(s\right)}\right)\hfill \\ \hfill & -{\kappa }_2{\tilde{T}}_2{\tilde{T}}_2^{*}\rho \left(\frac{{\tilde{T}}_2}{T\left(s\right)}\right)-{\kappa }_2{\tilde{T}}_2{\tilde{T}}_2^{*}\rho \left(\frac{T(s-{\tau }_1)T^{*}(s-{\tau }_1)}{{\tilde{T}}_2{T}^{*}\left(s\right)}\right)+\frac{{\kappa }_1{\tilde{T}}_2{\tilde{V}}_2}{p{\tilde{T}}_2^{*}}\mathsf{\Theta }(I,A)\left(2-\frac{V\left(s\right)}{{\tilde{V}}_2}-\frac{{\tilde{V}}_2}{V\left(s\right)}\right).\hfill \end{array}$

Thus, we have $\frac{\mathrm{d}{L}_2\left(s\right)}{\mathrm{d}s}\le 0$, and equality holds if and only if $T\left(s\right)={\tilde{T}}_2$, $T^{*}\left(s\right)={\tilde{T}}_2^{*}$, $V\left(s\right)={\tilde{V}}_2$, and $B\left(s\right)={\tilde{B}}_2$. Consequently, the singleton $\left\{{\tilde{U}}_2\right\}$ is the largest compact invariant set in $\{(T\left(s\right),T^{*}\left(s\right),B\left(s\right),V\left(s\right))\in {\mathbb{R}}_{+}^4:\frac{\mathrm{d}{L}_2\left(s\right)}{\mathrm{d}s}=0\}$. Therefore, by applying LaSalle’s invariance principle, the immunity-activated equilibrium, ${\tilde{U}}_2$, is globally asymptotically stable when ${\mathcal{R}}_w>1$ and ${\tau }_2=0$. □

Remark 4.

For Model (4) with $I=A=0$, it is easy to verify the global asymptotic stability of $U_0$, $U_1$, and $U_2$, respectively, by constructing

$\begin{array}{cc}\hfill L_0\left(s\right)=& T_0\rho \left(\frac{T\left(s\right)}{T_0}\right)+e^{{\mu }_c{\tau }_1}{T}^{*}\left(s\right)+\frac{{\kappa }_1{T}_0}{{\mu }_v}V\left(s\right)+\frac{q{\kappa }_1{T}_0}{\sigma {\mu }_v}B\left(s\right)+\frac{q{\kappa }_1{T}_0}{{\mu }_v}{\int }_{s-{\tau }_2}^{s}B\left(x_1\right)V\left(x_1\right)\mathrm{d}{x}_1\hfill \\ \hfill & +{\int }_{s-{\tau }_1}^s\left({\kappa }_{1}T\left(x_1\right)V\left(x_1\right)+{\kappa }_{2}T\left(x_1\right)T^{*}\left(x_1\right)\right)\mathrm{d}{x}_1\hfill \end{array}$

and same Lyapunov functionals as in Theorem 3. Thus we omit here.

Next, our focus is on examining the necessary and sufficient conditions for Hopf bifurcation to take place at ${\tilde{U}}_2$. The characteristic equation of (4) at ${\tilde{U}}_2$ is

(9)$\begin{array}{cc}\hfill & {\lambda }^4+{\overline{C}}_3{\lambda }^3+{\overline{C}}_2{\lambda }^2+{\overline{C}}_1\lambda +{\overline{C}}_0+({\overline{D}}_3{\lambda }^3+{\overline{D}}_2{\lambda }^2+{\overline{D}}_1\lambda +{\overline{D}}_0)e^{-({\mu }_c+\lambda ){\tau }_1}\hfill \\ \hfill +& ({\overline{E}}_3{\lambda }^3+{\overline{E}}_2{\lambda }^2+{\overline{E}}_1\lambda +{\overline{E}}_0)e^{-\lambda {\tau }_2}+({\overline{F}}_2{\lambda }^2+{\overline{F}}_1\lambda +{\overline{F}}_0)e^{-{\mu }_c{\tau }_1-({\tau }_1+{\tau }_2)\lambda }=0,\hfill \end{array}$

where

$\begin{array}{cc}\hfill {\overline{C}}_3& ={\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+2{\mu }_c+{\delta }_c+q{\tilde{B}}_2+{\mu }_v+\zeta ,{\overline{E}}_1=-\sigma {\tilde{V}}_2({\mu }_c+{\delta }_c)({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c)-\sigma {\tilde{V}}_2{\mu }_v({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+2{\mu }_c+{\delta }_c),\hfill \\ \hfill {\overline{C}}_2& =({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c)({\mu }_c+{\delta }_c)+({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+2{\mu }_c+{\delta }_c)(q{\tilde{B}}_2+{\mu }_v+\zeta )+\zeta (q{\tilde{B}}_2+{\mu }_v),\hfill \\ \hfill {\overline{C}}_1& =(q{\tilde{B}}_2+{\mu }_v+\zeta )({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c)({\mu }_c+{\delta }_c)+\zeta (q{\tilde{B}}_2+{\mu }_v)({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+2{\mu }_c+{\delta }_c),\hfill \\ \hfill {\overline{C}}_0& =\zeta (q{\tilde{B}}_2+{\mu }_v)({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c)({\mu }_c+{\delta }_c),{\overline{D}}_3=-{\kappa }_2{\tilde{T}}_2,{\overline{D}}_2=-{\kappa }_2{\tilde{T}}_2({\mu }_c+q{\tilde{B}}_2+{\mu }_v+\zeta )-p{\kappa }_1{\tilde{T}}_2,\hfill \\ \hfill {\overline{D}}_1& =-{\kappa }_2{\tilde{T}}_2{\mu }_c(q{\tilde{B}}_2+{\mu }_v+\zeta )-p{\kappa }_1{\tilde{T}}_2({\mu }_c+\zeta )-\zeta {\kappa }_2{\tilde{T}}_2(q{\tilde{B}}_2+{\mu }_v),{\overline{D}}_0=-\zeta {\kappa }_2{\tilde{T}}_2{\mu }_c(q{\tilde{B}}_2+{\mu }_v)-\zeta p{\kappa }_1{\tilde{T}}_2{\mu }_c,\hfill \\ \hfill {\overline{E}}_3& =-\sigma {\tilde{V}}_2,{\overline{E}}_2=-\sigma V_2({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+2{\mu }_c+{\delta }_c+{\mu }_v),{\overline{E}}_0=-\sigma {\tilde{V}}_2{\mu }_c({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c)({\mu }_c+{\delta }_c),\hfill \\ \hfill {\overline{F}}_2& =\sigma {\tilde{V}}_2{\kappa }_2{\tilde{T}}_2^{*},{\overline{F}}_1=\sigma {\tilde{V}}_2{\kappa }_2{\tilde{T}}_2({\mu }_c+{\mu }_v)+\sigma {\tilde{V}}_{2}p{\kappa }_1{\tilde{T}}_2,{\overline{F}}_0=\sigma {\tilde{V}}_2{\tilde{T}}_2{\mu }_c({\kappa }_2{\mu }_v+{\kappa }_{1}p).\hfill \end{array}$

For ${\tau }_2>0$, we will categorize the stability analysis of ${\tilde{U}}_2$ into two cases, respectively.

Case 1:

${\tau }_1=0$ and ${\tau }_2>0$.

Equation (9) is, in this case, reduced to

(10)$\begin{array}{cc}\hfill & {\lambda }^4+{\overline{P}}_3{\lambda }^3+{\overline{P}}_2{\lambda }^2+{\overline{P}}_1\lambda +{\overline{P}}_0+({\overline{Q}}_3{\lambda }^3+{\overline{Q}}_2{\lambda }^2+{\overline{Q}}_1\lambda +{\overline{Q}}_0)e^{-\lambda {\tau }_2}=0,\hfill \end{array}$

where ${\overline{P}}_3={\overline{C}}_3+{\overline{D}}_3$, ${\overline{P}}_2={\overline{C}}_2+{\overline{D}}_2$, ${\overline{P}}_1={\overline{C}}_1+{\overline{D}}_1$, ${\overline{P}}_0={\overline{C}}_0+{\overline{D}}_0$, ${\overline{Q}}_3={\overline{E}}_3$, ${\overline{Q}}_2={\overline{E}}_2+{\overline{F}}_2$, ${\overline{Q}}_1={\overline{E}}_1+{\overline{F}}_1$, ${\overline{Q}}_0={\overline{E}}_0+{\overline{F}}_0$. Let $\lambda =i\omega$ be a purely imaginary root of (10), and substituting $\lambda =i\omega$ into (10) and separating the real and imaginary parts, we obtain

(11)$\left\{\begin{array}{c}({\overline{Q}}_3{\omega }^3-{\overline{Q}}_1\omega )sin\left(\omega {\tau }_2\right)+({\overline{Q}}_2{\omega }^2-{\overline{Q}}_0)cos\left(\omega {\tau }_2\right)={\omega }^4-{\overline{P}}_2{\omega }^2+{\overline{P}}_0,\hfill \\ ({\overline{Q}}_3{\omega }^3-{\overline{Q}}_1\omega )cos\left(\omega {\tau }_2\right)-({\overline{Q}}_2{\omega }^2-{\overline{Q}}_0)sin\left(\omega {\tau }_2\right)=-{\overline{P}}_3{\omega }^3+{\overline{P}}_1\omega .\hfill \end{array}\right.$

Squaring and adding both above equations lead to the following equation

(12)$\begin{array}{c}\hfill {\omega }^8+{\overline{G}}_6{\omega }^6+{\overline{G}}_4{\omega }^4+{\overline{G}}_2{\omega }^2+{\overline{G}}_0=0,\end{array}$

where ${\overline{G}}_6={\overline{P}}_3^2-{\overline{Q}}_3^2-2{\overline{P}}_2$, ${\overline{G}}_4={\overline{P}}_2^2+2{\overline{P}}_0+2{\overline{Q}}_1{\overline{Q}}_3-2{\overline{P}}_1{\overline{P}}_3-{\overline{Q}}_2^2$, ${\overline{G}}_2={\overline{P}}_1^2+2{\overline{Q}}_0{\overline{Q}}_2-2{\overline{P}}_0{\overline{P}}_2-{\overline{Q}}_1^2$, ${\overline{G}}_0={\overline{P}}_0^2-{\overline{Q}}_0^2$. Let $u={\omega }^2$, then (12) is equivalent to

(13)$\begin{array}{c}\hfill f\left(u\right)=u^4+{\overline{G}}_6{u}^3+{\overline{G}}_4{u}^2+{\overline{G}}_{2}u+{\overline{G}}_0=0.\end{array}$

If ${\overline{G}}_0<0$, then it can be observed that (13) must have at least one positive root. This is because the function $f\left(u\right)$ satisfies $f\left(0\right)={\overline{G}}_0<0$ and ${lim}_{u\to \infty }f\left(u\right)=+\infty$. Therefore, there exists a value $u_0\in (0,+\infty )$ such that $f\left(u_0\right)=0$. Furthermore,

(14)$\begin{array}{c}\hfill f^{\prime }\left(u\right)=4u^3+3{\overline{G}}_6{u}^2+2{\overline{G}}_{4}u+{\overline{G}}_2.\end{array}$

Let $v=u+\frac{3{\overline{G}}_6}{4}$, then Equation (14) becomes $v^3+n_{1}v+n_2=0$, where $n_1=\frac{1}{2}{\overline{G}}_4^2-\frac{3}{16}{\overline{G}}_6^2$, $n_2=\frac{1}{32}{\overline{G}}_6^3-\frac{1}{8}{\overline{G}}_6{\overline{G}}_4+{\overline{G}}_2$. Let $\Delta ={\left(\frac{n_2}{2}\right)}^2+{\left(\frac{n_1}{3}\right)}^3$ and $j=\frac{-1+i\sqrt{3}}{2}$, then

$\begin{array}{cccc}\hfill & v_1=\sqrt[3]{-\frac{n_2}{2}+\sqrt{\Delta }}+\sqrt[3]{-\frac{n_2}{2}-\sqrt{\Delta }},\hfill & \hfill v_2& =j\sqrt[3]{-\frac{n_2}{2}+\sqrt{\Delta }}+j^2\sqrt[3]{-\frac{n_2}{2}-\sqrt{\Delta }},\hfill \\ \hfill & v_3=j^2\sqrt[3]{-\frac{n_2}{2}+\sqrt{\Delta }}+j\sqrt[3]{-\frac{n_2}{2}-\sqrt{\Delta }},\hfill & \hfill u_i& =v_i-\frac{3{\overline{G}}_6}{4},i=1,2,3.\hfill \end{array}$

Applying the method described in [36], we can obtain the following result regarding the root distribution of Equation (13).

Lemma 2.

For polynomial Equation (13),

• $\left(i\right)$

  if ${\overline{G}}_0<0$, then Equation (13) has at least one positive root;

• $\left(ii\right)$

  if ${\overline{G}}_0\ge 0$ and $\Delta \ge 0$, then Equation (13) has a positive root if and only if $u_1>0$ and $f\left(u_1\right)<0$;

• $\left(iii\right)$

  if ${\overline{G}}_0\ge 0$ and $\Delta <0$, then Equation (13) has a positive root if and only if there exists at least one $u^{*}\in \{u_1,u_2,u_3\}$, such that $u^{*}>0$ and $f\left(u^{*}\right)\le 0$.

If the conditions in Lemma 2 are invalid, then Equation (13) does not have any positive root. This implies that Equation (10) does not have purely imaginary roots. Furthermore, ${\tilde{U}}_2$ is locally asymptotically stable for all delay ${\tau }_2\ge 0$. As a result, the existence of Hopf bifurcation is not possible.

We now assume that one of the conditions in Lemma 2, which guarantees the existence of positive roots for Equation (13), is true and that Equation (13) has $k_0$ ($1\le k_0\le 4$) positive roots. Let $u_k$, $k=1,\dots ,k_0$ be the positive roots of (13), then Equation (12) has positive roots ${\omega }_k=\sqrt{u_k}$. For a given ${\omega }_k$, the corresponding critical value of delay ${\tau }_2$ is

$\begin{array}{c}\hfill {\tau }_{2,k}^{\left(j\right)}=\frac{1}{{\omega }_k}arccos\frac{({\omega }_k^4-{\overline{P}}_2{\omega }_k^2+{\overline{P}}_0)({\overline{Q}}_2{\omega }_k^2-{\overline{Q}}_0)+(-{\overline{P}}_3{\omega }_k^3+{\overline{P}}_1{\omega }_k)({\overline{Q}}_3{\omega }_k^3-{\overline{Q}}_1{\omega }_k)}{{({\overline{Q}}_3{\omega }_k^3-{\overline{Q}}_1{\omega }_k)}^2+{({\overline{Q}}_2{\omega }_k^2-{\overline{Q}}_0)}^2}+\frac{2j\pi }{{\omega }_k},\end{array}$

$j=0,1,2,\dots$. Then $\pm i{\omega }_k$ are a pair of purely imaginary roots of (10) with ${\tau }_2={\tau }_{2,k}^{\left(j\right)}$. Define

$\begin{array}{c}\hfill {\tau }_2^0={\tau }_{2,k_0}^{\left(0\right)}=min\left\{{\tau }_{2,k}^{\left(j\right)}|k=1,2,\dots ,k_0,j=0,1,\dots \right\},{\omega }_0={\omega }_{k_0},u_0=u_{k_0}.\end{array}$

Let $\lambda \left(\tau \right)=\alpha \left(\tau \right)+i\omega \left(\tau \right)$ be the root of Equation (10) at $\tau ={\tau }_2^0$ satisfying $\alpha \left({\tau }_2^0\right)=0$ and $\omega \left({\tau }_2^0\right)={\omega }_0$. Similar to the proof in [14], the transversality condition at ${\tau }_2={\tau }_2^0$ will be verified. Differentiating both sides of (10) with respect to ${\tau }_2$, we obtain

$\begin{array}{c}\hfill {\left(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2}\right)}^{-1}=-\frac{4{\lambda }^3+3{\overline{P}}_3{\lambda }^2+2{\overline{P}}_2\lambda +{\overline{P}}_1}{\lambda ({\lambda }^4+{\overline{P}}_3{\lambda }^3+{\overline{P}}_2{\lambda }^2+{\overline{P}}_1\lambda +{\overline{P}}_0)}+\frac{3{\overline{Q}}_3{\lambda }^2+2{\overline{Q}}_2\lambda +{\overline{Q}}_1}{\lambda ({\overline{Q}}_3{\lambda }^3+{\overline{Q}}_2{\lambda }^2+{\overline{Q}}_1\lambda +{\overline{Q}}_0)}-\frac{{\tau }_2}{\lambda }.\end{array}$

From Equation (11), one has

$\begin{array}{c}\hfill {({\overline{Q}}_3{\omega }_0^3-{\overline{Q}}_1{\omega }_0)}^2+{({\overline{Q}}_2{\omega }_0^2-{\overline{Q}}_0)}^2={({\omega }_0^4-{\overline{P}}_2{\omega }_0^2+{\overline{P}}_0)}^2+{({\overline{P}}_1{\omega }_0-{\overline{P}}_3{\omega }_0^3)}^2.\end{array}$

Then, evaluating ${\left(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2}\right)}^{-1}$ at ${\tau }_2={\tau }_2^0$$(\mathrm{i}.\mathrm{e}.$$\lambda =i{\omega }_0)$ and taking the real part, we obtain

$\begin{array}{c}\hfill \mathrm{Re}\left[{\left(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2}\right)}^{-1}{|}_{{\tau }_2={\tau }_2^0}\right]=\frac{4{\omega }_0^6+3{\overline{G}}_6{\omega }_0^4+2{\overline{G}}_4{\omega }_0^2+{\overline{G}}_2}{{({\overline{Q}}_3{\omega }_0^3-{\overline{Q}}_1{\omega }_0)}^2+{({\overline{Q}}_2{\omega }_0^2-{\overline{Q}}_0)}^2}=\frac{f^{\prime }\left({\omega }_0^2\right)}{{({\overline{Q}}_3{\omega }_0^3-{\overline{Q}}_1{\omega }_0)}^2+{({\overline{Q}}_2{\omega }_0^2-{\overline{Q}}_0)}^2}.\end{array}$

Therefore,

$\begin{array}{c}\hfill \mathrm{sign}\left\{\frac{\mathrm{d}\left(\mathrm{Re}\right(\lambda \left)\right)}{\mathrm{d}{\tau }_2}{|}_{{\tau }_2={\tau }_2^0}\right\}=\mathrm{sign}\left\{\mathrm{Re}\left[{\left(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2}\right)}^{-1}{|}_{{\tau }_2={\tau }_2^0}\right]\right\}=\mathrm{sign}\left\{f^{\prime }\left({\omega }_0^2\right)\right\}=\mathrm{sign}\left\{f^{\prime }\left(u_0\right)\right\}.\end{array}$

Based on the aforementioned analysis and the Hopf bifurcation theorem presented in [37], we can conclude the following result.

Theorem 4.

If all the conditions in Lemma 2 are not satisfied, then ${\tilde{U}}_2$ is locally asymptotically stable for all ${\tau }_2\ge 0$. If one of the conditions in Lemma 2 is satisfied, then ${\tilde{U}}_2$ is locally asymptotically stable for ${\tau }_2$ within the interval $[0,{\tau }_2^0)$. If one of the conditions in Lemma 2 is satisfied and $f^{\prime }\left(u_0\right)\ne 0$, then the transversality condition holds, and Model (4) will undergo a Hopf bifurcation at ${\tilde{U}}_2$ when ${\tau }_2={\tau }_2^0$.

Case 2:

${\tau }_1>0$ and ${\tau }_2>0$.

We take ${\tau }_2$ as the bifurcation parameter and assume the purely imaginary root of Equation (9) to be $\lambda =i\tilde{\omega }$ with $\tilde{\omega }>0$. By following the same computational process, it is clear that

(15)$\left\{\begin{array}{c}{\mathsf{\Pi }}_{1}cos\left(\tilde{\omega }{\tau }_2\right)-{\mathsf{\Pi }}_{2}sin\left(\tilde{\omega }{\tau }_2\right)={\mathsf{\Pi }}_3,\hfill \\ {\mathsf{\Pi }}_{1}sin\left(\tilde{\omega }{\tau }_2\right)+{\mathsf{\Pi }}_{2}cos\left(\tilde{\omega }{\tau }_2\right)={\mathsf{\Pi }}_4,\hfill \end{array}\right.$

where

$\begin{array}{cc}\hfill & {\mathsf{\Pi }}_1={\overline{E}}_2{\tilde{\omega }}^2-{\overline{E}}_0+({\overline{F}}_2{\tilde{\omega }}^2-{\overline{F}}_0)cos\left(\tilde{\omega }{\tau }_1\right)e^{-{\mu }_c{\tau }_1}-{\overline{F}}_1\tilde{\omega }sin\left(\tilde{\omega }{\tau }_1\right)e^{-{\mu }_c{\tau }_1},\hfill \\ \hfill & {\mathsf{\Pi }}_2={\overline{E}}_1\tilde{\omega }-{\overline{E}}_3{\tilde{\omega }}^3+({\overline{F}}_2{\tilde{\omega }}^2-{\overline{F}}_0)sin\left(\tilde{\omega }{\tau }_1\right)e^{-{\mu }_c{\tau }_1}+{\overline{F}}_1\tilde{\omega }cos\left(\tilde{\omega }{\tau }_1\right)e^{-{\mu }_c{\tau }_1},\hfill \\ \hfill & {\mathsf{\Pi }}_3={\tilde{\omega }}^4-{\overline{C}}_2{\tilde{\omega }}^2+{\overline{C}}_0+({\overline{D}}_0-{\overline{D}}_2{\tilde{\omega }}^2)cos\left(\tilde{\omega }{\tau }_1\right)e^{-{\mu }_c{\tau }_1}+({\overline{D}}_1\tilde{\omega }-{\overline{D}}_3{\tilde{\omega }}^3)sin\left(\tilde{\omega }{\tau }_1\right)e^{-{\mu }_c{\tau }_1},\hfill \\ \hfill & {\mathsf{\Pi }}_4={\overline{C}}_3{\tilde{\omega }}^3-{\overline{C}}_1\tilde{\omega }+({\overline{D}}_3{\tilde{\omega }}^3-{\overline{D}}_1\tilde{\omega })cos\left(\tilde{\omega }{\tau }_1\right)e^{-{\mu }_c{\tau }_1}+({\overline{D}}_0-{\overline{D}}_2{\tilde{\omega }}^2)sin\left(\tilde{\omega }{\tau }_1\right)e^{-{\mu }_c{\tau }_1}.\hfill \end{array}$

From (15), one has

(16)$\begin{array}{c}\hfill {\mathsf{Ϝ}}_1\left(\tilde{\omega }\right)+{\mathsf{Ϝ}}_2\left(\tilde{\omega }\right)cos\left(\tilde{\omega }{\tau }_1\right)+{\mathsf{Ϝ}}_3\left(\tilde{\omega }\right)sin\left(\tilde{\omega }{\tau }_1\right)=0,\end{array}$

where

$\begin{array}{cc}\hfill {\mathsf{Ϝ}}_1\left(\tilde{\omega }\right)=& \left(-{\overline{C}}_1^2+{\overline{E}}_1^2+2{\overline{C}}_0{\overline{C}}_2-2{\overline{E}}_0{\overline{E}}_2-{\overline{D}}_1^2{e}^{-2{\mu }_c{\tau }_1}+{\overline{F}}_1^2{e}^{-2{\mu }_c{\tau }_1}+2{\overline{D}}_0{\overline{D}}_2{e}^{-2{\mu }_c{\tau }_1}-2{\overline{F}}_0{\overline{F}}_2{e}^{-2{\mu }_c{\tau }_1}\right){\tilde{\omega }}^2\hfill \\ \hfill & +\left(-2{\overline{C}}_0-{\overline{C}}_2^2+{\overline{E}}_2^2+2{\overline{C}}_1{\overline{C}}_3-2{\overline{E}}_1{\overline{E}}_3-{\overline{D}}_2^2{e}^{-2{\mu }_c{\tau }_1}+{\overline{F}}_2^2{e}^{-2{\mu }_c{\tau }_1}+2{\overline{D}}_1{\overline{D}}_3{e}^{-2{\mu }_c{\tau }_1}\right){\tilde{\omega }}^4\hfill \\ \hfill & +{\tilde{\omega }}^8+\left(2{\overline{C}}_2+{\overline{E}}_3^2-{\overline{C}}_3^2-{\overline{D}}_3^2{e}^{-2{\mu }_c{\tau }_1}\right){\tilde{\omega }}^6+{\overline{E}}_0^2-{\overline{C}}_0^2+{\overline{F}}_0^2{e}^{-2{\mu }_c{\tau }_1}-{\overline{D}}_0^2{e}^{-2{\mu }_c{\tau }_1},\hfill \\ \hfill {\mathsf{Ϝ}}_2\left(\tilde{\omega }\right)=& (2{\overline{D}}_2-2{\overline{C}}_3{\overline{D}}_3)e^{-2{\mu }_c{\tau }_1}{\tilde{\omega }}^6+\left(-2{\overline{D}}_0-2{\overline{C}}_1{\overline{D}}_3-2{\overline{C}}_2{\overline{D}}_2+2{\overline{C}}_3{\overline{D}}_1+2{\overline{E}}_2{\overline{F}}_2-2{\overline{E}}_3{\overline{F}}_1\right)e^{-2{\mu }_c{\tau }_1}{\tilde{\omega }}^4\hfill \\ \hfill & +\left(2{\overline{C}}_0{\overline{D}}_2-2{\overline{C}}_1{\overline{D}}_1+2{\overline{C}}_2{\overline{D}}_0-2{\overline{E}}_0{\overline{F}}_2+2{\overline{E}}_1{\overline{F}}_1-2{\overline{E}}_2{\overline{F}}_0\right)e^{-2{\mu }_c{\tau }_1}{\tilde{\omega }}^2+\left(2{\overline{E}}_0{\overline{F}}_0-2{\overline{C}}_0{\overline{D}}_0\right)e^{-2{\mu }_c{\tau }_1},\hfill \\ \hfill {\mathsf{Ϝ}}_3\left(\tilde{\omega }\right)=& 2{\overline{D}}_3{e}^{-2{\mu }_c{\tau }_1}{\tilde{\omega }}^7+\left(2{\overline{C}}_0{\overline{D}}_3-2{\overline{C}}_1{\overline{D}}_2+2{\overline{C}}_2{\overline{D}}_1-2{\overline{C}}_3{\overline{D}}_0+2{\overline{E}}_1{\overline{F}}_2-2{\overline{E}}_2{\overline{F}}_1+2{\overline{E}}_3{\overline{F}}_0\right)e^{-2{\mu }_c{\tau }_1}{\tilde{\omega }}^3\hfill \\ \hfill & +\left(-2{\overline{D}}_1-2{\overline{C}}_2{\overline{D}}_3+2{\overline{C}}_3{\overline{D}}_2-2{\overline{E}}_3{\overline{F}}_2\right)e^{-2{\mu }_c{\tau }_1}{\tilde{\omega }}^5+\left(2{\overline{C}}_1{\overline{D}}_0-2{\overline{C}}_0{\overline{D}}_1+2{\overline{E}}_0{\overline{F}}_1-2{\overline{E}}_1{\overline{F}}_0\right)e^{-2{\mu }_c{\tau }_1}\tilde{\omega }.\hfill \end{array}$

Denote $\mathsf{\Gamma }\left(\tilde{\omega }\right)={\mathsf{Ϝ}}_1\left(\tilde{\omega }\right)+{\mathsf{Ϝ}}_2\left(\tilde{\omega }\right)cos\left(\tilde{\omega }{\tau }_1\right)+{\mathsf{Ϝ}}_3\left(\tilde{\omega }\right)sin\left(\tilde{\omega }{\tau }_1\right)$. Obviously, ${lim}_{\tilde{\omega }\to +\infty }\mathsf{\Gamma }\left(\tilde{\omega }\right)=+\infty$. Further, if $\left[{({\overline{E}}_0+{\overline{F}}_0)}^2-{({\overline{C}}_0+{\overline{D}}_0)}^2\right]e^{-2{\mu }_c{\tau }_1}+({\overline{E}}_0^2-{\overline{C}}_0^2)(1-e^{-2{\mu }_c{\tau }_1})<0$, that is, $\mathsf{\Gamma }\left(0\right)<0$, then Equation (16) has ${\tilde{k}}_0$ ($1\le {\tilde{k}}_0\le 8$) positive roots ${\tilde{\omega }}_k$. For each ${\tilde{\omega }}_k$, there is a corresponding critical value of delay ${\tau }_2$ shown as

${\tau }_{2,k}^{\left(j\right)}=\frac{1}{{\tilde{\omega }}_k}arccos\left(\frac{{\mathsf{\Pi }}_1{\mathsf{\Pi }}_3+{\mathsf{\Pi }}_2{\mathsf{\Pi }}_4}{{\mathsf{\Pi }}_1^2+{\mathsf{\Pi }}_2^2}\right)+\frac{2j\pi }{{\tilde{\omega }}_k},j=0,1,2,\dots .$

Let ${\tau }_2^{*}={\tau }_{2,{\tilde{k}}_0}^{\left(0\right)}=min\left\{{\tau }_{2,\tilde{k}}^{\left(j\right)}|\tilde{k}=1,2,\dots ,{\tilde{k}}_0,j=0,1,\dots \right\}$, and let ${\tilde{\omega }}_0={\tilde{\omega }}_{{\tilde{k}}_0}$ be the root of (9) with ${\tau }_2^{*}$. Differentiating both sides of (9) with respect to ${\tau }_2$, we obtain

$\begin{array}{cc}\hfill {\left(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2}\right)}^{-1}=& -\frac{4{\lambda }^3+3{\overline{C}}_3{\lambda }^2+2{\overline{C}}_2\lambda +{\overline{C}}_1+\left[(3{\overline{D}}_3{\lambda }^2+2{\overline{D}}_2\lambda +{\overline{D}}_1)-{\tau }_1({\overline{D}}_3{\lambda }^3+{\overline{D}}_2{\lambda }^2+{\overline{D}}_1\lambda +{\overline{D}}_0)\right]\mathbb{E}}{\lambda \left[{\lambda }^4+{\overline{C}}_3{\lambda }^3+{\overline{C}}_2{\lambda }^2+{\overline{C}}_1\lambda +{\overline{C}}_0+({\overline{D}}_3{\lambda }^3+{\overline{D}}_2{\lambda }^2+{\overline{D}}_1\lambda +{\overline{D}}_0)\mathbb{E}\right]}\hfill \\ \hfill & +\frac{2{\overline{E}}_3{\lambda }^2+2{\overline{E}}_2\lambda +{\overline{E}}_1+[2{\overline{F}}_2\lambda +{\overline{F}}_1-{\tau }_1({\overline{F}}_2{\lambda }^2+{\overline{F}}_1\lambda +{\overline{F}}_0)]\mathbb{E}}{\lambda [{\overline{E}}_3{\lambda }^3+{\overline{E}}_2{\lambda }^2+{\overline{E}}_1\lambda +{\overline{E}}_0+({\overline{F}}_2{\lambda }^2+{\overline{F}}_1\lambda +{\overline{F}}_0)\mathbb{E}]}-\frac{{\tau }_2}{\lambda },\hfill \end{array}$

where $\mathbb{E}=e^{-({\mu }_c+\lambda ){\tau }_1}$. Then evaluating ${\left(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2}\right)}^{-1}$ at ${\tau }_2={\tau }_2^{*}$$(\mathrm{i}.\mathrm{e}.$$\lambda =i{\tilde{\omega }}_0)$ and taking the real part, one has

$\begin{array}{c}\hfill \mathrm{Re}\left[{\left(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2}\right)}^{-1}{|}_{{\tau }_2={\tau }_2^{*}}\right]=\frac{C_{11}+D_{11}+E_{11}+F_{11}}{A_{11}^2+B_{11}^2},\end{array}$

where

$\begin{array}{cc}\hfill A_{11}=& {\overline{C}}_3{\tilde{\omega }}_0^4-{\overline{C}}_1{\tilde{\omega }}_0^2+({\overline{D}}_3{\tilde{\omega }}_0^4-{\overline{D}}_1{\tilde{\omega }}_0^2)e^{-{\mu }_c{\tau }_1}cos\left({\tilde{\omega }}_0{\tau }_1\right)+({\overline{D}}_0{\tilde{\omega }}_0-{\overline{D}}_2{\tilde{\omega }}_0^3)e^{-{\mu }_c{\tau }_1}sin\left({\tilde{\omega }}_0{\tau }_1\right),\hfill \\ \hfill B_{11}=& {\tilde{\omega }}_0^5-{\overline{C}}_2{\tilde{\omega }}_0^3+{\overline{C}}_0\tilde{\omega }-({\overline{D}}_3{\tilde{\omega }}_0^4-{\overline{D}}_1{\tilde{\omega }}_0^2)e^{-{\mu }_c{\tau }_1}sin\left({\tilde{\omega }}_0{\tau }_1\right)+({\overline{D}}_0{\tilde{\omega }}_0-{\overline{D}}_2{\tilde{\omega }}_0^3)e^{-{\mu }_c{\tau }_1}cos\left({\tilde{\omega }}_0{\tau }_1\right),\hfill \\ \hfill C_{11}=& \left\{-3{\overline{E}}_3{\tilde{\omega }}_0^2+{\overline{E}}_1+[{\overline{F}}_1-{\tau }_1({\overline{F}}_0-{\overline{F}}_2{\tilde{\omega }}_0^2)]e^{-{\mu }_c{\tau }_1}cos\left({\tilde{\omega }}_0{\tau }_1\right)+(2{\overline{F}}_2{\tilde{\omega }}_0-{\tau }_1{\overline{F}}_2{\tilde{\omega }}_0)e^{-{\mu }_c{\tau }_1}sin\left({\tilde{\omega }}_0{\tau }_1\right)\right\}\hfill \\ \hfill & \times \left[{\overline{E}}_3{\tilde{\omega }}_0^4-{\overline{E}}_1{\tilde{\omega }}_0^2-{\overline{F}}_1{\tilde{\omega }}_0^2{e}^{-{\mu }_c{\tau }_1}cos\left({\tilde{\omega }}_0{\tau }_1\right)+({\overline{F}}_0{\tilde{\omega }}_0-{\overline{F}}_2{\tilde{\omega }}_0^3)e^{-{\mu }_c{\tau }_1}sin\left({\tilde{\omega }}_0{\tau }_1\right)\right],\hfill \\ \hfill D_{11}=& \left\{2{\overline{E}}_2{\tilde{\omega }}_0+(2{\overline{F}}_2{\tilde{\omega }}_0-{\tau }_1{\overline{F}}_1{\tilde{\omega }}_0)e^{-{\mu }_c{\tau }_1}cos\left({\tilde{\omega }}_0{\tau }_1\right)+[{\tau }_1({\overline{F}}_0-{\overline{F}}_2{\tilde{\omega }}_0^2)-{\overline{F}}_1]e^{-{\mu }_c{\tau }_1}sin\left({\tilde{\omega }}_0{\tau }_1\right)\right\}\hfill \\ \hfill & \times \left[{\overline{E}}_0{\tilde{\omega }}_0-{\overline{E}}_2{\tilde{\omega }}_0^3+({\overline{F}}_0{\tilde{\omega }}_0-{\overline{F}}_2{\tilde{\omega }}_0^3)e^{-{\mu }_c{\tau }_1}cos\left({\tilde{\omega }}_0{\tau }_1\right)+{\overline{F}}_1{\tilde{\omega }}_0^2{e}^{-{\mu }_c{\tau }_1}sin\left({\tilde{\omega }}_0{\tau }_1\right)\right],\hfill \\ \hfill E_{11}=& \{3{\overline{E}}_3{\tilde{\omega }}_0^2-{\overline{C}}_1+[3{\overline{D}}_3{\tilde{\omega }}_0^2-{\overline{D}}_1+{\tau }_1({\overline{D}}_0-{\overline{D}}_2{\tilde{\omega }}_0^2)]e^{-{\mu }_c{\tau }_1}cos\left({\tilde{\omega }}_0{\tau }_1\right)+[{\tau }_1({\overline{D}}_1{\tilde{\omega }}_0-{\overline{D}}_3{\tilde{\omega }}_0^3)-2{\overline{D}}_2{\tilde{\omega }}_0]e^{-{\mu }_c{\tau }_1}\hfill \\ \hfill & \times sin\left({\tilde{\omega }}_0{\tau }_1\right)\}\left[{\overline{C}}_3{\tilde{\omega }}_0^4-{\overline{C}}_1{\tilde{\omega }}_0^2+({\overline{D}}_3{\tilde{\omega }}_0^4-{\overline{D}}_1{\tilde{\omega }}_0^2)e^{-{\mu }_c{\tau }_1}cos\left({\tilde{\omega }}_0{\tau }_1\right)+({\overline{D}}_0{\tilde{\omega }}_0-{\overline{D}}_2{\tilde{\omega }}_0^3)e^{-{\mu }_c{\tau }_1}sin\left({\tilde{\omega }}_0{\tau }_1\right)\right],\hfill \\ \hfill F_{11}=& \{4{\tilde{\omega }}_0^3-2{\overline{C}}_2{\tilde{\omega }}_0+[{\tau }_1({\overline{D}}_1{\tilde{\omega }}_0-{\overline{D}}_3{\tilde{\omega }}_0^3)-2{\overline{D}}_2{\tilde{\omega }}_0]e^{-{\mu }_c{\tau }_1}cos\left({\tilde{\omega }}_0{\tau }_1\right)+[{\tau }_1({\overline{D}}_2{\tilde{\omega }}_0^2-{\overline{D}}_0)+{\overline{D}}_1-3{\overline{D}}_3{\tilde{\omega }}_0^2]e^{-{\mu }_c{\tau }_1}\hfill \\ \hfill & \times sin\left({\tilde{\omega }}_0{\tau }_1\right)\}[{\tilde{\omega }}_0^5-{\overline{C}}_2{\tilde{\omega }}_0^3+{\overline{C}}_0\tilde{\omega }-({\overline{D}}_3{\tilde{\omega }}_0^4-{\overline{D}}_1{\tilde{\omega }}_0^2)e^{-{\mu }_c{\tau }_1}sin\left({\tilde{\omega }}_0{\tau }_1\right)+({\overline{D}}_0{\tilde{\omega }}_0-{\overline{D}}_2{\tilde{\omega }}_0^3)e^{-{\mu }_c{\tau }_1}cos\left({\tilde{\omega }}_0{\tau }_1\right)].\hfill \end{array}$

Therefore,

$\begin{array}{c}\hfill \mathrm{sign}\left\{\frac{\mathrm{d}\left(\mathrm{Re}\right(\lambda \left)\right)}{\mathrm{d}{\tau }_2}{|}_{{\tau }_2={\tau }_2^{*}}\right\}=\mathrm{sign}\left\{\mathrm{Re}\left[{\left(\frac{\mathrm{d}\lambda }{\mathrm{d}{\tau }_2}\right)}^{-1}{|}_{{\tau }_2={\tau }_2^{*}}\right]\right\}=\mathrm{sign}\left\{C_{11}+D_{11}+E_{11}+F_{11}\right\}.\end{array}$

Based on the above arguments, we can derive the following results.

Theorem 5.

For Model (4), given that ${\tau }_1>0$ and ${\tau }_2>0$, if $\mathsf{\Gamma }\left(0\right)<0$ and $\mathrm{sign}\{C_{11}+D_{11}+E_{11}+F_{11}\}\ne 0$, then ${\tilde{U}}_2$ is locally asymptotically stable for all ${\tau }_2\in [0,{\tau }_2^{*})$, and unstable for ${\tau }_2>{\tau }_2^{*}$. Additionally, Model (4) undergoes a Hopf bifurcation at ${\tilde{U}}_2$ when ${\tau }_2={\tau }_2^{*}$.

3.3. Direction and Stability of the Hopf Bifurcation

In this section, we will employ the standard form theory and the central manifold theorem presented in [38] to analyze the direction of the Hopf bifurcation and investigate the stability of the bifurcating periodic solution in Model (4) with respect to ${\tau }_2$. It is always assumed that Model (4) undergoes a Hopf bifurcation at ${\tilde{U}}_2$ when ${\tau }_1>0$ and ${\tau }_2={\tau }_2^{*}$. Based on the biological context, it is useful to assume that ${\tau }_1<{\tau }_2^{*}$, where ${\tau }_2\in [0,{\tau }_2^{*})$.

Let ${\tau }_2={\tau }_2^{*}+\tilde{\mu }$, $\tilde{\mu }\in \mathbb{R}$, then Model (4) will generate Hopf bifurcation at ${\tilde{U}}_2$ when $\tilde{\mu }=0$. Denote $\tilde{x}\left(s\right)={({\tilde{x}}_1\left(s\right),{\tilde{x}}_2\left(s\right),{\tilde{x}}_3\left(s\right),{\tilde{x}}_4\left(s\right))}^T\in {\mathbb{R}}^4$, ${\tilde{x}}_1\left(s\right)=T-{\tilde{T}}_2$, ${\tilde{x}}_2\left(s\right)=T^{*}-{\tilde{T}}_2^{*}$, ${\tilde{x}}_3\left(s\right)=V-{\tilde{V}}_2$, ${\tilde{x}}_4\left(s\right)=B-{\tilde{B}}_2$, normalizing the delay $s\to \left(\frac{s}{{\tau }_2}\right)$, then Model (4) can be transformed into

(17)$\begin{array}{c}\hfill \frac{\mathrm{d}\tilde{x}\left(s\right)}{\mathrm{d}s}=L_{\tilde{\mu }}\left({\tilde{x}}_s\right)-f(\tilde{\mu },{\tilde{x}}_s),\end{array}$

where maps $L_{\tilde{\mu }}:C\to {\mathbb{R}}^4$ and $f:\mathbb{R}\times C\to {\mathbb{R}}^4$, respectively, are given by

$\begin{array}{cc}\hfill L_{\tilde{\mu }}\left(\varphi \right)& =({\tau }_2^{*}+\tilde{\mu })\left(N_1\varphi \left(0\right)+N_2\varphi \left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right)+N_3\varphi (-1)\right),\hfill \\ \hfill f(\tilde{\mu },\varphi )& =({\tau }_2^{*}+\tilde{\mu })\left(\begin{array}{c}-{\kappa }_1{\varphi }_1\left(0\right){\varphi }_3\left(0\right)-{\kappa }_2{\varphi }_1\left(0\right){\varphi }_2\left(0\right)\\ {\kappa }_1{\mathrm{e}}^{-{\mu }_c{\tau }_1}{\varphi }_1\left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right){\varphi }_3\left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right)+{\kappa }_2{\mathrm{e}}^{-{\mu }_c{\tau }_1}{\varphi }_1\left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right){\varphi }_2\left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right)\\ -q{\varphi }_3\left(0\right){\varphi }_4\left(0\right)\\ \sigma {\varphi }_3(-1){\varphi }_4(-1)\end{array}\right),\hfill \end{array}$

where $\varphi \left(\theta \right)=({\varphi }_1\left(\theta \right),{\varphi }_2\left(\theta \right),{\varphi }_3\left(\theta \right),{\varphi }_4\left(\theta \right))\in C$, and

$\begin{array}{c}\hfill N_1=\left(\begin{array}{cccc}-{\kappa }_1{\tilde{V}}_2-{\kappa }_2{\tilde{T}}_2^{*}-{\mu }_c& -{\kappa }_2{\tilde{T}}_2^{*}& -{\kappa }_1{\tilde{V}}_2& 0\\ 0& -({\mu }_h+{\delta }_c)& 0& 0\\ 0& p& -q{\tilde{B}}_2-{\mu }_v& -q{\tilde{V}}_2\\ 0& 0& 0& -\zeta \end{array}\right),\end{array}$

$\begin{array}{c}\hfill N_2=\left(\begin{array}{cccc}0& 0& 0& 0\\ ({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}){\mathrm{e}}^{-{\mu }_c{\tau }_1}& {\kappa }_2{\tilde{T}}_2{\mathrm{e}}^{-{\mu }_c{\tau }_1}& {\kappa }_1{\tilde{T}}_2{\mathrm{e}}^{-{\mu }_c{\tau }_1}& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\end{array}\right),N_3=\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& \sigma {\tilde{B}}_2& \sigma {\tilde{V}}_2\end{array}\right).\end{array}$

Hence, from the Riesz representation theorem, there is a $4\times 4$ matrix function $\eta (\theta ,\tilde{\mu }):[-1,0]\to {\mathbb{R}}^3$ whose components are bounded, such that $L_{\tilde{\mu }}\left(\varphi \right)=$${\int }_{-1}^0\mathrm{d}\eta (\theta ,\tilde{\mu })\varphi \left(\theta \right)$ for $\varphi \in C$. In fact, we can choose

(18)$\eta (\theta ,\tilde{\mu })=\left\{\begin{array}{cc}({\tau }_2^{*}+\tilde{\mu })(N_1+N_2+N_3),\hfill & \hfill \theta =0,\\ ({\tau }_2^{*}+\tilde{\mu })(N_2+N_3),\hfill & \hfill \theta \in \left[\frac{-{\tau }_1}{{\tau }_2^{*}},0\right),\\ ({\tau }_2^{*}+\tilde{\mu })N_3,\hfill & \hfill \theta \in \left(-1,\frac{-{\tau }_1}{{\tau }_2^{*}}\right),\\ 0,\hfill & \hfill \theta =-1.\end{array}\right.$

For any $\varphi \in C([-1,0],{\mathbb{R}}^4)$, the adjoint operator, $\mathcal{A}$, of $\mathcal{A}\left(0\right)$ is defined by

$\mathcal{A}\left(\tilde{\mu }\right)\varphi =\left\{\begin{array}{cc}\frac{\mathrm{d}\varphi \left(\theta \right)}{\mathrm{d}\theta },\hfill & \hfill \theta \in [-1,0),\\ {\int }_{-1}^0\mathrm{d}\eta (\theta ,\tilde{\mu })\varphi \left(\theta \right),\hfill & \hfill \theta =0,\end{array}\right.\mathcal{R}\left(\tilde{\mu }\right)\varphi =\left\{\begin{array}{cc}0,\hfill & \hfill \theta \in [-1,0),\\ f(\tilde{\mu },\varphi ),\hfill & \hfill \theta =0.\end{array}\right.$

Then, Model (17) is equivalent to the following operator equation

$\begin{array}{c}\hfill \frac{\mathrm{d}{\tilde{x}}_s\left(\theta \right)}{\mathrm{d}s}=\mathcal{A}\left(\tilde{\mu }\right){\tilde{x}}_s+\mathcal{R}\left(\tilde{\mu }\right){\tilde{x}}_s,\end{array}$

where ${\tilde{x}}_s\left(\theta \right)=\tilde{x}(s+\theta )$ for $\theta \in [-1,0]$. For any $\psi \in C([0,1],{\left({\mathbb{R}}^4\right)}^{*})$, define

${\mathcal{A}}^{*}\psi \left(l\right)=\left\{\begin{array}{cc}-\frac{\mathrm{d}\psi \left(l\right)}{\mathrm{d}l},\hfill & \hfill l\in (0,1],\\ {\int }_{-1}^0\mathrm{d}{\eta }^T(l,0)\varphi (-l),\hfill & \hfill l=0,\end{array}\right.$

and a bilinear inner product is defined by

(19)$\begin{array}{c}\hfill ⟨\psi \left(l\right),\varphi \left(\theta \right)⟩=\overline{\psi }\left(0\right)\varphi \left(0\right)-{\int }_{-1}^0{\int }_{\xi =0}^{\theta }\overline{\psi }(\xi -\theta )\mathrm{d}\eta \left(\theta \right)\varphi \left(\xi \right)\mathrm{d}\xi ,\end{array}$

where $\eta \left(\theta \right)=\eta (\theta ,0)$. We assume that $\pm i{\omega }_2^{*}{\tau }_2^{*}$ are eigenvalues of $\mathcal{A}\left(0\right)$, and the other eigenvalues have strictly negative real parts. Hence, they are eigenvalues of ${\mathcal{A}}^{*}$ as well.

If $q\left(\theta \right)={(1,q_2,q_3,q_4)}^T{e}^{i\theta {\omega }_2^{*}{\tau }_2^{*}}$ is the eigenvector of $\mathcal{A}\left(0\right)$ corresponding to eigenvalue $i{\omega }_2^{*}{\tau }_2^{*}$, then $\mathcal{A}\left(0\right)q\left(\theta \right)=i{\omega }_2^{*}{\tau }_2^{*}q\left(\theta \right)$. From (18) and the definition of $\mathcal{A}\left(0\right)$, we have

${\tau }_2^{*}\left(\begin{array}{cccc}{a}_{11}& a_{12}& a_{13}& 0\\ a_{21}& a_{22}& a_{23}& 0\\ 0& a_{32}& a_{33}& a_{34}\\ 0& 0& a_{43}& a_{44}\end{array}\right)\left(\begin{array}{c}1\\ q_2\\ q_3\\ q_4\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\end{array}\right),$

where $a_{11}=i{\omega }_2^{*}+{\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c$, $a_{12}={\kappa }_2{\tilde{T}}_2$, $a_{13}={\kappa }_1{\tilde{T}}_2$, $a_{21}=-({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*})e^{-{\chi }_1}$, $a_{22}=i{\omega }_2^{*}-{\kappa }_2{\tilde{T}}_2{e}^{-{\chi }_1}+{\mu }_c+{\delta }_c$, $a_{23}=-{\kappa }_1{\tilde{T}}_2{e}^{-{\chi }_1}$, ${\chi }_1={\mu }_c{\tau }_1+i{\omega }_2^{*}{\tau }_1$, $a_{32}=-p$, $a_{33}=i{\omega }_2^{*}+q{\tilde{B}}_2+{\mu }_v$, $a_{34}=q{\tilde{V}}_2$, $a_{43}=-\sigma {\tilde{B}}_2{e}^{-i{\omega }_2^{*}{\tau }_2^{*}}$, $a_{44}=i{\omega }_2^{*}-\sigma {\tilde{V}}_2{e}^{-i{\omega }_2^{*}{\tau }_2^{*}}+\zeta$. Therefore, we can obtain $q{\left(0\right)}^T={(1,q_2,q_3,q_4)}^T$, where

$\begin{array}{cc}\hfill q_2=& \frac{-(i{\omega }_2^{*}+{\mu }_c)e^{-{\chi }_1}}{i{\omega }_2^{*}+{\mu }_c+{\delta }_c},q_3=\frac{{\kappa }_2(i{\omega }_2^{*}+{\mu }_c)e^{-{\chi }_1}}{{\kappa }_1(i{\omega }_2^{*}+{\mu }_c+{\delta }_c)}-\frac{i{\omega }_2^{*}+{\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c}{{\kappa }_1{\tilde{T}}_2},\hfill \\ \hfill q_4=& \frac{\sigma {\tilde{B}}_2{e}^{-i{\omega }_2^{*}{\tau }_2^{*}}[{\kappa }_2{\tilde{T}}_2(i{\omega }_2^{*}+{\mu }_c)e^{-{\chi }_1}-(i{\omega }_2^{*}+{\mu }_c+{\delta }_c)(i{\omega }_2^{*}+{\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c)]}{{\kappa }_1{\tilde{T}}_2(i{\omega }_2^{*}+{\mu }_c+{\delta }_c)(i{\omega }_2^{*}-\sigma {\tilde{V}}_2{e}^{-i{\omega }_2^{*}{\tau }_2^{*}}+\zeta )}.\hfill \end{array}$

Similarly, according to (18) and the definition of ${\mathcal{A}}^{*}$, if $q^{*}\left(l\right)=D(1,q_2^{*},q_3^{*},q_4^{*})e^{il{\omega }_2^{*}{\tau }_2^{*}}$ is the eigenvector of ${\mathcal{A}}^{*}$ corresponding to eigenvalue $-i{\omega }_2^{*}{\tau }_2^{*}$, we further have

$\begin{array}{cc}\hfill q_2^{*}& =\frac{{\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c-i{\omega }_2^{*}}{({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*})e^{-{\chi }_1}},q_3^{*}=\frac{({\mu }_c+{\delta }_c-i{\omega }_2^{*}-{\kappa }_2{\tilde{T}}_2{e}^{-{\chi }_1})({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c-i{\omega }_2^{*})}{p({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*})e^{-{\chi }_1}}+\frac{{\kappa }_2{\tilde{T}}_2}{p},\hfill \\ \hfill q_4^{*}& =\frac{q{\kappa }_2{\tilde{T}}_2{\tilde{V}}_2({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*})e^{-{\chi }_1}+q{\tilde{V}}_2({\mu }_c+{\delta }_c-i{\omega }_2^{*}-{\kappa }_2{\tilde{T}}_2{e}^{-{\chi }_1})({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c-i{\omega }_2^{*})}{p(\sigma {\tilde{V}}_2{e}^{-i{\omega }_2^{*}{\tau }_2^{*}}-\zeta +i{\omega }_2^{*})({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*})e^{-{\chi }_1}}.\hfill \end{array}$

In order to guarantee $⟨q^{*}\left(l\right),q\left(\theta \right)⟩=1$, the value of D needs to be determined. From (19), one has

$\begin{array}{cc}\hfill ⟨q^{*}\left(l\right),q\left(\theta \right)⟩=& \overline{D}(1,{\overline{q}}_2^{*},{\overline{q}}_3^{*},{\overline{q}}_4^{*}){(1,q_2,q_3,q_4)}^T-{\int }_{-1}^0{\int }_0^{\theta }\overline{D}(1,{\overline{q}}_2^{*},{\overline{q}}_3^{*},{\overline{q}}_4^{*})e^{-i(\xi -\theta ){\omega }_2^{*}{\tau }_2^{*}}\mathrm{d}\eta \left(\theta \right){(1,q_2,q_3,q_4)}^T{e}^{i\xi {\omega }_2^{*}{\tau }_2^{*}}\mathrm{d}\xi \hfill \\ \hfill =& \overline{D}\left[1+q_2{\overline{q}}_2^{*}+q_3{\overline{q}}_3^{*}+q_4{\overline{q}}_4^{*}+{\tau }_1{q}_2^{*}{e}^{-{\chi }_1}({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+q_2{\kappa }_2{\tilde{T}}_2+q_3{\kappa }_1{\tilde{T}}_2)+{\tau }_2^{*}{\overline{q}}_4^{*}{e}^{-i{\omega }_2^{*}{\tau }_2^{*}}(q_3\sigma {\tilde{B}}_2+q_4\sigma {\tilde{V}}_2)\right],\hfill \end{array}$

where $\overline{D}$ and ${\overline{q}}_k^{*}(k=2,3,4)$ are the conjugates of D and $q_k^{*}(k=2,3,4)$, respectively. Thus,

$\begin{array}{c}\hfill D={\left[1+\sum _{k=2}^4{\overline{q}}_k{q}_k^{*}+{\tau }_1{\overline{q}}_2^{*}{e}^{-{\mu }_c{\tau }_1+i{\omega }_2^{*}{\tau }_1}({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\overline{q}}_2{\kappa }_2{\tilde{T}}_2+{\overline{q}}_3{\kappa }_1{\tilde{T}}_2)+{\tau }_2^{*}{q}_4^{*}{e}^{i{\omega }_2^{*}{\tau }_2^{*}}({\overline{q}}_3\sigma {\tilde{B}}_2+{\overline{q}}_4\sigma {\tilde{V}}_2)\right]}^{-1},\end{array}$

where ${\overline{q}}_k(k=2,3,4)$ are the conjugates of $q_k(k=2,3,4)$, respectively.

Next, attention is focused on the coefficients for determining the Hopf bifurcation direction and the stability of the bifurcation period solution, which, according to the algorithm given in [39] and similar argumentation process in [40], are given by

$\begin{array}{cc}\hfill g_{20}=& 2{\tau }_2^{*}\overline{D}\left[-{\kappa }_1{q}_3-{\kappa }_2{q}_2+{\kappa }_1{\overline{q}}_2^{*}{q}_3{e}^{-{\chi }_2}+{\kappa }_2{q}_2{\overline{q}}_2^{*}{e}^{-{\chi }_2}-q{\overline{q}}_3^{*}{q}_3{q}_4+{\overline{q}}_4^{*}{q}_3{q}_4\sigma e^{-2i{\omega }_2^{*}{\tau }_2^{*}}\right],\hfill \\ \hfill g_{11}=& {\tau }_2^{*}\overline{D}[-{\kappa }_1{\overline{q}}_3-{\kappa }_1{q}_3-{\kappa }_2{\overline{q}}_2-{\kappa }_2{q}_2+({\kappa }_1{\overline{q}}_3{\overline{q}}_2^{*}+{\kappa }_1{q}_3{\overline{q}}_2^{*}+{\kappa }_2{\overline{q}}_2{\overline{q}}_2^{*}+{\kappa }_2{q}_2{\overline{q}}_2^{*})e^{-{\mu }_c{\tau }_1}\hfill \\ \hfill & -q{\overline{q}}_3^{*}{q}_3{\overline{q}}_4-q{\overline{q}}_3^{*}{\overline{q}}_3{q}_4+{\overline{q}}_4^{*}\sigma q_3{\overline{q}}_4+{\overline{q}}_4^{*}\sigma {\overline{q}}_3{q}_4],\hfill \\ \hfill g_{02}=& 2{\tau }_2^{*}\overline{D}\left[-{\kappa }_1{\overline{q}}_3-{\kappa }_2{\overline{q}}_2+{\kappa }_1{\overline{q}}_2^{*}{\overline{q}}_3{e}^{-{\mu }_c{\tau }_1+2i{\omega }_2^{*}{\tau }_1}+{\kappa }_2{\overline{q}}_2^{*}{\overline{q}}_2{e}^{-{\mu }_c{\tau }_1+2i{\omega }_2^{*}{\tau }_1}-q_3^{*}q{\overline{q}}_3{\overline{q}}_4+{\overline{q}}_4^{*}\sigma {\overline{q}}_3{\overline{q}}_4{e}^{2i{\omega }_2^{*}{\tau }_2^{*}}\right],\hfill \\ \hfill g_{21}=& {\tau }_2^{*}\overline{D}\{-2{\kappa }_1{W}_{11}^{\left(3\right)}\left(0\right)-{\kappa }_1{W}_{20}^{\left(3\right)}\left(0\right)-{\kappa }_1{\overline{q}}_3{W}_{20}^{\left(1\right)}\left(0\right)-2{\kappa }_1{q}_3{W}_{11}^{\left(1\right)}\left(0\right)-2{\kappa }_2{W}_{11}^{\left(2\right)}\left(0\right)-{\kappa }_2{W}_{20}^{\left(2\right)}\left(0\right)\hfill \\ \hfill & +{\overline{q}}_4^{*}\sigma \left[2q_3{W}_{11}^{\left(4\right)}(-1)e^{-i{\omega }_2^{*}{\tau }_2^{*}}+{\overline{q}}_3{W}_{20}^{\left(4\right)}(-1)e^{i{\omega }_2^{*}{\tau }_2^{*}}+{\overline{q}}_4{W}_{20}^{\left(3\right)}(-1)e^{i{\omega }_2^{*}{\tau }_2^{*}}+2q_4{W}_{11}^{\left(3\right)}(-1)e^{-i{\omega }_2^{*}{\tau }_2^{*}}\right]-{\kappa }_2{\overline{q}}_2{W}_{20}^{\left(1\right)}\left(0\right)\hfill \\ \hfill & -2{\kappa }_2{q}_2{W}_{11}^{\left(1\right)}\left(0\right)+{\overline{q}}_2^{*}{e}^{-{\mu }_c{\tau }_1}[2{\kappa }_1{W}_{11}^{\left(3\right)}\left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right)e^{-i{\omega }_2^{*}{\tau }_1}+{\kappa }_1{W}_{20}^{\left(3\right)}\left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right)e^{i{\omega }_2^{*}{\tau }_1}+{\overline{q}}_3{\kappa }_1{W}_{20}^{\left(1\right)}\left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right)e^{i{\omega }_2^{*}{\tau }_1}\hfill \\ \hfill & +2{\kappa }_1{q}_3{W}_{11}^{\left(1\right)}\left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right)e^{-i{\omega }_2^{*}{\tau }_1}+2{\kappa }_2{W}_{11}^{\left(2\right)}\left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right)e^{-i{\omega }_2^{*}{\tau }_1}+{\kappa }_2{W}_{20}^{\left(2\right)}\left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right)e^{i{\omega }_2^{*}{\tau }_1}+{\kappa }_2{\overline{q}}_2{W}_{20}^{\left(1\right)}\left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right)e^{i{\omega }_2^{*}{\tau }_1}\hfill \\ \hfill & +2{\kappa }_2{q}_2{W}_{11}^{\left(1\right)}\left(\frac{-{\tau }_1}{{\tau }_2^{*}}\right)e^{-i{\omega }_2^{*}{\tau }_1}]-{\overline{q}}_3^{*}q\left[2q_3{W}_{11}^{\left(4\right)}\left(0\right)+{\overline{q}}_3{W}_{20}^{\left(4\right)}\left(0\right)+{\overline{q}}_4{W}_{20}^{\left(3\right)}\left(0\right)+2q_4{W}_{11}^{\left(3\right)}\left(0\right)\right]\},\hfill \end{array}$

with

$\begin{array}{cc}\hfill W_{20}\left(\theta \right)=& \frac{i{g}_{20}}{{\omega }_2^{*}{\tau }_2^{*}}q\left(0\right)e^{i{\omega }_2^{*}{\tau }_2^{*}\theta }+\frac{i{\overline{g}}_{02}}{3{\omega }_2^{*}{\tau }_2^{*}}\overline{q}\left(0\right)e^{-i{\omega }_2^{*}{\tau }_2^{*}\theta }+E_1{e}^{2i{\omega }_2^{*}{\tau }_2^{*}\theta },\hfill \\ \hfill W_{11}\left(\theta \right)=& -\frac{i{g}_{11}}{{\omega }_2^{*}{\tau }_2^{*}}q\left(0\right)e^{i{\omega }_2^{*}{\tau }_2^{*}\theta }+\frac{i{\overline{g}}_{11}}{{\omega }_2^{*}{\tau }_2^{*}}\overline{q}\left(0\right)e^{-i{\omega }_2^{*}{\tau }_2^{*}\theta }+E_2,\hfill \end{array}$

where $E_1$ and $E_2$ can be calculated as follows, respectively

$\begin{array}{cc}\hfill E_2\left(\begin{array}{cccc}{c}_{11}& a_{12}& a_{13}& 0\\ c_{21}& c_{22}& c_{23}& 0\\ 0& a_{32}& c_{33}& a_{34}\\ 0& 0& c_{43}& c_{44}\end{array}\right)& =2\left(\begin{array}{c}-{\kappa }_1(q_3+{\overline{q}}_3)-{\kappa }_2(q_2+{\overline{q}}_2)\\ ({\kappa }_1(q_3+{\overline{q}}_3)+{\kappa }_2(q_2+{\overline{q}}_2))e^{-{\mu }_c{\tau }_1}\\ -q({\overline{q}}_3{q}_4+q_3{\overline{q}}_4)\\ \sigma ({\overline{q}}_3{q}_4+q_3{\overline{q}}_4)\end{array}\right),\hfill \\ \hfill E_1\left(\begin{array}{cccc}{b}_{11}& a_{12}& a_{13}& 0\\ b_{21}& b_{22}& b_{23}& 0\\ 0& a_{32}& b_{33}& a_{34}\\ 0& 0& b_{43}& b_{44}\end{array}\right)& =2\left(\begin{array}{c}-{\kappa }_1{q}_3-{\kappa }_2{q}_2\\ ({\kappa }_1{q}_3+{\kappa }_2{q}_2)e^{-{\chi }_1}\\ -q{q}_3{q}_4\\ \sigma q_3{q}_4{e}^{-2i{\omega }_2^{*}{\tau }_2^{*}}\end{array}\right),\hfill \end{array}$

where $b_{11}=2i{\omega }_2^{*}+{\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c$, $b_{21}=-({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*})e^{-{\chi }_2}$, $b_{22}=2i{\omega }_2^{*}-{\kappa }_2{\tilde{T}}_2{e}^{-{\chi }_2}+{\mu }_c+{\delta }_c$, $b_{23}=-{\kappa }_1{\tilde{T}}_2{e}^{-{\chi }_2}$, ${\chi }_2={\mu }_c{\tau }_1+2i{\omega }_2^{*}{\tau }_1$, $b_{33}=2i{\omega }_2^{*}+q{\tilde{B}}_2+{\mu }_v$, $b_{43}=-\sigma {\tilde{B}}_2{e}^{-2i{\omega }_2^{*}{\tau }_2^{*}}$, $b_{44}=2i{\omega }_2^{*}-\sigma {\tilde{V}}_2{e}^{-2i{\omega }_2^{*}{\tau }_2^{*}}+\zeta$, $c_{11}={\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*}+{\mu }_c$, $c_{21}=-({\kappa }_1{\tilde{V}}_2+{\kappa }_2{\tilde{T}}_2^{*})e^{-{\mu }_c{\tau }_1^{*}}$, $c_{22}=-{\kappa }_2{\tilde{T}}_2{e}^{-{\mu }_c{\tau }_1}+{\mu }_c+{\delta }_c$, $c_{23}=-{\kappa }_1{\tilde{T}}_2{e}^{-{\mu }_c{\tau }_1}$, $c_{33}=q{\tilde{B}}_2+{\mu }_v$, $a_{34}=q{\tilde{V}}_2$, $c_{43}=-\sigma {\tilde{B}}_2$, $c_{44}=-\sigma {\tilde{V}}_2+\zeta$. Therefore, we can calculate the following values:

$\begin{array}{cc}\hfill c_1\left(0\right)& =\frac{i}{2{\omega }_2^{*}{\tau }_2^{*}}(g_{11}{g}_{20}-2{|g_{11}|}^2-\frac{{|g_{02}|}^2}{3})+\frac{g_{21}}{2},{\mu }_2=-\frac{\mathrm{Re}\left\{c_1\left(0\right)\right\}}{\mathrm{Re}\left\{{\lambda }^{\prime }\left({\tau }_2^{*}\right)\right\}},\hfill \\ \hfill \breve{\beta }& =2\mathrm{Re}\left\{c_1\left(0\right)\right\},{\breve{T}}_2=-\frac{\mathrm{Im}\left\{c_1\left(0\right)\right\}+{\mu }_2\mathrm{Im}\left\{{\lambda }^{\prime }\left({\tau }_2^{*}\right)\right\}}{{\omega }_2^{*}{\tau }_2^{*}}.\hfill \end{array}$

Based on the preceding analysis, we can derive the following result.

Theorem 6.

If ${\mu }_2>0$$({\mu }_2<0)$, then the Hopf bifurcation is a supercritical bifurcation (subcritical bifurcation); if $\check{\beta }<0$$(\check{\beta }>0)$, then the bifurcated periodic solution is stable (unstable); if ${\breve{T}}_2>0$$({\breve{T}}_2<0)$, then the period of bifurcating periodic solutions increases (decreases).

4. Dynamics of Coupled Slow Time Model without Immune Response

In this section, we investigate the scenario where there is no immune response present in the host, that is, $B\left(s\right)=0$ in the within-host virus infection model (4). Further, in order to establish a meaningful coupling between the within-host model and the between-host model, we give attention to the case I, $A>0$ in Model (4). Therefore, it is assumed that ${\mathcal{R}}_w\le 1$ in the following discussion. Since the immunological process occurring within the host is much faster than the epidemiological process occurring between hosts, we always suppose that the state of the fast time model has reached its stable equilibrium while the state of the slow time model is without further change. Hence, based on the the stability of equilibrium ${\tilde{U}}_1$, one can assume that viral load $V\left(s\right)\to {\tilde{V}}_1$ as $s\to \infty$. Let $V={\tilde{V}}_1$ in Model (3), then we have the following model

(20)$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}& ={\mathsf{\Lambda }}_h-\frac{c_1{\beta }_1\left({\tilde{V}}_1\right)I\left(t\right)+c_2{\beta }_2\left({\tilde{V}}_1\right)A\left(t\right)}{1+\alpha \left(I\right(t)+A(t\left)\right)}S\left(t\right)-\mu S\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=& \frac{c_1{\beta }_1\left({\tilde{V}}_1\right)I(t-\tau )+c_2{\beta }_2\left({\tilde{V}}_1\right)A(t-\tau )}{1+\alpha \left(I\right(t-\tau )+A(t-\tau \left)\right)}S(t-\tau )e^{-\mu \tau }\hfill \\ \hfill & -(\mu +\xi +{\delta }_1\left({\tilde{V}}_1\right))I\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}=& \xi I\left(t\right)-(\mu +{\delta }_2\left({\tilde{V}}_1\right))A\left(t\right),\hfill \end{array}\right.$

with the initial condition

(21)$\begin{array}{cc}\hfill & S\left(\varrho \right)={\varphi }_1\left(\varrho \right)>0,I\left(\varrho \right)={\varphi }_3\left(\varrho \right)\ge 0,A\left(\varrho \right)={\varphi }_4\left(\varrho \right)\ge 0,\hfill \\ \hfill & \varrho \in [-\tau ,0],{\varphi }_1\left(0\right)>0,{\varphi }_i\left(0\right)\ge 0,i=3,4,\hfill \end{array}$

where $\varphi =({\varphi }_1,{\varphi }_3,{\varphi }_4)\in C([-\tau ,0],{\mathbb{R}}_{+}^3)$. Clearly, Model (20) has a unique solution $\left(S\right(t),I(t),A(t\left)\right)$ satisfying the initial condition (21), according to the fundamental theory of functional differential equations.

On the nonnegativity and ultimate boundedness of solutions for Model (4), the following result is straightforward.

Theorem 7.

Under the initial condition (21), any solution $\left(S\right(t),I(t),A(t\left)\right)$ of Model (20) is nonnegative and ultimately bounded for all $t\ge 0$. Particularly, ${lim\; sup}_{t\to \infty }(S(t-\tau )+I\left(t\right)e^{\mu \tau }+A\left(t\right)e^{\mu \tau })\le \frac{{\mathsf{\Lambda }}_h}{\mu }$, for all $t\ge 0$.

We denote the between-host reproduction number for Model (20) as follows

${\mathcal{R}}_0^h=\frac{(c_1{\beta }_1\left(V_1\right)(\mu +{\delta }_2\left(V_1\right))+c_2{\beta }_2\left(V_1\right)\xi ){\mathsf{\Lambda }}_h{e}^{-\mu \tau }}{\mu (\mu +{\delta }_2\left(V_1\right))(\mu +\xi +{\delta }_1\left(V_1\right))}.$

Denote $S_0=\frac{{\mathsf{\Lambda }}_h}{\mu }$ and let

$\begin{array}{c}\hfill {\mathcal{R}}_{01}=\frac{c_1{\beta }_1\left(V_1\right)S_0{e}^{-\mu \tau }}{\mu +\xi +{\delta }_1\left(V_1\right)},{\mathcal{R}}_{02}=\frac{c_2{\beta }_2\left(V_1\right)\xi S_0{e}^{-\mu \tau }}{(\mu +{\delta }_2\left(V_1\right))(\mu +\xi +{\delta }_1\left(V_1\right))}.\end{array}$

Remark 5.

Clearly, ${\mathcal{R}}_{01}$ is the basic reproduction number for the disease driven by the HIV-positive individual, and ${\mathcal{R}}_{02}$ is the basic reproduction number for the disease caused by the AIDS individual. Then we have ${\mathcal{R}}_0^h={\mathcal{R}}_{01}+{\mathcal{R}}_{02}$.

4.1. The Existence and Stability of the Disease-Free Equilibrium

Obviously, Model (20) always has a disease-free equilibrium $E_0=(S_0,0,0)$, and the following result is on the local stability of $E_0$.

Theorem 8.

For any $\tau \ge 0$, $E_0$ is locally asymptotically stable if ${\mathcal{R}}_0^h\le 1$ and unstable if ${\mathcal{R}}_0^h>1$.

Proof. 

The characteristic equation of Model (20) at $E_0$ is

(22)$\begin{array}{c}\hfill (\lambda +\mu )\left[(\lambda -c_1{\beta }_1\left(V_1\right)S_0{e}^{-(\mu +\lambda )\tau }+\mu +\xi +{\delta }_1\left(V_1\right))(\lambda +\mu +{\delta }_2\left(V_1\right))-c_2{\beta }_2\left(V_1\right)\xi S_0{e}^{-(\mu +\lambda )\tau }\right]=0.\end{array}$

Clearly, (22) has a negative real root, ${\lambda }_1=-\mu$, and other roots are determined by

$G\left(\lambda \right):=(\lambda -c_1{\beta }_1\left(V_1\right)S_0{e}^{-(\mu +\lambda )\tau }+\mu +\xi +{\delta }_1\left(V_1\right))(\lambda +\mu +{\delta }_2\left(V_1\right))-c_2{\beta }_2\left(V_1\right)\xi S_0{e}^{-(\mu +\lambda )\tau }=0.$

This yields from the above equation that

(23)$\begin{array}{c}\hfill \left(\frac{\lambda }{\mu +\xi +{\delta }_1\left(V_1\right)}+1\right)\left(\frac{\lambda }{\mu +{\delta }_2\left(V_1\right)}+1\right)=e^{-\lambda \tau }\left(\frac{\lambda }{\mu +{\delta }_2\left(V_1\right)}{\mathcal{R}}_{01}+{\mathcal{R}}_0^h\right).\end{array}$

Now, we claim that all roots of $G\left(\lambda \right)=0$ have negative real parts for ${\mathcal{R}}_0^h\le 1$. If it is invalid, then there exists a root ${\lambda }_2=x_2+i{y}_2$ with $x_2\ge 0$. For this case, it is obvious that

$\begin{array}{c}\hfill \left|\frac{{\lambda }_2}{\mu +\xi +{\delta }_1\left(V_1\right)}+1\right|>\left|e^{-{\lambda }_2\tau }\right|,\left|\frac{{\lambda }_2}{\mu +{\delta }_2\left(V_1\right)}+1\right|>\left|\frac{{\lambda }_2}{\mu +{\delta }_2\left(V_1\right)}{\mathcal{R}}_{01}+{\mathcal{R}}_0^h\right|.\end{array}$

It follows that

$\begin{array}{c}\hfill \left|\left(\frac{{\lambda }_2}{\mu +\xi +{\delta }_1\left(V_1\right)}+1\right)\left(\frac{{\lambda }_2}{\mu +{\delta }_2\left(V_1\right)}+1\right)\right|>\left|e^{-{\lambda }_2\tau }\left(\frac{{\lambda }_2}{\mu +{\delta }_2\left(V_1\right)}{\mathcal{R}}_{01}+{\mathcal{R}}_0^h\right)\right|,\end{array}$

due to the fact that ${\mathcal{R}}_0^h\le 1$. This contradicts to (23). Thus, if ${\mathcal{R}}_0^h\le 1$, all roots of Equation (22) have negative real parts, indicating that $E_0$ is locally asymptotically stable. Conversely, if ${\mathcal{R}}_0^h>1$, we have $G\left(0\right)=(\mu +\xi +{\delta }_1\left(V_1\right))(\mu +{\delta }_2\left(V_1\right))(1-{\mathcal{R}}_0^h)<0$ and $G\left(\lambda \right)\to +\infty$ as $\lambda \to +\infty$. It should be noted that $G\left(\lambda \right)$ is a continuous function with respect to $\lambda$. Consequently, Equation (22) has a positive real root, indicating instability of $E_0$. □

4.2. The Existence of the Positive Equilibrium

Due to the limitations of research methods, it is often challenging to perform a theoretical analysis of the existence of a positive equilibrium when the viral load-dependent transmission rate is subject to a saturated incidence rate, and disease-induced mortality depends on the viral load. Therefore, we consider the scenario where ${\beta }_i\left({\tilde{V}}_1\right)$, $i=1,2$, satisfying the conditions as outlined in Section 2, and the disease-induced mortality ${\delta }_i\left({\tilde{V}}_1\right)$ are replaced by positive constants ${\delta }_i$, $i=1,2$. Let ${\tilde{y}}_1=\frac{c_1{\eta }_1\tilde{I}+c_2{\eta }_2\tilde{A}}{1+\alpha (\tilde{I}+\tilde{A})}$, then the immunity-inactivated infection equilibrium, ${\tilde{U}}_1$, of Model (4) can be written as

$\begin{array}{cc}\hfill {\tilde{V}}_1& =\frac{{\tilde{y}}_1({\mu }_c+{\delta }_c)+p({\mathsf{\Lambda }}_c-{\mu }_c{\tilde{T}}_1)e^{-{\mu }_c{\tau }_1}}{{\mu }_v({\mu }_c+{\delta }_c)},{\tilde{T}}_1^{*}=\frac{({\mathsf{\Lambda }}_c-{\mu }_c{\tilde{T}}_1)e^{-{\mu }_c{\tau }_1}}{{\mu }_c+{\delta }_c},\hfill \\ \hfill {\tilde{T}}_1& =\frac{T_0}{2}\left(1+\frac{1}{{\mathcal{R}}_0^w}+\frac{{\kappa }_1{\tilde{y}}_1}{{\mu }_v{\mu }_c{\mathcal{R}}_0^w}-\sqrt{{\left[1+\frac{1}{{\mathcal{R}}_0^w}+\frac{{\kappa }_1{\tilde{y}}_1}{{\mu }_v{\mu }_c{\mathcal{R}}_0^w}\right]}^2-\frac{4}{{\mathcal{R}}_0^w}}\right).\hfill \end{array}$

Let ${\tilde{y}}_2=\frac{c_1{\beta }_1\left({\tilde{V}}_1\right)\tilde{I}+c_2{\beta }_2\left({\tilde{V}}_1\right)\tilde{A}}{1+\alpha (\tilde{I}+\tilde{A})}$, then the equilibrium, $\tilde{E}$, can be written as

$\begin{array}{c}\hfill \tilde{S}=\frac{{\mathsf{\Lambda }}_h-(\mu +\xi +{\delta }_1)e^{\mu \tau }\tilde{I}}{\mu },\tilde{A}=\frac{{\tilde{y}}_2{\mathsf{\Lambda }}_h\xi e^{-\mu \tau }}{(\mu +{\delta }_2)({\tilde{y}}_2+\mu )(\mu +\xi +{\delta }_1)},\tilde{I}=\frac{{\tilde{y}}_2{\mathsf{\Lambda }}_h{e}^{-\mu \tau }}{({\tilde{y}}_2+\mu )(\mu +\xi +{\delta }_1)}.\end{array}$

Further, ${\tilde{y}}_1$ can be described by ${\tilde{y}}_2$, that is

${\tilde{y}}_1=\frac{(c_1{\eta }_1(\mu +{\delta }_2)+c_2{\eta }_2\xi ){\mathsf{\Lambda }}_h{\tilde{y}}_2{e}^{-\mu \tau }}{({\tilde{y}}_2+\mu )(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)+\alpha {\mathsf{\Lambda }}_h{\tilde{y}}_2(\mu +\xi +{\delta }_2)e^{-\mu \tau }}.$

Thus, ${\tilde{V}}_1={\tilde{V}}_1\left({\tilde{y}}_2\right)$. According to ${\tilde{y}}_2=\frac{c_1{\beta }_1\left({\tilde{V}}_1\right)\tilde{I}+c_2{\beta }_2\left({\tilde{V}}_1\right)\tilde{A}}{1+\alpha (\tilde{I}+\tilde{A})}$, we have

${\tilde{y}}_2=\frac{(c_1{\beta }_1\left({\tilde{V}}_1\right)(\mu +{\delta }_2)+c_2{\beta }_2\left({\tilde{V}}_1\right)\xi ){\mathsf{\Lambda }}_h{\tilde{y}}_2{e}^{-\mu \tau }}{({\tilde{y}}_2+\mu )(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)+\alpha {\mathsf{\Lambda }}_h{\tilde{y}}_2(\mu +\xi +{\delta }_2)e^{-\mu \tau }}.$

Let

(24)$\begin{array}{c}\hfill H\left(y_2\right)=\frac{(c_1{\beta }_1\left({\tilde{V}}_1\left(y_2\right)\right)(\mu +{\delta }_2)+c_2{\beta }_2\left({\tilde{V}}_1\left(y_2\right)\right)\xi ){\mathsf{\Lambda }}_h{e}^{-\mu \tau }}{(y_2+\mu )(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)+\alpha {\mathsf{\Lambda }}_h{y}_2(\mu +\xi +{\delta }_2)e^{-\mu \tau }}-1,\end{array}$

then ${\tilde{y}}_2$ is a solution of $H\left(y_2\right)=0$. Subsequently, our objective is to determine the number of positive real roots of $H\left(y_2\right)=0$ by examining the characteristics of the function $H\left(y_2\right)$. It can be demonstrated that

$H\left(0\right)=\frac{(c_1{\beta }_1\left(V_1\right)(\mu +{\delta }_2)+c_2{\beta }_2\left(V_1\right)\xi ){\mathsf{\Lambda }}_h{e}^{-\mu \tau }}{\mu (\mu +{\delta }_2)(\mu +\xi +{\delta }_1)}-1={\mathcal{R}}_0^h-1\left\{\begin{array}{c}<0,{\mathcal{R}}_0^w>1,{\mathcal{R}}_0^h<1,\hfill \\ >0,{\mathcal{R}}_0^w>1,{\mathcal{R}}_0^h>1,\hfill \\ -1,{\mathcal{R}}_0^w<1,\hfill \end{array}\right.$

and $H(\infty )<0$. The derivative of $H\left(y_2\right)$ is

$\begin{array}{c}\hfill \frac{\mathrm{d}H\left(y_2\right)}{\mathrm{d}{y}_2}=\frac{{\mathsf{\Lambda }}_h{e}^{-\mu \tau }F\left(y_2\right)}{{[(y_2+\mu )(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)+\alpha {\mathsf{\Lambda }}_h{y}_2(\mu +\xi +{\delta }_2)e^{-\mu \tau }]}^2},\end{array}$

where

(25)$\begin{array}{cc}\hfill F\left(y_2\right)=& \left[(\mu +{\delta }_2)c_1\frac{\mathrm{d}{\beta }_1\left({\tilde{V}}_1\right)}{\mathrm{d}{V}_1}+\xi c_2\frac{\mathrm{d}{\beta }_2\left({\tilde{V}}_1\right)}{\mathrm{d}{V}_1}\right]\frac{\mathrm{d}{\tilde{V}}_1\left(y_2\right)}{\mathrm{d}{y}_2}[(y_2+\mu )(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)+\alpha {\mathsf{\Lambda }}_h{y}_2(\mu +\xi +{\delta }_2)e^{-\mu \tau }]\hfill \\ \hfill & -[(\mu +{\delta }_2)c_1{\beta }_1\left({\tilde{V}}_1\right)+\xi c_2{\beta }_2\left({\tilde{V}}_1\right)][(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)+\alpha {\mathsf{\Lambda }}_h(\mu +\xi +{\delta }_2)e^{-\mu \tau }].\hfill \end{array}$

It is clear that the sign of $\frac{\mathrm{d}H\left(y_2\right)}{\mathrm{d}{y}_2}$ can be determined by analyzing the properties of the function $F\left(y_2\right)$. So, we first examine the characteristics of the function $F\left(y_2\right)$. The derivative of $F\left(y_2\right)$ is

$\begin{array}{cc}\hfill \frac{\mathrm{d}F\left(y_2\right)}{\mathrm{d}{y}_2}=& \left\{\left[(\mu +{\delta }_2)c_1\frac{\mathrm{d}{\beta }_1\left({\tilde{V}}_1\right)}{\mathrm{d}{\tilde{V}}_1}+\xi c_2\frac{\mathrm{d}{\beta }_2\left({\tilde{V}}_1\right)}{\mathrm{d}{\tilde{V}}_1}\right]\frac{{\mathrm{d}}^2{\tilde{V}}_1\left(y_2\right)}{\mathrm{d}{y}_2^2}+\left[(\mu +{\delta }_2)c_1\frac{{\mathrm{d}}^2{\beta }_1\left({\tilde{V}}_1\right)}{\mathrm{d}{\tilde{V}}_1^2}+\xi c_2\frac{{\mathrm{d}}^2{\beta }_2\left({\tilde{V}}_1\right)}{\mathrm{d}{\tilde{V}}_1^2}\right]{\left[\frac{\mathrm{d}{\tilde{V}}_1\left(y_2\right)}{\mathrm{d}{y}_2}\right]}^2\right\}\hfill \\ \hfill & \times \left[(y_2+\mu )(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)+\alpha {\mathsf{\Lambda }}_h{y}_2(\mu +\xi +{\delta }_2)e^{-\mu \tau }\right]<0,\hfill \end{array}$

with

$\begin{array}{cc}\hfill \frac{\mathrm{d}{\tilde{V}}_1\left(y_1\right)}{\mathrm{d}{y}_1}& =\frac{1}{2{\mu }_v}\left[\frac{2{\kappa }_2{\mu }_v+{\kappa }_{1}p}{{\kappa }_2{\mu }_v+{\kappa }_{1}p}+\frac{{\kappa }_{1}p{T}_0\left(1+\frac{1}{{\mathcal{R}}_0^w}+\frac{{\kappa }_1{y}_1}{{\mu }_v{\mu }_c{\mathcal{R}}_0^w}\right)}{{\mu }_v({\mu }_c+{\delta }_c){\mathcal{R}}_0^w{e}^{{\mu }_c{\tau }_1}\sqrt{{\left(1+\frac{1}{{\mathcal{R}}_0^w}+\frac{{\kappa }_1{y}_1}{{\mu }_v{\mu }_c{\mathcal{R}}_0^w}\right)}^2-\frac{4}{{\mathcal{R}}_0^w}}}\right]>0,\hfill \\ \hfill \frac{\mathrm{d}{y}_1\left(y_2\right)}{\mathrm{d}{y}_2}& =\frac{{\mathsf{\Lambda }}_h\mu (\mu +{\delta }_2)(\mu +\xi +{\delta }_1)[c_1{\eta }_1(\mu +{\delta }_2)+c_2{\eta }_2\xi ]e^{-\mu \tau }}{{[(y_2+\mu )(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)+\alpha {\mathsf{\Lambda }}_h{y}_2(\mu +\xi +{\delta }_2)e^{-\mu \tau }]}^2}>0,\hfill \end{array}$

$\begin{array}{cc}\hfill \frac{\mathrm{d}{\tilde{V}}_1\left(y_2\right)}{\mathrm{d}{y}_2}& =\frac{\mathrm{d}{\tilde{V}}_1\left(y_1\right)}{\mathrm{d}{y}_1}\frac{\mathrm{d}{y}_1\left(y_2\right)}{\mathrm{d}{y}_2}>0,\frac{{\mathrm{d}}^2{\tilde{V}}_1\left(y_2\right)}{\mathrm{d}{y}_2^2}=\frac{{\mathrm{d}}^2{\tilde{V}}_1\left(y_1\right)}{\mathrm{d}{y}_1^2}{\left[\frac{\mathrm{d}{y}_1\left(y_2\right)}{\mathrm{d}{y}_2}\right]}^2+\frac{\mathrm{d}{\tilde{V}}_1\left(y_1\right)}{\mathrm{d}{y}_1}\frac{{\mathrm{d}}^2{y}_1\left(y_2\right)}{\mathrm{d}{y}_2^2}<0,\hfill \\ \hfill \frac{{\mathrm{d}}^2{\tilde{V}}_1\left(y_1\right)}{\mathrm{d}{y}_1^2}& =\frac{p{\mu }_c}{2{\mu }_v({\mu }_c+{\delta }_c)e^{{\mu }_c{\tau }_1}{T}_0\sqrt{{\left(1+\frac{1}{{\mathcal{R}}_0^w}+\frac{{\kappa }_1{y}_1}{{\mu }_v{\mu }_c{\mathcal{R}}_0^w}\right)}^2-\frac{4}{{\mathcal{R}}_0^w}}}\left[1-\frac{{\left(1+\frac{1}{{\mathcal{R}}_0^w}+\frac{{\kappa }_1{y}_1}{{\mu }_v{\mu }_c{\mathcal{R}}_0^w}\right)}^2}{{\left(1+\frac{1}{{\mathcal{R}}_0^w}+\frac{{\kappa }_1{y}_1}{{\mu }_v{\mu }_c{\mathcal{R}}_0^w}\right)}^2-\frac{4}{{\mathcal{R}}_0^w}}\right]<0,\hfill \\ \hfill \frac{{\mathrm{d}}^2{y}_1\left(y_2\right)}{\mathrm{d}{y}_2^2}& =\frac{-2[(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)e^{\mu \tau }+\alpha {\mathsf{\Lambda }}_h(\mu +\xi +{\delta }_2)]{\mathsf{\Lambda }}_h\mu (\mu +{\delta }_2)(\mu +\xi +{\delta }_1)[c_1{\eta }_1(\mu +{\delta }_2)+c_2{\eta }_2\xi ]}{e^{2\mu \tau }{[(y_2+\mu )(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)+\alpha {\mathsf{\Lambda }}_h{y}_2(\mu +\xi +{\delta }_2)e^{-\mu \tau }]}^3}<0.\hfill \end{array}$

It is evident that $F\left(y_2\right)$ is a monotonically decreasing function, and $F\left(y_2\right)\to -\infty$ as $y_2\to +\infty$. Consequently, if $F\left(0\right)>0$, as $y_2$ increases, the function $H\left(y_2\right)$ initially increases and then decreases. In this case, if $H\left(0\right)<0$ and $H_{max}>0$ hold, the equation $H\left(y_2\right)=0$ will have two distinct positive roots. If $H\left(0\right)<0$ and $H_{max}=0$ or only $H\left(0\right)>0$ hold, there will be a unique positive root. Additionally, if $F\left(0\right)<0$, the function $H\left(y_2\right)$ will always decrease. In such a scenario, if $H\left(0\right)<0$, there will be no positive roots for the equation $H\left(y_2\right)=0$. However, if $H\left(0\right)>0$, there will be only one positive root. Based on the aforementioned analysis, we can establish the following result regarding the existence of positive equilibria for Model (20).

Theorem 9.

Let $H=H\left(y_2\right)$, $F=F\left(y_2\right)$ are defined in (24) and (25), and  $H_{max}:={max}_{y_2\in [0,+\infty )}H\left(y_2\right)$,

• $\left(i\right)$

  for ${\mathcal{R}}_0^w>1$ and ${\mathcal{R}}_0^h<1$,

  • $\left(a\right)$

    if $F\left(0\right)>0$ and $H_{max}>0$, then two positive equilibria, ${\tilde{E}}_1$ and ${\tilde{E}}_2$, both exist;

  • $\left(b\right)$

    if $F\left(0\right)>0$ and $H_{max}=0$, then a unique positive equilibrium, $\tilde{E}$, exists;

  • $\left(c\right)$

    if $F\left(0\right)\le 0$ or $H_{max}<0$, then a positive equilibrium does not exist;

• $\left(ii\right)$

  for ${\mathcal{R}}_0^w>1$ and ${\mathcal{R}}_0^h>1$, a unique positive equilibrium, $\tilde{E}$, exists;

• $\left(iii\right)$

  for ${\mathcal{R}}_0^w<1$,

  • $\left(a\right)$

    if $F\left(0\right)>0$ and $H_{max}>0$, then two positive equilibria, ${\tilde{E}}_1$ and ${\tilde{E}}_2$, both exist;

  • $\left(b\right)$

    if $F\left(0\right)>0$ and $H_{max}=0$, then a unique positive equilibrium, $\tilde{E}$, exists;

  • $\left(c\right)$

    if $F\left(0\right)\le 0$ or $H_{max}<0$, then a positive equilibrium does not exist.

4.3. The Local Stability of the Positive Equilibrium

In this section, we are concerned with the local stability of equilibria provided in Theorem 9. Since it is hard to verify the stability evaluated at a positive equilibrium in Model (20) due to the complexity of the functions ${\beta }_i\left(V\right)$ and ${\delta }_i\left(V\right)$$(i=1,2)$, we therefore suppose that ${\beta }_i\left(V\right)={\overline{a}}_{i}V$, ${\delta }_i\left(V\right)={\delta }_i$, $i=1,2$, where ${\overline{a}}_i$ and ${\delta }_i$ are positive constants. The same assumptions are adopted as in [33].

Theorem 10.

Suppose that Model (20) admits that a unique positive equilibrium, $\tilde{E}$, exists, which is locally asymptotically stable if $\frac{\alpha {\tilde{V}}_1}{1+\alpha (\tilde{I}+\tilde{A})}>max\{{\tilde{V}}_{1_I}^{\prime },{\tilde{V}}_{1_A}^{\prime }\}$, where ${\tilde{V}}_{1_I}^{\prime }$ and ${\tilde{V}}_{1_A}^{\prime }$ are given below.

Proof. 

The characteristic equation of Model (20) at $\tilde{E}$ is

(26)$\begin{array}{cc}\hfill & (\lambda +m+\mu )(\lambda +\mu +\xi +{\delta }_1)(\lambda +\mu +{\delta }_2)-(\mu +\xi +{\delta }_1)(\lambda +\mu )e^{-\lambda \tau }{\left[\left(c_1{\overline{a}}_1{\tilde{V}}_1\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_1\tilde{A}\right)\left(1+\alpha (\tilde{I}+\tilde{A})\right)\right]}^{-1}\hfill \\ \hfill & \times \left\{(\lambda +\mu +{\delta }_2)\tilde{I}\left[c_1{\overline{a}}_1\left({\tilde{V}}_{1_I}^{\prime }\tilde{I}(1+\alpha (\tilde{I}+\tilde{A}))+{\tilde{V}}_1(1+\alpha \tilde{A})\right)+c_2{\overline{a}}_2\left({\tilde{V}}_{1_I}^{\prime }\tilde{A}(1+\alpha (\tilde{I}+\tilde{A}))-\alpha \tilde{A}{\tilde{V}}_1\right)\right]\right.\hfill \\ \hfill & \left.+\xi \tilde{I}\left[c_1{\overline{a}}_1\left({\tilde{V}}_{1_A}^{\prime }\tilde{I}(1+\alpha (\tilde{I}+\tilde{A}))-\alpha \tilde{I}{\tilde{V}}_1\right)+c_2{\overline{a}}_2\left({\tilde{V}}_{1_A}^{\prime }\tilde{A}(1+\alpha (\tilde{I}+\tilde{A}))+{\tilde{V}}_1(1+\alpha \tilde{I})\right)\right]\right\}=0,\hfill \end{array}$

where

$\begin{array}{cc}\hfill m& =\frac{c_1{\overline{a}}_1{\tilde{V}}_1\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_1\tilde{A}}{1+\alpha (\tilde{I}+\tilde{A})},{\tilde{V}}_{1_I}^{\prime }=\frac{{\rm Y}\left(c_1{\eta }_1(1+\alpha \tilde{A})-c_2{\eta }_2\alpha \tilde{A}\right)}{{(1+\alpha (\tilde{I}+\tilde{A}))}^2},{\tilde{V}}_{1_A}^{\prime }=\frac{{\rm Y}\left(c_2{\eta }_2(1+\alpha \tilde{I})-c_1{\eta }_1\alpha \tilde{I}\right)}{{(1+\alpha (\tilde{I}+\tilde{A}))}^2},\hfill \\ \hfill {\rm Y}& =\frac{1}{{\mu }_v}\left[1-\frac{{\kappa }_{1}p{T}_0}{2{\mu }_v({\mu }_c+{\delta }_c){\mathcal{R}}_0^w{e}^{{\mu }_c{\tau }_1}}\left(1-\frac{1+\frac{1}{{\mathcal{R}}_0^w}+\frac{{\kappa }_1{\tilde{y}}_1}{{\mu }_v{\mu }_c{\mathcal{R}}_0^w}}{\sqrt{{\left(1+\frac{1}{{\mathcal{R}}_0^w}+\frac{{\kappa }_1{\tilde{y}}_1}{{\mu }_v{\mu }_c{\mathcal{R}}_0^w}\right)}^2-\frac{4}{{\mathcal{R}}_0^w}}}\right)\right].\hfill \end{array}$

Due to the fact that $\xi \tilde{I}=(\mu +{\delta }_2)\tilde{A}$, Equation (26) becomes

(27)$\begin{array}{cc}\hfill & (\mu +\xi +{\delta }_1)(\lambda +\mu )e^{-\lambda \tau }\left[\frac{\alpha (c_1{\overline{a}}_1{\tilde{V}}_1\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_1\tilde{A})}{1+\alpha (\tilde{I}+\tilde{A})}-\left(c_1{\overline{a}}_1{\tilde{V}}_{1_I}^{\prime }\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_{1_I}^{\prime }\tilde{A}\right)\right]\frac{(\lambda +\mu +{\delta }_2)\tilde{I}}{c_1{\overline{a}}_1{\tilde{V}}_1\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_1\tilde{A}}\hfill \\ \hfill & +(\mu +\xi +{\delta }_1)(\lambda +\mu )e^{-\lambda \tau }\left[\frac{\alpha (c_1{\overline{a}}_1{\tilde{V}}_1\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_1\tilde{A})}{1+\alpha (\tilde{I}+\tilde{A})}-\left(c_1{\overline{a}}_1{\tilde{V}}_{1_A}^{\prime }\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_{1_A}^{\prime }\tilde{A}\right)\right]\frac{\xi \tilde{I}}{c_1{\overline{a}}_1{\tilde{V}}_1\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_1\tilde{A}}\hfill \\ \hfill & +(\lambda +\mu )(\mu +\xi +{\delta }_1)e^{-\lambda \tau }\frac{c_2{\overline{a}}_2{\tilde{V}}_1\lambda \tilde{A}}{c_1{\overline{a}}_1{\tilde{V}}_1\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_1\tilde{A}}+(\lambda +m+\mu )(\lambda +\mu +\xi +{\delta }_1)(\lambda +\mu +{\delta }_2)\hfill \\ \hfill & -(\lambda +\mu )(\mu +\xi +{\delta }_1)(\lambda +\mu +{\delta }_2)e^{-\lambda \tau }=0.\hfill \end{array}$

Now, we claim that all roots of (27) have negative real parts. If not, there is a root ${\lambda }_3=x_3+i{y}_3$ with $x_3\ge 0$. In this situation, it is obvious that

$\begin{array}{c}\hfill \left|{\lambda }_3+m+\mu \right|>|{\lambda }_3+\mu |,|{\lambda }_3+\mu +\xi +{\delta }_1|\ge |\mu +\xi +{\delta }_1|,1\ge \left|e^{-{\lambda }_3\tau }\right|.\end{array}$

If $\frac{\alpha {\tilde{V}}_1}{1+\alpha (\tilde{I}+\tilde{A})}>max\{{\tilde{V}}_{1_I}^{\prime },{\tilde{V}}_{1_A}^{\prime }\}$, it follows that

$\begin{array}{cc}\hfill & \left|(\mu +\xi +{\delta }_1)(\lambda +\mu )e^{-\lambda \tau }\left[\frac{\alpha (c_1{\overline{a}}_1{\tilde{V}}_1\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_1\tilde{A})}{1+\alpha (\tilde{I}+\tilde{A})}-\left(c_1{\overline{a}}_1{\tilde{V}}_{1_I}^{\prime }\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_{1_I}^{\prime }\tilde{A}\right)\right]\frac{(\lambda +\mu +{\delta }_2)\tilde{I}}{c_1{\overline{a}}_1{\tilde{V}}_1\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_1\tilde{A}}\right.\hfill \\ \hfill & +(\mu +\xi +{\delta }_1)(\lambda +\mu )e^{-\lambda \tau }\left[\frac{\alpha (c_1{\overline{a}}_1{\tilde{V}}_1\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_1\tilde{A})}{1+\alpha (\tilde{I}+\tilde{A})}-\left(c_1{\overline{a}}_1{\tilde{V}}_{1_A}^{\prime }\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_{1_A}^{\prime }\tilde{A}\right)\right]\frac{\xi \tilde{I}}{c_1{\overline{a}}_1{\tilde{V}}_1\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_1\tilde{A}}\hfill \\ \hfill & \left.+(\lambda +m+\mu )(\lambda +\mu +\xi +{\delta }_1)(\lambda +\mu +{\delta }_2)+(\mu +\xi +{\delta }_1)(\lambda +\mu )e^{-\lambda \tau }\frac{c_2{\overline{a}}_2{\tilde{V}}_1\lambda \tilde{A}}{c_1{\overline{a}}_1{\tilde{V}}_1\tilde{I}+c_2{\overline{a}}_2{\tilde{V}}_1\tilde{A}}\right|\hfill \\ \hfill >& \left|(\lambda +\mu )(\mu +\xi +{\delta }_1)(\lambda +\mu +{\delta }_2)e^{-\lambda \tau }\right|,\hfill \end{array}$

which contradicts to (27). Consequently, given that all roots of Equation (27) have negative real parts, it can be deduced that $\tilde{E}$ is locally asymptotically stable. □

Remark 6.

Theorem 9 indicates that Model (20) may have no, one, or two positive equilibria under different conditions, respectively. However, it is difficult to check analytically the stability of the coexistence of multiple positive equilibria, so some numerical simulations are conducted to confirm the corresponding stability for different cases. And we find that Model (20) may exhibit complex dynamic behavior, including backward bifurcation and Hopf bifurcation.

Remark 7.

It should be pointed out that if the disease-induced mortality, ${\delta }_i\left({\tilde{V}}_1\right)$, varies based on the viral load, rather than being replaced by positive constants, then the existence and stability of positive equilibria in Model (20) become more complicated, implying more complex dynamic behavior.

5. Dynamics of Coupled Slow Time Model with Immune Response

In this section, it is assumed that the immune response plays a role in the host and ${\tau }_2$, the delay of the immune response, in its stable interval $[0,{\tau }_2^{*})$. When $I,A>0$, the dynamics of the virus in the host satisfies $\frac{\mathrm{d}V\left(s\right)}{\mathrm{d}s}=\mathsf{\Theta }(I,A)+p{T}^{*}\left(s\right)-qB\left(s\right)V\left(s\right)-{\mu }_{v}V\left(s\right)$. And we always assume that ${\mathcal{R}}_0^w>1$ and ${\mathcal{R}}_w>1$ in the analysis below. Furthermore, considering that the viral infection process within the host occurs at a much faster rate compared with the disease transmission process between hosts, and based on the stability of the equilibrium, ${\tilde{U}}_2$, derived from Theorem 4, it can be assumed that the state $(T\left(s\right),T^{*}\left(s\right),V\left(s\right),B\left(s\right))$ of the fast time model has reached its stable equilibrium $({\tilde{T}}_2,{\tilde{T}}_2^{*},{\tilde{V}}_2,{\tilde{B}}_2)$. Therefore, the viral load, $V={\tilde{V}}_2$, in vivo is a positive constant, and Model (3) can be written as

(28)$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}=& {\mathsf{\Lambda }}_h-\frac{c_1{\beta }_1\left({\tilde{V}}_2\right)I\left(t\right)+c_2{\beta }_2\left({\tilde{V}}_2\right)A\left(t\right)}{1+\alpha \left(I\right(t)+A(t\left)\right)}S\left(t\right)-\mu S\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=& \frac{c_1{\beta }_1\left({\tilde{V}}_2\right)I(t-\tau )+c_2{\beta }_2\left({\tilde{V}}_2\right)A(t-\tau )}{1+\alpha \left(I\right(t-\tau )+A(t-\tau \left)\right)}S(t-\tau )e^{-\mu \tau }\hfill \\ \hfill & -(\mu +\xi +{\delta }_1\left({\tilde{V}}_2\right))I\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}=& \xi I\left(t\right)-(\mu +{\delta }_2\left({\tilde{V}}_2\right))A\left(t\right),\hfill \end{array}\right.$

with the initial condition

(29)$\begin{array}{cc}\hfill & S\left(\varrho \right)={\varphi }_1\left(\varrho \right)>0,I\left(\varrho \right)={\varphi }_3\left(\varrho \right)\ge 0,A\left(\varrho \right)={\varphi }_4\left(\varrho \right)\ge 0,\hfill \\ \hfill & \varrho \in [-\tau ,0],{\varphi }_1\left(0\right)>0,{\varphi }_i\left(0\right)\ge 0,i=3,4,\hfill \end{array}$

where $\varphi =({\varphi }_1,{\varphi }_3,{\varphi }_4)\in C([-\tau ,0],{\mathbb{R}}_{+}^3)$. It is obvious that Model (28) has a unique solution $\left(S\right(t),I(t),A(t\left)\right)$ satisfying the initial condition (29). Regarding the nonnegativity and boundedness of solutions for Model (28), we can draw a result similar to that of Model (20).

Theorem 11.

Under the initial condition (29), any solution $\left(S\right(t),I(t),A(t\left)\right)$ of Model (28) is nonnegative and ultimately bounded for all $t\ge 0$. Particularly, we also have ${lim}_{t\to \infty }(S(t-\tau )+I\left(t\right)e^{\mu \tau }+A\left(t\right)e^{\mu \tau })\le \frac{{\mathsf{\Lambda }}_h}{\mu }$, for all $t\ge 0$.

For the sake of brevity, we denote ${\beta }_1\left({\tilde{V}}_2\right)={\beta }_1$, ${\beta }_2\left({\tilde{V}}_2\right)={\beta }_2$, ${\delta }_1\left({\tilde{V}}_2\right)={\delta }_1$ and ${\delta }_2\left({\tilde{V}}_2\right)={\delta }_2$, and the between-host reproduction number for Model (28) is defined by

${\tilde{\mathcal{R}}}_0^h=\frac{(c_1{\beta }_1(\mu +{\delta }_2)+c_2{\beta }_2\xi )S_0{e}^{-\mu \tau }}{(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)}.$

5.1. The Existence and Stability of the Equilibria

Compared with Model (20), it is easy to show that the existence of the various equilibria of Model (28) is straightforward, and we can summarize the results as follows.

Theorem 12.

If ${\tilde{\mathcal{R}}}_0^h\le 1$, Model (28) only has a disease-free equilibrium, ${\tilde{E}}_0=(S_0,0,0)$; if ${\tilde{\mathcal{R}}}_0^h>1$, except for ${\tilde{E}}_0$, Model (28) also has a unique endemic equilibrium, $E^{*}=(S^{*},I^{*},A^{*})$, where

$\begin{array}{cc}\hfill S^{*}& =\frac{{\mathsf{\Lambda }}_h\alpha (\mu +\xi +{\delta }_2)+(\mu +\xi +{\delta }_1)(\mu +{\delta }_2)e^{\mu \tau }}{\mu \alpha (\mu +\xi +{\delta }_2)+c_1{\beta }_1(\mu +{\delta }_2)+c_2{\beta }_2\xi },I^{*}=\frac{\mu (\mu +{\delta }_2)({\tilde{\mathcal{R}}}_0^h-1)}{\mu \alpha (\mu +\xi +{\delta }_2)+c_1{\beta }_1(\mu +{\delta }_2)+c_2{\beta }_2\xi },\hfill \\ \hfill A^{*}& =\frac{\mu \xi ({\tilde{\mathcal{R}}}_0^h-1)}{\mu \alpha (\mu +\xi +{\delta }_2)+c_1{\beta }_1(\mu +{\delta }_2)+c_2{\beta }_2\xi }.\hfill \end{array}$

Since the proof of local asymptotic stability of ${\tilde{E}}_0$ is identical to Theorem 8, we omit it here. Then for the global stability of ${\tilde{E}}_0$, we have the following result.

Theorem 13.

For any $\tau \ge 0$, the disease-free equilibrium, ${\tilde{E}}_0$, of Model (28) is globally asymptotically stable if ${\mathcal{R}}_0^h\le 1$.

Proof. 

Define

$\begin{array}{c}\hfill {\tilde{L}}_0\left(t\right)=f_{1}I\left(t\right)+f_{2}A\left(t\right)+f_3{\int }_{t-\tau }^t\frac{c_1{\beta }_{1}I\left(\overline{s}\right)+c_2{\beta }_{2}A\left(\overline{s}\right)}{1+\alpha (I\left(\overline{s}\right)+A\left(\overline{s}\right))}S\left(\overline{s}\right)\mathrm{d}\overline{s},\end{array}$

where $f_1=\frac{(c_1{\beta }_1(\mu +{\delta }_2)+c_2{\beta }_2\xi )e^{\mu \tau }}{(\mu +\xi +{\delta }_1)(\mu +{\delta }_2)}$, $f_2=\frac{c_2{\beta }_2{e}^{\mu \tau }}{\mu +{\delta }_2}$, $f_3=\frac{c_1{\beta }_1(\mu +{\delta }_2)+c_2{\beta }_2\xi }{(\mu +\xi +{\delta }_1)(\mu +{\delta }_2)}$. Calculating the time derivative of ${\tilde{L}}_0\left(t\right)$ along the positive solutions of Model (28) yields that

$\begin{array}{c}\hfill \frac{\mathrm{d}{\tilde{L}}_0\left(t\right)}{\mathrm{d}t}=\frac{e^{\mu \tau }(c_1{\beta }_{1}I\left(t\right)+c_2{\beta }_{2}A\left(t\right))}{S_0}\left({\tilde{\mathcal{R}}}_0^h\frac{S\left(t\right)}{1+\alpha \left(I\right(t)+A(t\left)\right)}-S_0\right)\le e^{\mu \tau }(c_1{\beta }_{1}I\left(t\right)+c_2{\beta }_{2}A\left(t\right))({\tilde{\mathcal{R}}}_0^h-1).\end{array}$

Therefore, under the assumption ${\tilde{\mathcal{R}}}_0^h\le 1$, we have $\frac{\mathrm{d}{\tilde{L}}_0\left(t\right)}{\mathrm{d}t}\le 0$. Additionally, it can be easily verified that $S\left(t\right)=S_0$, $I\left(t\right)=0$, and $A\left(t\right)=0$ hold when $\frac{\mathrm{d}{\tilde{L}}_0\left(t\right)}{\mathrm{d}t}=0$. Hence, the singleton $\left\{{\tilde{E}}_0\right\}$ is the largest compact invariant set in $\{(S\left(t\right),I\left(t\right),A\left(t\right))\in {\mathbb{R}}_{+}^3:\frac{\mathrm{d}{\tilde{L}}_0\left(t\right)}{\mathrm{d}t}=0\}$. As a result, by employing LaSalle’s invariance principle, we can conclude that ${\tilde{E}}_0$ is globally asymptotically stable when ${\tilde{\mathcal{R}}}_0^h\le 1$. This completes the proof. □

Further, on the local asymptotic stability of $E^{*}$, the following result is established.

Theorem 14.

For any $\tau \ge 0$, the endemic equilibrium, $E^{*}$, is locally asymptotically stable if ${\tilde{\mathcal{R}}}_0^h>1$.

Proof. 

The characteristic equation of Model (28) at $E^{*}$ is

(30)$\begin{array}{cc}\hfill & \left(\lambda +\frac{c_1{\beta }_1{I}^{*}+c_2{\beta }_2{A}^{*}}{1+\alpha (I^{*}+A^{*})}+\mu \right)(\lambda +\mu +\xi +{\delta }_1)(\lambda +\mu +{\delta }_2)-(\lambda +\mu )(\mu +\xi +{\delta }_1)e^{-\lambda \tau }\hfill \\ \hfill & \times \frac{[c_2{\beta }_2(1+\alpha I^{*})-\alpha c_1{\beta }_1{I}^{*}]\xi I^{*}+[c_1{\beta }_1(1+\alpha A^{*})-\alpha c_2{\beta }_2{A}^{*}](\lambda +\mu +{\delta }_2)I^{*}}{(1+\alpha (I^{*}+A^{*}))(c_1{\beta }_1{I}^{*}+c_2{\beta }_2{A}^{*})}=0.\hfill \end{array}$

Note that $\xi I^{*}=(\mu +{\delta }_2)A^{*}$, Equation (30) becomes

(31)$\begin{array}{cc}\hfill & \left(\lambda +\frac{c_1{\beta }_1{I}^{*}+c_2{\beta }_2{A}^{*}}{1+\alpha (I^{*}+A^{*})}+\mu \right)(\lambda +\mu +\xi +{\delta }_1)(\lambda +\mu +{\delta }_2)-(\lambda +\mu )(\mu +\xi +{\delta }_1)(\lambda +\mu +{\delta }_2)e^{-\lambda \tau }\hfill \\ \hfill & +(\lambda +\mu )(\mu +\xi +{\delta }_1)e^{-\lambda \tau }\left(\frac{\alpha \xi I^{*}+\alpha (\lambda +\mu +{\delta }_2)I^{*}}{(1+\alpha (I^{*}+A^{*}))(c_1{\beta }_1{I}^{*}+c_2{\beta }_2{A}^{*})}+\frac{c_2{\beta }_2\lambda A^{*}}{c_1{\beta }_1{I}^{*}+c_2{\beta }_2{A}^{*}}\right)=0.\hfill \end{array}$

Now, we make the claim that all roots of Equation (31) have negative real parts. If this is not the case, then there exists a root ${\lambda }_4=x_4+i{y}_4$ with $x_4\ge 0$. In this situation, if ${\tilde{\mathcal{R}}}_0^h>1$, it is obvious that

$\begin{array}{c}\hfill \left|{\lambda }_4+\frac{c_1{\beta }_1{I}^{*}+c_2{\beta }_2{A}^{*}}{1+\alpha (I^{*}+A^{*})}+\mu \right|>|{\lambda }_4+\mu |,|{\lambda }_4+\mu +\xi +{\delta }_1|\ge |\mu +\xi +{\delta }_1|,1\ge \left|e^{-{\lambda }_4\tau }\right|.\end{array}$

This yields from the above inequalities that

$\begin{array}{cc}\hfill & \left|({\lambda }_4+\mu )(\mu +\xi +{\delta }_1)e^{-{\lambda }_4\tau }\left(\frac{\alpha \xi I^{*}+\alpha ({\lambda }_4+\mu +{\delta }_2)I^{*}}{(1+\alpha (I^{*}+A^{*}))(c_1{\beta }_1{I}^{*}+c_2{\beta }_2{A}^{*})}+\frac{c_2{\beta }_2{\lambda }_4{A}^{*}}{c_1{\beta }_1{I}^{*}+c_2{\beta }_2{A}^{*}}\right)\right.\hfill \\ \hfill & \left.+\left({\lambda }_4+\frac{c_1{\beta }_1{I}^{*}+c_2{\beta }_2{A}^{*}}{1+\alpha (I^{*}+A^{*})}+\mu \right)({\lambda }_4+\mu +\xi +{\delta }_1)({\lambda }_4+\mu +{\delta }_2)\right|\hfill \\ \hfill >& \left|({\lambda }_4+\mu )(\mu +\xi +{\delta }_1)({\lambda }_4+\mu +{\delta }_2)e^{-{\lambda }_4\tau }\right|.\hfill \end{array}$

This contradicts Equation (31). Therefore, we can conclude that all roots of Equation (30) have negative real parts, indicating that $E^{*}$ is locally asymptotically stable. The proof is completed. □

Remark 8.

In this section, almost complete conclusions about the dynamics of the coupled slow time model (28) have been established. Specifically, ${\tilde{E}}_0$ is globally asymptotically stable if ${\tilde{\mathcal{R}}}_0^h<1$, and, if ${\tilde{\mathcal{R}}}_0^h>1$, the unique endemic equilibrium, $E^{*}$, is locally asymptotically stable. In fact, $E^{*}$ is globally asymptotically stable if ${\tilde{\mathcal{R}}}_0^h>1$, which will be verified in numerical simulations. However, the construction of a suitable Lyapunov functional to prove the global asymptotic stability of $E^{*}$ remains an interesting open problem.

5.2. Uniform Persistence

We are now proceeding to analyze the uniform persistence of Model (28). To do so, let us denote $S(t,\varrho )$, $I(t,\varrho )$ and $A(t,\varrho )$ as the solutions of Model (28) with an initial function $\varrho =({\varrho }_1,{\varrho }_2,{\varrho }_3)\in C([-\tau ,0],{\mathbb{R}}_{+}^3)$. For any $t>0$, we define $(S_t\left(\varrho \right),I_t\left(\varrho \right),A_t\left(\varrho \right))=(S(t+\theta ,\varrho ),I(t+\theta ,\varrho ),A(t+\theta ,\varrho ))$, $\theta \in [-\tau ,0]$. Furthermore, let $X:=\left\{({\varrho }_1,{\varrho }_2,{\varrho }_3)\in C([-\tau ,0],{\mathbb{R}}_{+}^3):{\varrho }_2\ne 0,{\varrho }_3\ne 0\right\}$, $\partial X:=C([-\tau ,0],{\mathbb{R}}_{+}^3)/X=\{({\varrho }_1,{\varrho }_2,{\varrho }_3)\in C([-\tau ,0],{\mathbb{R}}_{+}^3):{\varrho }_2\equiv 0,\mathrm{or}{\varrho }_3\equiv 0\}$. Then, we can draw the following conclusion.

Theorem 15.

Assume that ${\tilde{\mathcal{R}}}_0^h>1$, then there exists a constant $\epsilon >0$ such that solution $S(t,\varrho )$, $I(t,\varrho )$, $A(t,\varrho )$ of Model (28) with initial function $\varrho =({\varrho }_1,{\varrho }_2,{\varrho }_3)\in X$ satisfies ${lim\; inf}_{t\to \infty }S(t,\varrho )\ge \epsilon$, ${lim\; inf}_{t\to \infty }I(t,\varrho )\ge \epsilon$, ${lim\; inf}_{t\to \infty }A(t,\varrho )\ge \epsilon$.

Proof. 

Considering the solution $u(t,\varrho )=\left(S\right(t,\varrho ),I(t,\varrho ),A(t,\varrho \left)\right)$ of Model (28) with the initial condition (29), we can derive the following expression from the first equation of Model (28)

$\begin{array}{c}\hfill \frac{\mathrm{d}S(t,\varrho )}{\mathrm{d}t}={\mathsf{\Lambda }}_h-({\beta }_{1}I(t,\varrho )+{\beta }_{2}A(t,\varrho )+\mu )S(t,\varrho )\ge {\mathsf{\Lambda }}_h-\left(({\beta }_1+{\beta }_2)\frac{{\mathsf{\Lambda }}_h}{\mu }+\mu \right)S(t,\varrho ).\end{array}$

According to the comparison principle, we have ${lim\; inf}_{t\to \infty }S(t,\varrho )\ge \frac{\mu {\mathsf{\Lambda }}_h}{({\beta }_1+{\beta }_2){\mathsf{\Lambda }}_h+{\mu }^2}$, which implies that $S\left(t\right)$ exhibits uniform persistence. It is evident that the set X serves as an invariant set for Model (28). Define $W_{\partial }=\{\varrho \in C([-\tau ,0],{\mathbb{R}}_{+}^3):(S_t\left(\varrho \right),I_t\left(\varrho \right),A_t\left(\varrho \right))\in \partial X,\forall t\ge 0\}$. Let $W_0=\left\{{\tilde{E}}_0\right\}$. It is obvious that $W_0\subset {\cup }_{{\varrho }_0\in W_{\partial }}\omega \left({\varrho }_0\right)$, where $\omega \left({\varrho }_0\right)$ represents the $\omega$-limit set of solution $(S_t\left({\varrho }_0\right),I_t\left({\varrho }_0\right),A_t\left({\varrho }_0\right))$ with the initial value $u\left(0\right)={\varrho }_0$. Given any $\varrho \in W_{\partial }$, it is observed that $(S_t\left(\varrho \right),I_t\left(\varrho \right),A_t\left(\varrho \right))\in W_{\partial }$ for all $t\ge 0$. Consequently, $I(t,\varrho )\equiv 0$ or $A(t,\varrho )\equiv 0$.

If $I(t,\varrho )\equiv 0$, then, we have $A(t,\varrho )\equiv 0$, by means of the second equation of Model (28). Therefore, Model (28) is reduced to

(32)$\begin{array}{c}\hfill \frac{\mathrm{d}S(t,\varrho )}{\mathrm{d}t}={\mathsf{\Lambda }}_h-\mu S(t,\varrho ).\end{array}$

From (32), we can obtain ${lim}_{t\to \infty }S(t,\varrho )=S_0$. Thus, $\omega \left({\varrho }_0\right)=\left\{{\tilde{E}}_0\right\}$.

If $A(t,\varrho )\equiv 0$, then by examining the third equation of Model (28), we can derive that $I(t,\varrho )\equiv 0$. In a similar manner, we can also obtain (32), thus ${lim}_{t\to \infty }S(t,\varrho )=S_0$, which further suggests that $\omega \left({\varrho }_0\right)=\left\{{\tilde{E}}_0\right\}$. From the foregoing discussion, we can easily see that $W_0\supset {\cup }_{{\varrho }_0\in W_{\partial }}\omega \left({\varrho }_0\right)$. Consequently, $W_0={\cup }_{{\varrho }_0\in W_{\partial }}\omega \left({\varrho }_0\right)$.

Next, we prove that $K^s\left({\tilde{E}}_0\right)\cap X=\varnothing$, where $K^s\left({\tilde{E}}_0\right)$ is the stable set of ${\tilde{E}}_0$. By contradiction, assume that there exists a $\varrho \in X$ such that ${lim}_{t\to \infty }u(t,\varrho )={\tilde{E}}_0$. Since ${\tilde{\mathcal{R}}}_0^h>1$, there is a sufficiently small $\epsilon >0$ such that

$\begin{array}{c}\hfill \frac{(c_1{\beta }_1(\mu +{\delta }_2)+c_2{\beta }_2\xi )(S_0-\epsilon )e^{-\mu \tau }}{(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)(1+2\alpha \epsilon )}-1>0.\end{array}$

Hence, with this chosen value of $\epsilon >0$, there exists a $t^{*}>0$ such that $S(t,\varrho )\ge S_0-\epsilon$, $I(t,\varrho )<\epsilon$, and $A(t,\varrho )<\epsilon$ hold for all $t\ge t^{*}$. Define

$\begin{array}{c}\hfill L(t,\varrho )=f_{1}I(t,\varrho )+f_{2}A(t,\varrho )+f_3{\int }_{t-\tau }^t\frac{c_1{\beta }_{1}I\left(\overline{s}\right)+c_2{\beta }_{2}A\left(\overline{s}\right)}{1+\alpha (I\left(\overline{s}\right)+A\left(\overline{s}\right))}S\left(\overline{s}\right)\mathrm{d}\overline{s}.\end{array}$

Then ${lim}_{t\to \infty }L(t,\varrho )=0$. When $t\ge t^{*}$,

$\begin{array}{cc}\hfill \frac{\mathrm{d}L(t,\varrho )}{\mathrm{d}t}=& \frac{(c_1{\beta }_1(\mu +{\delta }_2)+c_2{\beta }_2\xi )(c_1{\beta }_{1}I(t,\varrho )+c_2{\beta }_{2}A(t,\varrho ))S(t,\varrho )}{(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)(1+\alpha (I(t,\varrho )+A(t,\varrho )))}-(c_1{\beta }_{1}I(t,\varrho )+c_2{\beta }_{2}A(t,\varrho ))e^{\mu \tau },\hfill \\ \hfill \ge & \left(\frac{(c_1{\beta }_1(\mu +{\delta }_2)+c_2{\beta }_2\xi )(S_0-\epsilon )e^{-\mu \tau }}{(\mu +{\delta }_2)(\mu +\xi +{\delta }_1)(1+2\alpha \epsilon )}-1\right)(c_1{\beta }_{1}I(t,\varrho )+c_2{\beta }_{2}A(t,\varrho ))e^{\mu \tau }.\hfill \end{array}$

Obviously, $\frac{\mathrm{d}L(t,\varrho )}{\mathrm{d}t}>0$ for all $t\ge t^{*}$, which implies that $L(t,\varrho )$ is an increasing function on $t\ge t^{*}$, and hence ${lim}_{t\to \infty }L(t,\varrho )\ne 0$. This is a contradiction. Thus, we obtain $K^s\left({\tilde{E}}_0\right)\cap X=\varnothing$. According to the theory of persistence for dynamic systems, it is clear that for each $\varrho \in X$ there exists a positive constant, $\epsilon$, such that ${lim\; inf}_{t\to \infty }S(t,\varrho )\ge \epsilon ,{lim\; inf}_{t\to \infty }I(t,\varrho )\ge \epsilon ,{lim\; inf}_{t\to \infty }A(t,\varrho )\ge \epsilon$. Consequently, the positive solutions for Model (28) exhibit uniform persistence. The proof is finished. □

6. Dynamics of Coupled Slow Time Model with Immune Response and Stable Periodic Solution

In this section, we investigate the slow time model with an immune response in the host, where the delay, ${\tau }_2$, does not play a role in its stable interval until a stable periodic solution guaranteed by $\breve{\beta }<0$ occurs in Model (4). In this case, the state $(T\left(s\right),T^{*}\left(s\right),V\left(s\right),B\left(s\right))$ of the fast time model has reached its stable periodic solution. Therefore, ${\beta }_1\left(V\left(s\right)\right)={\overline{a}}_{1}V\left(s\right)={\overline{a}}_{1}V\left(ϵt\right):={\beta }_1\left(t\right)$ is the periodic function with period $\omega$. Similarly, ${\beta }_2\left(V\left(s\right)\right):={\beta }_2\left(t\right)$, ${\delta }_1\left(V\left(s\right)\right):={\delta }_1\left(t\right)$, ${\delta }_2\left(V\left(s\right)\right):={\delta }_2\left(t\right)$ are also periodic functions. And thus, Model (1) can be written as follows, a periodic system,

(33)$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}=& {\mathsf{\Lambda }}_h-\frac{c_1{\beta }_1\left(t\right)I\left(t\right)+c_2{\beta }_2\left(t\right)A\left(t\right)}{1+\alpha \left(I\right(t)+A(t\left)\right)}S\left(t\right)-\mu S\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}E\left(t\right)}{\mathrm{d}t}=& -\frac{c_1{\beta }_1(t-\tau )I(t-\tau )+c_2{\beta }_2(t-\tau )A(t-\tau )}{1+\alpha \left(I\right(t-\tau )+A(t-\tau \left)\right)}S(t-\tau )e^{-\mu \tau }-\mu E\left(t\right)\hfill \\ \hfill & +\frac{c_1{\beta }_1\left(t\right)I\left(t\right)+c_2{\beta }_2\left(t\right)A\left(t\right)}{1+\alpha \left(I\right(t)+A(t\left)\right)}S\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=& \frac{c_1{\beta }_1(t-\tau )I(t-\tau )+c_2{\beta }_2(t-\tau )A(t-\tau )}{1+\alpha \left(I\right(t-\tau )+A(t-\tau \left)\right)}S(t-\tau )e^{-\mu \tau }\hfill \\ \hfill & -(\mu +\xi +{\delta }_1\left(t\right))I\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}=& \xi I\left(t\right)-(\mu +{\delta }_2\left(t\right))A\left(t\right).\hfill \end{array}\right.$

Since the S, I, and A equations in Model (33) are independent of the variable E, it suffices to study the following model

(34)$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}=& {\mathsf{\Lambda }}_h-\frac{c_1{\beta }_1\left(t\right)I\left(t\right)+c_2{\beta }_2\left(t\right)A\left(t\right)}{1+\alpha \left(I\right(t)+A(t\left)\right)}S\left(t\right)-\mu S\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=& \frac{c_1{\beta }_1(t-\tau )I(t-\tau )+c_2{\beta }_2(t-\tau )A(t-\tau )}{1+\alpha \left(I\right(t-\tau )+A(t-\tau \left)\right)}S(t-\tau )e^{-\mu \tau }-(\mu +\xi +{\delta }_1\left(t\right))I\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}=& \xi I\left(t\right)-(\mu +{\delta }_2\left(t\right))A\left(t\right).\hfill \end{array}\right.$

Here, ${\beta }_1\left(t\right)$, ${\beta }_2\left(t\right)$, ${\delta }_1\left(t\right)$, and ${\delta }_2\left(t\right)$ are $\omega$-periodic, continuous, and positive functions.

6.1. Basic Reproduction Number

We now apply the method developed in [41] to derive the basic reproduction number, ${\mathcal{R}}_0$, for Model (34), and then give a characterization of ${\mathcal{R}}_0$. Letting $I=A=0$ in Model (34), one obtains

(35)$\begin{array}{c}\hfill \frac{\mathrm{d}S\left(t\right)}{\mathrm{d}t}={\mathsf{\Lambda }}_h-\mu S\left(t\right).\end{array}$

As a result, Model (35) admits a unique globally asymptotically stable equilibrium $\frac{{\mathsf{\Lambda }}_h}{\mu }$, that is, Model (34) exists as a disease-free periodic solution, $E_0\left(t\right)=(S^{*}\left(t\right),0,0)$, with periodic $\omega$, where $S^{*}\left(t\right)=\frac{{\mathsf{\Lambda }}_h}{\mu }$. After linearizing Model (34) at $E_0\left(t\right)$, we derive the periodic linear equations as follows

(36)$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}& =x_1\left(t\right)I(t-\tau )+x_2\left(t\right)A(t-\tau )-x_3\left(t\right)I\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}& =x_4\left(t\right)I\left(t\right)-x_5\left(t\right)A\left(t\right),\hfill \end{array}\right.$

where $x_1\left(t\right)=c_1{\beta }_1(t-\tau )S^{*}\left(t\right)e^{-\mu \tau }$, $x_2\left(t\right)=c_2{\beta }_2(t-\tau )S^{*}\left(t\right)e^{-\mu \tau }$, $x_3\left(t\right)=\mu +\xi +{\delta }_1\left(t\right)$, $x_4\left(t\right)=\xi$, $x_5\left(t\right)=\mu +{\delta }_2\left(t\right)$. Define $\mathcal{C}=C([-\tau ,0],{\mathbb{R}}^2)$ and ${\mathcal{C}}^{+}=C([-\tau ,0],{\mathbb{R}}_{+}^2)$. It can be established that $(\mathcal{C},{\mathcal{C}}^{+})$ forms an ordered Banach space, equipped with the maximum norm and the positive cone ${\mathcal{C}}^{+}$. Let $v=(v_1,v_2):[-\tau ,\check{\sigma })\to {\mathbb{R}}^2$ be a given continuous function with $\breve{\sigma }>0$. We define $v_t\in \mathcal{C}$ as $v_t\left(\theta \right)=(v_1(t+\theta ),v_2(t+\theta ))$, for all $t\in [0,\breve{\sigma })$, $\theta \in [-\tau ,0]$. Let

$\mathbb{F}\left(t\right)\psi =\left(\begin{array}{c}{x}_1\left(t\right){\psi }_1(-\tau )+x_2\left(t\right){\psi }_2(-\tau )\\ x_4\left(t\right){\psi }_1\left(0\right)\end{array}\right),\mathbb{V}\left(t\right)=\left(\begin{array}{cc}{x}_3\left(t\right)& 0\\ 0& x_5\left(t\right)\end{array}\right),$

where $\mathbb{F}:\mathbb{R}\to \mathcal{L}(\mathcal{C},{\mathbb{R}}^2)$ and $\mathbb{V}\left(t\right)$ is a continuous $2\times 2$ matrix function on $\mathbb{R}$.

Let ${\mathcal{C}}_{\omega }$ be the ordered Banach space of all continuous and $\omega$-periodic functions from $\mathbb{R}$ to ${\mathbb{R}}^2$, equipped with the maximum norm and the positive cone ${\mathcal{C}}_{\omega }^{+}:=\{v\in {\mathcal{C}}_{\omega }:v\left(t\right)\ge 0,\forall t\in \mathbb{R}\}$. According to [41] (Section 2), we define

$\left[Lv\right]\left(t\right)={\int }_0^{\infty }\mathsf{\Phi }(t,t-\overline{s})\mathbb{F}(t-\overline{s})v(t-\overline{s}+·)\mathrm{d}\overline{s},$

for all $t\in \mathbb{R}$, $v\in {\mathcal{C}}_{\omega }$, then $L:{\mathcal{C}}_{\omega }\to {\mathcal{C}}_{\omega }$ is a linear operator. And the evolution matrix $\mathsf{\Phi }(t,\overline{s}),t\ge \overline{s}$ satisfies

$\frac{\partial }{\partial t}\mathsf{\Phi }(t,\overline{s})=-\mathbb{V}\left(t\right)\mathsf{\Phi }(t,\overline{s}),\forall t\ge \overline{s},\mathrm{and}\mathsf{\Phi }(\overline{s},\overline{s})=\mathbb{I},\forall \overline{s}\in \mathbb{R},$

for all $\overline{s}\in \mathbb{R}$, where $\mathbb{I}$ is the $2\times 2$ identity matrix. Therefore, one can easily deduce that

$\mathsf{\Phi }(t,\overline{s})=\left(\begin{array}{cc}{e}^{-{\int }_{\overline{s}}^t{x}_3\left(r\right)dr}& 0\\ 0& e^{-{\int }_{\overline{s}}^t{x}_5\left(r\right)dr}\end{array}\right).$

And the exponential growth bound of $\mathsf{\Phi }(t,\overline{s})$ is defined as

$\widehat{\omega }\left(\mathsf{\Phi }\right)=\inf \left\{\tilde{\omega }:\parallel \mathsf{\Phi }(t+\overline{s},\overline{s})\parallel \le M{e}^{\tilde{\omega }t},\exists M\ge 1,\forall \overline{s}\in \mathbb{R},t\ge 0\right\}.$

Note that $\parallel \mathsf{\Phi }(t+\overline{s},\overline{s})\parallel \le max\{e^{-(\mu +\xi +\underline{{\delta }_1})t},e^{-(\mu +\underline{{\delta }_2})t}\}$, then

$\widehat{\omega }\left(\mathsf{\Phi }\right)\le -min\left\{\mu +\xi +\underline{{\delta }_1},\mu +\underline{{\delta }_2}\right\},$

where $\underline{{\delta }_1}={min}_{t\in [0,\omega ]}\left\{{\delta }_1\left(t\right)\right\}$, $\underline{{\delta }_2}={min}_{t\in [0,\omega ]}\left\{{\delta }_2\left(t\right)\right\}$. Therefore, $\mathbb{F}\left(t\right):\mathcal{C}\to {\mathbb{R}}^2$ is positive in a way that can be interpreted as $\mathbb{F}\left(t\right){\mathcal{C}}^{+}\subseteq {\mathbb{R}}_{+}^2$, $-\mathbb{V}\left(t\right)$ is cooperative, and $\widehat{\omega }\left(\mathsf{\Phi }\right)<0$. And one can calculate

$\begin{array}{cc}\hfill \left[Lv\right]\left(t\right)& ={\int }_0^{\infty }\left(\begin{array}{c}{e}^{-{\int }_{t-\overline{s}}^t{x}_3\left(r\right)\mathrm{d}r}(x_1(t-\overline{s})v_1(t-\overline{s}-\tau )+x_2(t-\overline{s})v_2(t-\overline{s}-\tau ))\\ e^{-{\int }_{t-\overline{s}}^t{x}_5\left(r\right)\mathrm{d}r}{x}_4(t-\overline{s})v_1(t-\overline{s})\end{array}\right)\mathrm{d}\overline{s}\hfill \\ \hfill & =\left(\begin{array}{c}{\int }_{\tau }^{\infty }{e}^{-{\int }_{t-\overline{s}+\tau }^t{x}_3\left(r\right)\mathrm{d}r}(x_1(t-\overline{s}+\tau )v_1(t-\overline{s})+x_2(t-\overline{s}+\tau )v_2(t-\overline{s}))\mathrm{d}\overline{s}\\ {\int }_0^{\infty }{e}^{-{\int }_{t-\overline{s}}^t{x}_5\left(r\right)\mathrm{d}r}{x}_4(t-\overline{s})v_1(t-\overline{s})\mathrm{d}\overline{s}\end{array}\right)\hfill \\ \hfill & ={\int }_0^{\infty }K(t,\overline{s})v(t-\overline{s})\mathrm{d}\overline{s},\forall t\in \mathbb{R},v={(v_1,v_2)}^{\mathrm{T}}\in {\mathcal{C}}_{\omega },\hfill \end{array}$

where, if $\overline{s}\ge \tau$, then

$\begin{array}{cc}\hfill K(t,\overline{s})& =\left(\begin{array}{cc}{e}^{-{\int }_{t-\overline{s}+\tau }^t{x}_3\left(r\right)\mathrm{d}r}{x}_1(t-\overline{s}+\tau )& e^{-{\int }_{t-\overline{s}+\tau }^t{x}_3\left(r\right)\mathrm{d}r}{x}_2(t-\overline{s}+\tau )\\ e^{-{\int }_{t-\overline{s}}^t{x}_5\left(r\right)\mathrm{d}r}{x}_4(t-\overline{s})& 0\end{array}\right);\hfill \end{array}$

and if $\overline{s}<\tau$, then

$\begin{array}{cc}\hfill K(t,\overline{s})& =\left(\begin{array}{cc}0& 0\\ e^{-{\int }_{t-\overline{s}}^t{x}_5\left(r\right)\mathrm{d}r}{x}_4(t-\overline{s})& 0\end{array}\right).\hfill \end{array}$

According to [41], we have ${\mathcal{R}}_0=r\left(L\right)$, where $r\left(L\right)$ represents the spectral radius of L.

Let $P\left(t\right)$ be the solution maps of (36), that is, $P\left(t\right)\psi =u_t\left(\psi \right)$, $t\ge 0$, where $u(t,\psi )$ is the unique solution of (36), with $u_0=\psi \in C([-\tau ,0],{\mathbb{R}}^2)$. Then $P:=P\left(\omega \right)$ is the Poincaré map associated with (36). Let $r\left(P\right)$ be the spectral radius of P. By the arguments shown in [41] (Theorem 2.1), we have the following conclusion.

Lemma 3.

The sign of ${\mathcal{R}}_0-1$ is the same as the sign of $r\left(P\right)-1$.

6.2. Threshold Dynamics

In this subsection, we aim to analyze the global dynamics of Model (34) with respect to ${\mathcal{R}}_0$. Let $\mathbb{Z}=C([-\tau ,0],{\mathbb{R}}^3)$, ${\mathbb{Z}}^{+}=C([-\tau ,0],{\mathbb{R}}_{+}^3)$, then the following Lemma 4 holds.

Lemma 4.

For any $\psi \in {\mathbb{Z}}^{+}$, there exists a unique nonnegative bounded solution, $u(t,\psi )$, of Model (34), with $u_0=\psi$ for all $t\ge 0$. Moreover, Model (34) gives rise to an ω-periodic semiflow $Q\left(t\right)=u_t:{\mathbb{Z}}^{+}\to {\mathbb{Z}}^{+}$, that is, $Q\left(t\right)\psi =u_t\left(\psi \right)$ for all $t\ge 0$, and the semiflow $Q:=Q\left(\omega \right)$ possesses a strong global attractor.

Proof. 

For any $\psi =({\psi }_1,{\psi }_2,{\psi }_3)\in {\mathbb{Z}}^{+}$, we define

$f(t,\psi )=\left(\begin{array}{c}{\mathsf{\Lambda }}_h-\frac{c_1{\beta }_1\left(t\right){\psi }_3\left(0\right)+c_2{\beta }_2\left(t\right){\psi }_4\left(0\right)}{1+\alpha ({\psi }_3\left(0\right)+{\psi }_4\left(0\right))}{\psi }_1\left(0\right)-\mu {\psi }_1\left(0\right)\\ \frac{c_1{\beta }_1(t-\tau ){\psi }_3(-\tau )+c_2{\beta }_2(t-\tau ){\psi }_4(-\tau )}{1+\alpha ({\psi }_3(-\tau )+{\psi }_4(-\tau ))}{\psi }_1(-\tau )e^{-\mu \tau }-(\mu +\xi +{\delta }_1\left(t\right)){\psi }_3\left(0\right)\\ \xi {\psi }_3\left(0\right)-(\mu +{\delta }_2\left(t\right)){\psi }_4\left(0\right)\end{array}\right).$

As $f(t,\psi )$ is continuous in $(t,\psi )\in {\mathbb{R}}_{+}\times {\mathbb{Z}}^{+}$ and Lipschitz in $\psi$ on all compact subsets of ${\mathbb{Z}}^{+}$, it can be inferred from [35] (Theorems 2.2.1 and 2.2.3) that Model (34) admits a unique solution, $u(t,\psi )$, on its maximal interval, $[0,{\sigma }_{\psi })$, of existence with initial condition $u_0=\psi$.

Let $\psi =({\psi }_1,{\psi }_2,{\psi }_3)$. If ${\psi }_i\left(0\right)=0$, $i\in \{1,2,3\}$, then $f_i(t,\psi )\ge 0$. By employing [42] (Remark 5.2.1 and Theorem 5.2.1), for all $t\in [0,{\sigma }_{\psi })$, the solution, $u(t,\psi )$, of Model (34) is unique with initial condition $u_0=\psi$ fulfilling $u_t\left(\psi \right)\in {\mathbb{Z}}^{+}$. Define

$\mathbb{X}=\left\{\varphi \in C([-\tau ,0],{\mathbb{R}}_{+}^4):{\varphi }_2\left(0\right)={\int }_{-\tau }^0\frac{c_1{\beta }_1\left(\overline{s}\right){\varphi }_3\left(\overline{s}\right)+c_2{\beta }_2\left(\overline{s}\right){\varphi }_4\left(\overline{s}\right)}{1+\alpha ({\varphi }_3\left(\overline{s}\right)+{\varphi }_4\left(\overline{s}\right))}{\varphi }_1\left(\overline{s}\right)e^{\mu \overline{s}}\mathrm{d}\overline{s}\right\}.$

It can be readily deduced that, for any $\varphi \in \mathbb{X}$, (33) has a unique nonnegative solution, $v(t,\varphi )$, which satisfies the initial condition, $v_0=\varphi$, for all $t\in [0,{\sigma }_{\psi })$. Let $N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+A\left(t\right)$. Consequently, $\frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}={\mathsf{\Lambda }}_h-\mu N\left(t\right)-{\delta }_1\left(t\right)I\left(t\right)-{\delta }_2\left(t\right)A\left(t\right)\le {\mathsf{\Lambda }}_h-\mu N\left(t\right)$, for all $t\in [0,{\sigma }_{\psi })$. Hence, $S\left(t\right)$, $E\left(t\right)$, $I\left(t\right)$, and $R\left(t\right)$ are bounded on $t\in [0,{\sigma }_{\psi })$. It implies that ${\sigma }_{\psi }=+\infty$ [35]. Then we have

(37)$\begin{array}{c}\hfill \frac{\mathrm{d}N\left(t\right)}{\mathrm{d}t}={\mathsf{\Lambda }}_h-\mu N\left(t\right)-{\delta }_1\left(t\right)I\left(t\right)-{\delta }_2\left(t\right)A\left(t\right)\le {\mathsf{\Lambda }}_h-\mu N\left(t\right),t\ge 0.\end{array}$

Thus, the global stability of $S^{*}\left(t\right)$ for (35) and the comparison argument imply that the solutions of Model (33) with initial values in $\mathbb{X}$, and hence Model (34) in ${\mathbb{Z}}^{+}$, are ultimately bounded and exist globally on $[0,+\infty )$. Consequently, $Q\left(t\right):{\mathbb{Z}}^{+}\to {\mathbb{Z}}^{+}$ is point dissipative. It is worth noting that, for all $t>\tau$, $Q\left(t\right)$ is compact [35]. By applying the results of [43], we can conclude that $Q:=Q\left(\omega \right)$ has a strong global attractor in ${\mathbb{Z}}^{+}$. □

The solution maps of (36) on $\mathbb{Y}:=C([-\tau ,0],{\mathbb{R}}_{+})\times {\mathbb{R}}_{+}$ are denoted as $\widehat{P}\left(t\right)$. The Poincaré map associated with (36) on $\mathbb{Y}$ is represented by $\widehat{P}:=\widehat{P}\left(\omega \right)$. The spectral radius of $\widehat{P}$ is denoted as $r\left(\widehat{P}\right)$. Based on similar arguments presented in [44] (Lemma 3.8), it can be concluded that $r\left(P\right)=r\left(\widehat{P}\right)$. Furthermore, it can be shown, in accordance with [35] (Theorem 3.6.1 and Lemma 5.3.2), that $\widehat{P}\left(t\right)$ is strongly positive and compact on $\mathbb{Y}$ for $t>3\tau$. Selecting an integer $n_0>0$, such that $n_0\omega >3\tau$, we can use [38] (Lemma 3.1) to deduce that $\lambda =r\left(\widehat{P}\right)>0$. Here, $\lambda$ represents a simple eigenvalue of $\widehat{P}$ with a strongly positive eigenvector, while the modulus of any other eigenvalue is smaller than $r\left(\widehat{P}\right)$. By referring to the proof of [45] (Lemma 1), we can conclude that the following result holds.

Lemma 5.

Let $\widehat{\mu }=\frac{lnr\left(P\right)}{\omega }=\frac{lnr\left(\widehat{P}\right)}{\omega }$, then there is a positive ω-periodic function $\overline{v}\left(t\right)={({\overline{v}}_1\left(t\right),{\overline{v}}_2\left(t\right))}^T$ such that $v^{*}\left(t\right)=e^{\widehat{\mu }t}\overline{v}\left(t\right)$ is a positive solution of (36).

Next, the threshold condition that determines disease extinction will be explored.

Theorem 16.

If ${\mathcal{R}}_0<1$, then $E_0\left(t\right)$ is globally attractive for Model (34) in ${\mathbb{Z}}^{+}$.

Proof. 

For the case of ${\mathcal{R}}_0<1$, one obtains $r\left(P\right)<1$. Let $P_{ϵ}$ denote the Poincaré map of the auxiliary model described below

(38)$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}& =(c_1{\beta }_1(t-\tau )I(t-\tau )+c_2{\beta }_2(t-\tau )A(t-\tau ))(S^{*}\left(t\right)+ϵ)e^{-\mu \tau },\hfill \\ \hfill \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}& =\xi I\left(t\right)-(\mu +{\delta }_2\left(t\right))A\left(t\right).\hfill \end{array}\right.$

Since ${lim}_{ϵ\to 0^{+}}r\left(P_{ϵ}\right)=r\left(P\right)<1$, we can fix a sufficiently small constant, $ϵ\in (0,\frac{{\mathsf{\Lambda }}_h}{\mu })$, such that $r\left(P_{ϵ}\right)<1$. According to Lemma 5, there is a positive $\omega$-periodic function ${\overline{v}}_{ϵ}\left(t\right)=({\overline{v}}_{1ϵ}\left(t\right),{\overline{v}}_{2ϵ}\left(t\right))$ such that $v_{ϵ}^{*}\left(t\right)=e^{{\widehat{\mu }}_{ϵ}t}{\overline{v}}_{ϵ}\left(t\right)$ is a solution of (38) on ${\mathcal{C}}^{+}$, where ${\widehat{\mu }}_{ϵ}=\frac{lnr\left(P_{ϵ}\right)}{\omega }<0$. For any given $\psi \in {\mathbb{Z}}^{+}$, let $v(t,\psi )=\left(S\right(t),I(t),A(t\left)\right)$. In view of (37) and the global stability of $S^{*}\left(t\right)$ for (35), it follows that there exists a sufficiently large integer, $n_1$, such that $n_1\omega \ge \tau$ and $S\left(t\right)\le S^{*}\left(t\right)+ϵ$, for $t\ge n_1\omega -\tau$. Thus, for $t\ge n_1\omega$, we have

$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}& \le (c_1{\beta }_1(t-\tau )I(t-\tau )+c_2{\beta }_2(t-\tau )A(t-\tau ))(S^{*}\left(t\right)+ϵ)e^{-\mu \tau },\hfill \\ \hfill \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}& \le \xi I\left(t\right)-(\mu +{\delta }_2\left(t\right))A\left(t\right).\hfill \end{array}\right.$

Therefore, according to the comparison theorem stated in [42] (Theorem 5.1.1), one has $(I\left(t\right),A\left(t\right))\le K{v}_{ϵ}^{*}\left(t\right)=K{e}^{{\widehat{\mu }}_{ϵ}t}{\overline{v}}_{ϵ}\left(t\right),$ for all $t\ge n_1\omega$, with constant $K>0$ satisfying $(I\left(t\right),A\left(t\right))\le K{v}_{ϵ}^{*}\left(t\right)$ for $t\in [n_1\omega -\tau ,n_1\omega ]$. And hence, ${lim}_{t\to \infty }(I\left(t\right),A\left(t\right))=(0,0)$. Moreover, by using the first equation of (34), we have that ${lim}_{t\to \infty }S\left(t\right)=\frac{{\mathsf{\Lambda }}_h}{\mu }$. This completes the proof. □

According to the epidemiological concept, if ${\mathcal{R}}_0<1$, the number of infected people tends to zero over time, and the disease is extinct. Conversely, the disease-free periodic solution is unstable, namely, no matter the initial state of infection or small the number is, the disease is persistent and eventually becomes endemic.

Let $X_0=\{\psi =({\psi }_1,{\psi }_2,{\psi }_3)\in {\mathbb{Z}}^{+}:{\psi }_2\left(0\right)>0$ and ${\psi }_3\left(0\right)>0\}$, $\partial X_0={\mathbb{Z}}^{+}\setminus X_0=\{\psi =({\psi }_1,{\psi }_2,{\psi }_3)\in {\mathbb{Z}}^{+}:{\psi }_2\left(0\right)=0$ or ${\psi }_3\left(0\right)=0\}$. Then the following conclusion holds.

Theorem 17.

If ${\mathcal{R}}_0>1$, then Model (34) admits a positive ω-periodic solution $\overline{E}\left(t\right)(\overline{S}\left(t\right),\overline{I}\left(t\right),\overline{A}\left(t\right))$. Further, there is a constant, $\tilde{\eta }>0$, such that the solution $u(t,\psi )=\left(S\right(t),I(t),A(t\left)\right)$ of Model (34) with $\psi \in X_0$ satisfies ${lim}_{t\to \infty }(I\left(t\right),A\left(t\right))\ge (\tilde{\eta },\tilde{\eta })$.

Proof. 

In the case where ${\mathcal{R}}_0>1$, we have $r\left(P\right)>1$. Recall that $Q\left(t\right)\psi =u_t\left(\psi \right)$ for $t\ge 0$, where $Q:=Q\left(\omega \right)$ represents the Poincaré map associated with (34). Considering the structure of (34), we can see that $Q\left(t\right)X_0\subset X_0$ holds for all $t\ge 0$. In order to establish the uniform persistence of Q in terms of $(X_0,\partial X_0)$, let $P_{\delta }$ be the Poincaré map of the auxiliary model given below

(39)$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}=& \frac{{\mathsf{\Lambda }}_h-\mu \delta }{\mu (1+2\alpha \delta )}(c_1{\beta }_1(t-\tau )I(t-\tau )+c_2{\beta }_2(t-\tau )A(t-\tau ))e^{-\mu \tau }-(\mu +\xi +{\delta }_1\left(t\right))I\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}=& \xi I\left(t\right)-(\mu +{\delta }_2\left(t\right))A\left(t\right).\hfill \end{array}\right.$

Since ${lim}_{\delta \to 0^{+}}r\left(P_{\delta }\right)=r\left(P\right)>1$, we can fix a constant, $\delta \in (0,\frac{{\mathsf{\Lambda }}_h}{\mu })$, such that $r\left(P_{\delta }\right)>1$. In addition, there exists a function, ${\overline{v}}_{\delta }\left(t\right)=({\overline{v}}_{1\delta }\left(t\right),{\overline{v}}_{2\delta }\left(t\right))$, so that $v_{\delta }^{*}\left(t\right)=e^{{\widehat{\mu }}_{\delta }t}{\overline{v}}_{\delta }\left(t\right)$ is a solution of (39) on ${\mathcal{C}}^{+}$, where ${\widehat{\mu }}_{\delta }=\frac{lnr\left(P_{\delta }\right)}{\omega }>0$.

Given $M_1=({\check{S}}^{*},\check{0},\check{0})$, where ${\check{S}}^{*}\left(\theta \right)=S^{*}$ and $\check{\phi }\left(\theta \right)=0$ for all $\theta \in [-\tau ,0]$, it follows that $Q\left(t\right)M_1=({\check{S}}^{*},\check{0},\check{0})$ holds for all $t\ge 0$, and $Q\left(M_1\right)=M_1$. Moreover, since ${lim}_{\psi \to M_1}(Q\left(t\right)\psi -Q\left(t\right)M_1)=0$ uniformly for $t\in [0,\omega )$, there is a function, $\overline{\delta }=\overline{\delta }\left(\delta \right)$, such that

$\parallel Q\left(t\right)\psi -Q\left(t\right)M_1\parallel \le \delta ,t\in [0,\omega ],\parallel \psi -M_1\parallel <\overline{\delta }.$

Claim 1.

lim ${\sup }_{n\to \infty }\parallel Q^n\left(\psi \right)-M_1\parallel \ge \overline{\delta }$, for all $\psi \in X_0$.

Assuming the claim is false, there exists $\phi \in X_0$ such that lim ${\sup }_{n\to \infty }\parallel Q^n\left(\phi \right)-M_1\parallel <\overline{\delta }$. This implies that there is an integer, $n_2\ge 1$, such that $\parallel Q^n\left(\phi \right)-M_1\parallel <\overline{\delta }$ for all $n\ge n_2$. Now, consider any $t\ge n_2\omega$ and write $t=n\omega +t^{\prime }$, where $n\ge n_2$ and $t^{\prime }\in [0,\omega )$, we obtain $\parallel Q\left(t\right)\phi -Q\left(t\right)M_1\parallel =\parallel Q\left(t^{\prime }\right)\left(Q^n\left(\phi \right)\right)-Q\left(t^{\prime }\right)M_1\parallel <\delta$. It is easy to see that

$\left|S(t,\phi )-\frac{{\mathsf{\Lambda }}_h}{\mu }\right|<\delta ,0<I(t,\phi ),A(t,\phi )<\delta ,t\ge n_2\omega .$

Therefore, for any $t\ge n_2\omega +\tau$, $I(t,\phi )$ and $A(t,\phi )$ satisfy

(40)$\left\{\begin{array}{cc}\hfill \frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}\ge & \frac{{\mathsf{\Lambda }}_h-\mu \delta }{\mu (1+2\alpha \delta )}(c_1{\beta }_1(t-\tau )I(t-\tau )+c_2{\beta }_2(t-\tau )A(t-\tau ))e^{-\mu \tau }-(\mu +\xi +{\delta }_1\left(t\right))I\left(t\right),\hfill \\ \hfill \frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}\ge & \xi I\left(t\right)-(\mu +{\delta }_2\left(t\right))A\left(t\right).\hfill \end{array}\right.$

where $I\left(t\right):=I(t,\phi )$ and $A\left(t\right):=A(t,\phi )$. Since $Q\left(t\right)X_0\subset X_0$ for all $t\ge 0$, we can choose a small enough $\nu >0$ such that $(I(t,\phi ),A(t,\phi ))\ge \nu v_{\delta }^{*}\left(t\right)=\nu e^{{\widehat{\mu }}_{\delta }t}{\overline{v}}_{\delta }\left(t\right),$ for all $t\in [n_2\omega ,n_2\omega +\tau ]$. According to comparison theorem [42] (Theorem 5.1.1), one has $(I(t,\phi ),A(t,\phi ))\ge \nu e^{{\widehat{\mu }}_{\delta }t}{\overline{v}}_{\delta }\left(t\right),$ for all $t\ge n_2\omega +\tau$. Thus, ${lim}_{t\to \infty }I(t,\phi )={lim}_{t\to \infty }A(t,\phi )=\infty$, a contradiction. Therefore, Claim 1 is valid.

This claim implies that $M_1$ is an isolated invariant set for Q in ${\mathbb{Z}}^{+}$ and the intersection $W^s\left(M_1\right)\cap X_0$ is empty, where $W^s\left(M_1\right)$ denotes the stable set of $M_1$ for Q. Let $M_{\partial }:=\{\psi \in \partial X_0:Q^n\left(\psi \right)\in \partial X_0,\forall n\ge 0\}.$ For any given $\phi \in M_{\partial }$, $Q^n\left(\phi \right)\in \partial X_0$ for all $n\ge 0$. Therefore, $I(n\omega ,\phi )=0$ or $A(n\omega ,\phi )=0$. Since, for all $t\ge 0$,

$\frac{\mathrm{d}I\left(t\right)}{\mathrm{d}t}\ge -(\mu +\xi +{\delta }_1\left(t\right))I\left(t\right),\frac{\mathrm{d}A\left(t\right)}{\mathrm{d}t}\ge -(\mu +{\delta }_2\left(t\right))A\left(t\right).$

We conclude that, if $I\left(t_0\right)>0$, or $A\left(t_0\right)>0$ for some $t_0\ge 0$, then $I\left(t\right)>0$ or $A\left(t\right)>0$ for all $t\ge t_0$. It implies that $I(t,\phi )\equiv 0$ or $A(t,\phi )\equiv 0$ for all $t\ge 0$. If $I(t,\phi )=0$ for each $t\ge 0$, according to the first and third equations of (34), we see that $A(t,\phi )\to 0$ as $t\to \infty$, and thus $S(t,\phi )-S^{*}\to 0$ as $t\to \infty$. If $I(t_0,\phi )>0$ for $t_0\ge 0$, the second equation of (34) implies

$I(t,\phi )\ge I(t_0,\phi )e^{-{\int }_{t_0}^t(\mu +\xi +{\delta }_1\left(\overline{s}\right))\mathrm{d}\overline{s}}>0,$

for all $t\ge t_0$. Hence, we can conclude that $A(t,\phi )=0$ for all $t\ge t_0$. Using the first and third equations in (34), there holds that $A(t,\phi )\to 0$ as $t\to \infty$, which results in $S(t,\phi )-S^{*}\to 0$ as $t\to \infty$. So one can establish that $\omega \left(\phi \right)=M_1$ for any $\phi \in M_{\partial }$, and $M_1$ cannot form a cycle for Q in $\partial X_0$. Based on the acyclicity theorem for uniform persistence of maps (refer to [43] (Theorem 1.3.1 and Remark 1.3.1)), it follows that $Q:{\mathbb{Z}}^{+}\to {\mathbb{Z}}^{+}$ exhibits uniform persistence regarding $X_0$.

Observe that, for any integer n with $n\omega >\tau$, $Q^n=Q\left(n\omega \right):{\mathbb{Z}}^{+}\to {\mathbb{Z}}^{+}$ is compact. By [43] (Theorem 3.5.1) and [46] (Theorem 4.1.1), there exists an equivalent norm for $\mathbb{Z}$ such that, for any $t>0$, the map $Q\left(t\right)$ acts as an $\alpha$-contraction on $\mathbb{Z}$. Therefore, it can be inferred from [43] (Theorem 1.3.10) that the map $Q:X_0\to X_0$ possesses a global attractor $Y_0$, and, within this attractor, Q has a fixed point denoted by ${\psi }^{★}\in Y_0$. Then, $u(t,{\psi }^{★})=(S(t,{\psi }^{★}),I(t,{\psi }^{★}),A(t,{\psi }^{★}))$ is an $\omega$-periodic solution of (34). And $S^{★}\left(t\right)\ge 0$, $I^{★}\left(t\right)\ge 0$, $A^{★}\left(t\right)\ge 0$.

We claim that there exists some $t^{★}\in [0,\omega ]$ such that $S^{★}\left(t^{★}\right)>0$. If this statement does not hold, it would imply that $S^{★}\left(t\right)\equiv 0$ for all $t\ge 0$, because of the periodicity of $S^{★}\left(t^{★}\right)>0$. Consequently, based on the first equation in (34), we would have $0=\frac{{\mathsf{\Lambda }}_h}{\mu }$, which leads to a contradiction. Since

$\begin{array}{c}\hfill \frac{\mathrm{d}{S}^{★}\left(t\right)}{\mathrm{d}t}\ge -\left(\frac{c_1{\beta }_1\left(t\right)I^{★}\left(t\right)+c_2{\beta }_2\left(t\right)A^{★}\left(t\right)}{1+\alpha (I^{★}\left(t\right)+A^{★}\left(t\right))}+\mu \right)S^{★}\left(t\right),\end{array}$

we can conclude that $S^{★}\left(t\right)>0$ for $t\ge t^{★}$. Due to the periodic nature of $S^{★}\left(t\right)$, it follows that $S^{★}\left(t\right)>0$ for all $t\ge 0$. By utilizing the third equation in (34) and the fact that $I^{★}\left(t\right)>0$ for $t\ge 0$, one can deduce that $A^{★}\left(t\right)>0$ for $t\ge 0$. Consequently, $(S^{★}\left(t\right),I^{★}\left(t\right),A^{★}\left(t\right))$ serves as a positive $\omega$-periodic solution of (34). Furthermore, by making an argument similar to the proof in [45] (Theorem 5), the practical uniform persistence can be obtained, i.e, there exists $\check{\eta }>0$ such that ${lim}_{t\to \infty }(I(t,\psi ),A(t,\psi ))\ge (\check{\eta },\check{\eta })$. The proof is finished. □

7. Numerical Simulation and Discussion

In this section, some numerical simulations are performed to explain our theoretical results and examine the effect of delay and immune response on the transmission of HIV-1 using the Runge-Kutta method in software MATLAB. More specifically, numerical simulations will be fulfilled for three examples: the first one considers the transmission dynamics of the coupled between-host model without an immune response, the second case considers the coupled slow time model with an immune response, in which the immune delay is in its stable interval, and the third example is dedicated to discuss the impact of the stable periodic solution formed in the host as a result of the immune delay on the disease transmission between hosts.

For the convenience of numerical simulations, we also take ${\beta }_i\left(V\right)={\overline{a}}_{i}V$, ${\delta }_i\left(V\right)={\delta }_i$$(i=1,2)$, as mentioned in Section 4.3 for the first example, and take ${\beta }_i\left(V\right)={\overline{a}}_{i}V$, ${\delta }_i\left(V\right)={\delta }_{i}V$$(i=1,2)$ for the last two examples. And the parameter values for the within-host model were obtained from [11,21,47]. For the between-host model, the parameter values were derived from [33].

Example 1.

The dynamics of the coupled slow time model, where the fast time model is without an immune response.

In Model (4), we choose parameters as follows: ${\mathsf{\Lambda }}_c=900$, ${\kappa }_1=8.75\times {10}^{-7}$, ${\kappa }_2=3.75\times {10}^{-7}$, ${\mu }_c=0.0129$, ${\delta }_c=0.52$, $p=80$, $q=0.053$, ${\mu }_v=0.4$, $\sigma =0.1\times {10}^{-5}$, $\zeta =0.68$, $c_1=10$, $c_2=5.3$, ${\eta }_1=0.4\times {10}^{-4}$, ${\eta }_2=0.4\times {10}^{-4}$, $\alpha =5.5\times {10}^{-5}$, ${\tau }_1=0.8$, ${\tau }_2=3.5$. By the direct calculation, one has ${\mathcal{R}}_w=0.4812<1$. As shown in Figure 2a,b, ${\tilde{U}}_1$ is globally asymptotically stable, which coincides with the conclusion $\left(i\right)$ of Theorem 3.

Further, taking the parameters ${\mathsf{\Lambda }}_h=2500$, $\mu =0.0045$, ${\overline{a}}_1=9\times {10}^{-9}$, ${\overline{a}}_2=9\times {10}^{-11}$, ${\delta }_1=0.2$, ${\delta }_2=0.25$, $\tau =5$, and $\xi =0.66$, we can calculate ${\mathcal{R}}_0^w=1.0001>1$ and ${\mathcal{R}}_0^h=0.4152<1$. Therefore, Model (20) has two positive equilibriums, ${\tilde{E}}_1$ and ${\tilde{E}}_2$, as stated in Theorem 9, where one is unstable, and the other is locally asymptotically stable. As depicted in Figure 3a, the solutions converge to either the disease-free equilibrium or the positive equilibrium, depending on initial conditions. However, when ${\mathsf{\Lambda }}_h=1900$ and other parameters are fixed, as in Figure 3a, we obtain ${\mathcal{R}}_0^w=1.0001>1$ and ${\mathcal{R}}_0^h=0.3155<1$. In this case, we find that the solutions with a lower initial value converge to the disease-free equilibrium, whereas the others with a higher initial value may oscillate periodically, which can be observed in Figure 3b. Next, we just change the values of parameters as ${\mathsf{\Lambda }}_h=2000$, ${\mathsf{\Lambda }}_c=998$, and other model parameters are chosen, as in Figure 3a, and we obtain ${\mathcal{R}}_0^w=0.9995<1$, and the dynamical behavior is shown in Figure 3c. What is more, if we tack ${\mathsf{\Lambda }}_c=990$, then ${\mathcal{R}}_0^w=0.9915<1$, which corresponds to Figure 3d. Clearly, the stability of equilibria given in Figure 3c,d are similar to those in Figure 3a and b, respectively. These show that the system may exhibit more diverse dynamics, except for the backward bifurcation, and the basic reproduction number, ${\mathcal{R}}_0^w$, no longer serves as a threshold condition determining disease persistence or extinction. Lastly, let ${\overline{a}}_2=9\times {10}^{-9}$, ${\mathsf{\Lambda }}_h=1800$, with other parameter values unchanged in Model (20). By numerical calculations, we obtain ${\mathcal{R}}_0^w=1.0001>1$ and ${\mathcal{R}}_0^h=11.8351>1$. Therefore, from Theorem 9, Model (20) has a unique positive equilibrium, $\tilde{E}$. And Figure 3e illustrate that $\tilde{E}$ is globally asymptotically stable. However, if we take ${\mathsf{\Lambda }}_h=1800$ and fix other model parameter values as above, then the basic reproduction numbers ${\mathcal{R}}_0^w=1.0001>1$ and ${\mathcal{R}}_0^h=5.0528>1$ can be obtained. In this situation, as shown in Figure 3f, Model (20) has a positive periodic solution.

It is noticed that the existence of equilibria concluded in Theorem 9 and their stability revealed in Figure 3a–f may indicate the existence of saddle-node bifurcation, backward bifurcation, and Hopf bifurcation for our coupled slow time model (20). That is, Model (20) has bistable attractors for ${\mathcal{R}}_0^h<1$ (or ${\mathcal{R}}_0^w<1$), and a unique stable positive equilibrium for ${\mathcal{R}}_0^w>1$ and ${\mathcal{R}}_0^h>1$. These are shown in Figure 4a,b, in which the number of infected individuals at positive equilibrium is expressed as a function of ${\mathcal{R}}_0^w$ and ${\mathcal{R}}_0^h$, with the black solid curve, red solid curve, and blue dashed curve representing the locally asymptotically stable, conditionally stable, and unstable equilibrium, respectively. Especially, Figure 4a shows that the model has two positive equilibria, ${\tilde{E}}_1$ and ${\tilde{E}}_2$, in which one is conditionally stable and the other is unstable for ${\mathcal{R}}_{hc}<{\mathcal{R}}_0^h<1$, and has no positive equilibrium for ${\mathcal{R}}_0^h<{\mathcal{R}}_{hc}$. This implies that the threshold value for disease eradication is not ${\mathcal{R}}_0^h=1$, but ${\mathcal{R}}_0^h={\mathcal{R}}_{hc}=0.252$. That is to say, the infectious disease may become extinct until ${\mathcal{R}}_0^h<{\mathcal{R}}_{hc}$. And starting from ${\mathcal{R}}_0^h={\mathcal{R}}_{hc}$, with the increase of ${\mathcal{R}}_0^h$, Model (20) undergoes saddle-node bifurcation at ${\mathcal{R}}_{hc}=0.252$, forward Hopf bifurcation that changes the stability of the positive equilibrium, ${\tilde{E}}_1$, from unstable to stable, and backward bifurcation at ${\mathcal{R}}_0^h=1$ in turn. When ${\mathcal{R}}_0^h>1$, Model (20) has a unique conditionally stable positive equilibrium, $\tilde{E}$, which exhibits different dynamical behavior under different values of ${\mathcal{R}}_0^h$. That is, it is globally asymptotically stable for larger ${\mathcal{R}}_0^h$, but loses stability for smaller ${\mathcal{R}}_0^h$, and the solution of the model oscillates periodically around it. These complex phenomena highlight the challenges in the control of disease. Similar to the situation revealed in Figure 4a, Figure 4b shows that complex dynamic behaviors can exist for ${\mathcal{R}}_0^w<1$, and the lower bound of the threshold condition for disease control is ${\mathcal{R}}_0^w={\mathcal{R}}_{wc}=0.99838$. This suggests that between-host disease may still persist, resulting from contact with infected individuals, even though the virus can be cleared within the host. The parameter values used in Figure 4a are ${\overline{a}}_1=3\times {10}^{-14}$, ${\overline{a}}_2=2\times {10}^{-14}$, ${\kappa }_1=8.75\times {10}^{-8}$, ${\kappa }_2=3.75\times {10}^{-8}$, ${\mu }_c=0.0129$, ${\tau }_1=0.8$, ${\mathsf{\Lambda }}_c=500$, ${\delta }_c=0.52$, $p=80$, $q=0.053$, ${\mu }_v=0.4$, $\sigma =0.003$, $\zeta =0.068$, $c_1=10$, $c_2=5.3$, ${\eta }_1=0.00004$, ${\eta }_2=0.00004$, $\mu =0.00009$, ${\mathsf{\Lambda }}_h=2000$, ${\delta }_1=0.00055$, ${\delta }_2=0.85$, $\tau =5$, $c_1=1000$, $c_2=80$, $\xi =0.2$, and $\alpha =0.000005$. In Figure 4b, parameters have the same values as in Figure 4a except that ${\overline{a}}_1=3\times {10}^{-14}$.

Example 2.

The dynamics of the coupled slow time model, where the fast time model has an immune response.

Now, we choose the values of parameters in Model (4) as follows: ${\mathsf{\Lambda }}_c=900$, ${\kappa }_1=8.75\times {10}^{-7}$, ${\kappa }_2=3.75\times {10}^{-7}$, ${\mu }_c=0.0129$, ${\delta }_c=0.52$, $p=80$, $q=0.053$, ${\mu }_v=0.5$, $\sigma =0.002$, $\zeta =0.68$, $c_1=10$, $c_2=5.3$, ${\eta }_1=0.00004$, ${\eta }_2=0.00004$, $\alpha =0.000055$, ${\tau }_1=0.8$, ${\tau }_2=2.8$, then the immune-activated reproduction number ${\mathcal{R}}_w=19.0010>1$ and ${\tau }_2^{*}=3.22$ can be calculated. Therefore, the fast time model (4) has a globally asymptotically stable immunity-activated infection equilibrium, ${\tilde{U}}_2$, from Theorem 3, and ${\tilde{V}}_2=\frac{\zeta }{\sigma }=340$. These are shown in Figure 5a,b.

Further, we let ${\mathsf{\Lambda }}_h=2500$, $\mu =0.0095$, ${\overline{a}}_1=2.07\times {10}^{-10}$, ${\overline{a}}_2=1.07\times {10}^{-10}$, ${\delta }_1=0.05$, ${\delta }_2=0.35$, $\tau =5$, $\xi =0.9$, $\alpha =0.55$, and replace the ${\tilde{V}}_2=340$ in the coupled slow time model (28), then ${\tilde{\mathcal{R}}}_0^h=0.3103<1$. And Theorem 12 suggests that the Model (28) has a unique equilibrium, $E_0$, which is globally asymptotically stable. Our numerical simulations confirm this result, as shown in Figure 6a,b. Next, if we take ${\mathsf{\Lambda }}_h=5000$, $\mu =0.009$, ${\overline{a}}_1=9.2\times {10}^{-10}$, ${\overline{a}}_2=9.2\times {10}^{-10}$, ${\delta }_1=0.05$, ${\delta }_2=0.22$, $\tau =5$, $\xi =0.15$, and $\alpha =0.02$, then ${\tilde{\mathcal{R}}}_0^h=10.7084>1$. As depicted in Figure 6c,d, the solutions originating from various initial values converge to the unique endemic equilibrium, $E^{*}$. Thus, the global asymptotic stability of $E^{*}$ stated in Remark 8 would be reasonable. Clearly, if immune cells are stimulated to produce, that is, humoral immunity is involved in host defense mechanisms, and the immune delay, ${\tau }_2$, is in its stable interval, then whether the infectious disease between hosts is extinct simply depends on whether ${\tilde{\mathcal{R}}}_0^h$ is less than 1.

Example 3.

The transmission dynamics of the coupled slow time model, where the fast time model has an immune response and stable periodic solution.

Taking the same parameter values as Figure 5 in the fast time model (4), except for ${\tau }_2=3.4>{\tau }_2^{*}=3.22$ (that is, the delay ${\tau }_2$ is not in its stability interval, for more details see Theorem 5), as shown in Figure 7a, the infection equilibrium, ${\tilde{U}}_2$, of Model (4) is unstable. And the plots in Figure 7b–d show that the solutions starting from a different initial value gradually move away from ${\tilde{U}}_2$ and eventually converge to a stable limit cycle.

In Section 6, we discuss in detail the dynamics of the coupled slow time model in the presence of a stable periodic solution for the fast time model (4). Since the exact expression for the periodic solution $(\tilde{T}\left(s\right),{\tilde{T}}^{*}\left(s\right),\tilde{V}\left(s\right),\tilde{B}\left(s\right))$ of the model is not available, we obtain an approximate expression for $\tilde{V}\left(s\right)$ by function fitting. This is shown in Figure 8a, where $V\left(s\right)$ exhibits a periodic oscillatory behavior when s is sufficiently large, which is exactly the starting point for our search for this periodic solution. Fortunately, one can observe in Figure 8b that the red curve described by function

$\begin{array}{cc}\hfill f\left(s\right)=& a_0+a_{1}cos\left(sw\right)+b_{1}sin\left(sw\right)+a_{2}cos\left(2sw\right)+b_{2}sin\left(2sw\right)+a_{3}cos\left(3sw\right)+b_{3}sin\left(3sw\right)\hfill \\ \hfill & +a_{4}cos\left(4sw\right)+b_{4}sin\left(4sw\right)+a_{5}cos\left(5sw\right)+b_{5}sin\left(5sw\right)+a_{6}cos\left(6sw\right)+b_{6}sin\left(6sw\right)\hfill \\ \hfill & +a_{7}cos\left(7sw\right)+b_{7}sin\left(7sw\right)+a_{8}cos\left(8sw\right)+b_{8}sin\left(8sw\right),\hfill \end{array}$

where $a_0=425.6$, $a_1=-253.2$, $b_1=-312.4$, $a_2=-40.24$, $b_2=128$, $a_3=36.29$, $b_3=-6.733$, $a_4=-5.112$, $b_4=-7.529$, $a_5=-1.208$, $b_5=1.705$, $a_6=0.441$, $b_6=0.1249$, $a_7=-0.00434$, $b_7=-0.1035$, $a_8=-0.02095$, $b_8=0.003557$, $w=0.3099$, is established through the software MATLAB, which fits almost exactly with $\tilde{V}\left(s\right)$.

In the following simulations, we take ${\mathsf{\Lambda }}_h=1000$, $\mu =9.0\times {10}^{-3}$, ${\overline{a}}_1=2.0\times {10}^{-10}$, ${\overline{a}}_2=1.0\times {10}^{-10}$, ${\delta }_1=3.5\times {10}^{-5}$, ${\delta }_2=8.5\times {10}^{-5}$, $\tau =5$, $\xi =0.6$, $\alpha =8.5\times {10}^{-3}$ and $ϵ=0.01$ in Model (34). And the estimations for the periodic parameters are given by ${\beta }_1\left(t\right)={\overline{a}}_{1}f\left(ϵt\right)$, ${\beta }_2={\overline{a}}_{2}f\left(ϵt\right)$, ${\delta }_1\left(t\right)={\delta }_{1}f\left(ϵt\right)$, and ${\delta }_2\left(t\right)={\delta }_{2}f\left(ϵt\right)$. Using the above parameter values and the direct method presented in [48], we can calculate ${\mathcal{R}}_0=0.6685<1$. In this case, the HIV-positive individuals, $I\left(t\right)$, and the AIDS individuals, $A\left(t\right)$, tend to 0, which can be noticed in Figure 9a,b. The plots mean that the disease eventually goes extinct, which is exactly what Theorem 16 describes. However, if we take ${\mathsf{\Lambda }}_h=2500$, ${\overline{a}}_1=2.0\times {10}^{-8}$, ${\overline{a}}_2=1.0\times {10}^{-8}$, ${\delta }_1=3.5\times {10}^{-3}$, ${\delta }_2=8.5\times {10}^{-3}$, with other parameters values unchanged, then ${\mathcal{R}}_0=7.4711>1$. As ensured by Theorem 17, a positive $\omega$-periodic solution, $\overline{E}\left(t\right)=(\overline{S}\left(t\right),\overline{I}\left(t\right),\overline{A}\left(t\right))$, exists for Model (34) under this scenario. And the plots in Figure 9c,d illustrate that the disease will persist and exhibit periodic fluctuation eventually.

What is more, with the effect of delay $\tau$, the contact rates, $c_1$ and $c_2$, on the positive periodic solutions in the case of ${\mathcal{R}}_0>1$ are demonstrated in Figure 10. The plots in Figure 10a,b suggest that the amplitude and period of the periodic solution of the slow time model (33) vary with delay. More specifically, in an infection cycle, for the number of HIV-positive, $I\left(t\right)$, as the latency delay increases, the time of the first peak is delayed and the intensity of the outbreak decreases, while the time of the second peak is earlier and the peak increases. However, for the number of AIDS patients, $A\left(t\right)$, with the increase of delay, the time when $A\left(t\right)$ reach the peak is advanced, and the peak is reduced. The contact rates, $c_1$ and $c_2$, rates have a more significant effect on the distribution of infectious individuals than latency delay. As the contact rates decrease, although the moment of peak arrival and the peak value of $I\left(t\right)$ do not change significantly, the duration of the high run of infected cases decreases dramatically, and the peak value of $A\left(t\right)$ decreases dramatically.

8. Conclusions

In this paper, a multiscale dynamical model based on the nested modeling approach to describe the transmission of HIV-1 is proposed. The model is described by a within-host virus infection fast time model with intracellular delay and immune delay, and a between-host disease transmission slow time model with latency delay. The two models are coupled bidirectionally by enhancing the generation rate of the virus in the host, which is based on the fact that the viral load will increase by direct contact with infected individuals, and introducing between-host parameters that explicitly depend on the within-host viral load. In addition, given the intricacy of the HIV-1 transmission process, in this model, we divided infected individuals into HIV-positive individuals and confirmed AIDS patients, and introduced saturation incidence, which makes theoretical analysis more challenging but realistic. This is one of the highlights of this article. In addition, the results suggest that the immune response in the host has a great influence on HIV-1 transmission between hosts. After the virus invades the human body, if humoral immunity does not work or has a weak effect in the host, the theoretical analysis and corresponding numerical simulation results in Section 4 show that backward bifurcation and Hopf bifurcation could occur. That is to say, even if the basic reproduction number is less than 1, the disease may persist or erupt periodically between hosts, which undoubtedly poses many challenges in controlling the spread of this epidemic and makes it difficult to treat the disease effectively. When hosts have an extensive immune response and the time required between viral appearance and the production of new immune particles is short, the results obtained in Section 5 indicate that the spread of disease between hosts will be easy to control. In the presence of an immune response, the basic reproduction number, ${\tilde{\mathcal{R}}}_0^h$, can be less than 1. If it is possible to reduce the contact rate with infected individuals, increase the production rate of B cells, and prolong the latent period of infected individuals, the epidemic can be effectively controlled. However, if the time required for new immune particles to emerge in response to antigenic stimulation exceeds a certain range, which is entirely possible for the mechanism of HIV infection, then the disease persists between hosts in periodic outbreaks when ${\mathcal{R}}_0$ is greater than 1, as can be concluded from Section 6. This is another highlight of this paper and something that has not yet been discussed in existing multiscale infectious disease models. In the event of a disease outbreak, it can be more difficult to control the spread of the disease once it has sustained periodic outbreaks, compared with a situation where the disease becomes endemic over time. Therefore, taking effective therapeutic measures to maximize the inhibition of viral replication, shorten the time that antigenic stimulation needs for generating B cells, and improve immune function is the treatment goal of AIDS. Globally, more than seventy percent of people living with HIV have undetectable or suppressed viral loads as a result of antiretroviral treatment, meaning they have a zero or negligible risk of transmitting HIV to their sexual partners [49].

In order to analyze the effects of immune response and delay on HIV-1 transmission, factors such as population heterogeneity, social networks, and population mobility influencing the HIV-1 transmission are ignored in our paper. For example, population heterogeneity is an important influencing factor in epidemics, since people of different ages have different behavioral capabilities; therefore, a multiscale model bidirectional coupling within- and between-host dynamics with age structure is worth investigating. Moreover, in reality, all individuals are not homogeneously mixed and the probability of an infected individual coming into contact with others is quite different. Thus, multiscale models for characterizing the communicable disease transmission on complex networks are also to be considered in future work.

Author Contributions

Conceptualization, M.W. and L.N.; methodology, M.W., Y.W., L.H. and L.N.; software, M.W.; validation, M.W., Y.W., L.H. and L.N.; formal analysis, M.W. and L.N.; writing—original draft preparation, M.W.; writing—review and editing M.W. and L.N. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Acknowledgments

The authors are thankful to the anonymous reviewers for their careful reading, useful comments, and constructive suggestions for the improvement of the present research work. We are also thankful to the editor for his speed and cooperation.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. A flowchart for the multiscale Model (1) and Model (2).

Figure 2. The existence and stability of the immunity-inactivated infection equilibrium, [Forumla omitted. See PDF.], of Model (4) when [Forumla omitted. See PDF.]: (a) the existence of [Forumla omitted. See PDF.]; (b) the stability of [Forumla omitted. See PDF.].

Figure 3. Dynamics of the full system with different cases, where the different colored lines represent the trajectories of [Forumla omitted. See PDF.] with different initial values: (a) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], where [Forumla omitted. See PDF.] may tend to zero or to a constant value, which depends on the state of the infected individuals at the initial moment; (b) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], where [Forumla omitted. See PDF.] may tend to zero or may show cyclic oscillatory behavior and the model admits a Hopf bifurcation; (c) [Forumla omitted. See PDF.], which is similar to (a); (d) [Forumla omitted. See PDF.], which is similar to (b); (e) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], the existence and global asymptotic stability of the endemic equilibrium; (f) [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], the existence of a positive periodic solution, which implies delay induces the Hopf bifurcation.

Figure 4. Illustration of the bifurcation diagrams, where [Forumla omitted. See PDF.] is a function of [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], with the black solid curve, red solid curve, and blue dashed curve representing the locally asymptotically stable, conditionally stable, and unstable equilibrium, respectively. (a) When [Forumla omitted. See PDF.], changing [Forumla omitted. See PDF.], the saddle-node bifurcation occurs at [Forumla omitted. See PDF.]. For [Forumla omitted. See PDF.], there are two positive equilibria, where one is conditionally stable and the other is unstable; for [Forumla omitted. See PDF.], there is a unique conditionally stable positive equilibrium. (b) Changing [Forumla omitted. See PDF.], the saddle-node bifurcation occurs at [Forumla omitted. See PDF.]. The existence and stability of the positive equilibrium for [Forumla omitted. See PDF.] is the same as that of the situation [Forumla omitted. See PDF.].

Figure 5. The existence and stability of immunity-activated infection equilibrium, [Forumla omitted. See PDF.], of the fast time model (4) when [Forumla omitted. See PDF.]: (a) the existence of [Forumla omitted. See PDF.]; (b) the stability of [Forumla omitted. See PDF.].

Figure 6. The dynamics of the coupled slow time model (28) with an immune response: (a,b) the existence and stability of [Forumla omitted. See PDF.] when [Forumla omitted. See PDF.]; (c,d) the existence and stability of [Forumla omitted. See PDF.] when [Forumla omitted. See PDF.].

Figure 7. The dynamics of the fast time model (4) with [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], where the infection equilibrium, [Forumla omitted. See PDF.], is unstable and the model (4) exists in a stable periodic solution: (a) the time series charts of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]; (b) three-dimensional phase diagram of [Forumla omitted. See PDF.]; (c,d) the stability of the periodic solution.

Figure 8. Periodic solution fitting: (a) the time series of [Forumla omitted. See PDF.] of the fast time model (4); (b) the function fitting of periodic solution [Forumla omitted. See PDF.], where the curve formed by the black dot represents the time series of [Forumla omitted. See PDF.], and the red represents the fitting curve.

Figure 9. Long-term behaviors of Model (34), where the blue, red and black lines represent the trajectories of [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.], respectively. (a,b) The existence and stability of the disease-free periodic solution [Forumla omitted. See PDF.] when [Forumla omitted. See PDF.]. (c,d) The existence of the positive periodic solution [Forumla omitted. See PDF.] when [Forumla omitted. See PDF.].

Figure 10. The influence of main parameters on the infected individuals [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] for Model (33) when [Forumla omitted. See PDF.]: (a,b) the effect of the latency period [Forumla omitted. See PDF.]; (c,d) the effect of contact rate [Forumla omitted. See PDF.]; (e,f) the effect of contact rate [Forumla omitted. See PDF.].

Table 1

The biological interpretations of parameters in Model (1).

Parameter	The Biological Meanings
${\mathsf{\Lambda }}_h$	Total birth/recruitment rate of hosts
$c_1$, $c_2$	The mean number of contacts with individuals who are HIV-positive and those with AIDS, respectively
${\beta }_1\left(V\right)$, ${\beta }_2\left(V\right)$	Transmission rates of HIV-positive individuals and AIDS patients, respectively
$\mu$	Natural mortality of hosts
$\tau$	Latent period of virus in vivo
$e^{-\mu \tau }$	The surviving rate of infected individuals during latent period $\tau$
$\xi$	Transfer rate from HIV-positive stage to the AIDS patients stage
${\delta }_1\left(V\right)$, ${\delta }_2\left(V\right)$	Disease-induced mortality of HIV-positive individuals and AIDS patients, respectively

Table 2

The biological interpretations of parameters in Model (2).

Parameter	The Biological Meanings
${\mathsf{\Lambda }}_c$	Recruitment rates of healthy cells
${\kappa }_1$, ${\kappa }_2$	Infection rates of virus-to-cell transmission and cell-to-cell transmission, respectively
${\mu }_c$	Natural mortality rate of cells (healthy cells and infected cells)
${\tau }_1$	Intracellular delay (see [21] for more detail)
$e^{-{\mu }_c{\tau }_1}$	The probability of surviving the time period from $s-{\tau }_1$ to s
${\delta }_c$	Disease-induced mortality of infected cells
${\eta }_1$, ${\eta }_2$	The quantity of virus that HIV-positive individuals and AIDS patients release
	during each contact
p	The quantity of viral particles generated by each infected cell
q	The number of viral particles killed by each B cell
${\mu }_v$	Virus clearance rate within the host
$\sigma$	Production rate of new B cells stimulated by antigen
${\tau }_2$	The time required for antigenic stimulation to generate B cells
$\zeta$	Death rate of B cells

References

1. Westem, K.A. National Institute of Allergy and Infectious Disease. Reseach Conf. Natl. Inst. Allergy Infect. Dis.; 2008; 1, pp. 3-8.

2. Whiteside, A. HIV and AIDS: A Very Short Introduction; Oxford University Press: Oxford, UK, 2016.

3. HIV Statistic Globally and by WHO Region, 2023. Available online: https://www.who.int/teams/global-hiv-hepatitis-and-stis-programmes/hiv/strategic-information/hiv-data-and-statistics (accessed on 18 July 2023).

4. Mccluskey, C.C. A model of HIV/AIDS with staged progression and amelioration. Math. Biosci.; 2003; 181, pp. 1-16. [DOI: https://dx.doi.org/10.1016/S0025-5564(02)00149-9]

5. Samanta, G.P. Permanence and extinction of a nonautonomous HIV/AIDS epidemic model with distributed time delay. Nonlinear Anal. Real; 2011; 12, pp. 1163-1177. [DOI: https://dx.doi.org/10.1016/j.nonrwa.2010.09.010]

6. Perelson, A.S.; Neumann, A.U.; Markowitz, M.; Leonard, J.M.; Ho, D.D. HIV-1 dynamics in vivo: Virion clearance rate, infected cell life-span, and viral generation time. Science; 1996; 271, pp. 1582-1586. [DOI: https://dx.doi.org/10.1126/science.271.5255.1582]

7. Zhu, H.Y.; Zou, X.F. Dynamics of a HIV-1 infection model with cell-mediated immune response and intracellular delay. Discrete Cont. Dyn. B; 2009; 12, pp. 511-524. [DOI: https://dx.doi.org/10.3934/dcdsb.2009.12.511]

8. Virgin, H.W.; Walker, B.D. Immunology and the elusive AIDS vaccine. Nature; 2010; 464, pp. 224-231. [DOI: https://dx.doi.org/10.1038/nature08898]

9. Roederer, M.; Keele, B.F.; Schmidt, S.D.; Mason, R.D.; Welles, H.C.; Fischer, W.; Labranche, C. Immunological and virological mechanisms of vaccine-mediated protection against SIV and HIV. Nature; 2014; 505, pp. 502-508. [DOI: https://dx.doi.org/10.1038/nature12893]

10. Wang, S.; Zou, D. Global stability of in-host viral models with humoral immunity and intracellular delays. Appl. Math. Model.; 2012; 36, pp. 1313-1322. [DOI: https://dx.doi.org/10.1016/j.apm.2011.07.086]

11. Murase, A.; Sasaki, T.; Kajiwara, T. Stability analysis of pathogen-immune interaction dynamics. J. Math. Biol.; 2005; 51, pp. 247-267. [DOI: https://dx.doi.org/10.1007/s00285-005-0321-y]

12. Huang, G.; Takeuchi, Y.; Ma, W.; Wei, D. Global stability for delay SIR and SEIR epidemic models with nonlinear incidence rate. B. Math. Biol.; 2010; 72, 1192. [DOI: https://dx.doi.org/10.1007/s11538-009-9487-6]

13. Xu, R. Global dynamics of an SEIS epidemiological model with time delay describing a latent period. Math. Comput. Simulat.; 2012; 785, pp. 90-102. [DOI: https://dx.doi.org/10.1016/j.matcom.2012.10.004]

14. Farah, E.M.; Amine, S.; Allali, K. Dynamics of a time-delayed two-strain epidemic model with general incidence rates. Chaos Soliton. Fract.; 2021; 153, 111527. [DOI: https://dx.doi.org/10.1016/j.chaos.2021.111527]

15. Tipsri, S.; Chinviriyasit, W. The effect of time delay on the dynamics of an SEIR model with nonlinear incidence. Chaos Soliton. Fract.; 2015; 75, pp. 153-172. [DOI: https://dx.doi.org/10.1016/j.chaos.2015.02.017]

16. Jiang, Z.; Ma, W.; Wei, J. Global hopf bifurcation and permanence of a delayed SEIRS epidemic model. Math. Comput. Simulat.; 2016; 122, pp. 35-54. [DOI: https://dx.doi.org/10.1016/j.matcom.2015.11.002]

17. Culshaw, R.V.; Ruan, S.; Webb, G. A mathematical model of cell-to-cell spread of HIV-1 that includes a time delay. J. Math. Biol.; 2003; 46, pp. 425-444. [DOI: https://dx.doi.org/10.1007/s00285-002-0191-5]

18. Xu, R. Global stability of an HIV-1 infection model with saturation infection and intracellular delay. Comput. Math. Appl.; 2011; 61, pp. 2799-2805. [DOI: https://dx.doi.org/10.1016/j.camwa.2011.03.050]

19. Liu, H.; Zhang, J.F. Dynamics of two time delays differential equation model to HIV latent infection. Phys. A; 2019; 541, pp. 384-395. [DOI: https://dx.doi.org/10.1016/j.physa.2018.09.087]

20. Pawelek, K.A.; Liu, S.; Pahlevani, I.F.; Rong, L.B. A model of HIV-1 infection with two time delays: Mathematical analysis and comparison with patient data. Math. Biosci.; 2012; 235, pp. 98-105. [DOI: https://dx.doi.org/10.1016/j.mbs.2011.11.002]

21. Lin, J.; Xu, R.; Tian, X. Threshold dynamics of an HIV-1 virus model with both virus-to-cell and cell-to-cell transmissions, tntracellular delay, and humoral immunity. Appl. Math. Comput.; 2017; 315, pp. 516-530.

22. Yang, Y.; Huang, G.; Dong, Y.P. Stability and Hopf bifurcation of an HIV infection model with two time delays. Math. Biosci. Eng.; 2023; 20, pp. 1938-1959. [DOI: https://dx.doi.org/10.3934/mbe.2023089]

23. Feng, Z.L.; Velasco-Hernandez, J.; Tapia-Santos, B.; Leite, M.C.A. A model for coupling within-host and between-host dynamics in an infectious disease. Nonlinear Dynam.; 2011; 68, pp. 401-411. [DOI: https://dx.doi.org/10.1007/s11071-011-0291-0]

24. Lingappa, J.R.; Hughes, J.P.; Wang, R.S.; Baeten, J.M.; Connie, C.; Gray, G.E.; Stevens, W.S.; Deborah, D.; Campbell, M.S.; Carey, F. Estimating the impact of plasma HIV-1 RNA reductions on heterosexual HIV-1 transmission risk. PLoS ONE; 2010; 5, e12598. [DOI: https://dx.doi.org/10.1371/journal.pone.0012598] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20856886]

25. Cohen, M.; Chen, Y.Q.; Mccauley, M.; Fleming, T.R. Prevention of HIV-1 infection with early antiretroviral therapy. New Engl. J. Med.; 2011; 365, pp. 493-505. [DOI: https://dx.doi.org/10.1056/NEJMoa1105243] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21767103]

26. Gilchrist, M.A.; Sasaki, A. Modeling host-parasite coevolution: A nested spproach based on mechanistic models. J. Theor. Biol.; 2002; 218, pp. 289-308. [DOI: https://dx.doi.org/10.1006/jtbi.2002.3076] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12381431]

27. Gilchrist, M.A.; Coombs, D. Evolution of virulence: Interdependence, constraints, and selection using nested models. Theor. Popul. Biol.; 2006; 69, pp. 145-153. [DOI: https://dx.doi.org/10.1016/j.tpb.2005.07.002] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16198387]

28. Boldin, B.; Diekmann, O. Superinfections can induce evolutionarily stable coexistence of pathogens. J. Math. Biol.; 2008; 56, pp. 635-672. [DOI: https://dx.doi.org/10.1007/s00285-007-0135-1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17924106]

29. Feng, Z.L.; Velasco-Hernandez, J.; Tapia-Santos, B. A mathematical model for coupling within-host and between-host dynamics in an environmentally infectious disease. Math. Biosci.; 2013; 241, pp. 49-55. [DOI: https://dx.doi.org/10.1016/j.mbs.2012.09.004] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/23041478]

30. Cen, X.L.; Feng, Z.L.; Zhao, Y.L. Emerging disease dynamics in a model coupling within-host and between-host systems. J. Theor. Biol.; 2014; 361, pp. 141-151. [DOI: https://dx.doi.org/10.1016/j.jtbi.2014.07.030]

31. Feng, Z.L.; Cen, X.L.; Zhao, Y.L.; Velasco-Hernandez, J. Coupled within-host and between-host dynamics and evolution of virulence. Math. Biosci.; 2015; 270, pp. 204-212. [DOI: https://dx.doi.org/10.1016/j.mbs.2015.02.012]

32. Wang, X.Y.; Wang, S.P.; Wang, J.; Rong, L.B. A multiscale model of COVID-19 dynamics. Bull. Math. Biol.; 2022; 84, 99. [DOI: https://dx.doi.org/10.1007/s11538-022-01058-8]

33. Xue, Y.Y.; Xiao, Y.N. Analysis of a multiscale HIV-1 model coupling within-host viral dynamics and between-host transmission dynamics. Math. Biosci. Eng.; 2020; 17, pp. 6720-6736. [DOI: https://dx.doi.org/10.3934/mbe.2020350]

34. Wang, J.L.; Zou, X.F.; Lang, J.Y. Analysis of an age structured HIV infection model with virus-to-cell infection and cell-to-cell transmission. Nonlinear Anal.-Real; 2017; 34, pp. 75-96. [DOI: https://dx.doi.org/10.1016/j.nonrwa.2016.08.001]

35. Hale, J.K.; Lunel, S.M.V. Introduction to Functional Differential Equations; Springer: New York, NY, USA, 1993.

36. Li, X.L.; Wei, J.J. On the zeros of a fourth degree exponential polynomial with applications to a neural network model with delays. Chaos Solitons Fractals; 2005; 26, pp. 519-526. [DOI: https://dx.doi.org/10.1016/j.chaos.2005.01.019]

37. Hale, J.K. Theory of Functional Differential Equations; Springer: New York, NY, USA, 1977.

38. Liang, X.; Zhao, X.Q. Asymptotic speeds of spread and traveling waves for monotone semiflows with applications. Commun. Pur. Appl. Math.; 2007; 60, pp. 1-40. [DOI: https://dx.doi.org/10.1002/cpa.20154]

39. Hassard, B.D.; Kazarinoff, N.D.; Wan, Y.H. Theory and Applications of Hopf Bifurcation; Cambridge University Press: Cambridge, UK, 1981.

40. Song, Y.L.; Wei, J.J. Bifurcation analysis for Chen’s system with delayed feedback and its application to control of chaos. Chaos Solitons Fractals; 2004; 22, pp. 75-91. [DOI: https://dx.doi.org/10.1016/j.chaos.2003.12.075]

41. Zhao, X.Q. Basic reproduction ratios for periodic compartmental models with time delay. J. Dyn. Differ. Equ.; 2017; 29, pp. 67-82. [DOI: https://dx.doi.org/10.1007/s10884-015-9425-2]

42. Smith, H.L. Monotone Dynamical Systems: An Introduction To the Theory of Competitive and Cooperative Systems; American Mathematical Society: Providence, RI, USA, 1995.

43. Zhao, X.Q. Dynamical Systems in Population Biology; Springer: Cham, Switzerland, 2003.

44. Wang, X.; Zhao, X.Q. A periodic vector-bias malaria model with incubation period. SIAM. J. Appl. Math.; 2017; 77, pp. 181-201. [DOI: https://dx.doi.org/10.1137/15M1046277]

45. Wang, X.; Zhao, X.Q. Dynamics of a time-delayed Lyme disease model with seasonality. SIAM. J. Appl. Dyn. Syst.; 2017; 16, pp. 853-881. [DOI: https://dx.doi.org/10.1137/16M1087916]

46. Hale, J.K. Asymptotic Behavior of Dissipative Systems; American Mathematical Society: Providence, RI, USA, 1988.

47. Nowak, M.A.; Bangham, C.R. Population dynamics of immune responses to persistent viruses. Science; 1996; 272, pp. 74-79. [DOI: https://dx.doi.org/10.1126/science.272.5258.74]

48. Yang, T.H.; Zhang, L. Remarks on basic reproduction ratios for periodic abstract functional differential equations. Discrete Cont. Dyn. B; 2019; 24, pp. 6771-6782. [DOI: https://dx.doi.org/10.3934/dcdsb.2019166]

49. World Health Organization. Available online: https://www.who.int/director-general/speeches/detail/who-director-general-s-opening-remarks-at-the-media-briefing---26-july-2023 (accessed on 26 July 2023).