The human immunodeficiency virus (HIV), a fast-replicating retrovirus within the lentivirus genus, specifically targets and disrupts the functionality of critical immune cells, with a strong preference for CD4⁺T cells. These immune cells play a vital role in regulating and coordinating immune defense mechanisms. Once HIV enters the body, it integrates its genetic material into the host genome, gradually depleting the number of CD4⁺T cells. In a healthy individual, the CD4⁺T cell count is typically around 1,000 cells per mm³ of blood. As the infection progresses, these levels decline steadily, often without immediate symptoms. When the count falls below 200 cells/mm³, the individual is considered to have developed acquired immunodeficiency syndrome (AIDS), a condition characterized by increased susceptibility to opportunistic infections and certain types of cancer [1, 2]. In the absence of appropriate antiretroviral therapy, the continuous weakening of the immune system can lead to severe complications and, ultimately, death. The adaptive immune response is vital in managing viral infections and operates primarily through two pathways: antibodies produced by B cells that neutralize HIV particles, and cytotoxic T lymphocytes (CTLs) that identify and destroy infected host cells. These coordinated responses work to inhibit viral replication and reduce the rate of disease advancement.

Mathematical modeling serves as an essential approach for exploring the complex interactions between viral infections and the host's immune defenses. By representing biological processes through mathematical equations, these models enable a deeper investigation into the factors that influence infection dynamics, the strength and limitations of immune responses, and the potential outcomes of various scenarios. This approach has greatly enhanced our comprehension of virus-immune system relationships and has become instrumental in guiding the design of effective treatments and shaping public health policies. A basic model for analyzing HIV dynamics was proposed by Nowak and Bangham [3], capturing the interactions between healthy CD4⁺T cells, infected cells, and free virus particles. Since its introduction, this foundational framework has been extended in various directions to incorporate the CTL response [4–9], the antibody-mediated immune response [10, 11], and more comprehensive formulations that include both CTL and antibody responses simultaneously [12–16].

Pyroptosis, unlike the controlled and non-inflammatory nature of apoptosis, is a pro-inflammatory programmed cell death mechanism closely linked to immune system activation and inflammatory responses. In HIV-infected individuals, Doitsh et al. [17] demonstrated that caspase-1, a cysteine-dependent protease, becomes activated and promotes the secretion of inflammatory mediators such as interleukin-1β (IL-1β). These cytokines sustain a chronic inflammatory environment that draws in uninfected CD4⁺T cells, rendering them vulnerable to death. As a result, a destructive feedback loop is created in which infected CD4⁺T cells undergoing pyroptosis release inflammatory factors that drive the death of surrounding uninfected CD4⁺T cells, accelerating the decline of immune function [17]. Studies indicate that merely about 5% of CD4⁺T cell death is attributed to apoptosis triggered by caspase-3, whereas the majority are lost through pyroptosis, driven by caspase-1 activation [18].

Wang et al. [19] developed a mathematical model of HIV infection that incorporates the contribution of pyroptosis to CD4⁺T cell depletion. The model is given by a system of ordinary differential equations (ODEs) which describes the interaction of healthy CD4⁺T cells, productively infected cells, abortively infected cells, inflammatory cytokines, and free HIV particles. In this framework, the impact of inflammatory cytokines on the basic reproduction number is not considered. As a result, the model overlooks the potential rise in infection rates driven by cytokine-induced recruitment of CD4⁺T cells to inflamed areas, which could lead to an underestimation of the basic reproduction number.

Recent research has highlighted the role of cytokine-driven mechanisms in the accumulation of infected CD4⁺T cells and their impact on HIV dynamics [20]. Failing to account for the rise in viral infection caused by the higher concentration of CD4⁺T cells, which are attracted to inflamed regions by cytokines, could lead to an underestimation of the basic reproduction number [21]. Cytokine-induced viral infections can disturb the balance between cell renewal and viral replication, particularly through mechanisms such as pyroptosis [18]. The growth of the infected CD4⁺T cell population is shaped by two main pathways: (i) direct infection of healthy CD4⁺T cells through viral contact, and (ii) an increase in infection rates linked to the elevated presence of CD4⁺T cells in inflamed tissues, where cytokines actively recruit them (a process known as cytokine-enhanced viral infection). To capture the combined effects of these mechanisms, the following system of equations models the dynamics of HIV infection under cytokine influence [20]:

