In-host modeling of Tfh–Tfr–fCD8+ T cell dynamics on HIV reservoir persistence

1 Introduction

Lymph nodes are a primary site for adaptive immune responses and contain regions that harbor T and B cells, enabling antigen-presenting cells to initiate immune responses [1]. In chronic, untreated HIV-1 infection, viral replication is concentrated in secondary lymphoid follicles [2]. The current Antiretroviral Therapy (ART) effectively suppresses HIV, but some CD4 T cells still harbor the virus in a quiescent state, allowing it to re-establish itself indefinitely. The situation can be dire, with suboptimal adherence or without ART. Thus, some T cells can serve as a reservoir for HIV in various sites, such as the lymph nodes and the central nervous system. T follicular helper cells (Tfh) and T follicular regulatory cells (Tfr) are both subsets of T cells in the secondary lymphoid follicles, a significant site of HIV replication [2–4]. They have a high affinity for HIV in the germinal centers (GCs) of the B-cell follicles and act as HIV reservoirs, which can lead to resurgence of HIV in the sub-optimal use of ART [5]. They are more susceptible to HIV than extrafollicular subsets of CD4⁺ T cells [6]. Although primarily antibody production centers, germinal centers constitute the largest reservoirs of HIV virions [1, 7]. The production of CD4+ T cells from naive CD4+ T cells primed by dendritic cells, and further differentiation yields Tfh cells that migrate to and remain resident in lymph nodes (LNs) and the spleen, where they help B cells. Their primary role is to protect against infectious diseases through priming the GC responses to infection [8]. The functions of Tfh cells must be regulated to prevent immunodeficiency and antibody-mediated autoimmune diseases. Tfr cells primarily prevent or eliminate autoreactive B cells generated during antigen-specific immune responses [8, 9]. The Tfr cells have been shown to directly suppress Tfh cells and to expand in the presence of HIV. Tfr cells are derived from conventional regulatory T cells [Foxp3+ Treg precursor cells [10]], extrafollicular regions, and subsequently migrate to and localize within the follicles [2]. Their regulatory mechanisms may impair Tfh cell function, and they are more permissive to HIV infection than Tfh cells [6].

Virus-specific follicular CD8 T (fCD8) cells play a significant role in controlling HIV infection by killing infected cells. They are present at significantly lower levels in the GCs of B follicles [4] because intact follicular and GC barriers keep fCD8 cells outside [1]. Their low levels are also attributed to a lack of CXCR5 expression, a follicular homing receptor. Although fCD8 cells, are effective at controlling HIV replication, they cannot suppress it entirely, as they are less cytotoxic than their extrafollicular counterparts [6]. The follicles are also highly permissive to free-circulating HIV, which finds a home there, attaches to follicular dendritic cells, and creates an extracellular reservoir of HIV [6].

The persistence of HIV, particularly in immune-privileged sites such as lymphoid tissues, has been a major obstacle to a complete cure. Tfh cells are preferentially infected and serve as a key reservoir for latent HIV infection. Tfh, Tfr, and fCD8 contribute to HIV persistence and impaired humoral immunity. Understanding their roles can facilitate the development of vaccines and strategies for an HIV cure. No therapies have been developed to target T follicular helper cells and directly address what is thought to be a significant source of HIV persistence [4].

Mathematical models have been used to represent HIV infection in vivo, and the analysis is therefore valuable for better understanding the prognosis of HIV infection. Some studies have focused extensively on extra-follicular HIV infection dynamics [11–13]. Still, recent developments point to follicular HIV infection as a source of HIV persistence with or without treatment. Popular treatment strategies involve the use of reverse transcriptase inhibitors (RTIs) and protease inhibitors (PIs), which have led to the suppression of HIV infections. Because HIV is harbored in the follicles, current strategies are inadequate, leaving room for recent discoveries of potential immunotherapies targeting Tfh and free virions that infect the machines under development. This calls for studies to better understand the immunopathogenesis of HIV in follicles. The study by Erwin and Cupe in 2017 [14] introduced a mathematical model of germinal center dynamics during acute and chronic HIV infection, focusing on B cells and Tfh cells to determine how B-cell selection and competition for T follicular helper cells and antigen influence the size and composition of germinal centers in both acute and chronic infections [14]. Their model captured both mutating and non-mutating antigens, as well as competition for Tfh cell stimulation between the B-cell clones. Their results showed that T follicular helper cells were limiting the emergence of high B-cell somatic hypermutation rounds for non-mutating and mutating viruses. Although the study is mechanistically coherent and biologically motivated, providing a solid foundation for reservoir-integrated models, the antigen is modeled phenomenologically rather than with equations for infected cells and viruses. Thus, there is no feedback between viral load and immune pressure, nor is there a reservoir structure (such as the Tfh infection), which limits the translational applicability to HIV cure modeling. The Tfr-mediated suppression of Tfh is not included, thereby omitting a key regulatory axis that influences the dynamics of HIV GC.

In this study, we develop and analyze a mechanistic in-host mathematical model of germinal center dynamics that quantitatively describes how the coupled interactions between T follicular helper (Tfh), Tfr, and follicular CD8+ (fCD8+) T cells control HIV replication, immune control, and the creation, upkeep, and stability of the viral reservoir.

2 Model formulation

We develop a modeling framework by introducing the following variables: Healthy T follicular helper cells, T_*fr*, infectious T follicular cells, I_*fh*, healthy T follicular regulatory cells, T_*fr*, infectious T follicular regulatory cells cells, I_*fr*, the HIV population, V, and the follicular CD8⁺ T cells, C_*f*.

The healthy T follicular helper cells are produced at a rate λ_*h* and die naturally at a rate μ_*h*. Free virions infect them at a rate β_*h*. They undergo clonal expansion [15] at a rate r_*h*, a process suppressed by the T follicular regulatory cells with a saturation constant α_*r*. They also become anergic due to the regulation by T follicular regulatory cells [16] at a rate α_*d*. The infected T follicular helper cells arise from the infection of healthy T follicular helper cells, and they die naturally at a rate μ_*h*. They also die due to the virus budding off or lysis at a rate δ_*h*. They are killed at a rate k_*h* by the follicular CD8⁺ T cells [4]. Healthy T follicular regulatory cells are produced at a rate λ_*r*. The study in Zhao et al. [17] and Lee et al. [18] showed that the T_*fr* cells expanded during HIV infection, which is dominated by the T_*fh* infection, which is the reservoir and main replication site of HIV in the GCs. Hence, we assume they expand due to the presence of the infected follicular helper cells [17], at a rate ρ_*r*. Here, ρ_*r* is the number of T_*fr* cells produced per unit time due to the generation of each infected T_*fh* cell. This growth term is also roughly proportional to the free virus, since infected T follicular helper cells are the sites of viral production. The T follicular regulatory cells die naturally at a rate μ_*r*. They move out through T_*fh*−induced instability that makes active T follicular helper cells transition into ex-T follicular regulator cells [18] at a rate σ_*r*. The infected T follicular regulatory cells are produced from the infection of healthy T follicular regulatory cells that die naturally at the rate μ_*r*. They also die due to the virus budding off or lysis at a rate δ_*r*. They are killed at a rate k_*r* by the follicular CD8⁺ T cells. The viral particles are produced by infected T follicular helper cells, with each actively infected cell producing an average of N_*h* viruses. They are also produced from infected T follicular regulatory cells, each actively infected cell producing N_*r* viruses on average [19, 20]. This production term excludes both viral influx into the germinal center (GC) from peripheral blood and viral efflux from the GC [3]. It allows the study to focus on the local immunological interactions within the GC microenvironment, where a significant fraction of HIV replication and immune dysregulation occurs. Similar approaches have been widely used in other HIV modeling studies that do not focus on a full systemic HIV model [11]. The virus dies naturally at a rate μ_*v*. The follicular CD8⁺ T cells grow in response to infected T follicular regulatory cells [21] at a rate ρ_*c*. They die naturally at a rate, μ_*c*, and their growth is suppressed by the healthy T follicular regulatory cells [22, 23].

The coupled T_*fh* − T_*fr* − fCD8 HIV infection dynamics is given by the following systems of Equations 1–6:

dTfhdt=λh-μhTfh-βhTfhV+rhTfhαrαr+Tfr-αdTfhTfr,(1)

dIfhdt=βhTfhV-(μh+δh)Ifh-khCfIfh,(2)

dTfrdt=λr+ρrIfh-βrTfrV-μrTfr-σrTfrTfh,(3)

dIfrdt=βrTfrV-(μr+δr)Ifr-krCfIfr,(4)

dVdt=NhδhIfh+NrδrIfr-μvV,(5)

dCfdt=ρcIfh-μcCf-δcTfrCf,(6)

subject to initial conditions T_*fh*(0) = T_*fh*0 > 0, I_*fh*(0) = I_*fh*0 ≥ 0, T_*fr*(0) = T_*fr*0 ≥ 0, I_*fr*(0) = I_*fr*0 ≥ 0, V(0) = V₀ ≥ 0, C_*f*(0) = C_*f*0 ≥ 0. In the following, we explain the formulation of each equation in the system (Equations 1–6):

