Two deadly viruses for the human body are the human immunodeficiency virus (HIV) and hepatitis B virus (HBV). In several patients, co-infection with these two viruses occurs because they share common transmission routes, such as blood transfusions, use of contaminated syringes, unprotected sexual intercourse, and others [1]. It is estimated that 8% to 10% of people diagnosed with HIV are also infected with HBV worldwide, making it important to study the co-infection of these two viruses. The co-infection with both viruses may increase the symptoms of each disease. This is due to the weakening and deterioration of the patient's immune system and the difficulty of determining effective treatments for both viruses [2, 3]. HIV is a single-stranded RNA virus that primarily infects CD4⁺ T cells, which are key components of the adaptive immune system. Over time, if left untreated, the number of CD4⁺ T cells declines, eventually leading to a disease called acquired immunodeficiency syndrome (AIDS) [4]. AIDS infection exposes the patient to opportunistic diseases and cancers. Besides infecting CD4⁺ T cells, HIV is also capable of invading other immune cells like dendritic cells and macrophages, resulting in their reduction and impaired function. The HIV also has the ability to infect different liver cells, including hepatocytes, Kupffer cells, and infiltrating T cells [5, 6]. HBV is a double-stranded DNA virus that primarily targets liver cells (hepatocytes) [7]. Over time, if left untreated, this infection can lead to chronic hepatitis, which may progress to cirrhosis and hepatocellular carcinoma.

Studying the interactions between viruses and target cells within the host, as well as immune responses, through laboratory experiments is not only difficult but also expensive. Therefore, mathematical modeling and simulation have become essential and important tools for understanding the dynamics of viruses within the host, as well as the role of the immune system in combating them. Furthermore, mathematical models may help design appropriate therapeutic strategies. In [8, 9], the basic models of HIV and hepatitis B virus (HBV) infection, respectively, have been formulated. These models consist of three main compartments: uninfected target cells (CD4⁺ T cells for HIV and hepatocytes for HBV), infected cells, and free viral particles. To take into account different biological phenomena not covered by these basic models, various extensions and generalizations have been made. One of the most important of these phenomena is the time delay during biological reactions. Examples of time delays include the time delay until a cell becomes infected, the time required for activation of the virus released from the cells, and the delay in the immune response. Time delays contribute to the development of accurate models that help understand the evolution of the virus within the host. Numerous studies have explored the incorporation of time delays into HIV single-infection models (see, e.g., [10–16]). On the other hand, HBV single-infection models with delay have also been formulated in some studies (e.g., [17–20]).

Mathematical modeling of HIV and hepatitis B virus co-infection at the population level has received considerable attention from researchers. These models play an important role in understanding patterns of infectious disease transmission within populations. These models may also play a significant role in formulating epidemic control strategies. Therefore, several mathematical models have been developed for the coinfection of HBV and HIV (see, e.g., [3, 21–26]). In contrast, within-host models that examine the biological interactions and co-dynamics of HIV and HBV at the cellular level are relatively scarce. To our knowledge, only a limited number of studies, notably [27, 28], have addressed the within-host interplay between these two viruses. This highlights a significant gap in the literature, underscoring the need for more comprehensive investigations into their intracellular co-infection dynamics.

The HIV-HBV co-infection model proposed by [28] is formulated as follows:

At time t, the variables x(t), y₁(t), y₂(t), v₁(t), and v₂(t) denotes the concentrations of uninfected hepatocytes, hepatocytes infected with HIV, hepatocytes infected with HBV, free HIV virions, and free HBV virions, respectively. τ₁ and τ₂ denote the respective time delays associated with the development of HIV-infected and HBV-infected hepatocytes. The factor e-niτi, i = 1, 2 represents the probability of survival of infected hepatocytes throughout the delay time [t − τ_i, t], and n_i > 0.

Notably, the model presented in [28] has certain limitations:

(i) it overlooks the dynamic evolution of both uninfected and HIV-infected CD4⁺ T cell populations, opting instead to assume a fixed rate of HIV production from infected CD4⁺ T cells, thereby ignoring their temporal variation;

(ii) it does not account for the maturation delay of newly produced virions;

(iii) Global stability was established exclusively for the infection-free equilibrium.

Limitation (i) was partially resolved in the model proposed by [27], which incorporated the dynamics of CD4⁺ T lymphocyte populations. However, the mathematical analysis in that study was largely confined to establishing the boundedness and positivity of the solutions. In addition, only the infection-free equilibrium was identified and analyzed for local stability, while a comprehensive investigation of other possible equilibria was not pursued.

