Delay-driven infection resurgence and intervention timing in within-host HIV infection

1 Introduction

Human immunodeficiency virus (HIV) infection remains a persistent global health challenge, with complex within-host dynamics involving susceptible CD4⁺ T cells, productively infected cells, free virus, and immune responses [1, 2]. Antiretroviral therapy (ART) has substantially improved the prognosis of people living with HIV by suppressing viral replication, reducing viral load, and supporting immune recovery [3, 4].

Nevertheless, durable viral control remains difficult because HIV infection is shaped by non-linear viral–immune interactions, latent reservoirs, heterogeneous treatment response, and the possibility of viral rebound after treatment interruption or insufficient suppression [5–7]. Understanding the mechanisms that determine suppression, persistence, immune-control failure, and post-control infection resurgence is therefore important for designing intervention strategies that remain effective under biological, immunological, and pharmacological delays.

Mathematical models provide a useful framework for studying HIV dynamics within the host [1, 2, 8]. Classical viral dynamic models describe the infection of target cells, the production of virus by infected cells, and the clearance of free virions [1, 3, 4]. Extensions of these models have incorporated immune responses, drug therapy, latent infection, immune escape, and heterogeneous disease progression [5, 6, 9–13]. These studies have clarified important threshold mechanisms associated with infection clearance, persistence, immune activation, and treatment response. However, many formulations assume that infection, immune activation, and treatment effects occur instantaneously. This assumption may obscure clinically relevant transient behavior, particularly when biological and pharmacological processes operate on different time scales.

Time delays are particularly important in HIV dynamics. After viral entry into a susceptible CD4⁺ T cell, there is an intracellular eclipse period before the cell becomes productively infected [14–16]. Similarly, cytotoxic T lymphocyte (CTL) activation does not occur immediately after antigenic stimulation, because immune recognition, activation, proliferation, and functional expansion require time [9–12]. Treatment effects may also be delayed because drug action, pharmacodynamic response, adherence patterns, and intervention uptake require time before altering viral replication or effective susceptibility [4, 17, 18]. Ignoring these delays may lead to an overly optimistic assessment of infection control, since delayed immune or treatment responses can allow infected-cell burden and viral load to increase before suppression becomes effective.

Delay differential equation models have therefore become important tools for representing biological lags in viral infection processes [14–16, 19]. In HIV modeling, delays have been used to describe intracellular infection stages, immune activation, treatment response, and intervention implementation [15, 16, 19]. Such models show that delays can change stability properties, generate sustained oscillations, and shift critical thresholds for infection persistence [8, 16, 19, 20]. From an infectious-disease-dynamics perspective, these oscillations should be interpreted as theoretical rebound-like or unstable immune-control regimes rather than direct patient-calibrated viral rebound, unless additional clinical mechanisms such as latent reservoirs, stochastic reactivation, adherence variability, and patient-specific viral-load data are explicitly included.

Despite this progress, several modeling issues remain important. First, delay terms should be biologically justified. If a delayed infection or immune-activation term is included, it should account for possible loss during the waiting period; otherwise, the model implicitly assumes that all infected cells or immune signals survive unchanged throughout the delay interval [14]. Second, prevention or protective effects should be distinguished from direct removal of susceptible cells. A protective intervention usually reduces effective susceptibility or the probability of successful infection rather than physically transferring CD4⁺ T cells out of the susceptible compartment [21, 22]. Third, immune activation, immune killing, prevention efficacy, and treatment action represent distinct biological mechanisms and should be parameterized separately. Finally, optimal-control results should be interpreted with respect to treatment-delay effects, robustness to cost weights, and post-intervention dynamics, because short-term suppression over a finite horizon does not necessarily imply long-term prevention of post-control resurgence [7, 18].

The present study addresses these issues by developing a delay-structured within-host HIV model that explicitly separates biological state delays from delayed intervention effects. The baseline model includes susceptible CD4⁺ T cells, productively infected cells, CTLs, free virus, and a normalized protective-response variable. The intracellular infection delay represents the time from viral entry to productive infection, whereas the immune-activation delay represents the time required for CTL expansion following antigenic stimulation. Survival factors are incorporated into the delayed terms to account for loss during these waiting periods. The protective variable reduces the effective infection rate rather than removing susceptible cells, thereby providing a clearer biological interpretation of prevention or prophylactic effects.

The analysis focuses on how delays reshape viral suppression, immune control, and delay-dependent infection dynamics. We derive biologically relevant equilibria and identify threshold conditions for infection clearance, persistent infection, and CTL-mediated control. Local stability and Hopf bifurcation analyses are used to determine when increasing delays destabilize steady states and generate oscillatory infection dynamics [8, 16, 19]. Numerical simulations then illustrate how short, intermediate, and long delays produce distinct within-host regimes, ranging from stable suppression to delay-dependent transient oscillations. A global sensitivity analysis is also performed to identify the parameters that most strongly influence peak infected-cell burden and infection-resurgence behavior [23, 24].

We further extend the baseline model to include delayed treatment and prevention controls. This optimal-control formulation is designed to evaluate how intervention timing affects viral suppression under pharmacological and implementation delays [17, 18]. In contrast to formulations that report only the optimal trajectory over a finite control horizon, we emphasize robustness with respect to cost weights and the post-control behavior after the terminal time. This is essential because a control strategy that appears effective during the prescribed treatment interval may still permit post-control resurgence once the intervention is relaxed or discontinued [7].

The main novelties of this study are as follows. First, the model separates intracellular infection delay, immune-activation delay, and treatment-response delay within a unified within-host HIV framework. Second, survival factors are incorporated into the delayed infection and immune-response terms, giving the delays a clearer biological interpretation. Third, the protective response is modeled as a reduction in effective susceptibility rather than as direct removal of susceptible cells. Fourth, the delay-dependent stability and Hopf bifurcation analysis is linked to rebound-like oscillatory regimes. Fifth, the delayed optimal-control formulation is supplemented with robustness and post-control resurgence analyses, allowing intervention timing to be evaluated beyond the finite treatment horizon.

The main contribution of this study is therefore not merely the introduction of another delayed HIV model, but the construction of a delay-aware framework linking biological lags, immune activation, intervention timing, and post-control resurgence risk in within-host HIV dynamics. By separating intracellular, immune, and treatment-response delays, the model provides a clearer basis for interpreting how delayed processes influence stability thresholds and viral resurgence. The results offer a quantitative theoretical framework for assessing when early intervention is necessary, when delayed immune feedback may destabilize infection control, and how robust treatment strategies should be evaluated under biologically plausible delay structures.

The remainder of the paper is organized as follows. Section 2 formulates the delay-structured within-host HIV model, states the modeling assumptions, establishes well-posedness, and presents the non-dimensional form. Section 3 derives the equilibria, threshold quantities, and local stability conditions. Section 4 analyzes delay-induced Hopf bifurcations and the emergence of oscillatory viral dynamics. Section 5 presents numerical simulations and global sensitivity results. Section 6 introduces delayed optimal controls for treatment and prevention, together with robustness and post-control resurgence analyses. Section 7 concludes with implications for intervention timing and future extensions.

2 Delay-structured within-host HIV model

We consider a within-host HIV model describing the interaction among susceptible CD4⁺ T cells, productively infected cells, cytotoxic T lymphocytes (CTLs), free virus, and a time-dependent protective-response variable. Let S(t), I(t), T(t), V(t), and P(t) denote, respectively, the densities of susceptible CD4⁺ T cells, productively infected cells, CTLs, free virus, and the normalized level of prevention or protective response at time t. The variables S(t), I(t), and V(t) follow the standard structure of within-host HIV dynamics, in which target cells are infected by virus, infected cells produce free virus, and free virus is cleared from the system [1–3]. The CTL compartment is included to represent cellular immune-mediated antiviral activity, which is known to play an important role in HIV pathogenesis and immune control [9–13].

The variable P(t) is interpreted as an aggregate protective factor that reduces effective susceptibility. It may represent prophylactic intervention, behavioral protection, or other prevention-related mechanisms that lower the probability of successful infection. Importantly, P(t) is not treated as a separate cell population. Instead, it acts by reducing the effective infection rate, consistent with prevention-oriented HIV modeling in which interventions modify transmission or susceptibility rather than physically removing susceptible cells from the system [21, 22]. In this study, P(t) is interpreted as a normalized aggregate protective-response variable rather than as a direct biological cell population. It may represent prevention-like or prophylaxis-like effects, adherence-supported protection, treatment-supported reduction in susceptibility, or other mechanisms that reduce the probability of successful viral entry. Since P(t) is normalized, 0 ≤ P(t) ≤ 1, where P(t) = 0 corresponds to the absence of a protective effect and P(t) = 1 corresponds to the maximal protective level represented in the model. Thus, the factor 1−εP(t) reduces the effective infection rate, rather than directly removing susceptible CD4⁺ T cells from the system. This interpretation provides a clearer biological meaning for the protective-response variable and avoids introducing an additional susceptible-cell compartment.

Before introducing the delayed system, we first present the corresponding instantaneous model. This non-delayed formulation serves as a baseline reference for identifying how biological delays alter the infection and immune-response dynamics

dS(t)dt=Λ-μSS(t)-β(1-εP(t))S(t)V(t),dI(t)dt=β(1-εP(t))S(t)V(t)-dI(t)-γT(t)I(t),dT(t)dt=ρV(t)T(t)-μTT(t),dV(t)dt=pI(t)-cV(t),dP(t)dt=q(1-P(t))-μPP(t).(1)

Model (1) assumes that infection, CTL stimulation, and protective effects act instantaneously. In contrast, the delayed model below incorporates two biological state delays: the intracellular infection delay τ_*s* and the immune-activation delay τ_*I*. These delays represent the time required for newly infected cells to become productively infected and for CTLs to respond to antigenic stimulation, respectively.

We now state the assumptions for the delayed model.

(A1) Susceptible CD4⁺ T cells are produced at a constant rate Λand die naturally at rate μ_*S*S(t).

(A2) Infection of susceptible cells by free virus follows a mass-action incidence term βS(t)V(t). The protective variable P(t) reduces the effective infection rate through the factor 1−εP(t), where 0 ≤ εP(t) ≤ 1. Thus, prevention reduces the probability of successful infection but does not remove susceptible cells directly.

(A3) Newly infected cells require an intracellular eclipse period before becoming productively infected. This biological delay is denoted by τ_*s*. During this waiting period, infected cells may die or fail to survive to productive infection; therefore, the delayed infection term is multiplied by the survival factor e-mτs, where m>0 is the loss rate during the eclipse phase [14–16].

(A4) Productively infected cells die at rate dI(t) and are additionally removed by CTL-mediated killing at rate γT(t)I(t).

(A5) Free virus is produced by productively infected cells at rate pI(t) and cleared at rate cV(t).

(A6) CTLs are stimulated by viral antigen after an immune-activation delay τ_*I*. The CTL expansion term is represented by ρV(t-τI)T(t-τI)e-mIτI, where m_*I*>0 accounts for loss during the activation interval. CTLs decay naturally at rate μ_*T*T(t). The use of a delayed immune-response term is motivated by the fact that CTL recognition, activation, proliferation, and functional expansion are not instantaneous biological processes [9–12].