3 Mathematical analysis

3.1 Positivity and boundedness

For the model system (Equations 1–6) to be biologically feasible, we require that the cell densities be non-negative and bounded. Hence, we state the following proposition.

Proposition 1. *For any initial conditions T_*fh*0 > 0, I_*fh*0 ≥ 0, T_*fr*0 > 0, I_*fr*0 ≥ 0, V₀ ≥ 0, C_*f*0 ≥ 0, model system (1) has a unique solution that is nonnegative and bounded for all t > 0.*

Proof. We use classical differential equations theory [24] to demonstrate that all solutions to the model system (Equations 1–E6) are nonnegative and bounded. This confirms the existence of a unique local solution

T_*fh*, I_*fh*, T_*fr*, I_*fr*, V, C_*f*, for all t ∈ [0, t_*). Hence, we have

dTfhdt∣Tfh=0=λh≥0, dIfhdt∣Ifh=0=βhTfhV≥0, dTfrdt∣Tfr=0=λr+ρrIfh≥0, dIfrdt∣Ifr=0=βrTfrV≥0, dVdt∣V=0=NhδhIfh+NrδrIfr≥0,dCfdt∣Cf=0=ρcIfh,≥0.(7)

This shows that all solutions are nonnegative for t ∈ [0, t_*). We state the following theorem:

Theorem 1. *The feasible region of systems Equations (1) - (6) are such T_*h* = T_*fh* + I_*fh*, T_*r* = T_*fr* + I_*fr*, X = (T_*fh*, I_*fh*, T_*fr*, I_*fr*, V, C_*f*) and Ω = Ω₁ × Ω₂ × Ω₃ × Ω₄ where*

Ω1={(Tfh,Ifh)∈ℝ+2|0<Th<T0, μh>rh},Ω2={(Tfr,Ifr)∈ℝ+2|0<Tr<T1,μr>rh},Ω3={V∈ℝ+|0<V<T2,μh>rh},Ω4={Cf∈ℝ+|0<Cf<T3,μh>rh}.

and

Ω={(Tfh,Ifh,Tfr,Ifr,V,Cf)∈ℝ+6|0<Th≤T0,0<Tr<T1,0≤V≤T2,0≤Cf≤T3,μh>rh}.

The regions Ω_i*, i = 1, ⋯ , 4 are positively invariant, attracting all positive orbits in ℝ₊ and bounded, and every solution in Ω_*i*, i = 1, ⋯ , 4 is in Ω*.

For T_*h* = T_*fh* + I_*fh*, we have

dThdt≤λh-μhTh+rhThαrαr+Tfh≤λh-μhTh+rhTh=λh+(rh-μh)Th,

provided, μ_*h* > r_*h*. Solving the differential inequality gives,

Th≤λhrh-μh+(Th0-λhrh-μh)=K0+(Th0-K0)e(rh-μh)t,K0=λhrh-μh.

If T_*h*0 ≥ K₀, then, (Th0-K0)e(rh-μh)t≤Th0-K0 and T_*h* ≤ K₀ + (T_*h*0 − K₀) = T_*h*0. If T_*h*0 < K₀, then, (Th0-K0)e(rh-μh)t≤0 and T_*h* ≤ K₀. Hence,

Th≤max{Th0,λhμh-rh}:=T0,μh>rh.(8)

Similarly, for T_*r* = T_*fr* + I_*fr*, adding Equations 3, 4, we have

dTrdt=λr+ρrIfh-μr(Tfh+Ifr)≤λr+ρrIfh-μrTr≤λr+ρrT0-μrTr=A-μrTr,(9)

where A = λ_*r* + ρ_*r*T₀. Solving the differential inequality yields

Tr(t)≤max{Tr(0),Aμr}:=T1 μh>rh.

Considering Equation 5, we get

dVdt=NhδhIfh+NrδrIfr-μvV≤NhδhT0+NrδrT1-μvV=B-μvV,μh>rh,B=NhδhT0+NrδrT1.(10)

It follows that

V≤max{V(0),Bμv}:=T2,μh>rh.(11)

Equation 6 yields,

Cf≤{Cf(0),Cμc}:=T3,C=ρcT1,μh>rh.(12)

□

3.2 Disease-free equilibrium and the reproduction number

The disease-free equilibrium (DFE) is obtained by setting the model system (Equations 1–6) equal to zero, and solving for the cases where the infected compartments are equal to zero; we obtain

E0=(Tfh0,0,Tfr0,0,0,0), where Tfr0=λrμr,Tfh0=λhμh

We compute the basic reproduction number R0 for system (Equation 1) using the next generation matrix method as described in Van den Driessche and Watmough [25] to obtain

R0=βhTfh0Nhδhμv(μh+δh)+βrTfr0Nrδrμv(μr+δr), where Tfr0=λrμr,Tfh0=λhμh

Theorem 2. The disease-free equilibrium of system (Equation 1) is locally asymptotically stable if R0<1, and unstable if R0>1.

Proof. The result follows from the next-generation matrix method in Van den Driessche and Watmough [25]. When R0<1, the spectral radius of the next-generation matrix is less than one, implying that the infection cannot invade the population and that the DFE is locally asymptotically stable. □

3.2.1 Bifurcation analysis

We proceed using center manifold theory for bifurcation analysis. We choose the bifurcation βh=βh* for which R0=1. We let the state vector X = (S, Z) where S = (T_*fh*, T_*fr*) are the susceptible variables and Z=(Ifh,Ifr,V,Cf)⊤ are the infected and immune response variables. We write the system (Equation 1) in the form