Here, N(t), U(t), S(t) and B(t) represent the time-dependent concentrations of healthy CD4⁺T cells, CD4⁺T cells infected with HIV, inflammatory cytokines, and free HIV virions, respectively. The parameter δ represents the production rate of uninfected CD4⁺T cells. The term ω₁NB corresponds to the rate at which CD4⁺T cells become infected through direct interaction with the virus, while ω₂NS captures the cytokine-induced enhancement of infection. Cytokine and viral particle release from infected cells are modeled by the terms β₂U and μU, respectively. The expression φ_ηη denotes the natural death rate of a given compartment η, and β₁U specifically accounts for the loss of infected CD4⁺T cells through pyroptotic cell death. Various versions of this model have been developed to incorporate key biological aspects, including (i) time delays [8, 20, 22], (ii) cell-to-cell transmission [23–25], (iii) CTL response [26, 27], (iv) antibody response [28, 29], and (v) the combined effects of CTL and antibody responses [25, 30].

Model (Equation 1) and its aforementioned extension presume a homogeneous distribution of both cells and viruses, thereby ignoring their spatial movement and localized interactions. While this assumption simplifies analysis, it fails to account for the critical influence of spatial organization on viral behavior. In actual biological systems, the uneven distribution of infected and healthy cells, along with constraints on viral mobility, often creates spatial variability. This heterogeneity can significantly shape how infections emerge, spread, and are contained, influencing immune response and treatment outcomes. Phenomena such as clustering of infections, formation of infection hotspots, and restricted viral spread are examples of effects that uniform mixing models cannot represent. Therefore, integrating spatial dynamics into modeling frameworks offers a more realistic and nuanced understanding of viral infections and their progression. Brainard et al. [31] and Tattermusch and Bangham [32] highlighted that T cells are capable of moving along concentration gradients, typically from regions of higher density to lower density. More recent research has extended this idea, proposing that immune cells, infected cells, and viral particles may also undergo migration from areas of abundance to those with lower presence [see, for example, [33–41]].

Wang and Zhang [42] proposed a nonlocal partial differential equations (PDEs) model incorporating time delays to analyze the role of pyroptosis in infection dynamics. Their framework includes both latent and actively infected cell populations. The infection process between healthy cells and the virus is modeled using a saturating function. Similarly, the influence of inflammatory cytokines on healthy CD4⁺T cells via pyroptosis is described by saturating function. In Wang et al. [43], the authors formulated a reaction-diffusion model that includes time-periodic parameters, spatial heterogeneity, and a latent infection stage to explore the role of pyroptosis in CD4⁺T cell loss during HIV infection. The model employs a Beddington-DeAngelis type functional response to characterize the formation rate of newly infected cells. The analysis focuses on threshold behavior governed by the basic reproduction number, providing insights into conditions for disease persistence or clearance. Wang et al. [44] introduced a partial differential equation model to investigate the effects of gasdermin D inhibition on pyroptosis in settings characterized by spatial and temporal heterogeneity. The framework includes both productively and abortively infected cell populations and accounts for the role of CTL response. The interaction between the virus and host cells is described using a general incidence function. Wang et al. [45] developed a periodic partial differential equation model to explore the dynamics of infected cell production, incorporating the heightened infection risk caused by cytokine-driven T cell migration to inflamed tissues. Their model also included the effect of necrosulfonamide, a pharmacological agent that inhibits pyroptosis. Wang and Feng [21] developed a partial differential equation model that incorporates spatial heterogeneity by introducing general functional forms to represent the dynamics of healthy CD4⁺T cell regeneration, virus-cell interactions, and cytokine-enhanced viral infection. However, in Wang and Feng [21], while the spatial dynamics of inflammatory cytokines and the virus are included, the movement of both healthy and infected CD4⁺T cells is omitted. Additionally, the immune response is not taken into account. In Chen et al. [46], the authors introduced a reaction-diffusion framework to study HIV dynamics, employing a general incidence function while accounting for cytokine-mediated enhancement of infection, CTL response, and the influence of time delay.

The models developed in Chen et al. [25] and Dahy et al. [30] assume that the presence of HIV and infected cells solely triggers CTL and antibody responses, overlooking the potential for immune suppression, commonly referred to as immune impairment. As noted in Lydyard et al. [47], HIV can weaken the immune system's functionality. Several studies have incorporated immune impairment into viral dynamics, focusing either on CTL dysfunction [e.g., [48–54]] or B-cell impairment [e.g., [54–57]]. More recently, AlShamrani et al. [58, 59] examined HIV dynamics by considering the impairment of both CTL and antibody responses. However, these studies did not consider the role of inflammatory cytokines. Although Song et al. [60] introduced CTL impairment alongside cytokine effects, it did not account for target cell population dynamics or the spatial mobility of viruses, infected cells and immune cells.