(A7) The protective variable P(t) increases through a bounded adoption or activation process with maximum normalized level one, represented by q(1−P(t)), and wanes at rate μ_*P*P(t).

The modeling assumptions are summarized in Table 1 to improve readability and to clarify the biological role of each assumption.

Table 1. Summary of modeling assumptions and their biological interpretation.

Table 1 shows that the delays and protective response are introduced for biological reasons, not only for mathematical generality. In particular, τ_*s* represents the intracellular waiting time before productive infection, τ_*I* represents the immune-activation waiting time before CTL expansion, and P(t) represents a normalized protection level that reduces effective susceptibility.

Under these assumptions, the state-delayed within-host HIV model is given by

dS(t)dt=Λ-μSS(t)-β(1-εP(t))S(t)V(t),dI(t)dt=βe-mτs(1-εP(t-τs))S(t-τs)V(t-τs)-dI(t) -γT(t)I(t),dT(t)dt=ρe-mIτIV(t-τI)T(t-τI)-μTT(t),dV(t)dt=pI(t)-cV(t),dP(t)dt=q(1-P(t))-μPP(t).(2)

This formulation differs from the instantaneous model (1) in two important respects. First, the infection process includes the intracellular delay τ_*s*, which represents the time between viral entry and the appearance of productively infected cells. Second, CTL expansion is not assumed to occur immediately after antigenic stimulation; instead, the immune-activation delay τ_*I* represents the time required for immune recognition and expansion. The survival factors e-mτs and e-mIτI account for losses during the corresponding delay intervals. This avoids the unrealistic assumption that all delayed infected cells or immune signals survive unchanged throughout the waiting period [14–16].

The parameters ε, ρ, and γ represent distinct mechanisms. The parameter ε measures the efficacy of protection in reducing effective susceptibility, ρ measures antigen-driven CTL stimulation, and γ measures CTL-mediated killing of productively infected cells. This separation prevents the same parameter from being used for biologically different processes.

2.1 Initial history and well-posedness

Let

τmax=max{τs,τI}.

The initial history functions are prescribed on the interval [−τ_max, 0 ] as given in Equation 3.

S(θ)=ϕ1(θ), I(θ)=ϕ2(θ), T(θ)=ϕ3(θ), V(θ)=ϕ4(θ), P(θ)=ϕ5(θ),(3)

where

ϕi∈C([-τmax,0],ℝ+), i=1,…,5,

and

0≤ϕ5(θ)≤1, θ∈[-τmax,0].

Since the right-hand side of system (2) is locally Lipschitz continuous with respect to the state variables on the phase space C([-τmax,0],ℝ+5), standard existence and uniqueness results for delay differential equations ensure the existence of a unique local solution for each admissible initial history.

The non-negative orthant is positively invariant. Indeed, if S(t) = 0, then

S′(t)=Λ>0.

If I(t) = 0, then

I′(t)=βe-mτs(1-εP(t-τs))S(t-τs)V(t-τs)≥0.

If T(t) = 0, then

T′(t)=ρe-mIτIV(t-τI)T(t-τI)≥0.

If V(t) = 0, then

V′(t)=pI(t)≥0,

and if P(t) = 0, then

P′(t)=q>0.

Therefore, solutions starting from non-negative initial histories remain non-negative for all t≥0 on their interval of existence.

The protective variable is also bounded above. If P(t) = 1, then

P′(t)=q(1−1)−μP=−μP<0.

Hence, the interval 0 ≤ P(t) ≤ 1 is positively invariant. The biologically feasible region is therefore defined in Equation 4.

Ω={(S,I,T,V,P)∈ℝ+5: 0≤P≤1}.(4)

2.2 Baseline model without control delay

System (2) describes the baseline state-delayed infection dynamics. At this stage, no treatment-control delay is included. The biological delays τ_*s* and τ_*I* represent intrinsic within-host processes, namely intracellular infection progression and immune activation. The pharmacological or implementation delay associated with treatment response is introduced later in the optimal-control formulation.

This separation is important. The biological delays τ_*s* and τ_*I* determine the intrinsic stability, persistence, and rebound-like oscillatory structure of the uncontrolled infection process. By contrast, the treatment-response delay determines how quickly an intervention can modify viral replication or effective susceptibility. Thus, throughout the baseline analysis, the delays are biological state delays only. The control delay, denoted later by τ_*c*, is reserved for the delayed intervention model.

2.3 Model schematic and biological interpretation

Figure 1 summarizes the structure of system (2). Susceptible CD4⁺ T cells are recruited at rate Λ, die naturally at rate μ_*S*S, and become infected through effective viral contact at rate β(1−εP)SV. The protective response P reduces the probability of successful infection but does not remove susceptible cells from the system. After an intracellular delay τ_*s*, a surviving fraction e-mτs of newly infected cells enters the productively infected class I. Productively infected cells are removed by natural infected-cell death at rate dI and by CTL-mediated killing at rate γTI.

Productively infected cells produce free virus at rate pI, while free virus is cleared at rate cV. Viral antigen stimulates CTL expansion after an immune-activation delay τ_*I*, with survival factor e-mIτI. CTLs decay at rate μ_*T*T. Finally, the protective variable P increases through the bounded adoption term q(1−P) and wanes at rate μ_*P*P.

This schematic is placed before the stability analysis so that the biological role of each model term is clear before the equilibria, thresholds, and characteristic equations are derived.

2.4 Non-dimensional form

For analytical convenience, we introduce the dimensionless variables

s=SS0, i=II0, z=TT0, v=VV0, π=P, t^=dt,(5)

where

S0=ΛμS, I0=Λd, V0=pI0c, T0=dγ.(6)

The rescaled delays are

τ^s=dτs, τ^I=dτI.(7)

After substituting (5)–(7) into system (2), rescaling time by t^=dt, dropping hats for notational simplicity, and defining the dimensionless parameters

a=μSd,b=βV0d,σ=md,σI=mId,r=ρV0d,ν=μTd,κ=pI0dV0,χ=cd,ω=qd,ξ=μPd,(8)

we obtain

ds(t)dt=a-as(t)-b(1-επ(t))s(t)v(t),di(t)dt=be-στs(1-επ(t-τs))s(t-τs)v(t-τs)-i(t)-z(t)i(t),dz(t)dt=re-σIτIv(t-τI)z(t-τI)-νz(t),dv(t)dt=κi(t)-χv(t),dπ(t)dt=ω(1-π(t))-ξπ(t).(9)

The non-dimensionalization reduces the number of independent quantities and highlights the mechanisms most directly associated with viral persistence, immune control, and delay-dependent infection dynamics: effective infectivity b, protection efficacy ε, CTL activation r, CTL decay ν, viral clearance χ, and the delay-survival products e-στs and e-σIτI. In particular, the products στ_*s* and σ_*I*τ_*I* quantify the effective loss accumulated during the intracellular infection and immune-activation waiting periods. In the subsequent analysis, we use the dimensionless system (9) unless otherwise stated.

3 Equilibria and threshold quantities

In this section, we derive the biologically relevant equilibria and threshold quantities of the dimensionless system (9). Since the subsequent stability and bifurcation analysis is based on the non-dimensional model, all quantities in this section are expressed in terms of system (9) unless otherwise stated. For convenience, we define the auxiliary quantities in Equations 10 and 11.

π-=ωω+ξ, Θ=1-επ-, Es=e-στs, EI=e-σIτI.(10)

Here, π- is the equilibrium level of the protective response, Θ is the remaining effective susceptibility under protection, and E_*s* and E_*I* are the survival factors associated with the intracellular infection delay and immune-activation delay, respectively. Throughout the analysis, we assume

0<Θ≤1,(11)

so that the protective response reduces but does not completely eliminate effective infection.

To make the derivation of the characteristic equations transparent, we first state the general linearization structure. Let

X(t)=(s(t),i(t),z(t),v(t),π(t))T

and let

E=(s,i,z,v*,π^)

be any equilibrium of system (9). For a perturbation Y(t) = X(t)−E^*, the linearized delay system can be written in the form

Y′(t)=A0(E)Y(t)+As(E)Y(t-τs)+AI(E*)Y(t-τI),

where A0(E), As(E), and AI(E*) are the Jacobian matrices associated with the instantaneous terms, the intracellular-delay terms, and the immune-activation-delay terms, respectively. They are given by

A0(E)=(-a-bΘv00-bΘsεbsv0-(1+z)-i0000-ν000κ0-χ00000-(ω+ξ)),As(E)=(00000bEsΘv00bEsΘs-εbEssv000000000000000),

and

AI(E)=(000000000000rEIvrEIz*00000000000).

Substituting a trial solution Y(t) = ηe^λ*t*, where η≠0, into the linearized delay system gives

[λI-A0(E)-As(E)e-λτs-AI(E*)e-λτI]η=0.

Therefore, the characteristic equation at E^* is obtained from

det[λI-A0(E)-As(E)e-λτs-AI(E*)e-λτI]=0.

This matrix formulation explains how the characteristic equations below are obtained from the linearized delay system.

3.1 Infection-free equilibrium

The infection-free equilibrium is obtained by setting

i=v=z=0.

From the first and fifth equations of system (9), we obtain

s=1, π=π-.

Therefore, the infection-free equilibrium is given in Equation 12.

E0=(1,0,0,0,π-).(12)

This equilibrium represents a state in which no productively infected cells, free virus, or CTL response are present, while the susceptible CD4⁺ T-cell population and protective response remain at their baseline equilibrium levels.

At E₀, the general Jacobian structure above reduces to a block form. The variables s, z, and π contribute the eigenvalues −a, −ν, and −(ω+ξ), respectively. The remaining stability condition is determined by the infected-cell/free-virus subsystem. Linearizing the i-and v-equations at E₀ gives the subsystem in Equation 13.

di(t)dt=bEsΘv(t-τs)-i(t), dv(t)dt=κi(t)-χv(t).(13)

Using the trial solutions i(t)=I0eλt and v(t)=V0eλt, we obtain

(λ+1)I0-bEsΘe-λτsV0=0,-κI0+(λ+χ)V0=0.

The determinant condition for non-trivial (I₀, V₀) yields the characteristic equation

(λ+1)(λ+χ)-bκEsΘe-λτs=0.(14)

This motivates the basic infection threshold

R0=bκEsΘχ=bκe-στs(1-επ-)χ.(15)

The quantity R0 measures the expected infection-generating potential of the virus in the presence of protection and intracellular delay. The factor e-στs reduces the effective production of productively infected cells by accounting for loss during the intracellular waiting period. Thus, increasing τ_*s*, increasing σ, or increasing protection efficacy ε decreases R0.

Theorem 3.1.

If R0<1, then the infection-free equilibrium E₀ is locally asymptotically stable for all τ_*s*≥0. If R0>1, then E₀ is unstable.

Proof.

The remaining eigenvalues associated with s, z, and π at E₀ are negative, namely −a, −ν, and −(ω+ξ). Hence, the stability of E₀ is determined by the roots of Equation 14.