This study builds upon and extends the models presented in [27, 28] by addressing key limitations. We explicitly model the dynamics of both uninfected and HIV-infected CD4⁺ T lymphocytes, replace the assumption of a constant HIV production rate with a more realistic approach, and include a biologically relevant time delay to reflect delays in viral production and virion maturation. A rigorous mathematical analysis, including equilibrium identification and global stability examination, provides deeper insights into the system's long-term dynamics. These findings are supported by numerical simulations and sensitivity analyses that highlight the influence of critical parameters.

The main contributions of this study are:

In this section formulate our proposed model on based on the next Hypotheses:

Based on Hypotheses (H1)-(H8), we formulate an HIV-HBV co-dynamics model with delays as follows:

The model incorporates five distinct time delays:

The initial conditions for model (Equations 1–7) are given as follows:

where τm=max{τ1,τ2,⋯,τ5},ϕi∈C([-τm,0],ℝ≥0), and C:[-τm,0]→ℝ≥0 is the Banach space of continuous functions. The norm is defined as ∥ϕi∥=sup-τm≤θ≤0|ϕi(θ)| for ϕ_i ∈ C, and i = 1, 2, ⋯ , 7. The standard theory of functional differential equations shows that, model (Equations 1–7) with initial conditions (Equation 2) has a unique solution [30, 31].

The model's variables and parameters are defined in Table 1.

Let Fi=e-niτi for i = 1, 2, ⋯ , 5, σ₁ = min{d₁, a₁, a₂} and σ₂ = min{d₂, a₃}.

Lemma 1. The dynamics (Equations 1–7) with the initial states (Equation 2) admits nonnegative and ultimately bounded solutions.

Proof. Clearly, Equations 1, 4 of model (Equations 1–7) give ẋ_i∣_{xi = 0} = λ_i > 0, which implies that x_i(t) > 0 for i = 1, 2 and all t ≥ 0. Moreover, for t ∈ [0, τ^m] we have

which gives (x₁, y₁, y₂, x₂, z₁, v₁, v₂)(t) ≥ 0. Hence, by recursive argumentation, we obtain that (x₁, y₁, y₂, x₂, z₁, v₁, v₂)(t) ≥ 0 for any t ≥ 0.

Next we prove the ultimately boundedness of the model's solution (x₁, y₁, y₂, x₂, z₁, v₁, v₂). From Equation 1, we have lim sup_{t → ∞} x₁(t) ≤ Γ₁, where Γ1=λ1d1. Concerning, y₁(t) and y₂(t), we define

Then, we get

Thus lim sup_{t → ∞} ψ₁(t) ≤ Γ₂, and hence lim sup_{t → ∞} y₁(t) ≤ Γ₂ and lim sup_{t → ∞} y₂(t) ≤ Γ₂, where Γ2=λ1σ1. From Equation 4, we have lim sup_{t → ∞} x₂(t) ≤ Γ₃, where Γ3=λ2d2. Define

Then, we get

This gives lim sup_{t → ∞} ψ₂(t) ≤ Γ₄, and hence lim sup_{t → ∞} z₁(t) ≤ Γ₄ where Γ4=λ2σ2.

Note that

Consequently, lim sup_{t → ∞} v₁(t) ≤ Γ₅, where Γ5=1c1[k1a1Γ2+k3a3Γ4].

Similarly, we have

This yields lim sup_{t → ∞} v₂(t) ≤ Γ₆, where Γ6=k2a2Γ2c2.

Based on Lemma 1, one can show

is positively invariant w.r.t. system (Equation 1–7).

The global asymptotic stability of all equilibrium points in the system (Equations 1–7) is investigated using Lyapunov's direct method in this section. The development of Lyapunov functions is based on the approach outlined in [33]. Let F represent the candidate Lyapunov function, and let Ω′ denote the largest invariant subset of Ω={(x1,y1,y2,x2,z1,v1,v2):dFdt=0}. Consider the function Φ(X) = X − ln X − 1. Let us denote by (x₁, y₁, y₂, x₂, z₁, v₁, v₂) = (x₁, y₁, y₂, x₂, z₁, v₁, v₂)(t), and (x1τ,y1τ,y2,x2τ,z1τ,v1τ,v2τ)=(x1,y1,y2,x2,z1,v1,v2)(t-τ).

Theorem 1. The infection-free equilibrium E0 is globally asymptotically stable (GAS) when RC≤1.

Proof. Consider F0(x1,y1,y2,x2,z1,v1,v2)