Accordingly, the objective of this work is to construct a comprehensive HIV infection model that captures several critical aspects of within-host viral dynamics. Specifically, the model aims to incorporate:

By integrating these features, the proposed model seeks to provide deeper insights into HIV pathogenesis and to support the development of more effective therapeutic strategies.

We establish the well-posedness of the model by proving the existence and boundedness of global solutions, identify the equilibria and derive the basic reproduction number, and analyze the global asymptotic stability of the equilibria using appropriate Lyapunov functions and LaSalle's invariance principle. The theoretical findings are supported by detailed numerical simulations, and a sensitivity analysis of the basic reproduction number is conducted to assess the influence of key parameters.

We formulate a reaction-diffusion model based on partial differential equations to describe the variation of the concentrations of six compartments with respect to spatial location ϰ and time t; healthy CD4⁺T cells N(ϰ, t), HIV-infected CD4⁺T cells U(ϰ, t), inflammatory cytokines S(ϰ, t), free HIV particles B(ϰ, t), Cytotoxic T Lymphocytes (CTLs) M(ϰ, t), and antibodies H(ϰ, t).

where ϰ = (ϰ₁, ϰ₂, ..., ϰ_ℓ) ∈ Λ, and t > 0. The spatial region Λ ⊂ ℝ^ℓ, is a connected, bounded domain, and has a smooth boundary ∂Λ, where ℓ is an integer such that ℓ ≥ 1. The diffusion coefficient DG is positive for each G in the set {N, U, S, B, M, H}. The Laplace operator, represented by Δ, is defined as ∂2∂ϰ2. The infection rate resulting from cellular infection is ω₃NU. The terms ξU and ψB represent the rates at which CTLs and antibodies proliferate, respectively, from infected cells and free HIV particles. The rate at which infected cells are killed by CTLs is denoted as ϱ₁UM, whereas free HIV particles are neutralized by antibodies at rate ϱ₂BH. The rates at which CTL and antibody immunities are impaired are labeled as λ₁UM and λ₂BH, respectively. The natural death rates associated with CTLs and antibodies are represented by φ_M and φ_H, respectively. As previously outlined in Section 1, all other parameters share an identical biological interpretation.

In the following, the initial conditions, as well as the homogeneous Neumann boundary conditions adopted for system (Equation 2), are given as:

Here, the functions Ij(ϰ), for j = 1, ..., 6, are both continuous and non-negative. Meanwhile, ∂∂Z→ represents the outward normal derivative on the boundary ∂Λ. The boundary conditions (Equation 4) ensure that all populations are prohibited from traversing the isolated boundary ∂Λ [61].

Remark 1. The mathematical model developed in Hattaf [6] incorporates CTL immune responses and considers the spatial mobility of both immune cells and viruses. However, it omits several key immunological factors, such as the contribution of antibody-mediated immunity, the regulatory role of pro-inflammatory cytokines, and the potential dysfunction or exhaustion of CTL response that may occur during chronic infection. Conversely, the study by Hajhouji et al. [62] focuses on the antibody response in the context of HIV dynamics but does not include CTL-mediated immunity, cytokine-driven inflammation, or the spatial migration of immune cells and viral particles. Together, these limitations highlight the need for more comprehensive models that integrate multiple layers of the immune system and spatial effects to better capture the complexity of HIV pathogenesis. These limitations have been addressed in the present model, which incorporates both CTL and antibody-mediated immune responses, accounts for the effects of pro-inflammatory cytokines, and includes the mobility of immune cells and viruses. By integrating these critical biological factors, our model offers a more comprehensive framework for analyzing HIV infection dynamics and immune system interactions.

The following result addresses the existence, uniqueness, non-negativity, and boundedness of solutions for model (Equation 2), which describe the densities of healthy CD4⁺T cells, HIV-infected CD4⁺T cells, inflammatory cytokines, free HIV particles, CTLs, and antibodies.

Lemma 1. Let the assumption DN=DU=DS=DB=DM=DH=D^ hold true. For any initial function I=(I1,I2,I3,I4,I5,I6)T∈Y+ satisfying the initial conditions (Equation 3), model (Equation 2) has a unique, non-negative and bounded solution (N(ϰ, t), U(ϰ, t), S(ϰ, t), B(ϰ, t), M(ϰ, t), H(ϰ, t)) defined on Λ̄×[0,+∞).