Suppose that R0<1 and that Equation 14 has a root with Re(λ)≥0. Taking moduli in (14) gives

∣(λ+1)(λ+χ)∣=bκEsΘ∣e-λτs∣.(16)

Since Re(λ)≥0, we have

∣λ+1∣≥1,∣λ+χ∣≥χ,∣e-λτs∣≤1.

Therefore,

χ≤∣(λ+1)(λ+χ)∣≤bκEsΘ=χR0<χ,

which is impossible. Hence all roots of Equation 14 have negative real parts, and E₀ is locally asymptotically stable.

If R0>1, then evaluating the left-hand side of Equation 14 at λ = 0 gives

χ-bκEsΘ=χ(1-R0)<0,

whereas the characteristic function tends to +∞ as λarrow+∞. Hence, by continuity, Equation 14 has a positive real root. Therefore, E₀ is unstable.

3.2 CTL-free endemic equilibrium

When R0>1, the infection may persist even in the absence of a sustained CTL response. Setting z = 0, π=π-, and solving the steady-state equations of system (9) yields the CTL-free endemic equilibrium

E1=(s1,i1,0,v1,π-),(17)

where

s1=1R0,(18)

v1=a(R0-1)bΘ,(19)

and

i1=χv1κ=aχ(R0-1)bκΘ.(20)

Thus, E₁ exists biologically if and only if

R0>1.(21)

The equilibrium E₁ represents persistent infection without sustained CTL activation. In this regime, viral replication is strong enough to maintain infected cells and free virus, but the CTL response does not persist at a positive equilibrium level.

At E₁, the characteristic equation obtained from the general Jacobian matrices factorizes into a CTL-direction component and a remaining infection-virus-protection component. In particular, the CTL-direction perturbation satisfies

dz(t)dt=rEIv1z(t-τI)-νz(t),

which gives the characteristic equation

λ+ν-rEIv1e-λτI=0.

The remaining characteristic component is obtained from the (s, i, v, π)-subsystem. Its determinant form is

(λ+ω+ξ)[(λ+a+bΘv1)(λ+1)(λ+χ)−bκEsΘ{s1(λ+a+bΘv1)−bΘs1v1}e−λτs]=0.

This factorization shows explicitly how the CTL-invasion condition is separated from the infection-virus-protection subsystem near the CTL-free endemic equilibrium. It also makes the characteristic component consistent with the Hopf analysis in Section 4.2.

3.3 CTL-activation threshold

To determine whether the CTL population can invade the CTL-free endemic equilibrium, we linearize the z-equation around E₁. Since v(t) = v₁ at equilibrium, the linearized CTL equation is

dz(t)dt=rEIv1z(t-τI)-νz(t).(22)

The corresponding characteristic equation is

λ+ν-rEIv1e-λτI=0.(23)

This motivates the CTL-activation threshold

RT=rEIv1ν=re-σIτIν·a(R0-1)bΘ.(24)

If RT<1, then the CTL response cannot persist near E₁. If RT>1, then CTL activation can invade the CTL-free endemic state. Thus, RT measures the ability of immune activation to overcome CTL decay in the persistent-infection regime.

Theorem 3.2.

The CTL-free endemic equilibrium E₁ is locally stable in the CTL direction if RT<1 and unstable in the CTL direction if RT> 1.

Proof.

The CTL characteristic equation is given by (23). Suppose that RT<1 and that (23) has a root with Reλ≥0. Taking moduli gives

∣λ+ν∣=rEIv1∣e-λτI∣≤rEIv1.(25)

Since Re(λ)≥0, we have

∣λ+ν∣≥ν.

ν≤rEIv1,

which contradicts RT<1. Hence all CTL-direction roots have negative real parts when RT< 1.

If RT>1, then the left-hand side of Equation 23 evaluated at λ = 0 is

ν-rEIv1=ν(1-RT)<0,

while it tends to +∞ as λarrow+∞. Therefore, Equation 23 has a positive real root, and the CTL-free equilibrium is unstable in the CTL direction.

3.4 CTL-activated endemic equilibrium

When both infection and CTL activation persist, system (9) admits a CTL-activated endemic equilibrium

E2=(s2,i2,z2,v2,π-),(26)

with all components positive. From the CTL equation, positivity of z₂ implies

rEIv2=ν,

and hence

v2=νrEI=νeσIτIr.(27)

v2=κi2χ,

so that

i2=χv2κ=χνeσIτIκr.(28)

The susceptible-cell equation gives

a-as2-bΘs2v2=0,

and therefore

s2=aa+bΘv2.(29)

Finally, the infected-cell equation gives

bEsΘs2v2-(1+z2)i2=0.

Using v₂ = κi₂/χ, we obtain

z2=R0s2-1.(30)

Therefore, E₂ is biologically feasible if

z2>0,

R0s2>1.(31)

This condition is equivalent to

v2<v1,

which is precisely

RT>1.(32)

Thus, the CTL-activated endemic equilibrium exists if

R0>1 and RT>1.(33)

The equilibrium E₂ represents a persistent infection state in which viral replication is balanced by CTL-mediated immune activity. The condition RT>1 indicates that the viral level at the CTL-free endemic state is sufficiently high to sustain CTL activation.

Biologically, this means that infection persists, but the immune response is strong enough to remain active and exert continuous pressure on productively infected cells.

In contrast to the CTL-free endemic equilibrium E₁, the equilibrium E₂ includes a positive CTL component. Hence, perturbations in infected cells, free virus, and CTLs are dynamically coupled near E₂. This coupling is important for the subsequent delay-dependent stability and Hopf bifurcation analysis, because both the intracellular infection delay τ_*s* and the immune-activation delay τ_*I* can influence the local eigenvalue structure.

3.5 Characteristic equation at a general equilibrium

Let

E=(s,i,z,v*,π-)

be any biologically feasible equilibrium of system (9). Consider perturbations

s(t)=s+x1(t), i(t)=i+x2(t),z(t)=z+x3(t),v(t) =v+x4(t),π(t)=π-+x5(t).

The linearized system can be written in the form

x′(t)=A0x(t)+Asx(t-τs)+AIx(t-τI),(34)

where A₀ contains the instantaneous terms, A_*s* contains the delayed infection terms, and A_*I* contains the delayed CTL-activation terms.

More explicitly, the Jacobian matrix corresponding to the instantaneous terms is

A0=(-a-bΘv00-bΘsεbsv0-(1+z)-i0000-ν000κ0-χ00000-(ω+ξ)).

The Jacobian matrix associated with the intracellular infection delay τ_*s* is

As=(00000bEsΘv00bEsΘs-εbEssv000000000000000).

The Jacobian matrix associated with the immune-activation delay τ_*I* is

AI=(000000000000rEIvrEIz00000000000).

These matrices are obtained by taking partial derivatives of the right-hand side of system (9) with respect to the current state x(t), the τ_*s*-delayed state x(t−τ_*s*), and the τ_*I*-delayed state x(t−τ_*I*), respectively. Thus, A₀ describes the instantaneous local interactions, A_*s* describes the delayed infection feedback, and A_*I* describes the delayed CTL-activation feedback.

To derive the characteristic equation, we substitute

x(t)=ηeλt,η≠0,

into the linearized system (34). This gives

ληeλt=A0ηeλt+Asηeλ(t-τs)+AIηeλ(t-τI).

Dividing by e^λ*t*, we obtain

[λI-A0-Ase-λτs-AIe-λτI]η=0.

For a non-trivial perturbation vector η, the determinant must vanish. Hence, the corresponding characteristic equation is

Δ(λ;τs,τI)=det[λI-A0-Ase-λτs-AIe-λτI]=0.(35)

Equation 35 provides the basis for the delay-dependent stability and Hopf bifurcation analysis in the next section. In particular, purely imaginary roots λ = iω, ω>0, identify possible stability-switching points at which increasing biological delays may destabilize an equilibrium and generate delay-dependent oscillatory infection dynamics. Thus, the characteristic equation links the threshold analysis in Section 3 with the bifurcation analysis in Section 4.

This derivation also clarifies why different equilibria lead to different characteristic structures. At the infection-free equilibrium E₀, the infected-cell and viral-load variables form the critical subsystem. At the CTL-free endemic equilibrium E₁, the CTL direction can be separated from the remaining infection-virus-protection dynamics. At the CTL-activated endemic equilibrium E₂, infected cells, free virus, and CTLs are coupled, and the characteristic equation generally contains both delayed factors e-λτs and e- λτI.

4 Delay-dependent stability and Hopf bifurcation analysis

In this section, we analyze how the biological delays τ_*s* and τ_*I* affect the local stability of the biologically relevant equilibria derived in Section 3. The analysis is based on the characteristic Equation 35. Since delay differential equations may lose stability through the crossing of characteristic roots across the imaginary axis, we focus on the existence of purely imaginary roots and the associated transversality condition. The resulting oscillations are interpreted as delay-dependent infection oscillations or rebound-like fluctuations within the theoretical model, rather than as patient-calibrated viral rebound.

In response to the need for a more complete Hopf bifurcation analysis, we explicitly verify the eigenvalue crossing condition and state the normal-form quantities used to determine the direction, stability, and period of the bifurcating periodic solutions. The Hopf analysis is carried out by treating the delay as the bifurcation parameter and by following the standard center-manifold and normal-form framework for retarded functional differential equations.

4.1 Stability of the infection-free equilibrium

At the infection-free equilibrium

E0=(1,0,0,0,π-),

the characteristic equation decomposes into stable scalar factors associated with s, z, and π, together with the infected-cell and viral-load subsystem. The corresponding characteristic equation is

(λ+a)(λ+ν)(λ+ω+ξ)[(λ+1)(λ+χ)-bκEsΘe-λτs]=0.(36)

The roots

λ=-a,λ=-ν,λ=-(ω+ξ)(37)

are always negative. Therefore, the stability of E₀ is determined by

(λ+1)(λ+χ)-bκEsΘe-λτs=0.(38)

Using the threshold quantity R0 defined in (15), Equation 38 can be written as

(λ+1)(λ+χ)-χR0e-λτs=0.(39)

As shown in Theorem 3.1, if R0<1, then E₀ is locally asymptotically stable for all τ_*s*≥0. If R0>1, then E₀ is unstable. Hence, in this formulation, the intracellular delay modifies the value of R0 through the survival factor Es=e-στs, but once R0<1, it does not generate a delay-induced Hopf bifurcation around E₀.

This result has a useful biological interpretation. When the infection-generating potential is below one, viral invasion is suppressed regardless of the intracellular waiting time. However, when R0>1, infection can invade the disease-free state, and the dynamics move toward an endemic regime. Thus, delay alone cannot destabilize the infection-free state when the infection threshold remains below unity.

4.2 Stability of the CTL-free endemic equilibrium

The CTL-free endemic equilibrium is

E1=(s1,i1,0,v1,π-),

where s₁, i₁, and v₁ are given in Equations 18–20. The CTL direction is governed by