{Ṡ=G(S,Z,βh),Ż=H(S,Z,βh).

At the DFE, we have Z = 0 and S=S0=(Tfh0,Tfr0). By hypothesis, the Jacobian of the full system evaluated at the DFE with βh=βh* has a simple eigenvalue 0 corresponding to the bifurcation direction (center subspace), and all other eigenvalues have nonzero real parts. Center manifold reduction [26–28] yields a smooth mapping

S=Φ(Z), Φ(0)=S0, DΦ(0)=0,

so that the reduced ODE gives the dynamics restricted to the center manifold

Ż=H~(Z,βh):=H(Φ(Z),Z,βh).

We expand Φ and H~ up to second order in Z. Write

Φ(Z)=S0+12Q(Z,Z)+O(||Z||3),

where Q is a symmetric bilinear map Q:ℝ⁴×ℝ⁴ → ℝ². (The linear term vanishes because S is non-centered at the DFE). Substituting into H and Taylor expanding around (S0,0,βh*), we obtain

H˜(Z,βh)=AZ+12B(Z,Z)+(βh−βh)C(Z)+O(‖Z‖3,(βh−βh)‖Z‖2),

where

• A=DZH(S0,0,βh*) is the linearization of the infected subsystem at DFE,
• B(Z, Z) collects the quadratic terms in Z coming both from the original quadratic terms in H and from the substitution, S = Φ(Z) which introduces quadratic dependence of H on Z,
• C(Z)=∂βhH(S0,0,βh*)[Z] is the linear variation with respect to the parameter.

Let w ∈ ℝ⁴ (column) and v ∈ ℝ⁴ (row) be the right and left eigenvectors of A associated with the simple eigenvalue 0:

Aw=0, v⊤A=0, v⊤w=1.

Under the standard normalizations (see Hassard and Wan [28]), the first Lyapunov coefficient a and the parameter derivative coefficient b are given by

a=12∑i,j,k=14viwjwk∂2H~i∂Zj∂Zk(0,βh*),(13)

b=∑i,j=14viwj∂2H~i∂Zj∂βh(0,βh*).(14)

Using the expansion H˜(Z,βh)=AZ+12B(Z,Z)+(βh-βh*)C(Z)+⋯, we identify

∂2H˜i∂Zj∂Zk(0,βh)=Bi(ej,ek), ∂2H˜i∂Zj∂βh(0,βh)=∂βhCi(ej),

where e_*j* are the coordinate vectors in ℝ⁴. Hence

a=12∑i,j,kviwjwkBi(ej,ek).

3.2.1.1 Computation of the quadratic form B

The quadratic map B receives contributions from two sources:

• (I) the intrinsic quadratic terms already present in the infected equations H(S, Z, β_*h*) when expanding around (S0,0,βh*);
• (II) the induced quadratic terms coming from substituting S=Φ(Z)=S0+12Q(Z,Z)+⋯ into the linear parts of H that depend on S. Concretely, if H depends on S linearly by terms of the form ℓ(S)·Z, then replacing S=S0+12Q(Z,Z) yields an extra quadratic contribution 12ℓ(Q(Z,Z)).

In our model the infected block Z = (I_*fh*, I_*fr*, V, C_*f*) and the original right-hand sides contain:

H1=βhTfhV-(μh+δh)Ifh-khCfIfh,H2=βrTfrV-(μr+δr)Ifr-krCfIfr,H3=NhδhIfh+NrδrIfr-μvV,H4=ρcIfh-μcCf-δcTfrCf.

We note that

• The intrinsic quadratic terms present in H (source I) are e.g., −k_*h*C_*f*I_*fh*, −k_*r*C_*f*I_*fr*, −δ_*c*T_*fr*C_*f* these yield explicit bilinear contributions to B whose coefficients are products of positive parameters.
• The induced contributions (source II) arise because H contains terms linear in S multiplied by Z (for example, the infection terms β_*h*T_*fh*V and β_*r*T_*fr*V are linear in T_*fh*, T_*fr*). Since Φ has a quadratic part Q(Z, Z), substituting Tfh=Tfh0+12Q1(Z,Z)+⋯ etc. produces quadratic terms of the form βh12Q1(Z,Z)V, etc.

We therefore compute Q (the quadratic dependence of S on Z) by solving the invariance equation for the center manifold up to second order. The invariance equation reads

DSG(S0,0,βh)Q(Z,Z)+2DZG(S0,0,βh)[Z,Q(Z,Z)]=2DZG(S0,0,βh*)[Z,Z]+⋯ ,

and, at quadratic order, it is a linear system for the components of Q because the matrix DSG(S0,0,βh*) is invertible on the noncenter directions and has a unique solution with components expressed as rational combinations of positive parameters. In particular, the entries of Q are linear combinations with positive coefficients of quadratic monomials in the entries of Z (since all biological parameters are positive and the source terms driving S from Z are nonnegative). Thus, the components Q₁, Q₂ that give the quadratic variation of T_*fh* and T_*fr* have the form

Q1(Z,Z)=∑j,kα1,jkZjZk, Q2(Z,Z)=∑j,kα2,jkZjZk,

with coefficients α_*m,jk* ≥ 0 (and in fact are strictly positive for the relevant index pairs that contribute to infection dynamics, because the driving linear terms and the inverse of D_*S*G have strictly positive entries).

Now substitute S=S0+12Q(Z,Z) into the linear-in-S infection terms of H. For example, the β_*h*T_*fh*V term becomes

βh(Tfh0+12Q1(Z,Z))V=βhTfh0V+12βhQ1(Z,Z)V.

The second term is quadratic in Z and contributes to B. Because β_*h* > 0, Tfh0>0, and the coefficients α_1,*jk* are nonnegative (indeed positive for relevant pairs), the induced contribution 12βhQ1(Z,Z)V produces positive coefficients in B for the quadratic monomials involving V and those Z_*j*, Z_*k* present in Q₁.

Collecting all contributions (intrinsic and induced), we obtain that each coordinate B_*i*(e_*j*, e_*k*) is a sum of terms, each product of model parameters (all positive) and possibly elements of the inverse of D_*S*G (which are finite and, under the biological parameter regime, nonnegative). Therefore, each relevant B_*i*(e_*j*, e_*k*) that appears in the sum defining a is nonnegative, and at least one of them is strictly positive (because infection and feedback terms are present). Concretely, there are contributions to B of the sign: B₁(e₃, e₃) contains a positive term from β_*h* · (coef of V² in Q₁), and similarly other B_*i*(e_*j*, e_*k*) > 0 for indices corresponding to I_*fh*, I_*fr*, V.

3.2.1.2 Sign of a

Using Equation 13,

a=12∑i,j,kviwjwkBi(ej,ek).

By the Perron–Frobenius type structure [29, 30] of the positive infection block at DFE (and the usual choice of normalization), we may choose the right eigenvector w and the left eigenvector v associated with the simple zero eigenvalue to have strictly positive components. We apply the theory to the system (Equation 1), as follows:

From the full model (Equation 1), consider the infected subsystem

Z=(Ifh,Ifr,V,Cf)⊤.

The Jacobian of the infected subsystem, evaluated at the DFE (S0,0,βh*), is

A=DZH(S0,0,βh*).

Using the ordering Z = (I_*fh*, I_*fr*, V, C_*f*), the Jacobian matrix is

A=(-(μh+δh)0βh*Tfh000-(μr+δr)βrTfr00NhδhNrδr-μv0ρc00-(μc+δcTfr0)).

We notice that nonlinear terms such as −k_*h*C_*f*I_*fh* are quadratic in Z and therefore vanish in the linearization. Assume that the matrix A has a simple zero eigenvalue when βh=βh*; this corresponds to the condition R0=1. Let

w=(w1,w2,w3,w4)⊤, v=(v1,v2,v3,v4)⊤

be the right and left eigenvectors, respectively, satisfying

Aw=0, v⊤A=0, v⊤w=1.

3.2.1.3 Right eigenvector w

Solving the equation Aw = 0 gives

w1=βh*Tfh0μh+δhw3, w2=βrTfr0μr+δrw3.

Hence, we choose w₃ to be positive (w₃ > 0). Then

w4=ρcμc+δcTfr0w1.

3.2.1.4 Left eigenvector v

The equation v^⊤A = 0 yields

v3=Nhδhμh+δhv1 and v3=Nrδrμr+δrv2.

Hence, v₁ and v₂ are proportional to v₃ and we choose v₃ > 0. Then

v1=Nhδhμh+δhv3, v2=Nrδrμr+δrv3, v4=0.

3.2.1.5 Normalization

We compute the scalar product v^⊤w. Since v₄ = 0, we obtain,

v⊤w=v1w1+v2w2+v3w3=v3w3S whereS:=1+Nhδhβh*Tfh0(μh+δh)2+NrδrβrTfr0(μr+δr)2>0.

To enforce the normalization v^⊤w = 1, we fix the free scalings in many equivalent ways. Two convenient choices are

• **We choose v₃ = 1**. Then set w3=1S. All components are determined as
• We choose w₃ = 1. Then set v3=1S. This yields the alternative normalized vectors

We proceed to compute the nonzero second derivatives D²H at the DFE by using the infected-variable ordering

Z=(Z1,Z2,Z3,Z4)⊤=(Ifh,Ifr,V,Cf)⊤.

Denoting by H=(H1,H2,H3,H4)⊤, the infected-subsystem vector field such that Ż = H(S, Z, β_*h*), we evaluate all second partial derivatives of H at the disease-free equilibrium (S0,0,βh), where S0=(Tfh0,Tfr0). The only quadratic terms in the infected equations that involve products of infected variables are −k_hC_fI_fh* in H₁ and −k_*r*C_*f*I_*fr* in H₂. Consequently, the only nonzero second partial derivatives of H with respect to the infected variables Z_*j*, Z_*k* at the DFE are the following mixed derivatives (symmetric in indices):

∂2H1∂Z1∂Z4|(S0,0,βh)=∂2H1∂Z4∂Z1|(S0,0,βh)=−kh,∂2H2∂Z2∂Z4|(S0,0,βh) =∂2H2∂Z4∂Z2|(S0,0,βh)=−kr.(15)

All other second derivatives ∂2Hi/(∂Zj∂Zk) vanish at the DFE because the corresponding terms in the model are at most linear in the infected variables when evaluated at S = S⁰.

The only dependence of H on β_*h* that yields a non-zero mixed second derivative at the DFE is from the term β_*h*T_*fh*V appearing in H₁. Evaluating at the DFE gives the single nonzero mixed derivative

∂2H1∂Z3∂βh|(S0,0,βh)=∂2H1∂βh∂Z3|(S0,0,βh)=Tfh0.(16)

All other mixed derivatives ∂2Hi/(∂Zj∂βh) vanish at the DFE.

We now compute the quadratic part Q of the center-manifold function

S=Φ(Z)=S0+12Q(Z,Z)+O(||Z||3),

which expresses the (non-center) variables S=(Tfh,Tfr)⊤ as a quadratic function of the infected variables Z=(Ifh,Ifr,V,Cf)⊤. The coefficients of Q are determined by the invariance (or tangency) equation for the center manifold. We write the S-block of the system as

Ṡ=G(S,Z,βh),

with components (recalling S₁ = T_*fh*, S₂ = T_*fr*)

G1(S,Z,βh)=λh-μhS1-βhS1V+rhαrS1αr+S2-αdS1S2,G2(S,Z,βh)=λr+ρrIfh-βrS2V-μrS2-σrS2S1.

Evaluate the Jacobian of G with respect to S at the DFE (S0,0,βh*):

L:=DSG(S0,0,βh*)=(-μh+rhαrαr+Tfr0-αdTfr0-rhαrTfh0(αr+Tfr0)2-αdTfh0-σrTfr0-μr-σrTfh0).

Remark 1. At the disease-free equilibrium, the susceptible compartments are strictly stable (their linearization has negative eigenvalues). Therefore, their 2 × 2 matrix L is invertible, and they do not affect the critical bifurcation dynamics. Only the infected subsystem drives the center manifold.

The invariance equation for the center manifold of indices 1 ≤ j ≤ k ≤ 4,

LQ(ej,ek)=Rjk,(17)

where Q(ej,ek)∈ℝ2 denotes the column vector of coefficients of the bilinear form Q for the pair (e_*j*, e_*k*), and Rjk∈ℝ2 collects the forcing terms coming from the second-order expansion of G with respect to the infected variables Z and from the interaction of linear S-dependence with the quadratic part of Φ (see e.g., Carr [26]). Concretely,

Rjk=DZ,Z2G(S0,0,βh*)[ej,ek]+T

where T are the terms involving D_*Z*G and Q, and will be determined as follows: We seek the quadratic part Q(Z, Z) of the center manifold

S=Φ(Z)=S0+12Q(Z,Z)+O(||Z||3),

and the invariance equation,

DΦ(Z)H~(Z,βh)=G(Φ(Z),Z,βh),(18)

where H~ is the reduced infected field (whose linear part at the DFE is A). Expanding (Equation 18) to second order in Z yields, for each pair 1 ≤ j ≤ k ≤ 4, a linear equation of the form

LQ(ej,ek)=Rjk,(19)

where L=DSG(S0,0,βh*) and R_*jk* is the forcing vector depending on second derivatives of G, on D_*Z*G, on the quadratic form Q itself, and on the infected linear operator A.

Writing the Taylor expansions about (S0,0,βh*), for G, we have

G(S,Z,βh)=G(S0,0,βh)+DSG(S0,0,βh)(S-S0)+DZG(S0,0,βh)[Z] +12DS,S2G(S0,0,βh)[S-S0,S-S0] +DS,Z2G(S0,0,βh)[S-S0,Z]+12DZ,Z2G(S0,0,βh*)[Z,Z]+O(||S-S0||3,||Z||3).(20)

At the DFE the constant term G(S0,0,βh*)=0 and for the center-manifold map Φ, we write

Φ(Z)=S0+12Q(Z,Z)+O(||Z||3),

so S-S0=12Q(Z,Z)+O(||Z||3). For the reduced infected dynamics, we write

H~(Z,βh*)=AZ+12B(Z,Z)+O(||Z||3),(21)

where A=DZH(S0,0,βh*) and B is the intrinsic bilinear form. Finally, the derivative of Φ is

DΦ(Z)=D(12Q(Z,Z))+O(||Z||2)=Q(Z,·)+O(||Z||2),

where Q(Z, ·) denotes the linear map u ↦ Q(Z, u). We note that DΦ(0) = 0. Expanding the left-hand side of Equation 18, we have,

LHS=DΦ(Z)H~(Z,βh*)=(Q(Z,·)+O(||Z||2))(AZ+12B(Z,Z)+O(||Z||3)).(22)

We collect quadratic terms in Z: the only quadratic contributions from Equation 22 are

Q(Z,AZ)

(since Q(Z,·) is linear in the second slot and AZ is linear in Z).

Evaluated on a coordinate pair, (e_*j*, e_*k*) the contribution is

Q(ej,Aek)+Q(ek,Aej),

because the bilinear form is symmetric in its two arguments when we collect the Z_*j*Z_*k* coefficient. Using Equation 20 and substituting S-S0=12Q(Z,Z), the quadratic part of the right-hand side is

RHS=DSG(S0,0,βh)[12Q(Z,Z)]+DZG(S0,0,βh)[Z] +12DZ,Z2G(S0,0,βh)[Z,Z] +DS,Z2G(S0,0,βh)[12Q(Z,Z),Z]+O(||Z||3).(23)

We equate the quadratic coefficients on the LHS and the RHS for each pair (j, k), to obtain

Q(ej,Aek)+Q(ek,Aej)=DSG(S0,0,βh)[Q(ej,ek)]+DZ,Z2G(S0,0,βh)[ej,ek]+Tjk,

where Tjk collects the remaining contributions that involve the linear map D_*Z*G and possibly mixed second derivatives DS,Z2G combined with components of Q. Rearranging,

LQ(ej,ek)=DZ,Z2G(S0,0,βh*)[ej,ek]+{Tjk-(Q(ej,Aek)+Q(ek,Aej))}︸=:T^jk,(24)

with L=DSG(S0,0,βh*). Thus, the forcing term appearing in Equation 19 is

Rjk=DZ,Z2G(S0,0,βh*)[ej,ek]+T^jk.

we write a more explicit formula for the vector R_*jk*:

Rjk=DZ,Z2G(S0,0,βh)[ej,ek]+DZG(S0,0,βh)[ej](Q(ek,·)) +DZG(S0,0,βh)[ek](Q(ej,·)) -(Q(ej,Aek)+Q(ek,Aej)) +DS,Z2G(S0,0,βh)[Q(ej,ek),·]+⋯

We let T to be

T=DZG(S0,0,βh)[ej](Q(ek,·))+DZG(S0,0,βh)[ek](Q(ej,·))-(Q(ej,Aek)+Q(ek,Aej))+DS,Z2G(S0,0,βh*)[Q(ej,ek),·]+⋯

Here “·” denotes the slot where the remaining vector argument is inserted, and “⋯ ” denotes any other terms of the same quadratic order.

Remark 2. In many compartmental epidemic models

• DZ,Z2G(S0,0,βh)≡0 because the S-equations are linear in the infected variables at the DFE (they contain terms like S · V but not V² or Ifh2). Hence, the first term is zero*.
• *DZG(S0,0) is nonzero in general, but if Q ≡ 0 (which can happen when the invariance equation admits only the trivial quadratic solution), then the middle terms vanish*.
• *The terms Q(e_*j*, Ae_*k*)+Q(e_*k*, Ae_*j*) depend on Q and A and can thus be moved to the left to form LQ*.

*Thus, in our model, the right-hand side often simplifies to LQ(ej,ek)=0, which implies Q ≡ 0.*

Vanishing of the forcing terms R_*jk*, at the DFE (S = S⁰, Z = 0), we have

• G₁ contains the term −β_*h*S₁V; as a function of Z this is linear in V (i.e., linear in Z₃), hence DZ,Z2G1(S0,0)=0.
• G₂ contains the term ρ_*r*I_*fh* (linear in Z₁) and −β_*r*S₂V (linear in V at DFE); thus DZ,Z2G2(S0,0)=0.
• There are no terms in G that are quadratic purely in the infected variables Z (no V², Ifh2, Ifr2, IV-type intrinsic quadratic) when evaluated at S = S⁰.

Therefore, the direct second derivatives of G with respect to infected variables vanish, and

DZ,Z2G(S0,0,βh*)≡0.

Consequently, the right-hand side R_*jk* in Equation 17 reduces to terms that are linear in Q (coming from the coupling D_*Z*G and the substitution), but these terms are proportional to Q. At quadratic order, the invariance equation therefore takes the homogeneous form

LQ(ej,ek)=0 for all 1≤j≤k≤4.

Hence, since L is invertible, the unique solution of LY=0 is Y = 0. Thus, for every pair (j, k) we obtain Q(e_*j*, e_*k*) = 0 and the quadratic part of the center manifold vanishes:

Q(Z,Z)≡0.

Considering that the quadratic part of the center manifold Q(Z, Z) vanishes, the bilinear map B:ℝ⁴×ℝ⁴ → ℝ⁴ that appears in the reduced vector field (Equation 21) is exactly the intrinsic quadratic part of the infected subsystem. The only nonzero components of B are those coming from the quadratic terms identified and from Equation 15, we have

∂2H1∂Z1∂Z4|DFE=∂2H1∂Z4∂Z1|DFE=−kh,∂2H2∂Z2∂Z4|DFE=∂2H2∂Z4∂Z2|DFE=−kr.

Hence, the bilinear form B has the only nonzero entries

B1(e1,e4)=B1(e4,e1)=-kh,B2(e2,e4)=B2(e4,e2)=-kr,

and all other components B_*i*(e_*j*, e_*k*) vanish. Hence, using the nonzero entries of B above with the contributions involving the pairs (j, k) = (1, 4) and (2, 4):

a=12[v1(w1w4B1(e1,e4)+w4w1B1(e4,e1))+v2(w2w4B2(e2,e4)+w4w2B2(e4,e2))]=12[v1(2w1w4(-kh))+v2(2w2w4(-kr))]=-(khv1w1w4+krv2w2w4).

3.2.1.6 Sign of a

Under the biological sign convention, all parameters k_*h*, k_*r* are strictly positive, and the components of the left and right eigenvectors associated with the infection subspace can be chosen strictly positive (cf. Step 1). Therefore, each term in the parentheses is purely positive, and we conclude a < 0

3.2.1.7 Sign of b

We found that the only nonzero mixed second derivative of the infected population with respect to the parameter β_*h* and an infected variable is

∂2H1∂V∂βh|(S0,0,βh*)=Tfh0,

where Z = (Z₁, Z₂, Z₃, Z₄) = (I_*fh*, I_*fr*, V, C_*f*) and all other mixed derivatives ∂2Hi/(∂Zj∂βh) vanish at the DFE. Therefore the sum defining b reduces to a single term (with i = 1, j = 3):

b=∑i,j=14viwj∂2H~i∂Zj∂βh|(0,βh)=v1w3∂2H~1∂Z3∂βh|(0,βh)=v1w3Tfh0

Since Tfh0=λhμh>0 and, with the chosen normalization, the eigenvector components v₁ > 0 and w₃ > 0 (cf. Step 1), we conclude

b=v1w3Tfh0>0.

Combined with the previously computed sign b > 0 (see Step 2), the sign pair (a < 0, b > 0) corresponds to a forward supercritical bifurcation at R0=1. Hence, a stable endemic branch emerges for R0>1.

3.3 Sensitivity analysis of the basic reproduction number

Consider the basic reproduction number:

R0=NhδhβhTfh0μv(μh+δh)+NrδrβrTfr0μv(μr+δr)

The normalized forward sensitivity index of R0 with respect to a parameter p is defined by:

ΥpR0=∂R0∂p×pR0

and we have the following Table 1 using the parameters values of Table 2.

Table 1. Sensitivity indices of the basic reproduction number R0 with respect to model parameters.

Table 2. Model parameters with biological interpretation and sources/justifications.

Parameters are classified as experimentally supported, literature-calibrated, or phenomenological assumptions [12, 31–39].

From Table 1, the infection is driven mainly by the T_*fh* pathway, ≈ 89%. β_*h* and N_*h* increasing has the strongest impact. The T_*fr* pathway contributes about 11%. Viral clearance and has a perfect negative sensitivity −1. Increasing δ_*h*, increases R0 but weakly with 10%. The sensitivity of R0 to δ_*r* is roughly 3%. This confirms that the T_*fr* pathway plays a relatively minor role compared to the T_*fh* pathway in driving infection dynamics.

3.4 Endemic equilibrium

In this section, we study the endemic equilibrium

Proposition 2. Let

Ω:={X=(Tfh,Ifh,Tfr,Ifr,V,Cf)∈ℝ6:Tfh≥0,Ifh≥0,Tfr≥0,Ifr≥0,V≥0,Cf≥0}.

*Define F:Ω → ℝ⁶ by the right-hand side of system (Equation 1), i.e.,*

F(X)=(F1,…,F6)T

with

F1=λh-μhTfh-βhTfhV+rhαrTfhαr+Tfr-αdTfhTfr,F2=βhTfhV-(μh+δh)Ifh-khCfIfhF3=λr+ρrIfh-βrTfrV-μrTfr-σrTfrTfh, F4=βrTfrV-(μr+δr)Ifr-krCfIfr,F5=NhδhIfh+NrδrIfr-μvV, F6=ρcIfh-μcCf-δcTfrCf.

*Then F ∈ C^∞(Ω). In particular, F is locally Lipschitz on Ω, that is, for every compact set K ⊂ Ω there exists a constant L_*K* > 0 such that*

where ||·|| is any norm on ℝ⁶. One may take, for instance,

LK=max1≤i≤6∑j=16supZ∈K|∂xjFi(Z)|.

Proof. Each F_*i* is a finite sum of polynomial terms in the coordinates and for the rational term rhαrTfhαr+Tfr. On Ω, we have T_*fr* ≥ 0 and α_*r* > 0, hence, α_*r* + T_*fr* ≥ α_*r* > 0 and the denominator never vanishes on Ω. Therefore, each F_*i* is composed of smooth functions (polynomials and division by a strictly positive affine function), so F_*i* is C^∞ on Ω. Consequently F ∈ C^∞(Ω).Since F is C¹ on Ω, its Jacobian matrix

DF(X)=(∂xjFi(X))i,j=16

is continuous on Ω. Let K ⊂ Ω be compact. By continuity of DF on K, each entry ∂_{x_*j*}F_*i* attains a finite maximum on K. Define

Mij:=supZ∈K|∂xjFi(Z)|<∞, LK:=max1≤i≤6∑j=16Mij.

Then, for any X, Y ∈ K the mean-value form gives

F(X)-F(Y)=(∫01DF(Y+s(X-Y))ds)(X-Y),

and taking the ℓ^∞-operator norm yields

By the equivalence of norms on finite-dimensional spaces, the same inequality holds for any norms, possibly with a different constant. This establishes the local Lipschitz property on K.

For completeness, one can compute the nonzero partial derivatives. Writing X = (x₁, x₂, x₃, x₄, x₅, x₆) = (T_*fh*, I_*fh*, T_*fr*, I_*fr*, V, C_*f*), the non-zero partials are

∂x1F1=-μh-βhx5+rhαrαr+x3-αdx3,∂x3F1= -rhαrx1(αr+x3)2-αdx1,∂x5F1=-βhx1,∂x1F2=βhx5,∂x2F2=-(μh+δh)-khx6,∂x5F2=βhx1,∂x6F2=-khx2,∂x1F3=-σrx3,∂x2F3=ρr,∂x3F3=-βrx5-μr-σrx1,∂x5F3=-βrx3,∂x3F4=βrx5,∂x4F4=-(μr+δr)-krx6,∂x5F4=βrx3,∂x6F4=-krx4,∂x2F5=Nhδh,∂x4F5=Nrδr,∂x5F5=-μv,∂x2F6=ρc,∂x3F6=-δcx6,∂x6F6=-μc-δcx3.

Each of these entries is continuous on Ω and therefore bounded on any compact K ⊂ Ω; these bounds are precisely the M_*ij* used to form L_*K*.

Thus, F is C^∞(Ω) and locally Lipschitz on Ω, which proves the proposition. Hence, by the Picard–Lindelöf (Cauchy–Lipschitz) theorem [40–42], the initial value problem

dFdt=F(t,X(t)), F(0)=F0,

has a unique local solution F(t) passing through F0=(Tfh0,Ifh0,Tfr0,Ifr0,V0,Cf0)T. The system (Equation 1) thus, admits a unique local solution in ℝ+6. □

To qualitatively understand the unique solution, we proceed to analyze the equilibrium point EE=(Tfh,Tfr,Ifh,Ifr,V,Cf), associated with the solution by setting dFdt=0, and obtain

Ifh=βhTfhVμh+δh+khCf, Ifr=βrTfrVμr+δr+krCf,Cf=ρcβhTfhV(μh+δh+khCf)(μc+δcTfr)V=NhδhβhTfhVμv(μh+δh+khCf)+NrδrβrTfrVμv(μr+δr+krCf),Tfr=λrβrV+μr+σrTfh+ρrβhTfhV(μh+δh+khCf)(βrV+μr+σrTfh).(25)

Substituting the expressions for Tfr*, V^*, Cf* into the equilibrium equation for Tfh*, we obtain a nonlinear algebraic equation of the form:

P(Tfh)=a4(Tfh)4+a3(Tfh)3+a2(Tfh)2+a1Tfh*+a0=0

where the coefficients a₀, a₁, a₂, a₃, a₄ are explicit functions of the system parameters defined by

ai=fi(λh,λr,μh,μr,μv,δh,δr,δc,ρr,ρc,Nh,Nr,βh,βr,kh,kr,αr,αd,rh,σr).

We can show that the unique equilibrium of model (Equation 1) is globally asymptotically stable by using a Volterra-type function:

H(Tfh,Ifh,Tfr,Ifr,V,Cf)=∑X∈{Tfh,Ifh,Tfr,Ifr,V,Cf}(Xi-Xi-Xiln XiXi*)=H1+H2+H3+H4+H5+H6,

where H1=Tfh-Tfh-Tfhln(TfhTfh), H2=Ifh-Ifh-Ifhln(IfhIfh), H3=Tfr-Tfr-Tfrln(TfrTfr), H4=Ifr-Ifr-Ifrln(IfrIfr), H5=V-V-Vln(VV), H6=Cf-Cf-Cfln(CfCf). We have

H˙=∑X∈{Tfh,Ifh,Tfr,Ifr,V,Cf}(1-Xi*Xi)Ẋi=H˙1+H˙2+H˙3+H˙4+H˙5+H˙6.

Let Tfh=xTfh*, V = yV^*, Ifh=mIfh, Ifr=pIfr, Tfr=zTfr* and Cf=uCf. Using the arithmetic-geometric mean and the inequality ln(w*) ≤ w − 1, for any w > 0, we have for H˙i, i = 1, ⋯ , 6

ℋ.1=−μhTfh(x−1)2x−βhTfhV(1−x)(1−xy)x+(rhαrTfr)x· αrTfr(x−1)2+(x−1)(x−z)(αrTfr+z)(αr+Tfr) −αdTfhTfr(1−x)(1−xz)x ≤−μhTfh(x−1)2x−βhTfhVx[(−ln(x))(−ln(xy))] −αdTfhTfrx[(−ln(x))(−ln(xz))] −(rhαrTfr)x·ln(αrTfr(x−1)2+x2+z22+2)(αrTfr+z)(αr+Tfr)≤0.

ℋ.2=m−1m [βhTfhV(xy−1)−(μh+δh)Ifh(m−1)−khCfIfh(um−1)] ≤−12βhTfhVln(m+(xy)2m)−(μh+δh)Ifhm(m−1)2− khCfIfhm[−ln(m)(−ln(um))]≤0.

ℋ.3=z−1z[ρrIfh(m−1)−βrTfrV(zy−1)−μrTfr(z−1) −σrTfrTfh(zx−1)] ≤1z[−ρrIfh2ln(m2+z2+3)−βrTfrV[−ln(z)×−ln(zy)] −μrTfr[ln(z)]2] −1z[σrTfrTfh[−ln(z)×−ln(zx)]]≤0.

ℋ.4=p−1p[βrTfrV(zy−1)−(μr+δr)Ifr(p−1) −krCfIfr(up−1)] ≤1p[−βrTfrV2ln(p2+(zy)2+3)−(μr+δr)Ifr(p−1)2 −krCfIfr((−ln(p)×−ln(up))]≤0.

ℋ.5=y−1y[NhδhIfh(m−1)+NrδrIfr(p−1)−μvV(y−1)] ≤−NhδhIfh2yln(m2+y2+3)−NrδrIfr2yln(p2+y2+3)− μvV(y−1)2y≤0.

ℋ.6=u−1u[ρcIfh(m−1)−μcCf(u−1)−δcTfrCf(zu−1)] ≤−ρcIfh2uln(m2+u2+3)−μcCfu(u−1)2− δcTfrCfu[−ln(u)×−ln(zu)]≤0.

Since H≥0 and H˙≤0 with equality only at EE, by the Lyapunov-LaSalle Theorem [43, 44], the endemic equilibrium EE is globally asymptotically stable in the positive orthant.

3.5 Parameter estimation

The parameter values used in this study (see Table 2) are obtained from a combination of experimental studies, HIV within-host modeling literature, and biologically motivated assumptions where direct estimates are unavailable. For instance, using an initial value of T_fh₀ = 10, 000 and the natural death rate μ_*h* = 1/30, to ensure population conservation, we calculate the value of λ_*h* = μ_*h*T_fh₀ = 10, 000/30. The value of μ_*h* is nominal as the lifespan of the T_*fh* cells is dynamic, and we used a value within the lifespan related to CD4 ⁺ T cells [34, 35]. The estimation of λ_*r* follows a similar approach used for λ_*h*. The studies in Rychagov [36] and Lou et al. [37] estimates δ_*h* = 0.24. The studies [32, 36–38] show variability in viral clearance rates depending on the site of infection and the clearance rate. The values range from 2.4 to 500 per day, with lower values indicating slow clearance and higher values indicating rapid clearance. We used μ_*v* = 2.6, which is indicative of low viral clearance in the germinal center (GC). The study in Samaddar et al. [39] indicates that HIV infection of CD4 ⁺ T cells results in immune system impairment. The remaining parameters, which are less directly measurable and not direct experimental measurements, were estimated phenomenologically and mechanistically via model calibration, yielding baseline simulations that represent an ideal establishment of infection from the infection-free equilibrium while accounting for population conservation. Given this heterogeneity, the results should be interpreted qualitatively and structurally, with regard to threshold dynamics and sensitivity rankings. However, the simulations would provide an idea of the effects of changing these parameters.

4 Numerical methods and numerical simulations

In this section, we use Nonstandard Finite Difference (NSFD). The NSFD was early used by the pioneer Mickens [45] and mathematically improved by Anguelov and Lubuma [46]. For all about the nonstandard finite difference method and its construction, we refer the reader to the following Anguelov and Lubuma [46] (Theorem 12, 13) [45–47].

The choice of the NSFD scheme is motivated by the need to preserve the qualitative properties of the epidemiological model, such as the positivity and boundedness of solutions, as well as the preservation of equilibrium stability properties. Standard adaptive stiff solvers, while efficient and widely used, do not, in general, guarantee the preservation of these structural properties and sometimes introduce infeasible negative solutions or distort long-term dynamics, especially in strongly nonlinear epidemiological systems.

4.1 NSFD discretization of the model

In this section, we discretize our model (Equation 1) and investigate its dynamical consistency (positivity, boundedness). In the following, we propose the Theorem (Equation 26)

Theorem 3. The constructed NSFD following the scheme (Equation 26)

{Tfhk+1=[1+ζ(ℓ)rhαrαr+Tfrk]Tfhk+ζ(ℓ)λh1+ζ(ℓ)(μh+βhVk+αdTfrk),Ifhk+1=Ifhk+ζ(ℓ)βhTfhk+1Vk1+ζ(ℓ)(μh+δh+khCfk),Tfrk+1=Tfrk+ζ(ℓ)(λr+ρrIfhk)1+ζ(ℓ)(βrVk+μr+σrTfhk+1),Ifrk+1=Ifrk+ζ(ℓ)βrTfrk+1Vk1+ζ(ℓ)(μr+δr+krCfk),Vk+1=Vk+ζ(ℓ)(NhδhIfhk+1+NrδrIfrk+1)1+ζ(ℓ)μv,Cfk+1=Cfk+ζ(ℓ)ρcIfhk+11+ζ(ℓ)(μc+δcTfrk+1).(26)

is dynamically consistent with the original model (refer to Equation 1), provided it is well posed (please refer to Section 3) and furthermore preserves its positivity and the boundedness of the solutions with respect to Theorem 1.

Proof. Let N_*t* be a positive integer denoting the number of subintervals in the time direction. Define the time step size by ℓ=Δt=TNt. The time interval [0, T] is then uniformly partitioned as follows

t0=0<t1<⋯<tNt-1<tNt=T,

where each node is given by t_*k*+1 = t_*k* + ℓ, for k = 0, 1, …, N_*t* − 1. At each grid point t_*k*, we denote the numerical approximations to the functions T_*fh*(t_*k*), I_*fh*(t_*k*), T_*fr*(t_*k*), I_*fr*(t_*k*), V(t_*k*), and C_*f*(t_*k*) by

Tfhk, Ifhk, Tfrk, Ifrk, Vk, and Cfk,

respectively.

To construct the NSFD, we will discretize the system (Equation 1) by approximating the temporal derivatives using a generalized first-order method. For example, the discrete approximation of the derivative of T_*fh*(t) is given by

dTfhdt=Tfhk+1-Tfhkζ(ℓ)+O(ζ(ℓ)), as ℓ→0,

where ζ(ℓ) is a function of the time step size ℓ and Tfhk denotes the approximation of T_*fh*(t_*k*), which is a real-valued function satisfying

ζ(ℓ)=ℓ+O(ℓ2), for all ℓ>0.

Here, ℓ denotes the time-step size, and Tfhk is the numerical approximation of T_*fh*(t_*k*). Applying this to the evolution equation of T_*fh*, we obtain the discrete form

Tfhk+1-Tfhkζ(ℓ)=λh-μhTfhk+1-βhTfhk+1Vk+rhTfhkαr+Tfrkαr-αdTfhk+1Tfrk.

The same for the rest Ifhk+1,Tfrk+1,Ifrk+1,Vk+1,Cfk+1 and after arranging, we have

{Tfhk+1=[1+ζ(ℓ)rhαrαr+Tfrk]Tfhk+ζ(ℓ)λh1+ζ(ℓ)(μh+βhVk+αdTfrk),Ifhk+1=Ifhk+ζ(ℓ)βhTfhk+1Vk1+ζ(ℓ)(μh+δh+khCfk),Tfrk+1=Tfrk+ζ(ℓ)(λr+ρrIfhk)1+ζ(ℓ)(βrVk+μr+σrTfhk+1),Ifrk+1=Ifrk+ζ(ℓ)βrTfrk+1Vk1+ζ(ℓ)(μr+δr+krCfk),Vk+1=Vk+ζ(ℓ)(NhδhIfhk+1+NrδrIfrk+1)1+ζ(ℓ)μv,Cfk+1=Cfk+ζ(ℓ)ρcIfhk+11+ζ(ℓ)(μc+δcTfrk+1).(27)

where

q=λh+μh+βh+rh+αr+αd+δh+kh+λr+ρr+βr+μr+σr+βr+σr+kr+δr+μv+ρc+μc+δc

and ζ(ℓ)=1-e-qℓq. The constructed NSFD (Equation 26) is positive since all variables and parameters are positive. For boundedness, we just need to prove that

Thk+1=Tfhk+1+Ifhk+1≤λhμh-rhby assuming Thk=Tfhk+Ifhk≤λhμh-rh.

We have

{[1+ζ(ℓ)(μh+βhVk+αdTfrk)]Tfhk+1=[1+ζ(ℓ)rhαrαr+Tfrk]Tfhk+ζ(ℓ)λh,[1+ζ(ℓ)(μh+δh+khCfk)]Ifhk+1=Ifhk+ζ(ℓ)βhTfhk+1Vk(28)

By adding and arranging, we have,

[Tfhk+1+Ifhk+1](1+ζ(ℓ)μh)=Tfhk+ζ(ℓ)rhαrαr+TfrkTfhk+ζ(ℓ)λh−ζ(ℓ) (βhVk+αdTfrk)Tfhk+1 +Ifhk+ζ(ℓ)βhTfhk+1Vk−ζ(ℓ)(δh+khCfk))Ifhk+1

And since αrαr+Tfrk<1, Tfhk≤Thk and Thk≤λhμh-rh, we have (1+ζ(ℓ)μh)Thk+1≤[λhμh-rh](1+ζ(ℓ)μh). Hence Thk+1≤λhμh-rh. With the same reasoning for Trk+1,Vk+1,Cfk+1, ∀μ_*h* > r_*h*, we have

Trk+1=Tfrk+1+Ifrk+1≤Aμr,A=λr+ρrλhμh-rh,Vk+1≤Bμv,B=Nhδhλhμh-rh+NrδrAμr,Cfk+1≤Cμc,C=ρcλhμh-rh.

For consistency, consider Xk+1=(Tfhk+1,Ifhk+1,Tfrk+1,Ifrk+1,Vk+1,Cfk+1). Taylor expansion about the point (t_*k* + ζ(ℓ)) gives

Xk+1=X(tk+ζ(ℓ))=X(tk)+dX(tk)dtζ(ℓ)+d2Xdt2ζ(ℓ)22+O(ζ(ℓ)k3).

Applying this to each variable Tfhk+1,Ifhk+1,Tfrk+1,Ifrk+1,Vk+1,Cfk+1, we have

Tfhk+1=Tfh+Tfhdtζ(ℓ)+d2Tfhdt2ζ(ℓ)22+O(ζ(ℓ)3)=[1+ζ(ℓ)rhαrαr+Tfr]Tfh+ζ(ℓ)λh1+ζ(ℓ)(μh+βhV+αdTfr)

After simplification for ζ(ℓ), we have

(μh+βhV+αdTfr)Tfh+Tfhdt+(μh+βhV+αdTfr)Tfhdtζ(ℓ)+d2Tfhdt2ζ(ℓ)2+(μh+βhV+αdTfr)d2Tfhdt2ζ(ℓ)2+O(ζ(ℓ)2)=rhαrαr+TfrTfh+λh

As ζ(ℓ) → 0, we have

(μh+βhV+αdTfr)Tfh+Tfhdt=rhαrαr+TfrTfh+λh

Hence the consistency for Tfhk+1. It is easy to prove with the same reasoning for Ifhk+1,Tfrk+1,Ifrk+1,Vk+1,Cfk+1 respectively and ζ(ℓ) → 0, we have the consistency. □

To estimate the error for the constructed NSFD, let X(t)=(Tfh(t),Ifh(t),Tfr(t),Ifr(t),V(t),Cf(t))T∈ℝ6 be the solution of the corresponding continuous system, and define the discrete approximations Xk≈X(tk) with t_*k* = kζ(ℓ).

Theorem 4 (First-order convergence in time). *Assume the exact solution X(t) is sufficiently smooth and the nonlinear functions in the update scheme are Lipschitz continuous. Then the global error satisfies*

*for 0 ≤ k ≤ N, where Nζ(ℓ) = T is the final time, and C is a constant independent of ζ(ℓ)*.

Proof. It is easy to prove that for 0 ≤ k ≤ N, where Nζ(ℓ) = T is the final time, and C is a constant independent of ζ(ℓ), and ek=X(tk)-Xk, we have ||X(tk)-Xk||=||ek||≤Cζ(ℓ).

□

Remark 3. The NSFD (Equation 26) scheme constructed imposed static caps such as

Thmax=λhμh-rh,

which becomes singular or negative when μ_h* ≤ r_*h*. To remove this restriction, we replace these static caps with time-dependent comparison bounds obtained from the differential inequalities satisfied by the total population. Specifically, letting r = r_*h* − μ_*h* and denoting the initial total host population by T_*h*(0), we set*

[Thcap]k={λhμh-rh,μh>rh,Th0+λhk,μh=rh,(Th0+λhκ)erk-λhκ,μh<rh,

*which provides a non-singular upper bound that is valid across all parameter regimes. The corresponding bounds for T_*r*, V, and C_*f* are then propagated as*

[Trcap]k=λr+ρr[Thcap]kμr, [Vcap]k=Nhδh[Thcap]k+Nrδrμr[Trcap]kμv, [Cfcap]k=ρc[Thcap]kμc.

These time-dependent caps replace the singular static expressions and ensure that the projection steps in the numerical scheme (Equation 26) remain well-defined and biologically consistent even when μ_h* ≤ r_*h**.

4.2 Numerical simulations

In this section, we simulate the model (Equation 1) using the constructed NSFD scheme (Equation 26) and the parameter values in Table 2. We present the baseline curves in Figure 1, representing the ideal scenario for the establishment of infection. We then provide a case-scenario analysis as follows: (i) the effects of the Tfh parameters on all infection variables, (ii) the effects of the Tfr parameters on all infection variables, and (iii) the effects of the CD8+ T cell parameters on all infection variables.

Figure 1 shows the baseline graphs of the T_*fh* follicular helper cells, infected I_*fh* follicular helper cells, T_*fr* follicular regulatory cells, infected I_*fr* follicular regulatory cells, viral cells V produced from infected T follicular helper cells, and follicular CD8⁺ T cells over time in days (500 days). The figure presents an ideal scenario for the establishment of infection, in which the infection-free variables decrease, and the infected population increases over time. We use the baseline graphs to explore the cases (i)–(iii) on the persistence of infection.

4.2.1 The effects of T_*fh* parameters on I_*fh*, I_*fr*, and V

Figure 2 illustrates the influence of the T_*fh* conversion to the anergic state rate, α_*d*, on the within-host dynamics of HIV infection. The simulations show that the increase in α_*d* is associated with a reduction in the population of infected follicular helper T cells (I_*fh*), infected follicular regulatory T cells (I_*fr*), and viral (V). This trend suggests that higher levels of T_*fh* conversion to an anergic state by the T_*fr* cells hinder the ability of the virus to establish and maintain infection within the germinal center microenvironment. These results imply that targeting the enhancement of T_*fr* effects on regulating the T_*fh* cells to their anergic state may be one of the ways to suppress viral replication and promote immune-mediated control of the infection. Therefore, an increase in α_*d* creates conditions under which complete viral clearance becomes possible.

Figure 3 demonstrates the impact of the killing rate, k_*h*, of follicular CD8⁺ T cells on the progression of HIV infection within the germinal center. The simulations reveal that increasing k_*h* enhances the cytotoxic activity of follicular CD8⁺ T cells against infected populations, resulting in a decline in the populations of infected follicular helper T cells (I_*fh*), infected follicular regulatory T cells (I_*fr*), and free viral (V). This reduction reflects the ability of the immune response to limit viral replication and suppress the spread of infection. However, despite the decrease in both infected cell populations and viral load, neither I_*fh* nor I_*fr* is reduced to zero, converging to positive steady-state levels, indicating the continued presence of infected cells within the follicular microenvironment. Similarly, the viral population is sustained above zero, suggesting that viral replication persists even in the presence of enhanced follicular CD8⁺ T-cell-mediated killing. Thus, in this case, the immune response system contributes to improved immune control and reduces the infection burden, but is insufficient to eliminate the infection. The system approaches an endemic state in which the virus and infected cell populations are reduced but at lower persistent levels in the germinal centers. These findings suggest that strengthening the cytotoxic efficacy of follicular CD8⁺ T cells alone may not be adequate for viral clearance and that additional immune mechanisms or therapeutic interventions may be required to achieve complete infection clearance.

Figure 4 illustrates the effect of the lysis-induced death rate, δ_*h*, of infected follicular helper T cells (I_*fh*) on the dynamics of the infected compartments and viral population. The simulations show that increasing δ_*h* produces different effects among the infected populations. Specifically, a higher lysis death rate leads to a reduction in the population of infected follicular regulatory T cells (I_*fr*), indicating that the infection becomes less established within this compartment. However, the same increase in δ_*h* is accompanied by an increase in both the infected follicular helper T-cell population (I_*fh*) and the free virus concentration (V). This suggests that increased lysis of infected follicular helper cells creates conditions that favor increased viral replication and persistence, potentially shifting the infection to T_*fh* cells. Furthermore, Figure 4 also shows a backward shift in the population peaks of I_*fh*, I_*fr*, and V with higher and earlier shifts occurring with increased values of the lysis rate. This alters the timing of disease prognosis. The elevated levels of I_*fh* and free virus imply that the infection is established and sustained, early ensuring prolonged viral persistence and continued circulation of infected follicular helper cells. The results highlight the complex role of the lysis death rate δ_*h* in simultaneously promoting viral maintenance within the germinal centers.

In Figure 5, the variation of the saturation constant governing the regulatory activity of T follicular regulatory cells, α_*r*, is investigated for its impact on the infected compartments. As α_*r* increases, the regulatory feedback is associated with a delayed onset of growth in the infected follicular helper cell population, the infected follicular regulatory cell population, and the viral load. This delay indicates a slower early-stage reduction of infection under stronger saturation effects in the regulatory pathway. However, despite these temporal shifts and the corresponding increase in peak amplitudes in these infected populations, the long-term dynamics remain unchanged: the system converges to the same equilibrium states, indicating that α_*r* primarily influences transient behavior rather than steady-state infection levels.

4.2.2 The effects of T_*fr* parameters on I_*fh*, I_*fr*, and V

In Figure 6, the increase in the conversion of T_*fr* cells to ex-T_*fr* rate, σ_*r*, is associated with a decrease in both the infected T_*fr* cell population and the viral load V while the infected T_*fh* cell population exhibits a transient increase, indicating a redistribution of infection dynamics toward the T follicular helper cell population. The dynamics reveal that the peaks of all infected populations occur earlier at larger values of σ_*r*, suggesting accelerated disease progression under stronger conversion. The highest peak values are observed for the infected T_*fh* cells and viral population when σ_*r* is large, whereas the infected T_*fr* cells display a lower peak. Despite these shifts in infection dynamics, the system does not achieve full viral clearance.

In Figure 7, the increase in the killing rate k_*r* of follicular CD8⁺ T cells is associated with a reduction in all infected populations, indicating an enhanced immune-mediated clearance mechanism. In particular, a higher decrease is observed in the infected T_*fr* cell population, reflecting its higher susceptibility to cytotoxic elimination, whereas only a relatively small reduction is observed in the infected T_*fh* cells and the viral load. The location of the peaks remains unchanged, indicating that the killing rate primarily affects the magnitude of infection rather than its timing. This suggests that k_*r* modulates the intensity of the immune response without altering the intrinsic progression speed of the infection within the system.

In Figure 8, the increase in the lysis rate δ_*r* of T follicular regulatory (T_*fr*) cells is associated with a growth of infection dynamics in the infected follicular helper cell population I_*fh* as well as the viral load. This behavior suggests that increased viral production by regulatory cells facilitates viral replication and spread. The growth in the infections is more pronounced in the viral population than the infected T_*fh* cell population. The higher infection peaks are observed in the infected, T_*fh*, and viral populations. The peak reduction is observed in the infected T_*fh* population. These results highlight a shift in infection growth toward regulatory cells during lysis of infected helper cells, overall amplifying viral persistence in the germinal center.

In Figure 9, the increase in the expansion rate ρ_*r* of healthy T follicular regulatory (T_*fr*) cells is associated with a reduction in all infected classes, indicating that enhanced proliferation of regulatory cells strengthens the control and reduces infection propagation across the system. The peak locations of the infected populations remain unchanged as ρ_*r* increases, indicating that the intrinsic timing of the infection dynamics is unaffected by changes in the expansion rate. However, the amplitudes of these peaks decrease as ρ_*r* increases, reflecting reduced infection severity rather than a shift in the disease progression timeline. Overall, these results highlight the protective role of T_*fr* expansion in reducing infection burden while preserving the underlying dynamical structure of the system.

4.2.3 The effect of follicular CD8⁺ T cells parameters on I_*fh*, I_*fr*, and V

In Figure 10, the increase in the average number N_*h* of viral particles produced per infected T follicular helper (T_*fh*) cell is associated with an increase in the infection dynamics, as shown in the increase in the infected follicular helper cell population I_*fh*, the viral load, and, the infected follicular regulatory cell population I_*fr*. Thus, a higher viral production rate from infected T_*fh* cells enhances the infectivity and accelerates disease progression within the system. There is also a backward shift in the peaks of all infected classes as N_*h* increases, indicating an earlier onset of high infection levels. The corresponding increase in the peaks reflects a higher infection burden and intensified viral replication. A similar qualitative trend is observed in Figure 11 as the average number N_*r* of viral particles produced by infected T follicular regulatory (T_*fr*) cells increases. In this case as well, higher values of N_*r* increase the infected T_*fh* and T_*fr* cell populations and the viral load, thereby reinforcing the role of infected T_*fr* cells as an additional source of viral amplification. The results confirm that enhanced viral production from infected T_*fh*, T_*fr*, or both cells accelerates infection dynamics and intensifies the overall disease burden.

In Figure 12, the increase in the effect induced by the presence of infected T follicular regulatory (T_*fr*) cells, quantified through the parameter δ_*c*. The increase in δ_*c* is associated with an increase in the infected T follicular helper cell population I_*fh*, the infected T follicular regulatory cell population I_*fr*, as well as the viral load. This suggests that the increased effects induced by the infected T_*fr* on the immune response cytotoxic cells enhance viral persistence and growth of infected compartments across the system. The peaks of infection occur at the same time but rise in magnitude as δ_*c* increases, reflecting a higher infection burden. These results indicate that reducing immune cell growth alters their killing capacity and, consequently, contributes to the persistence and spread of infection in the germinal centers.

Remark 4. For each consecutive pair of parameter values a1=αd(i) and a2=αd(i+1), the corresponding numerical solutions (I_fh*, I_*fr*, V) are extracted from the simulation results. The sensitivity of the system with respect to variations in α_*d* is then quantified by computing the maximum absolute differences*

max|Ifh(a1)-Ifh(a2)|, max|Ifr(a1)-Ifr(a2)|, max|V(a1)-V(a2)|.

These differences are small (e.g. 10⁻¹⁰ or smaller), hence the green and red curves are literally on top of each other.

5 Discussion and conclusion

We developed and analyzed a mechanistic in-host mathematical model of germinal center dynamics to quantitatively explore how the coupled interactions among T follicular helper, T follicular regulatory (Tfr), and follicular CD8+ (fCD8+) T cells control HIV replication, immune control, and the creation, maintenance, and stability of the viral reservoir. We conducted mathematical analyses of the model and case-scenario analyses using numerical simulations to better understand this interplay.

T follicular helper (Tfh), T follicular regulatory (Tfr), and follicular CD8+ (fCD8+) T cells interact closely to sustain HIV persistence in germinal centers, as our mechanistic model illustrates. As the structural backbone of the viral reservoir, Tfh cells can collapse the linked system toward viral elimination if their support is compromised. When fCD8+ cells eliminate infected cells by cytotoxic means, the infectious burden is reduced. However, infection is not entirely cleared, underscoring the nonlinear feedback loop that sustains infection and the germinal center's limited immunological privilege. The number of infected helper cells and viral load might increase as the lysis rate of infected Tfh cells increases. This illustrates the amplification brought by viral release and the compartmental redistribution of infection.

Tfr cells primarily regulate the amplitude and timing of infection rather than structural persistence. Increased Tfr deregulation reduces the number of infected Tfr cells and viral load while slowly increasing Tfh infection, resulting in earlier infection peaks. Improved healthy Tfr growth, on the other hand, reduces all infected populations without altering peak infection levels. Tfr infection is significantly decreased by cytotoxic targeting of infected Tfr cells. The effects of the Tfr saturation constant result in transient peaks without altering equilibrium levels, indicating a loss of early regulatory efficiency and little effect on Tfh-driven viral production. Together, the findings show that the germinal center viral pool is a nonlinear, coupled network in which fCD8+ cells regulate viral burden, Tfr cells influence infection growth, and Tfh cells maintain persistence. The complex interaction presents a challenge: achieving successful eradication requires simultaneously disrupting the regulatory balance and the helper-driven viral support, rather than focusing on a single immunological component in isolation.

Tfh cells form the structural backbone of the germinal center viral reservoir, which can sustain HIV replication and infection. Cytotoxic killing alone reduced infection without completely clearing the reservoir, indicating partial immune function in follicular cells. Compromising the Tfh function can destabilize the linked system, resulting in decreased viral load and possible clearance. The model results indicated that even minor increases in cell death can paradoxically accelerate short-term viral production rather than enhance control, as increasing the lysis rate of infected Tfh cells may unexpectedly increase viral load. Regulatory feedback from Tfr cells and cytotoxic pressure from fCD8+ cells modulate the amplitude and timing of infection; however, persistence largely depends on Tfh-driven viral production. To successfully disrupt the germinal center HIV reservoir, one requires strategies that target the structural support provided by Tfh cells within the coupled immune network rather than focusing solely on isolated single-cell populations.

The model results support the interconnected interplay among immune cells that underlies HIV persistence despite robust immune responses and successful systemic antiretroviral therapy. Tfh-driven viral generation, Tfr-mediated control, and fCD8+-mediated cytotoxic pressure interact nonlinearly within the germinal center reservoir, which functions as a dynamically stabilized network. Successful reservoir eradication will likely require targeted interventions that simultaneously address helper-driven viral amplification, regulatory feedback, and cytotoxicity.

Overall, this study shows that linked immunological dynamics across several T-cell subsets influence germinal center HIV persistence. The model expands understanding of how the immunological network topology regulates viral replication, immune containment, and long-term reservoir stability in the follicular setting through quantitative analysis of these interactions.

A quantitative, mechanistic model of HIV persistence in the germinal center, driven by a coupled Tfh-Tfr-fCD8+ system, demonstrated through parameter exploration that no single infected-cell population is responsible for the viral burden, nor can the burden be eliminated by modifying a single immunological pathway. Instead, reservoir stability arose from nonlinear feedback among helper-driven viral production, regulatory stabilization, and cytotoxic immune pressure. Hence, mechanistic immunology and mathematical analysis can be used to improve our understanding of how the immune network architecture governs HIV replication, immunological containment, and long-term reservoir structure. Therapeutic interventions can be strategically designed using insights from this theoretical foundation.

Future work will focus on incorporating uncertainty quantification to assess how parameter variability and biological heterogeneity affect model predictions. These analyses will improve the model's robustness, reliability, and clinical relevance for studying HIV persistence and informing therapeutic strategies.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Feng Rao, Nanjing Tech University, China

Reviewed by: Ryan Zurakowski, University of Delaware, United States

Changjing Zhuge, Beijing University of Technology, China