Proof. Let Y=BUC(Λ̄,ℝ6) be defined as the set of all functions from Λ̄ to ℝ⁶ that are both bounded and uniformly continuous, with the norm ∥ϕ∥Y=supϰ∈Λ̄|ϕ(ϰ)|. We denote the positive cone Y+=BUC(Λ̄,ℝ+6)⊂Y which establishes a partial order on Y. This characterization demonstrates that the space (Y,∥·∥Y) forms a Banach lattice [63, 64].

Concerning every initial data I=(I1,I2,I3,I4,I5,I6)T∈Y+, we define K=(K1,K2,K3,K4,K5,K6)T:Y+→Y as follows:

We observe that K is locally Lipschitz on Y+, which is a fact straightforward to verify (see Corollary 4 in Martin and Smith [65]). System (Equation 2) subject to initial conditions (Equation 3) and boundary conditions (Equation 4) can be reformulated as the following abstract functional differential equation:

where L=(N,U,S,B,M,H)T and ΓL=(DNΔN,DUΔU,DSΔS,DBΔB,DMΔM,DHΔH)T. One can prove that

Follows the work of Xu and Xu [63], Zhang and Xu [64], and Smith [66], we deduce that for any I∈Y+, system (Equation 2) subject to (Equations 3, 4) possesses a unique non-negative mild solution (N(ϰ, t), U(ϰ, t), S(ϰ, t), B(ϰ, t), M(ϰ, t), H(ϰ, t)). This solution is defined on Λ̄×[0,TM), where [0,TM) is the maximal time interval over which the solution remains in existence. Furthermore, this solution constitutes a classical solution to the problem at hand.

To confirm that the solutions have a bounded nature, we introduce

Based on the fact that DN=DU=DS=DB=DM=DH=D^, using system (Equation 2) we derive

where ϖ=min{φN,β1+φU8,φS,φB2,φM,φH}. Consequently, P(ϰ,t) fulfills the subsequent system

Consider a solution, P~(t), for the following ordinary differential equation:

This results in P~(t)≤max{δϖ,maxϰ∈Λ̄P(ϰ,0)}. With reference to the comparison principle (refer to Protter and Weinberger [67]), we find that P(ϰ,t)≤P~(t). From this, we derive

which demonstrates that N(ϰ, t), U(ϰ, t), S(ϰ, t), B(ϰ, t), M(ϰ, t), and H(ϰ, t) are bounded on Λ̄×[0,TM). According to the standard theory for semi-linear parabolic systems, we infer that TM=+∞ [68]. The solution (N(ϰ, t), U(ϰ, t), S(ϰ, t), B(ϰ, t), M(ϰ, t), H(ϰ, t)) exists and is uniquely determined and non-negative for all ϰ ∈ Λ, t > 0.

In this section, we assess the equilibria and identify the threshold parameter necessary to confirm their existence. The results are outlined in the subsequent lemma:

Lemma 2. Considering system (Equation 2), a basic reproduction number

can be identified, which fulfills the following statements:

Proof. The basic reproduction number, R0, is computed through the next-generation matrix technique described in van den Driessche and Watmough [69]. To accomplish this, we can represent the right-hand side of system (Equation 2) as J1-J2 with

System (Equation 2) consistently exhibits an HIV-free equilibrium FE=(N0,0,0,0,0,0), where N0=δφN.

Upon computing the Jacobian matrices, J1 and J2, at the equilibrium FE, we find

Note that, the next generation matrix is in the following form:

The basic reproduction number R0 is determined by the spectral radius of the matrix product J1J2-1, expressed as:

where

To clarify, the contributions of viral and cellular infections are represented, respectively, by R0B and R0U, whereas R0S denotes the influence of inflammatory cytokines.

To identify the additional equilibrium beyond FE, we assume (N, U, S, B, M, H) represents any equilibrium that fulfills the following equations:

Referring to Equations 10, 11, we derive

Replacing the values from Equation 12 in Equation 9, we obtain

By substituting the expression from Equation 13 into Equation 12, we yield

Replacing the values from Equation 13 in Equation 8 gives

From Equations 6, 7, we get

Substituting from Equations 13, 14 into Equation 16, we get

Substituting from Equations 13–15, 17 into Equation 7, we get

where