λ+ν-rEIv1e-λτI=0.(40)

Using the CTL-activation threshold RT defined in Equation 24, this becomes

λ+ν-νRTe-λτI=0.(41)

Thus, the CTL-free equilibrium is stable in the CTL direction if RT<1 and unstable in the CTL direction if RT>1. The remaining infection subsystem determines whether delay-induced oscillations arise within the CTL-free endemic state. Therefore, the stability of E₁ must be interpreted in two parts: CTL invasion through Equation 41, and infection-subsystem stability through the characteristic equation below.

Linearizing the (s, i, v, π)-subsystem around E₁ gives a transcendental characteristic equation of the form

Δ1(λ;τs)=P1(λ)+Q1(λ)e-λτs=0,(42)

where P₁ contains the instantaneous terms and Q₁ contains the delayed infection terms. Explicitly,

P1(λ)=(λ+ω+ξ)[(λ+a+bΘv1)(λ+1)(λ+χ)],(43)

and

Q1(λ)=-(λ+ω+ξ)[bEsΘ{s1(λ+a+bΘv1)-bΘs1v1}κ].(44)

Therefore,

Δ1(λ;τs)=(λ+ω+ξ)[(λ+a+bΘv1)(λ+1)(λ+χ) −bκEsΘ{s1(λ+a+bΘv1)−bΘs1v1}e−λτs].(45)

Since λ = −(ω+ξ) is always stable, the relevant characteristic equation is

(λ+a+bΘv1)(λ+1)(λ+χ)−bκEsΘ{s1(λ+a+bΘv1)−bΘs1v1}e−λτs=0.(46)

Equation 46 shows explicitly how the intracellular delay affects the CTL-free endemic equilibrium through two mechanisms: the survival factor Es=e-στs and the delayed feedback term e-λτs.

4.3 Hopf bifurcation at the CTL-free endemic equilibrium

To locate possible Hopf bifurcations of E₁, we set

λ=iω1,ω1>0.(47)

Substituting Equation 47 into Equation 46 gives

P1(iω1)+Q1(iω1)e-iω1τs=0,(48)

where

P1*(λ)=(λ+a+bΘv1)(λ+1)(λ+χ),(49)

and

Q1*(λ)=-bκEsΘ{s1(λ+a+bΘv1)-bΘs1v1}.(50)

Taking moduli in (48), we obtain

∣P1(iω1)∣2=∣Q1(iω1)∣2.(51)

Equation 51 determines the admissible Hopf frequencies ω₁>0. If such a frequency exists, then the corresponding critical delays are given by

τs,n(1)=1ω1[arg(-P1(iω1)Q1(iω1))+2nπ],n=0,1,2,….(52)

The smallest positive value is denoted by

τs,0(1)=minn≥0{τs,n(1):τs,n(1)>0}.(53)

Since Es=e-στs, the modulus condition (51) and the phase condition (52) are solved consistently at the candidate critical delay. This ensures that both the delay-induced phase lag and the survival loss during the intracellular waiting period are included in the Hopf calculation.

To verify the Hopf crossing explicitly, define

F1(λ,τs)=P1(λ)+Q1(λ)e-λτs.

Differentiating F₁(λ, τ_*s*) = 0 with respect to τ_*s* gives

dλdτs=-∂F1/∂τs∂F1/∂λ.(54)

Since

∂F1∂τs=-λQ1(λ)e-λτs+∂Q1∂τse-λτs,

and

∂F1∂λ=dP1dλ+dQ1dλe-λτs-τsQ1*(λ)e-λτs,

the transversality condition at λ = iω₁ and τs=τs,0(1) is

Re[dλdτs]λ=iω1, τs=τs,0(1)≠0.(55)

This condition confirms that the conjugate pair of characteristic roots crosses the imaginary axis with non-zero speed as τ_*s* passes through the critical value. Therefore, the occurrence of purely imaginary roots is not merely tangential but corresponds to a genuine Hopf bifurcation.

Theorem 4.1.

Assume that R0>1, RT<1, and that the CTL-free endemic equilibrium E₁ is locally asymptotically stable when τ_*s* = 0. If Equation 51 admits a positive root ω₁, and the transversality condition

Re[dλdτs]λ=iω1, τs=τs,0(1)≠0(56)

holds, then E₁ undergoes a Hopf bifurcation at

τs=τs,0(1).(57)

Consequently, a branch of small-amplitude periodic solutions bifurcates from E₁ as the intracellular delay crosses the critical value. These periodic solutions represent delay-induced oscillatory infection dynamics in the CTL-free endemic regime.

Direction, stability, and period.

To determine the direction and stability of the bifurcating periodic solutions, we use the standard center-manifold and normal-form reduction for retarded functional differential equations. Let c₁(0) denote the first Lyapunov coefficient obtained from the reduced normal form at the critical eigenvalues ±iω₁. Define

μ2=-Re{c1(0)}Re(dλdτs)τs=τs,0(1),(58)

β2=2Re{c1(0)},(59)

and

T2=-Im{c1(0)}+μ2Im(dλ/dτs)τs=τs,0(1)ω1.(60)

The sign of μ₂ determines the direction of the Hopf bifurcation. If μ₂>0, the Hopf bifurcation is supercritical and the periodic solutions exist for τs>τs,0(1). If μ₂ < 0, the Hopf bifurcation is subcritical and the periodic solutions exist for τs<τs,0(1). The sign of β₂ determines the orbital stability of the bifurcating periodic solutions. If β₂ < 0, the bifurcating periodic solutions are orbitally stable, whereas if β₂>0, they are unstable. The period of the bifurcating periodic solutions is approximated by

THopf=2πω1(1+T2(τs−τs,0(1))+O((τs−τs,0(1))2)).(61)

In particular, the leading-order period is 2π/ω₁. Thus, the critical Hopf frequency provides a direct estimate of the oscillation period near the bifurcation point.

Proof. The characteristic equation near E₁ is given by (46). Setting λ = iω₁ and separating the exponential term yields (48). The modulus condition (51) guarantees the existence of purely imaginary roots. The phase condition then gives the critical delay sequence (52). If τs,0(1) is the smallest positive critical delay and the transversality condition (56) holds, then a conjugate pair of characteristic roots crosses the imaginary axis with non-zero speed. The Hopf bifurcation theorem for retarded functional differential equations implies the existence of a local branch of periodic solutions near E₁. The normal-form coefficients in (58)–(60) determine the direction, orbital stability, and period correction of the bifurcating periodic solutions.

4.4 Stability of the CTL-activated endemic equilibrium

The CTL-activated endemic equilibrium is

E2=(s2,i2,z2,v2,π-),

where the components are given in Equations 27–30. At this equilibrium, both the intracellular infection delay and the immune-activation delay enter the characteristic equation. The linearized system has the form

x′(t)=A0(2)x(t)+As(2)x(t-τs)+AI(2)x(t-τI),(62)

and the characteristic equation is

Δ2(λ;τs,τI)=det[λI-A0(2)-As(2)e-λτs-AI(2)e-λτI]=0.(63)

For analytical clarity, the characteristic equation can be written in the form

Δ2(λ;τs,τI)=P2(λ)+Q2(λ)e-λτs+R2(λ)e-λτI

+H2(λ)e-λ(τs+τI)=0,(64)

where P₂ contains the instantaneous contribution, Q₂ represents delayed infection, R₂ represents delayed CTL activation, and H₂ contains mixed delayed terms arising from the coupling between infection and immune activation. This form shows that the CTL-activated equilibrium can be affected not only by each biological delay separately, but also by their interaction through the combined delay τ_*s*+ τ_*I*.

Compared with the CTL-free endemic equilibrium, the stability of E₂ is more complex because infected cells, free virus, and CTLs are all positive and dynamically coupled. Therefore, intracellular infection delay and immune-activation delay may jointly alter the eigenvalue distribution. A Hopf bifurcation may occur when a conjugate pair of characteristic roots crosses the imaginary axis as either τ_*s* or τ_*I* is varied.

4.5 Hopf bifurcation at the CTL-activated endemic equilibrium

To analyze delay-dependent oscillatory infection dynamics around E₂, we consider two cases.

First, if τ_*I* is fixed and τ_*s* is treated as the bifurcation parameter, we set

λ=iω2, ω2>0,(65)

in Equation 64. The characteristic equation becomes

P2(iω2)+R2(iω2)e-iω2τI+[Q2(iω2)+H2(iω2)e-iω2τI]e-iω2τs =0.(66)

Solving for e-iω2τs, we obtain

e-iω2τs=-P2(iω2)+R2(iω2)e-iω2τIQ2(iω2)+H2(iω2)e-iω2τI.(67)

Thus, admissible Hopf frequencies satisfy

∣P2(iω2)+R2(iω2)e−iω2τI∣2=∣Q2(iω2)+H2(iω2)e−iω2τI∣2.(68)

If ω₂>0 satisfies (68), the critical values of τ_*s* are

τs,n(2)=1ω2[arg(-P2(iω2)+R2(iω2)e-iω2τIQ2(iω2)+H2(iω2)e-iω2τI)+2nπ], n=0,1,2,….(69)

Second, if τ_*s* is fixed and τ_*I* is treated as the bifurcation parameter, an analogous calculation gives

e-iω2τI=-P2(iω2)+Q2(iω2)e-iω2τsR2(iω2)+H2(iω2)e-iω2τs.(70)

The corresponding frequency condition is

∣P2(iω2)+Q2(iω2)e-iω2τs∣2=∣R2(iω2)+H2(iω2)e-iω2τs∣2.(71)

The critical immune-activation delays are

τI,n(2)=1ω2[arg(-P2(iω2)+Q2(iω2)e-iω2τsR2(iω2)+H2(iω2)e-iω2τs)+2nπ],n=0,1,2,….(72)

Equations 68, 71 show that the Hopf frequencies at the CTL-activated endemic equilibrium depend on both biological delays. Even when one delay is treated as the bifurcation parameter, the other delay still enters the frequency condition through the phase terms.

To verify that the purely imaginary roots produce a genuine Hopf bifurcation, the transversality condition must be checked explicitly. Differentiating the characteristic equation Δ₂(λ; τ_*s*, τ_*I*) = 0 with respect to the bifurcation delay gives

dλdτs=-∂Δ2/∂τs∂Δ2/∂λ,dλdτI=-∂Δ2/∂τI∂Δ2/∂λ.(73)

Thus, if τ_*s* is the bifurcation parameter, the transversality condition is

Re[dλdτs]λ=iω2, τs=τs,0(2)≠0.(74)

If τ_*I* is the bifurcation parameter, the transversality condition is

Re[dλdτI]λ=iω2, τI=τI,0(2)≠0.(75)

Conditions (74–75) ensure that the conjugate pair of characteristic roots crosses the imaginary axis with non-zero speed as the selected delay passes through its critical value. Hence, the crossing is not tangential and corresponds to a Hopf bifurcation.

Theorem 4.2.

Assume that E₂ exists, namely

R0>1,RT>1,

and that E₂ is locally asymptotically stable when the relevant delay is below its first critical value. If either Equation 68 or Equation 71 admits a positive Hopf frequency ω₂, and the corresponding transversality condition (Equation 74) or (Equation 75) holds at the critical delay, then E₂ undergoes a Hopf bifurcation.

Consequently, small-amplitude periodic solutions emerge near the CTL-activated endemic equilibrium as either the intracellular delay or the immune-activation delay crosses its critical value. These periodic solutions represent delay-induced oscillatory infection dynamics under CTL-mediated immune control.

Direction, stability, and period.

The direction and stability of the bifurcating periodic solutions are determined using the standard center-manifold and normal-form theory for retarded functional differential equations. Let c1(2)(0) denote the first Lyapunov coefficient at the Hopf point of E₂. When τ_*s* is the bifurcation parameter, define

μs(2)=-Re{c1(2)(0)}Re(dλ/dτs), βs(2)=2Re{c1(2)(0)}.(76)

If μs(2)>0, then the Hopf bifurcation is supercritical and the periodic solutions exist for τs>τs,0(2). If μs(2)<0, then the Hopf bifurcation is subcritical and the periodic solutions exist for τs<τs,0(2). If βs(2)<0, the bifurcating periodic solutions are orbitally stable; if βs(2)>0, they are unstable.

Similarly, when τ_*I* is the bifurcation parameter, the corresponding quantities are

μI(2)=-Re{c1(2)(0)}Re(dλ/dτI), βI(2)=2Re{c1(2)(0)}.(77)

The signs of μI(2) and βI(2) determine the direction and orbital stability of the bifurcating periodic solutions when the immune-activation delay is varied.

The leading-order period of the bifurcating periodic solutions is determined by the critical Hopf frequency and is given by

THopf(2)=2πω2.(78)

Thus, the Hopf frequency provides a direct estimate of the oscillation period near the critical delay. Higher-order corrections to the period can be obtained from the imaginary part of the normal-form coefficient, but the leading-order period 2π/ω₂ is sufficient for identifying the dominant oscillatory time scale near the Hopf point.

Proof. The characteristic equation at E₂ is given by Equation 64. When one delay is fixed and the other is varied, the equation can be rearranged into the form (67) or (70). The modulus conditions (68) and (71) determine the existence of purely imaginary characteristic roots, while the phase conditions yield the critical delay sequences (69) and (72). If the crossing speed is non-zero, as required by the transversality conditions (74) and (75), then a conjugate pair of characteristic roots crosses the imaginary axis transversely. The Hopf bifurcation theorem for retarded functional differential equations then implies the emergence of small-amplitude periodic solutions near E₂. The first Lyapunov coefficient c1(2)(0) determines the direction and stability of these periodic solutions, while the critical frequency ω₂ determines their leading-order period.

4.6 Biological interpretation of delay-induced oscillations

The above results show that biological delays influence HIV dynamics in two distinct ways. First, the intracellular delay τ_*s* may also destabilize the CTL-free endemic state and generate delay-dependent oscillatory infection dynamics. Second, the immune-activation delay τ_*I* reduces the effective CTL stimulation through e-σIτI, thereby decreasing RT. At the same time, delayed CTL feedback can destabilize the CTL-activated endemic equilibrium, leading to recurrent fluctuations in infected cells, free virus, and immune response.

Mathematically, this occurs because the delayed terms introduce phase lags into the infection and immune-feedback loops. If the phase lag is sufficiently large, the system may lose steady-state stability through Hopf bifurcation and exhibit periodic infection dynamics. The intracellular delay mainly affects the timing of productive infection, whereas the immune-activation delay affects the timing of CTL expansion. Their interaction can therefore produce oscillatory behavior even when the equilibrium is biologically feasible.

From a biological perspective, these oscillations may be interpreted as theoretical delay-induced infection fluctuations or unstable immune-mediated control within the model. In the absence of sufficiently rapid treatment or immune activation, viral replication may increase before the delayed suppressive mechanisms become effective. This provides a mechanistic explanation for why intervention timing and delay-aware control strategies are important in within-host HIV dynamics.

However, these oscillations should not be interpreted as direct clinical predictions of patient-specific viral rebound. The present model does not explicitly include latent reservoirs, stochastic viral reactivation, adherence variability, pharmacokinetic–pharmacodynamic drug exposure, or patient-specific calibration. Therefore, the Hopf bifurcation results should be viewed as theoretical evidence that biological delays can create oscillatory infection regimes and reduce the robustness of immune-mediated control.

Overall, the Hopf analysis shows that the stability of within-host HIV dynamics depends not only on threshold quantities such as R0 and RT, but also on the timing of intracellular infection, immune activation, and intervention response. This supports the need for delay-aware modeling when evaluating immune control and treatment timing in HIV infection.

5 Numerical simulations and sensitivity analysis

In this section, numerical simulations are performed to illustrate the theoretical results obtained in Sections 3 and 4. The objectives are threefold: first, to examine how the intracellular infection delay τ_*s* and the immune-activation delay τ_*I* influence transient and long-term infection dynamics; second, to identify parameter regimes associated with stable suppression, persistent infection, and delay-dependent oscillatory behavior; and third, to quantify the relative influence of model parameters on peak infected-cell burden using global sensitivity analysis.

All simulations are based on the dimensionless system (9). The delay differential equations are solved using MATLAB's built-in solver dde23. Unless otherwise stated, the relative and absolute tolerances are set to standard solver tolerances, and robustness with respect to tighter tolerances is examined in Section 5.6. The initial history functions are taken to be constant on [−τ_max, 0], where

τmax=max{τs,τI}.(79)

The baseline history is chosen as

s(θ)=s0, i(θ)=i0, z(θ)=z0, v(θ)=v0, π(θ)=π0, θ∈[-τmax,0].(80)

The numerical values used in the baseline simulations are selected to represent biologically plausible within-host dynamics and to illustrate the stability transitions predicted by the theoretical analysis. The baseline parameter set is

a=0.10,b=2.00,ε=0.40,σ=0.20,σI=0.15,r=1.20,ν=0.40,κ=1.50,χ=0.80,ω=0.30,ξ=0.20.(81)

The parameter values are not intended to represent a patient-specific calibrated HIV data set. Instead, they form a dimensionless, literature-motivated baseline regime that is consistent with standard within-host HIV modeling assumptions and is suitable for illustrating infection persistence, immune response, delay-dependent stability switching, and post-control resurgence mechanisms [1–4, 8, 14–18]. The selected values produce biologically feasible equilibria and allow the effects of intracellular delay, immune-activation delay, and protection efficacy to be examined clearly.

Table 2 summarizes the parameter meanings, baseline values, units in the dimensionless formulation, and the justification for their selection.

Table 2. Baseline parameter values used in the numerical simulations.

The table clarifies that the simulations are designed to explore theoretical delay-driven mechanisms rather than to estimate patient-specific disease trajectories. Parameters associated with viral production, viral clearance, immune activation, and protection are varied later in the global sensitivity analysis to assess the robustness of the qualitative conclusions.

The corresponding protective-response equilibrium is

π-=ωω+ξ=0.60,(82)

and hence

Θ=1-επ-=0.76.(83)

5.1 Numerical threshold evaluation

Using the baseline parameter values, the basic infection threshold is

R0=bκe-στs(1-επ-)χ.(84)

Thus, for fixed parameters except τ_*s*, increasing the intracellular delay reduces R0 through the survival factor e-στs. Similarly, the CTL-activation threshold is

RT=re-σIτIν·a(R0-1)bΘ.(85)

Equation 85 shows that τ_*I* reduces the effective CTL stimulation, while τ_*s* indirectly affects CTL activation through R0. These two delays therefore act through different but coupled mechanisms.

For illustration, we consider three representative delay regimes

(τs,τI)=(0.5,0.5), (τs,τI)=(1.5,1.5), (τs,τI)=(3.0,3.0).(86)

These regimes are used to compare short-delay, intermediate-delay, and long-delay dynamics. They are not intended as patient-specific delay estimates, but as representative values for illustrating the qualitative delay-dependent behavior of the model.

In addition to these three representative cases, the two-dimensional delay-response analysis in Section 5.3 and the global sensitivity analysis in Section 5.4 provide additional simulations under broader delay and parameter variations. These simulations were added to strengthen the numerical support for the delay-dependent stability and bifurcation results.

5.2 Time-series dynamics under increasing delays

Figure 2 presents the time-series trajectories of s(t), i(t), z(t), and v(t) under the three delay regimes in (86). For short delays, the trajectories approach a stable endemic equilibrium with mild transient variation. In this regime, infection persists but remains controlled, and the CTL response converges smoothly to its equilibrium level.

For intermediate delays, the convergence becomes slower and more oscillatory. The infected-cell and viral-load trajectories exhibit transient infection peaks before approaching equilibrium. This behavior reflects the fact that delayed infection progression and delayed immune activation create a phase lag between viral growth and immune-mediated suppression.

For long delays, more pronounced delay-dependent oscillatory behavior may occur. In this case, viral replication increases before the delayed CTL response becomes effective, leading to recurrent peaks in infected cells and free virus. These oscillations are consistent with the Hopf bifurcation mechanism described in Section 4. From a biological perspective, they should be interpreted as theoretical delay-induced infection fluctuations driven by delayed immune feedback, rather than direct clinical viral rebound episodes.

The peak infected-cell burden over the simulation interval is defined by

Imax=max0≤t≤Tsimi(t),(87)

where T_*sim* is the final simulation time. In the numerical experiments, T_*sim* = 200 is used to capture both transient and long-term behavior.

Figure 2 therefore provides a direct numerical illustration of how increasing biological delays can change the qualitative infection trajectory from smooth convergence to more pronounced transient or oscillatory behavior. This supports the theoretical prediction that delay-induced phase lags can alter local stability and generate oscillatory infection dynamics.

5.3 Delay-response patterns of peak infected-cell burden

To quantify how delays influence the peak infected-cell burden, we compute I_max over a two-dimensional delay grid

(τs,τI)∈[0,τsmax]×[0,τImax].(88)

For each pair (τ_*s*, τ_*I*), system (9) is simulated and I_max is recorded. This produces a delay-response surface, as shown in Figure 3.

(τs,τI)↦Imax(τs,τI).(89)

This analysis provides a numerical visualization of how intracellular infection delay and immune-activation delay jointly shape the maximum infected-cell burden. It also complements the Hopf analysis by identifying regions in the delay plane where transient peaks become larger or smaller.

This additional delay-grid simulation directly addresses the need to investigate how different delay values affect model behavior. Rather than relying only on three representative time-series simulations, the delay-response surface summarizes the influence of a broad range of (τ_*s*, τ_*I*) combinations on the peak infected-cell burden.

The resulting delay-response surface shows that the peak infected-cell burden is strongly shaped by the intracellular infection delay τ_*s*. In the present parameter setting, increasing τ_*s* tends to reduce I_max, reflecting attenuation of effective infection through the survival factor e-στs. The effect of the immune-activation delay τ_*I* is comparatively weaker over the explored range, but it still contributes to the delay-dependent transient response. These simulations therefore support the analytical finding that biological delays can reshape infection dynamics by balancing delay-induced attenuation against delay-induced destabilization of persistent infection.

The decrease of I_max with increasing τ_*s* in this parameter regime should not be interpreted as a universal effect of intracellular delay. Rather, it reflects the combined influence of delayed phase lag and survival loss e-στs. In other parameter regimes, stronger delay-induced phase lag may increase transient oscillations even when survival loss reduces the effective infection rate. This motivates the global sensitivity analysis below.

5.4 Global sensitivity analysis

A global sensitivity analysis is performed to identify which model parameters most strongly influence the peak infected-cell burden. We use Latin hypercube sampling because it provides efficient stratified coverage of multidimensional parameter spaces [24]. The response variable is I_max, defined in (87). A total of N = 1, 000 parameter sets are generated, and partial rank correlation coefficients (PRCCs) are computed to quantify the monotonic association between each parameter and I_max. The resulting PRCC values are shown in Figure 4.

Let

p=(a,b,ε,σ,σI,r,ν,κ,χ,ω,ξ,τs,τI)(90)

denote the vector of sampled parameters. Each parameter is varied independently within a prescribed range around its baseline value. Specifically, for a baseline value pj0, the sampled interval is

pj∈[0.7pj0,1.3pj0],(91)

unless biological feasibility requires a narrower range. For example, parameter ranges are restricted when necessary to preserve non-negative effective susceptibility and biologically meaningful delay values. For each Latin hypercube sample, the system is simulated over [0, T_*sim*], and I_max is recorded.

The ±30% sampling range was chosen to test whether the qualitative conclusions remain stable under moderate parameter uncertainty. This is appropriate because the baseline parameter set is dimensionless and illustrative rather than patient-calibrated. The use of PRCCs allows the analysis to identify which parameters most consistently increase or decrease peak infected-cell burden across the sampled parameter space.

The monotonic relationship between each parameter and I_max was quantified using partial rank correlation coefficients. If R(p_*j*) denotes the rank-transformed values of parameter p_*j*, and R(I_max) denotes the rank-transformed response values, then the PRCC for p_*j* was computed after removing the linear rank effects of all other parameters. Parameters with large positive PRCC values increase the peak infected-cell burden, whereas parameters with large negative PRCC values decrease the peak infected-cell burden. In this analysis, viral infectivity band viral production κ show strong positive associations with I_max, while viral clearance χ, intracellular delay τ_*s*, and loss rate σ show strong negative associations. These results indicate that both viral kinetic parameters and delay-survival mechanisms play major roles in determining the maximum infected-cell burden.

This parameter-variation analysis strengthens the numerical conclusions by showing that the observed infection peaks are not controlled by a single parameter alone. Instead, peak infected-cell burden is shaped jointly by infection rate, viral production, viral clearance, protection efficacy, immune activation, and delay-survival effects.

5.5 Interpretation of sensitivity results

The sensitivity results identify the effective infection parameter b, viral production factor κ, viral clearance rate χ, CTL activation rate r, CTL decay rate ν, protection efficacy ε, and the delays τ_*s* and τ_*I* as important drivers of the peak infected-cell burden.

A positive PRCC for b indicates that stronger effective infection increases peak infected-cell burden. A positive PRCC for κ indicates that higher viral production amplifies infection. A negative PRCC for χ indicates that stronger viral clearance reduces infection peaks. Similarly, a negative PRCC for ε indicates that stronger protection lowers the effective infection rate and suppresses peak infection.

The immune-related parameters have a more nuanced interpretation. A larger CTL activation rate r may reduce infection by strengthening immune control, whereas a larger CTL decay rate ν weakens immune persistence and may increase the infected-cell burden. The immune-activation delay τ_*I* is expected to have a positive association with the peak infected-cell burden, because delayed CTL expansion allows the infected-cell and viral-load levels to increase before immune suppression becomes effective. The intracellular delay τ_*s* may have mixed effects, because it reduces productive infection through the survival factor e-στs, but it may also contribute to delay-dependent oscillatory behavior when infection persists. These results support the conclusion that both viral kinetic parameters and biological delays jointly shape the magnitude and timing of infection peaks.

Overall, the sensitivity analysis supports the theoretical conclusion that peak infected-cell burden is governed by a balance among viral infectivity, viral production, viral clearance, immune activation, protection efficacy, and delay-survival effects. This also shows that the numerical outcomes are not driven by a single arbitrarily chosen parameter, but by interacting mechanisms that are biologically meaningful in within-host HIV dynamics.

The strong influence of b and κ indicates that infection establishment and viral amplification remain dominant drivers of infected-cell burden. In contrast, the negative influence of χ, ε, and τ_*s* in the present parameter regime suggests that viral clearance, protective effects, and loss during the intracellular waiting period can reduce the maximum infected-cell level. These findings are consistent with the threshold analysis, where R0 decreases when protection increases or when the survival factor e-στs becomes smaller.

The role of τ_*I* is particularly important from an immune-control perspective. A longer immune-activation delay can allow infected cells and free virus to increase before CTL expansion becomes effective. Thus, even if the CTL response is eventually activated, delayed immune feedback may produce larger transient infection peaks. This supports the need to include immune-response timing when evaluating within-host HIV persistence and delayed intervention effects.

5.6 Numerical robustness checks

To ensure that the numerical conclusions are not artifacts of a particular solver configuration, tolerance setting, or simulation horizon, numerical robustness checks were performed. The simulations were repeated using different relative and absolute tolerances in MATLAB dde23. The qualitative regimes were considered robust when the same stability pattern, transient response, or delay-dependent oscillatory behavior was observed under tighter tolerance settings.

In addition, the final simulation time T_*sim* is increased to verify whether apparent oscillations persist or eventually decay. This distinction is important because long transient oscillations can be mistaken for stable periodic behavior if the simulation horizon is too short. Therefore, a trajectory is classified as sustained oscillatory only if oscillations persist over an extended simulation interval after the initial transient has passed.

The numerical results are summarized using the following outputs, defined in Equations 92–94.

time - series trajectories of s(t),i(t),z(t),v(t),(92)

delay - response surfaces for Imax(τs,τI),(93)

and

PRCC values for the sampled parameter set p.(94)

Together, these outputs connect the theoretical stability thresholds to observable infection patterns, including transient infection peaks, delay-dependent oscillatory behavior, and parameter-driven changes in the peak infected-cell burden. This provides a numerical consistency check for the analytical results and supports the robustness of the main conclusions.

These robustness checks also address the concern that the simulations merely illustrate one selected parameter set. By varying solver tolerances, extending the simulation horizon, exploring delay-response surfaces, and performing global sensitivity analysis, the numerical section provides broader support for the analytical delay-dependent stability and bifurcation results.

Nevertheless, the simulations remain theoretical and dimensionless. They are intended to demonstrate qualitative mechanisms rather than to provide patient-specific predictions. Future work should calibrate the model to longitudinal viral-load, CD4⁺ T-cell, and immune-response data to assess whether the predicted delay-dependent infection patterns occur in clinically observed HIV dynamics.

6 Delayed optimal control formulation and post-control resurgence analysis

In this section, the baseline delay-structured HIV model is extended to include time-dependent treatment and prevention controls. The aim is to evaluate how delayed intervention effects influence viral suppression, immune control, treatment effort, and post-control resurgence. In contrast to the baseline model in Sections 2–5, where only biological state delays are considered, this section introduces a pharmacological or implementation delay associated with the response of the system to control actions. This distinction is important because treatment or prevention may be initiated at a given time, but its effect on viral replication or susceptibility may occur only after a biologically or operationally meaningful response delay.

Let u₁(t) denote the intensity of antiretroviral treatment, and let u₂(t) denote the intensity of prevention or protective intervention. The control u₁(t) reduces viral production by productively infected cells, while u₂(t) strengthens the protective effect that reduces effective susceptibility. The admissible controls are assumed to satisfy

0≤u1(t)≤1, 0≤u2(t)≤1, t∈[0,Tf].(95)

Here, u_*i*(t) = 0 represents no intervention of type i, whereas u_*i*(t) = 1 represents the maximum admissible intervention intensity. The terminal time T_*f* denotes the end of the prescribed control horizon. After T_*f*, the controls are withdrawn in the post-control simulations in order to evaluate whether infection remains suppressed or whether post-control resurgence occurs. Thus, the control problem is used to study not only suppression during the treatment interval, but also the persistence of suppression after intervention withdrawal.

6.1 Controlled model with pharmacological delay

We introduce a control-response delay τ_*c*≥0, representing the time required for treatment or prevention actions to affect viral replication and effective susceptibility. The controlled dimensionless system is given by

ds(t)dt=a−as(t)−b[1−επ(t)−η2u2(t−τc)]s(t)v(t), di(t)dt=be−στs[1−επ(t−τs) −η2u2(t−τs−τc)]s(t−τs)v(t−τs) −i(t)−z(t)i(t),dz(t)dt=re−σIτIv(t−τI)z(t−τI)−νz(t),dv(t)dt=κ[1−η1u1(t−τc)]i(t)−χv(t),dπ(t)dt=ω(1−π(t))−ξπ(t).(96)

The parameters η₁ and η₂ satisfy

0≤η1≤1, 0≤η2≤1,(97)

and represent the maximum efficacy of treatment and prevention controls, respectively. To preserve non-negative effective infection, we assume

1-επ(t)-η2u2(t-τc)≥0(98)

for admissible trajectories, or equivalently impose the truncated effective-susceptibility factor

[1-επ(t)-η2u2(t-τc)]+=max{0,1-επ(t)-η2u2(t-τc)}.(99)

In the numerical implementation, the admissible control range is chosen so that the effective-susceptibility factor remains non-negative. This avoids biologically unrealistic negative infection rates. The control history is prescribed on the interval [−τ_*c*, 0] as

u1(θ)=u1,0, u2(θ)=u2,0, θ∈[-τc,0].(100)

For delayed infection terms involving u₂(t−τ_*s*−τ_*c*), the control history is extended on

[-(τs+τc),0].(101)

This formulation separates the biological delays τ_*s*, τ_*I* from the pharmacological or implementation delay τ_*c*. The biological delays govern intrinsic infection and immune-response processes, whereas τ_*c* determines how quickly control actions influence the within-host system. This distinction is important because even a strong intervention may have limited short-term effect if its pharmacological response is delayed.

Biologically, u₁(t) may be interpreted as the relative intensity of antiretroviral treatment that reduces viral production from infected cells. The term 1−η₁u₁(t−τ_*c*) therefore represents the remaining viral production capacity after delayed treatment response. Similarly, u₂(t) may be interpreted as the relative intensity of prevention or protective support that reduces effective susceptibility. The delayed terms involving u₂ indicate that the protective effect does not act instantaneously on either current infection or delayed infection production.

6.2 Objective functional

The objective is to reduce infected-cell burden and viral load while penalizing excessive intervention effort. We define the performance index

J(u1,u2)=∫0Tf[Ai(t)+Bv(t)+C12u12(t)+C22u22(t)]dt,(102)

where

A>0, B>0, C1>0, C2>0.(103)

The weights A and B measure the relative importance assigned to infected-cell burden and viral load, respectively. The weights C₁ and C₂ penalize treatment and prevention effort. Quadratic control costs are used because they provide a convex penalty for high intervention intensity and are standard in biological optimal-control problems [23]. This formulation also permits the derivation of explicit projection formulas for the optimal controls.

The biological interpretation of these weights is as follows. The coefficient A represents the priority assigned to reducing productively infected cells, which reflects the importance of limiting infected-cell burden and preventing further viral production. The coefficient B represents the priority assigned to reducing free virus, which reflects the importance of lowering viral load and reducing ongoing infection pressure. The coefficient C₁ represents the relative cost, toxicity, adherence burden, or feasibility constraint associated with antiretroviral treatment intensity. The coefficient C₂ represents the relative cost or implementation burden associated with prevention or protective intervention. Thus, larger values of A or B favor stronger suppression of infection, whereas larger values of C₁ or C₂ discourage intensive control use.

These weights are not intended to prescribe a patient-specific clinical regimen. Instead, they allow the model to explore theoretical trade-offs between infection suppression and intervention effort. To avoid conclusions depending on a single arbitrary choice of weights, robustness with respect to different control-cost regimes is examined in the numerical simulations below.

The admissible control set is

U={(u1,u2)∈L∞(0,Tf)2:0≤u1(t)≤1,0≤u2(t)≤1}.(104)

An optimal control pair (u1∗,u2∗)∈U is sought such that

J(u1∗,u2∗)=min(u1,u2)∈UJ(u1,u2).(105)

The optimal-control problem is therefore formulated over a finite intervention horizon [0, T_*f*]. The post-control period is not included in the objective functional, but is simulated afterward to assess whether the controlled system remains suppressed or exhibits post-control resurgence after withdrawal of control.

6.3 Hamiltonian and necessary optimality conditions

Let

λ1(t), λ2(t),λ3(t),λ4(t),λ5(t)

be the adjoint variables associated with s(t), i(t), z(t), v(t), π(t), respectively. For compact notation, the Hamiltonian is first written using the delayed state terms while suppressing the explicit time-shift notation in the controls. In the full delayed-control formulation, u₁ and u₂ are replaced by their delayed values u₁(t−τ_*c*) and u₂(t−τ_*c*), as specified in Equation 96. The Hamiltonian is

H=Ai+Bv+C12u12+C22u22 +λ1[a-as-b(1-επ-η2u2)sv] +λ2[be-στs(1-επτs-η2u2,τs)sτsvτs-i-zi] +λ3[re-σIτIvτIzτI-νz] +λ4[κ(1-η1u1)i-χv] +λ5[ω(1-π)-ξπ],(106)

where the shorthand notation

sτs=s(t-τs), vτs=v(t-τs), πτs=π(t-τs),(107)

and similarly, for τ_*I*, is used. In the full delayed-control setting, u₁ and u₂ appearing in the state equations are replaced by their delayed values u₁(t−τ_*c*) and u₂(t−τ_*c*). This introduces advanced terms in the adjoint equations, as is standard in optimal-control problems with delays.

The adjoint equations satisfy

−λj'(t)=∂H∂xj(t)+1[0,Tf−τs](t)∂H(t+τs)∂xj(t) +1[0,Tf−τI](t)∂H(t+τI)∂xj(t),(108)

where x_*j*∈{s, i, z, v, π}, and 1 denotes the indicator function. The terminal transversality conditions are

λj(Tf)=0, j=1,…,5.(109)

The delayed terms in (108) arise because a state variable at time t may influence the system at future times t+τ_*s* or t+τ_*I*. This is a key distinction between instantaneous optimal control and delayed optimal control. For the pharmacological delay τ_*c*, the same principle applies to the delayed controls; in numerical implementation, the corresponding delayed control values are evaluated by interpolation from the current control iterate.

The necessary optimality conditions are used here as a theoretical and computational framework for evaluating delay-aware control strategies. Because the model is dimensionless and not patient-calibrated, the resulting control profiles should be interpreted as illustrative intervention patterns rather than direct clinical dosing schedules.

6.4 Characterization of the optimal controls

Differentiating the Hamiltonian with respect to the controls gives the stationarity conditions. For the treatment control u₁, the control enters through the viral-production term. The formal stationarity condition is

∂H∂u1=C1u1-η1κλ4i=0.(110)

u1(t)=η1κλ4(t)i(t)C1.(111)

For the prevention control u₂, the control enters through the effective infection term. The formal stationarity condition is

∂H∂u2=C2u2+η2bλ1sv-η2be-στsλ2sτsvτs=0.(112)

u2(t)=η2be-στsλ2(t)s(t-τs)v(t-τs)-η2bλ1(t)s(t)v(t)C2.(113)

Projecting onto the admissible interval [0, 1], the optimal controls are characterized by

u1*(t)=min{1,max{0,η1κλ4(t)i(t)C1}},(114)

and

u2∗(t)=min{1,max{0,η2be−στsλ2(t)s(t−τs)v(t−τs)−η2bλ1(t)s(t)v(t)C2}}.(115)

Equations 114, 115 show how treatment and prevention intensities depend on the current state, delayed infection terms, adjoint variables, and relative cost weights. Because the state equations depend on delayed controls u₁(t−τ_*c*) and u₂(t−τ_*c*), the projection formulas are implemented with the appropriate time shift in the forward–backward sweep algorithm. Specifically, a control value updated at time t affects the state equations at time t+τ_*c*. Therefore, in the numerical implementation, delayed control terms are obtained by interpolation from the current control iterate.

In the simulations reported below, the delayed optimal-control framework is used to motivate the treatment and prevention scenarios, while the post-control analysis focuses on illustrative intervention profiles and their robustness after control withdrawal. This distinction prevents overinterpretation of the numerical curves as patient-calibrated optimal treatment schedules and keeps the results within the scope of a theoretical delay-based model.

6.5 Numerical solution by forward–backward sweep

The optimality system is solved numerically using a forward–backward sweep method for delayed optimal-control problems [18]. The interval [0, T_*f*] is discretized into N time steps. Starting from an initial guess for u₁(t) and u₂(t), the state system is integrated forward in time using the prescribed histories. The adjoint system is then integrated backward from the terminal conditions (Equation 109). The controls are updated using the projection formulas (Equations 114, 115). To improve convergence, a relaxation update is used

ui(k+1)(t)=θui(k)(t)+(1-θ)ũi(k+1)(t), i=1,2,(116)

where 0 < θ < 1, and ũi(k+1) denotes the newly computed projected control. The iteration is stopped when

maxt∈[0,Tf](∣u1(k+1)(t)-u1(k)(t)∣+∣u2(k+1)(t)-u2(k)(t)∣)<tol.(117)

Because the controlled system contains delayed controls, the delayed terms u₁(t−τ_*c*) and u₂(t−τ_*c*) are evaluated by interpolation from the current control iterate. This implementation ensures that the pharmacological delay is consistently represented in both the forward state equations and the backward update of the control profiles.

The convergence criterion in Equation 117 ensures that the control profiles no longer change appreciably between successive iterations. The relaxation parameter θ is used to reduce numerical oscillations in the forward–backward sweep method, which can occur in delayed optimal-control systems because current control updates affect future state dynamics through τ_*c*.

6.6 Numerical control and post-control simulation setup

To connect the delayed optimal-control formulation with observable infection outcomes, we conduct numerical control simulations using the same baseline parameter set as in Section 5. In the simulations reported here, the control profiles are used as illustrative delayed-intervention scenarios motivated by the optimal-control formulation, rather than as patient-calibrated optimal treatment schedules. Unless otherwise stated, the control horizon and post-control simulation time are set to

Tf=100, Tpost=200,(118)

so that the controlled dynamics are simulated on [0, T_*f*], followed by an uncontrolled post-control phase on [T_*f*, T_*post*]. The efficacy parameters are chosen as

η1=0.75, η2=0.60,(119)

where η₁ represents the maximum reduction in viral production due to antiretroviral treatment and η₂ represents the maximum strengthening of prevention-induced susceptibility reduction. The values of η₁ and η₂ represent strong but incomplete intervention efficacy. They do not imply complete viral suppression or perfect prevention, which would be biologically unrealistic in a theoretical within-host model with delayed response.

The baseline objective weights are

A=1, B=1, C1=1, C2=1,(120)

and the baseline control-response delay is

τc=1.(121)

The weights A and B represent the relative priority assigned to reducing infected cells and free virus, respectively. The weights C₁ and C₂ represent the relative cost, toxicity, adherence burden, or implementation difficulty associated with treatment and prevention controls. Thus, larger A and B promote stronger infection suppression, whereas larger C₁ and C₂ discourage intensive intervention. The baseline choice A = B = C₁ = C₂ = 1 is used as a balanced illustrative case rather than as a patient-specific clinical cost estimate.

The regularization constant used in the post-control resurgence index is

δ=10-6.(122)

Three intervention scenarios are compared

No control: u1(t)=u2(t)=0,(123)

Treatment only: u1(t)>0,u2(t)=0,(124)

and

Combined treatment and prevention: u1(t)>0, u2(t)>0.(125)

For each scenario, the trajectories of infected cells i(t), free virus v(t), CTLs z(t), and the controls u₁(t), u₂(t) are recorded. The controls are applied only during [0, T_*f*]. After T_*f*, both controls are withdrawn and the state system is simulated without intervention to assess whether infection remains suppressed or exhibits post-control resurgence.

The numerical results are summarized by comparing the infected-cell burden during the control phase, the viral-load trajectory, the total control effort, and the post-control resurgence index. This allows the effectiveness of early intensive treatment, prevention support, and delayed intervention response to be evaluated within the same theoretical framework. The delayed-intervention and post-control resurgence scenarios are illustrated in Figure 5.

6.7 Robustness with respect to cost weights

The cost weights A, B, C₁, C₂ determine the trade-off between viral suppression and intervention effort. To avoid overinterpreting a single optimal-control solution, robustness analyses are performed by varying these weights across multiple scenarios. In particular, we compare

(C1,C2)=(0.1,0.1),(1,1),(10,10),(126)

representing low, moderate, and high intervention-cost regimes. When C₁ and C₂ are small, aggressive controls are less costly and the optimal strategy is expected to use stronger intervention intensities. When C₁ and C₂ are large, the optimal strategy is expected to reduce intervention effort, possibly at the cost of higher infected-cell and viral-load levels.

The sensitivity of the optimal trajectories to A and B is also examined. Larger A places greater emphasis on reducing infected cells, whereas larger B places greater emphasis on reducing free virus. These comparisons determine whether the qualitative conclusions are robust to alternative cost assumptions. In the numerical experiments, increasing the control-cost weights generally reduces the magnitude and duration of the optimal controls, whereas increasing the infection or viral-load weights promotes stronger suppression. Thus, the robustness analysis prevents the conclusions from depending on one arbitrary choice of objective-function weights.

Biologically, this means that the model can distinguish between settings in which strong intervention is acceptable and settings in which intervention intensity is constrained by toxicity, adherence burden, cost, or implementation feasibility. Therefore, the weight parameters should be interpreted as relative priorities in a theoretical control framework, not as fixed clinical values.

6.8 Post-control resurgence analysis

A key limitation of finite-horizon optimal-control simulations is that they may show effective suppression during the control interval but fail to determine what happens after the intervention ends. To address this issue, the state system is simulated beyond the terminal time T_*f* after controls are withdrawn

u1(t)=0, u2(t)=0, t>Tf.(127)

The post-control simulation is performed on

[Tf,Tpost],Tpost>Tf.(128)

The post-control resurgence magnitude is quantified by

Rpost=maxTf≤t≤Tposti(t),(129)

i(Tf).(130)

A post-control resurgence index is defined as

Bpost=Rposti(Tf)+δ,(131)

where δ>0 is a small regularization constant used to avoid division by zero. If

Bpost>1,(132)

then the infected-cell burden increases after the control horizon, indicating post-control post-control resurgence. This analysis is motivated by clinical evidence that viral rebound after treatment interruption is common and that post-treatment control occurs only in a subset of individuals [7]. However, in the present model, this quantity should be interpreted as a theoretical rebound-like index rather than a patient-calibrated viral rebound measure.

In the numerical experiments, controls reduce i(t) and v(t) during [0, T_*f*]. After the controls are withdrawn, the trajectories may either remain near a low-infection state or rebound toward the uncontrolled endemic regime, depending on the control-response delay τ_*c* and the cost weights C₁, C₂. Larger values of τ_*c* and larger control-cost weights generally produce weaker suppression and larger values of Bpost. Thus, post-control resurgence analysis is essential for distinguishing temporary suppression from sustained control.

This analysis extends the finite-horizon optimal-control result by testing whether infection remains controlled after intervention withdrawal. Therefore, it helps evaluate the robustness of the intervention strategy beyond the prescribed control interval.

6.9 Interpretation of delayed intervention effects

The control-response delay τ_*c* affects the timing and efficacy of intervention. When τ_*c* is small, treatment and prevention controls rapidly reduce viral replication and effective susceptibility. When τ_*c* is large, the system continues to evolve under near-uncontrolled dynamics before the intervention becomes effective. This delay can allow infected cells and viral load to increase substantially before suppression begins.

Therefore, even when the same total intervention effort is applied, a longer control-response delay can produce larger transient peaks and greater post-control resurgence risk. This emphasizes that intervention timing is not simply a numerical feature of the optimal-control problem but a biologically meaningful determinant of viral suppression and post-control resurgence.

Overall, the delayed optimal-control analysis shows that intervention intensity, intervention timing, control-response delay, and post-control behavior must be evaluated together. A strategy that appears effective during [0, T_*f*] may still be insufficient if delayed response or high control costs permit post-control resurgence after treatment withdrawal.

These findings support the need to evaluate delayed intervention effects together with post-control dynamics. In particular, treatment strategies should not be assessed only by their ability to reduce infection during the active control interval, but also by whether they prevent resurgence after controls are removed. This distinction is important in theoretical within-host HIV modeling because finite-horizon suppression does not necessarily imply long-term stability or sustained immune control.

7 Discussion and conclusion

This study developed a delay-structured within-host HIV model to investigate how biological delays, immune feedback, protective response, and delayed intervention effects shape viral persistence, immune control, delay-dependent oscillatory behavior, and post-control resurgence. The model separates three delay mechanisms: the intracellular infection delay τ_*s*, the immune-activation delay τ_*I*, and the pharmacological or implementation delay τ_*c* associated with treatment and prevention controls. This separation is important because these delays have different biological meanings and affect the system through different mechanisms. The intracellular delay modifies the transition from viral entry to productive infection, the immune-activation delay governs delayed CTL response, and the control-response delay determines how rapidly interventions influence viral replication and effective susceptibility.

Compared with classical within-host HIV models, which mainly describe target-cell infection, infected-cell production, and viral clearance [1–4], the present framework incorporates additional biological timing mechanisms. Earlier HIV delay models have shown that intracellular infection delays can alter stability and viral persistence [14–16], while immune-response models have emphasized the role of CTL-mediated control and immune escape [9–13]. However, many existing models focus on either intracellular delay or immune feedback separately. The present study extends this line of work by combining intracellular infection delay, immune-activation delay, and delayed treatment/prevention response within a single framework. This allows infection persistence, immune-control failure, delay-induced oscillations, and post-control resurgence to be studied together.

The baseline model extends standard within-host HIV dynamics by incorporating a normalized protective-response variable P(t), which reduces effective susceptibility through the factor 1−εP(t). This formulation avoids treating prevention as a direct loss of susceptible CD4⁺ T cells and provides a clearer biological interpretation of protective effects. In addition, survival factors e-mτs and e-mIτI are included in the delayed infection and immune-activation terms. These factors account for loss during the waiting periods and prevent the unrealistic assumption that all infected cells or immune signals survive unchanged until the delayed effect occurs.

The threshold analysis identified two main quantities: the basic infection threshold R0 and the CTL-activation threshold RT. The infection-free equilibrium is locally asymptotically stable when R0<1, whereas infection can persist when R0>1. When infection persists, the CTL-free endemic equilibrium exists, and the ability of the CTL response to invade this state is determined by RT. If RT>1, the CTL-activated endemic equilibrium becomes biologically feasible. These thresholds show how protection efficacy, viral clearance, immune activation, and delay-survival effects jointly determine whether infection is suppressed, persists without sustained CTL response, or persists under immune-mediated control.

The delay-dependent stability and Hopf bifurcation analysis showed that biological delays can destabilize endemic equilibria and generate delay-dependent oscillatory infection dynamics. Such oscillations should be interpreted as theoretical delay-induced infection fluctuations or unstable immune-control regimes, rather than direct patient-calibrated viral rebound. The intracellular delay τ_*s* has a dual effect: it reduces effective infection through the survival factor e-στs, but it may also contribute to delay-induced destabilization when infection persists. Similarly, the immune-activation delay τ_*I* reduces effective CTL stimulation through e-σIτI and creates a phase lag between viral expansion and immune suppression. This phase lag provides a mechanistic explanation for infection peaks when immune feedback is delayed.

The numerical simulations supported the analytical results. Short biological delays produced relatively smooth convergence to equilibrium, whereas longer delays generated stronger transient oscillations and rebound-like peaks in infected cells and free virus. The delay-response surface I_max(τ_*s*, τ_*I*) further illustrated how peak infected-cell burden changes across the delay plane. This analysis helps distinguish between regimes in which delays attenuate infection potential and regimes in which delayed feedback destabilizes persistent infection. The global sensitivity analysis identified the parameters most strongly associated with peak infected-cell burden, including infection-related, immune-related, clearance-related, and delay-related parameters. These results indicate that peak infected-cell burden and rebound-like behavior are controlled not only by viral transmission and production but also by immune timing and delay-dependent survival effects.

The delayed optimal-control formulation extended the baseline model by introducing treatment and prevention controls with a control-response delay τ_*c*. This extension addresses a practical limitation of finite-horizon control studies: suppression during the control interval does not necessarily imply sustained suppression after intervention withdrawal. The post-control resurgence analysis therefore evaluates whether infected cells remain suppressed after T_*f* or increase again once treatment and prevention controls are removed. The post-control resurgence index Bpost provides a quantitative measure for distinguishing temporary suppression from sustained control. The formulation emphasizes that intervention timing, control intensity, control cost, and post-control behavior must be evaluated together.

The results highlight the importance of early and sufficiently strong intervention. When the control-response delay is small, treatment and prevention can reduce viral replication and effective susceptibility before high infection peaks occur. In contrast, when τ_*c* is large, the infection may continue to evolve under near-uncontrolled dynamics before the intervention becomes effective. Thus, even if the same total intervention effort is applied, delayed response can lead to larger transient peaks and greater post-control resurgence risk. This finding supports the broader conclusion that intervention timing is a central determinant of within-host HIV control.

Several limitations should be acknowledged. First, the model is theoretical and uses representative parameter values rather than patient-specific calibration. Although the simulations illustrate biologically meaningful mechanisms, future work should incorporate longitudinal viral-load, CD4⁺ T-cell, and immune-response data to estimate patient-specific parameters and validate the predicted delay-dependent infection patterns. Earlier viral-kinetic studies have emphasized rapid viral turnover and clearance in HIV infection, which should be considered when calibrating future patient-specific extensions [25]. Second, the protective-response variable P(t) is represented in aggregate form. In practice, prevention, adherence, drug exposure, and immune status may involve heterogeneous mechanisms that could be modeled separately, including prevention- or vaccination-related mechanisms considered in earlier epidemic-control studies [26]. Third, the CTL response is represented by a single immune compartment. More detailed models could include immune exhaustion, antigenic escape, latent reservoirs, or multiple infected-cell stages. In particular, antigenic diversity and immune escape, as well as HIV latency and reactivation, are important mechanisms that may alter long-term persistence and post-control resurgence [27, 28]. Fourth, the optimal-control formulation is deterministic and finite-horizon. Future extensions should incorporate stochastic variability, adherence fluctuations, pharmacokinetic–pharmacodynamic drug exposure, and patient-specific treatment response. Finally, the resurgence quantities reported here should be interpreted as theoretical indicators of delayed infection resurgence, not as direct clinical predictions of viral rebound in individual patients.

Future research should extend the present framework in several directions. A first direction is to calibrate the model using patient-specific or experimental data so that delay parameters and immune-response parameters can be estimated more accurately. A second direction is to include latent reservoirs, stochastic reactivation, immune exhaustion, and adherence variability to obtain a more clinically realistic description of post-treatment infection resurgence. A third direction is to couple the control formulation with pharmacokinetic–pharmacodynamic drug models so that treatment delays can be linked directly to drug exposure and antiviral effect. Finally, stochastic or data-informed versions of the model could be developed to quantify uncertainty in delay thresholds, Hopf bifurcation points, and post-control outcomes.

In conclusion, this study provides a delay-aware mathematical framework for analyzing within-host HIV dynamics, immune control, and intervention timing. By distinguishing intracellular infection delay, immune-activation delay, and treatment-response delay, the model clarifies how delayed biological and pharmacological processes influence threshold behavior, stability switching, oscillatory infection patterns, and control robustness. The findings suggest that effective HIV intervention strategies should be assessed not only by their ability to reduce infection during treatment, but also by their robustness to delayed response and their ability to prevent post-control resurgence.

Generative AI statement

The author(s) declared that Generative AI was used in the creation of this manuscript. Generative AI was used to assist with language editing, grammar checking, and improving clarity and readability of the manuscript. It was not used to generate original research data or perform the mathematical analysis. The author(s) reviewed, edited, and verified all AI-assisted content and take full responsibility for the final version of the 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: Alain Miranville, University of Le Havre, France

Reviewed by: Mukhtiar Khan, National Taiwan University, Taiwan

Zhouhong Li, Yuxi Normal University, China

Krithika Balu, Amrita Vishwa Vidyapeetham University, India