Bifurcation analysis of HIV infection model with cell-to-cell transmission and non-cytolytic cure

Surya Prakash; Prashant K. Srivastava; Anuj Kumar Umrao

doi:10.1515/cmb-2023-0111

Open Access Published by De Gruyter Open Access December 31, 2023

Bifurcation analysis of HIV infection model with cell-to-cell transmission and non-cytolytic cure

Surya Prakash , Prashant K. Srivastava and Anuj Kumar Umrao

From the journal Computational and Mathematical Biophysics

https://doi.org/10.1515/cmb-2023-0111

Abstract

A mathematical model is proposed and discussed to study the effect of cell-to-cell transmission, the non-cytolytic process, and the effect of logistic growth on the dynamics of HIV in vivo. The model system consists of one disease-free steady state and another endemic steady state. The disease-free steady state is globally asymptotically stable and the disease eradicated if the basic reproduction number is smaller than one. However, the endemic steady state is globally stable under specific parametric conditions, when it exists. At R 0 = 1 , the forward transcritical bifurcation is obtained. Also, by considering proliferation rate as bifurcation parameter, we get Hopf and Hopf–Hopf bifurcations. We have performed numerical simulations using MATLAB to support our analytical results and show the effects of cell-to-cell infection, proliferation rate, and non-cytolytic cure on all three populations. In the end, we have performed data fitting and note the same behaviour of observed data with predicted data.

Keywords: cell-to-cell transmission; non-cytolytic cure; global stability; basic reproduction number

MSC 2010: 34D20; 34D23; 34C23; 92C37; 92D30

1 Introduction

AIDS caused by human immunodeficiency virus (HIV) is a devastating disease affecting millions globally. HIV is a single-stranded RNA virus that attacks CD4+ T cells, macrophages, and dendritic cells, which leads to HIV infection. AIDS is the last phase of HIV infection. The primary infection with HIV starts via blood, semen, vaginal fluids, and breast milk when these fluids transfer from one HIV-infected person to another healthy person. Once a virus attacks the host cell, it inserts into it and releases its genetic material, i.e., single standard RNA and other enzymes, and then reverse transcription, integration, transcription, and assembly occur. A new virus comes out from the host cells [18,33]. HIV infection occurs in three phases, acute stage, asymptomatic stage, and symptomatic HIV infection. In the first stage, a person can feel flu-like symptoms. This time cytotoxic T lymphocytes (CTLs) control the viral load, which peaks and then declines. In the second stage, a person may have no symptoms from 2 to 10 years or more. In the third stage, person’s immune system is damaged and leads to AIDS (CD4+ T cells count below less than 200). Once the immune response knows the penetration of HIV, then CD4+ T cells become active and send signals to CD8+ T cells to destroy the free virus and infected cells Nowak and coworkers [33].

The following three-dimensional ordinary differential equation system gives the basic viral infection model studied by Nowak and coworkers [18,19].

(1.1) d T d t = Λ − k V T − μ T , d T ⋆ d t = k V T − δ T * , d V d t = N δ T * − c V ,

where T , T * , and V represent the number of uninfected CD4+ T cells, infected CD4+ T cells, and virus, respectively, at any time t . Λ , k , and N are the production rate of uninfected CD4+ T cells, the rate of infection of CD4+ T cells with the virus, and the burst size, respectively. The mortality rates of uninfected cells, infected cells, and free virus are μ , δ , and c , respectively.

Since the development of the model system (1.1), there has been a significant focus on understanding and representing the dynamic behaviours and mechanisms involved in the viral infection process. Several modifications/generalizations of mathematical model (1.1) have been proposed to study the dynamics of HIV [2–4,25]. These models have been used to explain various aspects of HIV dynamics as disease progression, viral replication, eclipse phase, excess drug use, delayed model, and mode of transmission. CD4+ T cells are continuously generated from precursor cells at a constant rate and have a natural death rate. When antigen is stimulated, CD4+ T cells undergo proliferation with a rate r following logistic growth. To make the model more realistic, we need to incorporate logistic proliferation terms for the uninfected target cells [1,12,30]. There are two types of immune response, cytolytic immune response and non-cytolytic immune response. In the cytolytic immune response, immune response destroys infected cells, whereas infected cells can be cured in the non-cytolytic immune response [7,8,20,35]. The dynamics of HIV infection within host has been investigated using a variety of mathematical models [21,24,25]. The virus-to-cell mode of infection was the main focus of all of these models. Viral infection, however, can also happen directly between cells. The processes underlying cell-to-cell spread one yet to be fully understood. Virological synapses allow viruses to spread from infected cells to other infected cells [10]. Comparatively to cell-to-virus infection, cell-to-cell transmission may provide more favourable conditions for viral infection. For example, viruses transmitted from cell to cell may be less likely to be eliminated by CD8+ T cells or neutralized by neutralising antibodies [16]. By incorporating cell-to-cell transmission mechanisms, we can better understand HIV infection in vivo. In studies by Wang et al. [29,32], the mathematical study of virus models with cell-to-cell transmission has been carried out. It was established by Sigal et al. that the key mode of HIV infection is cell-to-cell transmission [23]. The global dynamics of an HIV infection that includes cell-to-cell transmission have been addressed by Li and Wang [12]. In the meanwhile, Lin et al. [14] investigated how HIV dynamics was impacted by cell-to-cell transmission. More recently, age-structured infection models with viruses that contained cell-to-cell transmission have been studied in the study by Wang et al. [31] and established the threshold condition for global behaviour.

Motivated by the aforementioned work, we are concerned about the combined impact of virus-to-cell and cell-to-cell transmission on the dynamics of HIV infection. Here, we study a primary infection model with logistic growth of CD4+ T cell, two modes of transmission (cell-to-cell [36] and virus-to-cell), and infected cells that can be cured by a non-cytolytic process [34] and revert to the uninfected cell population. We establish local and global behaviour of the disease-free steady-state and infected steady state, establish transcritical bifurcation at R 0 = 1 , and also show the effect of different parameters on uninfected CD4+ T cell population, infected CD4+ T cell population, and virus population. In the end, we estimate model parameters and fit the observed data with CD4+ T cell data.

This article is structured as follows: The model that includes cell-to-cell infection and non-cytolytic cure is described in the next section, and the analysis of local and global stability of disease-free and endemic steady states of the proposed model is expressed in the following section. Conditions for transcritical and Hopf bifurcation are carried out in Section 3. The numerical simulation and parameter estimation are established in Sections 4 and 5, respectively. Finally, we discuss our obtained results in Section 6.

2 Mathematical model

In this section, we formulate the HIV infection model with CD4+ T cells and constant CTL immune response. Let T ( t ) , T * ( t ) , and V ( t ) are the uninfected CD4+ T cells, infected CD4+ T cells, and the virus population at any instant of time t > 0 , respectively. We construct our mathematical model based on the following assumptions:

The growth of CD4+ T cells is logistic with r as its growth rate [1] and carrying capacity K ; the uninfected CD4+ T cells dies naturally at the rate μ ; the natural inflow of CD4+ T cells is represented by Λ .
The bilinear interaction, i.e., k V ( t ) T ( t ) and α T ( t ) T * ( t ) , represents the infection of CD4+ T cells with virus and cells, respectively, where k and α represent the rate of infection of CD4+ T cells with virus and cells, respectively.
The infected CD4+ T cells are cured with rate q due to constant CTL immune response E , i.e., q E T * ( t ) ; infected CD4+ T cells dies naturally at the rate δ . It is natural to assume μ ≤ δ .
As more CD4+ T cells become infected, more free virus is produced, N δ T * ( t ) , where N is the average number of free virus produced by an infected CD4+ T cells; free virus dies naturally at the rate c .

Combining all the aforementioned said assumptions lead to the following model:

(2.1) d T d t = Λ − k V T − μ T + r T 1 − T + T * K − α T T * + q E T * , d T ⋆ d t = k V T − δ T * + α T T * − q E T * , d V d t = N δ T * − c V

with initial conditions T ( 0 ) = T 0 > 0 , T * ( 0 ) = 0 , and V ( 0 ) = V 0 ≥ 0 . All parameters are positive and explained in Table 1.

Table 1

List of default values for parameters with sources

Parameter	Description	Value and unit	Source
Λ	Production rate of uninfected cells	10 mm ‒ 3 day ‒ 1	[26]
k	Virus-to-cell infection rate	0.0000024 mm 3 virion ‒ 1 day ‒ 1	[22]
α	Infection rate from cell to cell	0.0000006 mm 3 cell ‒ 1 day ‒ 1	Assumed
μ	Rate of mortality of uninfected cells	0.04 day ‒ 1	Assumed
δ	Rate of mortality of infected cells	1 day ‒ 1	[1]
q	Non-cytolytic cure rate	0.01 day ‒ 1	Assumed
K	Carrying capacity	5,000–10,000 cells mm ‒ 3	Assumed
N	Burst size	1,000	[27]
r	Proliferation rate of CD4+ T-cell	0.03 – 3 day ‒ 1	[22]
c	Clearance of free virus	0.40 day ‒ 1	Asuumed
E	Constant CTLs response	0.01 day ‒ 1	Asuumed

2.1 Positivity and boundedness

The populations should always be non-negative, much like in biological models. Therefore, for the population’s positivity in our model, we assume that the starting population size is selected in such a way that all population components remain positive over all future time. From the model equation (2.1), we derive the following:

d T d t { T = 0 , T * > 0 , V > 0 } = Λ + q E T * > 0 , d T * d t { T > 0 , T * = 0 , V > 0 } = k T V > 0 , d V d t { T > 0 , T * > 0 , V = 0 } = N δ T * ≥ 0 .

On the boundary plane of the cone R + 3 , the aforementioned rates are non-negative. As the vector fields on the boundary planes are oriented inward, we can say that if we begin in the interior of the cone, the solutions will always be contained within that positive cone of R + 3 . So long as the initial conditions are positive, the solutions to the model system (2.1) remain positive and corresponding the biologically feasible region of the model system (2.1) is Γ .

Γ = ( T , T * , V ) ∈ R + 3 : 0 ≤ T + T * ≤ K , V ≤ N δ K c .

2.2 Basic reproduction number calculation

By using Vanden Driessche and Watmough proposed work on next-generation matrix approach [28], we determine the basic reproduction number, for the model system (2.1), which is given as follows:

R 0 = α c + k N δ c ( δ + q E ) T 0 .

2.3 Feasibility and stability (local and global) of system’s steady states

In this section, we discuss the existence of steady states of the model system (2.1) and their stability. The model system (2.1) has the following two steady states:

The disease-free steady state
E 0 ≔ ( T 0 , 0 , 0 ) = K ( r − μ ) + ( μ − r ) 2 + 4 r Λ K 2 r , 0 , 0 ,
which always exists, provided r > μ . This steady state represents a situation where infection is not present. It is to be noted that, we assume r > μ , throughout our analysis, where r and μ are defined in Table 1.
The endemic steady state E 1 ≔ ( T ¯ , T ⋆ ¯ , V ¯ ) , where T ¯ , T ⋆ ¯ , and V ¯ are the positive solutions of the following isoclines:
Λ − k V ¯ T ¯ − μ T ¯ + r T ¯ 1 − T ¯ + T * ¯ K − α T ¯ T * ¯ + q E T * ¯ = 0 , k V ¯ T ¯ − δ T * ¯ + α T ¯ T * ¯ − q E T * ¯ = 0 , N δ T * ¯ − c V ¯ = 0 .
From the third and second equations, we have
V ¯ = N δ T * ¯ c , T ¯ = c ( δ + q E ) k N δ + c α .
Further, we put the value of V ¯ and T ¯ in the first equation, so that
T * ¯ = 1 − 1 R 0 ( Λ R 0 K + r T 0 2 ) δ K R 0 + r T 0 .
Therefore, the unique endemic steady state is given by
(2.2) E 1 ≔ ( T ¯ , T * ¯ , V ¯ ) = T 0 R 0 , 1 − 1 R 0 ( Λ R 0 K + r T 0 2 ) δ K R 0 + r T 0 , N δ c 1 − 1 R 0 ( Λ R 0 K + r T 0 2 ) δ K R 0 + r T 0 ,
provided R 0 > 1 .

The following theorems address the local and global stability of the two steady states present in the model system (2.1).

Theorem 2.1

The disease-free steady state E 0 is locally stable if R 0 < 1 and unstable if R 0 > 1 .
The endemic steady state E 1 exists when R 0 > 1 and is locally asymptotically stable if the conditions A > 0 , B > 0 , C > 0 , and A B − C > 0 holds, where A , B , and C are defined in the proof.

Proof

The Jacobian matrix at a steady state is given by,

J = − k V + r − μ − 2 r T K − r T * K − α T * − r T K − α T + q E − k T k V + α T * − δ + α T − q E k T 0 N δ − c .

The characteristic equation of the Jacobian matrix at E 0 is given by
r − μ + 2 r T 0 K − λ ( λ 2 + ( δ + c + q E − α T 0 ) λ + ( δ c + q E c ) ( 1 − R 0 ) ) = 0 .
Clearly, all roots of the characteristic equation of J ( E 0 ) are with negative real part when R 0 < 1 , Hence, E 0 is locally asymptotically stable when R 0 < 1 .
The characteristic equation of J ( E 1 ) is,

(2.3) λ 3 + A λ 2 + B λ + C = 0 ,
where
A = μ + c + δ + k V ¯ + α T * ¯ + ( 2 T ¯ + T ¯ * ) r K + q E ¯ − α T ¯ − r , B = c ( δ + q E ) + ( μ − r ) ( c + δ + q E ) + α ( r − μ − c ) − k N δ + 2 r K ( c + δ + α + q E ) T ¯ + k ( c + δ ) V ¯ + α c + r K ( δ + q E − α ) T * ¯ + α r K ( T ¯ T ¯ * ) + k r V ¯ T ¯ K , C = k N δ r + α k N δ + 2 r q E c K − k N δ μ − α r c − α μ c T ¯ − k r c V ¯ T ¯ − k r N δ K + α 2 c T ¯ T * ¯ + − k N δ α + r c δ K − r α c k + r q E c K + q E c α T * ¯ + ( c δ k + q E c k ) V ¯ + c δ μ + q E μ c − q E r c .
If A > 0 , B > 0 , C > 0 , and A B − C > 0 , then all inequalities of the Routh-Hurwitz criterion are satisfied. Therefore, the infected steady state E 1 is locally asymptotically stable for R 0 > 1 .□

Theorem 2.2

The disease-free steady state E 0 is globally asymptotically stable in Γ if R 0 < 1 .

Proof

Global stability of disease-free steady state obtained by constructing the Lyapunov function. A Lyapunov function of the model system (2.1) is defined as follows:

L = α c + k N δ k ( δ + q E ) T * + V .

Differentiating L with respect to t along the solution of model system (2.1), we obtain

L ˙ = α c + N δ k c ( δ + q E ) T − 1 α c T * k + c V .

Clearly, L ˙ ≤ 0 in Γ , since R 0 ≤ 1 and T ( t ) ≤ T 0 . Furthermore, if N is the set of solutions of the model system (2.1), where L ˙ = 0 if and only if T * = 0 , V = 0 , T = T 0 on R 0 = 1 . In both cases, N = { E 0 } . Thus, the Lyapunov–LaSalle principle [11] dictates all paths in Γ approach E 0 , and hence the global stability of disease-free steady state E 0 .□

Next, in the following theorem, we discuss the global stability of E 1 for the model system (2.1) using the method developed by Li and Muldowney [13].

Theorem 2.3

The unique endemic steady state E 1 of the model system (2.1), exists if R 0 > 1 and is globally asymptotically stable when Δ < 0 , where Δ is defined in the proof.

Proof

The Jacobian matrix J of the model system (2.1) is given as follows:

J = − k V + r − μ − 2 r T k − r T * K − α T * − r T K − α T + q E − k T k V + α T * − δ + α T − q E k T 0 N δ − c .

The second additive compound matrix that corresponds to it [17] J [ 2 ] is expressed as follows:

J [ 2 ] = L − δ + α T − q E − α T * k T k T N δ A − α T * − C − r T K − α T + q E 0 k V + α T * − δ + α T − q E − c .

where L = − k V + r − μ − 2 r T K − r T * K . Now, consider the function

Q = Q ( T , T * , V ) = 1 0 0 0 T * V 0 0 0 T * V , and Q f = 0 0 0 0 T ˙ * V − T * V ˙ V 2 0 0 0 T ˙ * V − T * V ˙ V 2 ,

we obtain

Q f Q − 1 = 0 0 0 0 T ˙ * T * − V ˙ V 0 0 0 T ˙ * T * − V ˙ V ,

and

B = Q f Q − 1 + Q J [ 2 ] Q − 1 = L − δ + α T − q E − α T * k T V T * k T V T * N δ T * V T ˙ * T * − V ˙ V + L − α T * − c − r T K − α T + q E 0 k V + α T * T ˙ * T * − V ˙ V − δ + α T − q E − c = B 11 B 12 B 21 B 22 ,

where

B 11 = L − δ + α T − q E − α T * , B 12 = k T V T * k T V T * , B 21 = N δ T * V 0 , and

B 22 = T ˙ * T * − V ˙ V + L − α T * − c − r T K − α T + q E k V + α T * T ˙ * T * − V ˙ V − δ + α T − q E − c .

Further, we define Lozinskii measure corresponding to matrix B as follows:

ν ( B ) ≤ max { g 1 , g 2 } ,

where g 1 = ν ( B 11 ) + ∥ B 12 ∥ and g 2 = ∥ B 21 ∥ + ν ( B 22 ) ( ∥ B 12 ∥ and ∥ B 21 ∥ are vector norm). Hence, we have

ν ( B 11 ) = − k V + r − μ − 2 r T K − r T * K − δ + α T − q E − α T * , ν ( B 22 ) = max T ˙ * T * − V ˙ V + L − α T * − c + k V + α T * , M + T ˙ * T * − V ˙ V − δ + α T − q E − c = T ˙ * T * − V ˙ V + r − μ − 2 r T K − r T * K − c , since μ ≤ δ ,

where M = − r T K − α T + q E . Now, from the differential equations in T * and V of the model system (2.1), we have

T ˙ * T * = k T V T * − δ − α T − q E , N δ T * V = V ˙ V + c .

Hence,

g 1 = T * ˙ T * − k V + r − μ − 2 r T K − r T * K − α T * , g 2 = T ˙ * T * + r − μ − 2 r T K − r T * K .

Therefore,

μ ( B ) ≤ max ( g 1 , g 2 ) , μ ( B ) ≤ T ˙ * T * + Δ ,

where Δ = r − μ . Furthermore, let ( T ( t ) , T * ( t ) , V ( t ) ) be any solution starting in Γ and let t ¯ be the consistent time as such ( T ( t ) , T * ( t ) , V ( t ) ) belongs to any compact subset of Γ for all t ≥ t ¯ , and then for t > t ¯ , we have

1 t ∫ 0 t ν ( B ) d s ≤ 1 t ∫ 0 t ν ( B ) d s + 1 t ln T * ( t ) T * ( t ¯ ) + t − t ¯ t Δ .

In the end, q ¯ 2 = limsup t → ∞ sup x 0 ∈ E 1 t ∫ 0 t ν ( B ) d s ≤ Δ < 0 is obtained from the definition of q ¯ 2 and the boundedness of T * ( t ) . This proves the global stability of E 1 .□

3 Bifurcation analysis

Now, we explore the emergence and disappearance of steady states and the changes in their stability through various types of local bifurcations. Notably, the model system (2.1) demonstrates local bifurcations, including transcritical bifurcation and Hopf bifurcation. We present the conditions of these local bifurcations in the following theorems.

3.1 Transcritical bifurcation

Theorem 3.1

The model system (2.1) undergoes a forward transcritical bifurcation, as R 0 crosses 1.

Proof

We prove bifurcation at R 0 = 1 with the help of the approach that Castillo-Chavez and Song [5] described. we choose the virus to cell transmission rate k as the bifurcation parameter and by R 0 = 1 , we have

k = k * = c ( δ + q E ) − α c T 0 N δ T 0 > 0 .

Thus, the model system (2.1) can be written as follows:

d T d t = Λ − k V T − μ T + r T 1 − T + T ⋆ K − α T T ⋆ + q E T ⋆ ≔ f 1 , d T ⋆ d t = k V T − δ T ⋆ + α T T ⋆ − q E T ⋆ ≔ f 2 , d V d t = N δ T ⋆ − c V ≔ f 3 .

Now, the Jacobian matrix of the model system (2.1) evaluated at E 0 and k * is given by

J ( E 0 , k * ) = r − μ − 2 r T 0 K − r T 0 K − α T 0 + q E − k * T 0 0 − δ + α T 0 − q E k * T 0 0 N δ − c .

This Jacobian matrix J ( E 0 , k * ) has simple zero eigenvalue at R 0 = 1 , and all other eigenvalue of J ( E 0 , k * ) has negative real parts. Further,

w = ( w 1 , w 2 , w 3 ) T = k * T 0 ( δ K + r T 0 ) ( − δ + α T 0 − q E ) ( 2 r T 0 + K μ − r K ) , − k * T 0 − δ + α T 0 − q E , 1 T ,

and v = ( v 1 , v 2 , v 3 ) = 0 , 1 , k * T 0 c are the right and left eigenvectors of J ( E 0 , k * ) corresponding to the eigenvalue 0, respectively. The coefficients a and b are defined in the study by Castillo-Chavez and Song [5] are given as follows:

a = ∑ i , j , l = 1 3 v l w i w j ∂ 2 f l ∂ x i ∂ x j ( E 0 , K * ) , b = ∑ i , l = 1 3 v l w i ∂ 2 f l ∂ x i ∂ K ( E 0 , k * ) ,

where the second-order partial derivatives of f 1 , f 2 , and f 3 which are nonzero at J ( E 0 , k * ) are given by

∂ 2 f 1 ∂ T 2 = − 2 r K , ∂ 2 f 1 ∂ T ∂ T * = − r K − α , ∂ 2 f 1 ∂ T ∂ V = − k * , ∂ 2 f 2 ∂ T ∂ T * = α , ∂ 2 f 2 ∂ T ∂ V = k * , ∂ 2 f 2 ∂ V ∂ k = T 0 .

Hence, we obtain

a = − ( c ( δ + q E ) − α c T 0 ) 2 ( δ k * + r T 0 ) ( δ + q E ) N 2 δ 2 T 0 ( − δ + α T 0 − q E ) 2 ( 2 r T 0 + K μ − r K ) , and b = T 0 .

Clearly, a < 0 and b > 0 . Thus, the model system (2.1) undergoes a forward transcritical bifurcation at the disease-free equilibrium E 0 .□

3.2 Existence of Hopf bifurcation at endemic steady state

Here, we examine the occurrence of Hopf bifurcation around the endemic steady-state E 1 . Specifically, we consider r (proliferation rate of CD4+ T cell) as the bifurcation parameter while keeping the other parameters constant, and the coefficients of the characteristic equation ∣ J ( E 1 ) − λ I 3 × 3 ∣ = 0 is a continuously differentiable function of r . If the following criteria are met, Hopf bifurcation exists around the endemic steady state E 1 :

At the endemic steady state E 1 , the Jacobian matrix J ( E 1 ) has two purely imaginary eigenvalues and one negative real eigenvalue
d ( Re ( λ ) ) d r r = r crit ≠ 0 , where r c r i t is the equivalent bifurcation threshold value of r and λ denotes the eigenvalue of J ( E 1 ) at r = r crit .

We apply the criterion developed by Liu [15] to achieve oscillatory behaviour via Hopf bifurcation. If the coefficients of ∣ J ( E 1 ) − λ I 3 × 3 ∣ = 0 satisfy the conditions such as: B ( r crit ) > 0 , C ( r crit ) > 0 , and Δ = ( A B − C ) ( r crit ) = 0 , we obtain a pair of purely imaginary roots of ∣ J ( E 1 ) − λ I 3 × 3 ∣ = 0 . Further, to calculate d ( Re ( λ ) ) d r r = r crit , we consider

(3.1) P ( λ ) ≔ λ 3 + A λ 2 + B λ + C = 0 ,

and λ = ± ι ω is a pair of purely imaginary roots of ∣ J ( E 1 ) − λ I 3 × 3 ∣ = 0 . Further, we differentiate equation (3.1) with respect to r , and we obtain

d P d r = 3 λ 2 d λ d r + 2 λ A d λ d r + B d λ d r + λ 2 d A d r + λ d B d r + d C d r = 0 .

We have

d λ d r = − λ 2 d A d r + λ d B d r + d C d r 3 λ 2 + 2 A λ + B .

By substituting λ = ι ω , we obtain

d λ d r λ = ι ω = − − ω 2 d A d r + ι ω d B d r + d C d r ( − 3 ω 2 + B − 2 A ι ω ) ( − 3 ω 2 + B ) 2 + ( 2 A ω ) 2 ,

Re d λ d r λ = ι ω = − 3 ω 4 d A d r + ω 2 2 A d B d r − 3 d C d r − B d A d r + B d C d r ( B − 3 ω 2 ) 2 + ( 2 A ω ) 2 .

Therefore,

Re d λ d r λ = ι ω ≠ 0 , if 3 ω 4 d A d r + ω 2 2 A d B d r − 3 d C d r − B d A d r + B d C d r ≠ 0 .

The following theorem summarizes the aforementioned discussion.

Theorem 3.2

The Hopf bifurcation through the endemic steady state E 1 occurs under the necessary and sufficient conditions that a critical threshold value of r = r crit exists, such that following specified condition:

B ( r crit ) > 0 , C ( r crit ) > 0 , and Δ = ( A B − C ) ( r crit ) = 0 ,
3 ω 4 d A d r + ω 2 2 A d B d r − 3 d C d r − B d A d r + B d C d r ≠ 0 holds.

3.3 Direction and stability of Hopf bifurcation

Here, we employ normal form theory to determine the direction of Hopf-bifurcation and investigate the stability characteristics of the bifurcating periodic solution in the system described by equation (2.1). This analysis is detailed in the literature [9]. First, we find eigenvector V 1 and V 3 corresponding to the eigen values λ 1 = ι ω and λ 3 = − A , respectively, at r = r crit , where ω = B . Here,

V 1 = b 11 − i b 12 b 21 − i b 22 b 31 − i b 32 , V 3 = b 13 b 23 b 33 ,

with

b 11 = r T ̄ K + δ c − ω 2 r − μ − 2 r T ̄ K − r T * ̄ K − r T ̄ K + δ + c ω 2 N δ r − μ − 2 r T ̄ K − r T * ̄ K 2 + ω 2 b 12 = − r T ̄ K + δ c ω − ω 3 + r − μ − 2 r T ̄ K − r T * ̄ K r T ̄ K + δ + c ω 2 N δ r − μ − 2 r T ̄ K − r T * ̄ K 2 + ω 2 b 21 = c N δ , b 22 = − ω N δ , b 31 = 1 , b 32 = 0 b 13 = − r T ̄ K − δ + A − c + A N δ r − μ − 2 r T ̄ K − r T * ̄ K + A , b 23 = c − A N δ , b 33 = 1 .

We use the following transformation

x = x * + b 11 x 1 + b 12 y 1 + b 13 z 1 , y = y * + b 21 x 1 + b 22 y 1 + b 23 z 1 , z = z * + b 31 x 1 + b 32 y 1 + b 33 z 1

to reduce the given system as

(3.2) d x 1 d t = H 1 , d y 1 d t = H 2 , d z 1 d t = H 3 ,

where

H 1 = b 22 E 1 − b 12 E 2 + ( b 12 b 23 − b 22 b 13 ) E 3 D , H 2 = ( − b 21 + b 23 ) E 1 + ( b 11 − b 13 ) E 2 + ( − b 11 b 23 + b 13 b 21 ) E 3 D , H 3 = − b 22 E 1 + b 12 E 2 + ( b 11 b 22 − b 12 b 21 ) E 3 D ,

with

D = ( b 11 b 22 − b 21 b 12 + b 12 b 23 − b 22 b 13 ) , E 1 = Λ − k ( x * + b 11 x 1 + b 12 y 1 + b 13 z 1 ) ( z * + b 31 x 1 + b 32 y 1 + b 33 z 1 ) − μ ( x * + b 11 x 1 + b 12 y 1 + b 13 z 1 ) + r ( x * + b 11 x 1 + b 12 y 1 + b 13 z 1 ) 1 − x * + b 11 x 1 + b 12 y 1 + b 13 z 1 + z * + b 31 x 1 + b 32 y 1 + b 33 z 1 K − α ( x * + b 11 x 1 + b 12 y 1 + b 13 z 1 ) ( z * + b 31 x 1 + b 32 y 1 + b 33 z 1 ) + q E ( y * + b 21 x 1 + b 22 y 1 + b 23 z 1 ) E 2 = k ( z * + b 31 x 1 + b 32 y 1 + b 33 z 1 ) ( x * + b 11 x 1 + b 12 y 1 + b 13 z 1 ) − δ ( y * + b 21 x 1 + b 22 y 1 + b 23 z 1 ) + α ( x * + b 11 x 1 + b 12 y 1 + b 13 z 1 ) ( z * + b 31 x 1 + b 32 y 1 + b 33 z 1 ) − q E ( z * + b 31 x 1 + b 32 y 1 + b 33 z 1 ) E 3 = N δ ( y * + b 21 x 1 + b 22 y 1 + b 23 z 1 ) − c ( z * + b 31 x 1 + b 32 y 1 + b 33 z 1 ) .

Clearly one equilibrium point of system (3.2) is ( 0 , 0 , 0 ) . So Jacobian matrix of new system (3.2) becomes

J ( E ) = ∂ H 1 ∂ x 1 ∂ H 1 ∂ y 1 ∂ H 1 ∂ z 1 ∂ H 2 ∂ x 1 ∂ H 2 ∂ y 1 ∂ H 2 ∂ z 1 ∂ H 3 ∂ x 1 ∂ H 3 ∂ y 1 ∂ H 3 ∂ z 1

with ∂ H 1 ∂ x 1 = ∂ H 2 ∂ y 1 = ∂ H 1 ∂ z 1 = ∂ H 3 ∂ x 1 = ∂ H 3 ∂ y 1 = ∂ H 2 ∂ z 1 = 0 , − ∂ H 1 ∂ y 1 = ∂ H 2 ∂ x 1 = ω , and ∂ H 3 ∂ z 1 = A . Now we calculate the value of g 11 , g 02 , g 20 , G 101 , G 110 , G 21 , h 11 , h 20 , w , w 20 , w 11 , and g 21 by using the following relations:

g 11 = 1 4 ∂ 2 H 1 ∂ x 1 2 + ∂ 2 H 2 ∂ y 1 2 + i ∂ 2 H 2 ∂ x 1 2 + ∂ 2 H 1 ∂ y 1 2 , g 02 = 1 4 ∂ 2 H 1 ∂ x 1 2 − ∂ 2 H 1 ∂ y 1 2 − 2 ∂ 2 H 2 ∂ x 1 ∂ y 1 + i ∂ 2 H 2 ∂ x 1 2 − ∂ 2 H 2 ∂ y 1 2 + 2 ∂ 2 H 1 ∂ x 1 ∂ y 1 , g 20 = 1 4 ∂ 2 H 1 ∂ x 1 2 − ∂ 2 H 1 ∂ y 1 2 + 2 ∂ 2 H 2 ∂ x 1 ∂ y 1 + i ∂ 2 H 2 ∂ x 1 2 − ∂ 2 H 2 ∂ y 1 2 − 2 ∂ 2 H 1 ∂ x 1 ∂ y 1 , G 21 = 1 8 ∂ 3 H 1 ∂ x 1 3 + ∂ 3 H 1 ∂ x 1 ∂ y 1 2 + ∂ 3 H 2 ∂ x 1 2 ∂ y 1 + ∂ 3 H 2 ∂ y 1 3 + i ∂ 3 H 2 ∂ x 1 3 + ∂ 3 H 2 ∂ x 1 ∂ y 1 2 − ∂ 3 H 1 ∂ x 1 2 ∂ y 1 − ∂ 3 H 1 ∂ y 1 3 , ω = − ∂ H 1 ∂ y 1 , h 11 = 1 4 ∂ 2 H 3 ∂ x 1 2 + ∂ 2 H 3 ∂ y 1 2 , h 20 = 1 4 ∂ 2 H 3 ∂ x 1 2 − ∂ 2 H 3 ∂ y 1 2 − 2 i ∂ 2 H 3 ∂ x 1 ∂ y 1 .

By solving the following equations, we have to find w 11 and w 20

A w 11 = − h 11 , ( D − 2 i ω ) w 20 = − h 20 ,

where

G 110 = 1 2 ∂ 2 H 1 ∂ x 1 ∂ z 1 + ∂ 2 H 2 ∂ y 1 ∂ z 1 + i ∂ 2 H 2 ∂ x 1 ∂ z 1 − ∂ 2 H 1 ∂ y 1 ∂ z 1 , G 101 = 1 2 ∂ 2 H 1 ∂ x 1 ∂ z 1 − ∂ 2 H 2 ∂ y 1 ∂ z 1 + i ∂ 2 H 2 ∂ x 1 ∂ z 1 + ∂ 2 H 1 ∂ y 1 ∂ z 1 , g 21 = G 21 + 2 G 110 w 11 + G 101 w 20 .

We compute the following quantities by using the aforementioned values:

C 1 ( 0 ) = i 2 ω g 20 g 11 − 2 ∣ g 11 ∣ 2 − 1 3 ∣ g 02 ∣ 2 + 1 2 g 21 , μ 2 = − Re { C 1 ( 0 ) } α ′ ( 0 ) β 2 = 2 Re { C 1 ( 0 ) } T 2 = − Im { C 1 ( 0 ) } + μ 2 ω ′ ( 0 ) ω ,

where α ′ ( 0 ) = d d k ( Re { λ 1 ( r ) } ) ∣ r = r crit and ω ′ ( 0 ) = d d r ( Im { λ 1 ( r ) } ) ∣ r = r crit . Here, μ 2 decides the direction of the Hopf-bifurcation. The Hopf-bifurcation is supercritical (sub-critical) if μ 2 > 0 ( < 0 ) ; β 2 decides the stability of the periodic solutions: If β 2 < 0 ( > 0 ) , then the solutions are stable (unstable), and T 2 determines the period of the bifurcating periodic solutions: the period increases (decreases) if T 2 > 0 ( < 0 ) .

4 Numerical simulation

This section uses MATLAB for numerical simulations to back the analytical findings. A hypothetical parameter set meeting conditions and demonstrating dynamic properties is chosen for the simulations. We fix the parameter values as Λ = 10 , μ = 0.04 , δ = 1 , q = 0.01 , N = 500 , c = 0.4 , and E = 0.01 , and other parameter values are varied in the following examples.

Example 1

In this illustration, we numerically validate the existence of forward transcritical bifurcation through which endemic steady state E 1 appears and disease-free steady state E 0 exchanges its stability with E 1 . For this, we set α = 0.0000006 , k = 0.0000024 , K = 8,000 , and r ∈ ( 0.001 , 0.0195 ) so that R 0 ∈ ( 0.76867648327 , 1.38728978183 ) . In this case, disease-free steady state E 0 always exists and loses its stability as R 0 crosses unity from the left, and a unique steady state E 1 appears, which is stable for R 0 > 1 . Figure 1(a) demonstrates the forward transcritical bifurcation as R 0 crosses the threshold value from the left in the considered range of r . Further, we plot phase diagram (see Figure 1(b)) for fixed r = 0.01820240481 (so that R 0 = 1.3161680312 > 1 ). We notice that the disease-free steady state E 0 = ( 438.678814 , 0 , 0 ) is unstable and the solutions curve initiating from different initial conditions ( 76 , 44 , 200 ) , ( 800 , 20 , 100 ) , ( 400 , 50 , 150 ) , and ( 700 , 30 , 50 ) converge to unique steady state E 1 = ( 333.300007 , 2.4802199 , 3100.27482 ) , it suggests that E 1 is locally stable.

Example 2

In this illustration, we validate the existence of Hopf bifurcation through which endemic steady state E 1 loses its stability and periodic oscillation appears. We set α = 0.000001 , k = 0.0000035 , K = 8,000 , and r ∈ ( 0.001 , 0.5 ) . For this set of parameter values, R 0 > 1 , and a unique endemic steady state E 1 exists for the considered range of r . We observe that E 1 loses its stability through Hopf bifurcation at threshold value r crit = 0.2739245332 . At this critical value, there exists a unique endemic steady state E 1 = ( 228.5420 , 61.19429 , 76492.858 ) and the characteristic equation (3.1) is given by

λ 3 + A λ 2 + B λ + C = 0 ,

where A = 1.451479267541721 , B = 0.074374258011956 , and C = 0.107952693543153 . Certainly, A > 0 , B > 0 , C > 0 , and A B − C = − 1.11022 × 1 0 − 16 ≅ 0 . Therefore, the first condition of Theorem 3.2 is satisfied. Now, from the Figure 2(b), we can observe that d ( Re ( λ ) ) d r r = r crit > 0 . Also the eigenvalues of the characteristic equation for r = r crit are − 1.45147927 , and ± 0.2727164 ι . Hence, the transversality condition of Theorem 3.2 is satisfied. Thus, the model system (2.1) undergoes the Hopf bifurcation at the endemic steady state E 1 for the critical value r = r crit . In Figure 2(a), we plot the Hopf bifurcation diagram that depicts the emergence of the limit cycle at the critical value r = r crit . For determining the stability, direction, and nature of periodic oscillation at the Hopf-bifurcation point, we evaluate C 1 ( 0 ) = − 0.385682 × 1 0 − 9 + ι 0.33708 × 1 0 − 10 , μ 2 = 0.2322 × 1 0 − 9 , β 2 = − 0.7714 × 1 0 − 9 , and T 2 = 0.2885 × 1 0 − 9 using algorithms outlined in the direction and stability of Hopf bifurcation (Subsection 3.3). It confirms the occurrence of bifurcating periodic solution. The oscillation exhibits supercritical behaviour and remains stable, indicated by β 2 < 0 and μ 2 > 0 . Further, time series solutions of T , T * , and V at r = 0.3 are plotted in Figure 2(c)–(e). Also, we plot phase diagram to show the stable limit cycle numerically by choosing two different initial conditions: ( 270 , 83 , 81,000 ) and ( 100 , 60 , 60,000 ) (see Figure 2(f)).

$Figure 1 (a) Existence of forward transcritical bifurcation and (b) phase portrait for different initial conditions at r = 0.0182024048 ( R 0 = 1.3161680312 ) r=0.0182024048\hspace{0.33em}\left({R}_{0}=1.3161680312) showing stability of E 1 {E}_{1} .$

Figure 1

(a) Existence of forward transcritical bifurcation and (b) phase portrait for different initial conditions at r = 0.0182024048 ( R 0 = 1.3161680312 ) showing stability of E 1 .

$Figure 2 (a) A bifurcation diagram with a single parameter is used to illustrate how a Hopf bifurcation occurs when the parameter r r crosses r = r crit = 0.2739245 r={r}_{{\rm{crit}}}=0.2739245 in the r T * V r{T}^{* }V -plane, (b) plot for the real part of the roots of the characteristic equation corresponding to E 1 {E}_{1} , (c)–(e) time series plots at r = 0.3 r=0.3 for uninfected CD4+ T cells, infected CD4+ T cells, and virus, respectively, and (f) phase plot shows the existence of stable limit cycle for two different initial conditions: ( 270 , 83 , 81,000 ) \left(270,83,\hspace{0.1em}\text{81,000}\hspace{0.1em}) and ( 100 , 60 , 60,000 ) \left(100,60,\hspace{0.1em}\text{60,000}\hspace{0.1em}) .$

Figure 2

(a) A bifurcation diagram with a single parameter is used to illustrate how a Hopf bifurcation occurs when the parameter r crosses r = r crit = 0.2739245 in the r T * V -plane, (b) plot for the real part of the roots of the characteristic equation corresponding to E 1 , (c)–(e) time series plots at r = 0.3 for uninfected CD4+ T cells, infected CD4+ T cells, and virus, respectively, and (f) phase plot shows the existence of stable limit cycle for two different initial conditions: ( 270 , 83 , 81,000 ) and ( 100 , 60 , 60,000 ) .

4.1 Hopf–Hopf bifurcation and multiple stability switch

In the next example, we numerically validate the sequential changes in the system’s behaviour as a parameter varies, transitioning from a stable steady state to an oscillatory state after losing stability via the first Hopf bifurcation and regaining it through the second Hopf bifurcation.

Example 3

In this illustration, we use α = 0.0000006 , k = 0.0000024 , K = 5,000 , and r ∈ ( 0.00599 , 4 ) so that R 0 ∈ ( 0.87323 , 14.85906 ) . For this set of parameter values, we notice that unique endemic steady state E 1 loses its stability at the through Hopf bifurcation at the threshold value r = r crit 1 = 0.30690043495 . At this critical value of r , the characteristic equation (3.1) is given by

λ 3 + 1.450388093430468 λ 2 + 0.076248032383685 λ + 0.110589238316796 = 0 .

Clearly, every coefficient in the aforementioned equation is positive and also A B − C = 1.45716772 × 1 0 − 15 ≅ 0 . Therefore, at r crit 1 , the characteristic equation has a pair of purely imaginary eigenvalues, i.e., ± 0.276130462614476 ι and one other eigenvalues is − 1.4503881 . From Figure 3(b), we notice that d ( Re ( λ ) ) d r r = r crit 1 > 0 . Thus, both conditions of Theorem 3.2 are satisfied. Hence, the model system (2.1) experiences a Hopf bifurcation at r = r crit 1 around E 1 .

$Figure 3 (a) Diagram illustrating the presence of the Hopf–Hopf bifurcation in the r T * V r{T}^{* }V -plane, (b) plot for the real part of the roots of the characteristic equation corresponding to E 1 {E}_{1} , (c) phase diagram at r = 0.25 < r crit 1 r=0.25\lt {r}_{{{\rm{crit}}}_{1}} , (d) phase diagram at r = r crit 1 r={r}_{{{\rm{crit}}}_{1}} , (e) phase diagram at r = r crit 2 r={r}_{{{\rm{crit}}}_{2}} , and (f) phase diagram at r = 2.4 > r crit 2 r=2.4\gt {r}_{{{\rm{crit}}}_{2}} .$

Figure 3

(a) Diagram illustrating the presence of the Hopf–Hopf bifurcation in the r T * V -plane, (b) plot for the real part of the roots of the characteristic equation corresponding to E 1 , (c) phase diagram at r = 0.25 < r crit 1 , (d) phase diagram at r = r crit 1 , (e) phase diagram at r = r crit 2 , and (f) phase diagram at r = 2.4 > r crit 2 .

Further, if we increase the parameter r , then we observe another Hopf bifurcation at the threshold value r = r crit 2 = 2.288072354 through which endemic steady state E 1 become stable. At this critical value, the coefficient of characteristic equation (3.1) are A = 1.582610352 , B = 0.5372807 , and C = 0.850306001 , which are positive; moreover, A B − C = − 8.8817842 × 10 − 16 ≅ 0 . Thus, the characteristic equation has a pair of purely imaginary roots, i.e., ± 0.73299434 ι , and one other eigenvalue is − 1.582610352 . Also, d ( Re ( λ ) ) d r r = r crit 1 < 0 . Therefore, from Theorem 3.2, we can conclude that the model system (2.1) goes through Hopf bifurcation at the critical value r = r crit 2 around E 1 . Now, to observe this multiple stability switch, we plot the bifurcation diagram in Figure 3(a), which depicts the stability switch from stable state to unstable state through Hopf bifurcation at threshold value r = r crit 1 and again become stable via another Hopf bifurcation as the parameter r crosses the threshold value r = r crit 2 . Figure 3(b) shows the maximum real part of the complex eigenvalues of the Jacobian matrix J ( E 1 ) by varying r ∈ ( 0.00599 , 3.5 ) . Phase portraits at r = 0.25 < r crit 1 , r = r crit 1 , r = r crit 2 , and r = 2.4 > r crit 2 are plotted in Figure 3(c)–(f), respectively.

4.2 Combined effect of cell to cell infection rate and proliferation rate of CD4+ T cell

Now, we are interested in exploring how parameters α and r simultaneously impact the dynamics of the model system (2.1). To do so, we have given the following examples.

Example 4

In this example, we choose K = 5,000 , α ∈ ( 0.000001 , 0.001 ) , and r ∈ ( 0.001 , 3 ) . For this set of parameter values, we notice in Figure 4(a) that α r -parametric space is divided into three different colours (red, yellow, and blue) regions. In the red and yellow colour region disease-free steady state E 0 and endemic steady state E 1 always exist, in which E 0 is always unstable. Endemic steady state E 1 is stable and surrounded by a stable limit cycle in the yellow and red regions, respectively. In the blue colour region, only a globally stable disease-free steady state exists. Further, to observe the dynamics on each different colour region, we plot phase diagram in Figure 4c–e, for fixed value of ( α , r ) = ( 0.000102939 , 1.2508 ) , ( α , r ) = ( 0.000408755 , 1.83712 ) , and ( α , r ) = ( 0.00063302 , 0.001 ) , respectively. Clearly, we can see the oscillatory coexistence of all populations in the red region, stable E 1 in the yellow region, and stable E 0 in the blue region.

Example 5

Finally, in this example, we choose all parameter values the same as in Example 4 except K = 8,000 . We notice that for this set of parameter values instability region (red region) of unique endemic steady state E 1 increases (Figure 4(b)). Therefore, this example exhibits a very interesting phenomenon, i.e., the paradox of enrichment. As the carrying capacity increases, there are higher numbers of uninfected CD4+ T cells. However, this also increases infected CD4+ T cells and free HIV, causing instability. The amplifying process intensifies the spread, potentially leading to population crashes and increased disease severity. To effectively combat HIV, strategies must target both viral replication and infection rate to mitigate unintended consequences and achieve better disease control.

$Figure 4 (a) and (b) Two-parameter stability region in α r \alpha r -plane for K = 5,000 K=\hspace{0.1em}\text{5,000}\hspace{0.1em} and K = 8,000 K=\hspace{0.1em}\text{8,000}\hspace{0.1em} , respectively. The description of all different colour regions is described in the text, (c) phase plot for a fixed value of ( α , r ) = ( 0.000102939 , 1.2508 ) \left(\alpha ,r)=\left(0.000102939,1.2508) in the red colour region, (d) phase plot for a fixed value of ( α , r ) = ( 0.000408755 , 1.83712 ) \left(\alpha ,r)=\left(0.000408755,1.83712) in the yellow colour region, and (e) phase plot for a fixed value of ( α , r ) = ( 0.00063302 , 0.001 ) \left(\alpha ,r)=\left(0.00063302,0.001) in the blue colour region.$

Figure 4

(a) and (b) Two-parameter stability region in α r -plane for K = 5,000 and K = 8,000 , respectively. The description of all different colour regions is described in the text, (c) phase plot for a fixed value of ( α , r ) = ( 0.000102939 , 1.2508 ) in the red colour region, (d) phase plot for a fixed value of ( α , r ) = ( 0.000408755 , 1.83712 ) in the yellow colour region, and (e) phase plot for a fixed value of ( α , r ) = ( 0.00063302 , 0.001 ) in the blue colour region.

4.3 Effect of non-cytolytic cure rate ( q )

In this section, we analyze the effect of non-cytolytic cure rate ( q ). We plot two-parameter stability region by varying parameters r ∈ ( 0.0000001 , 1.5 ) and q ∈ ( 0 , 4 ) (Figure 5(a)) and keeping other parameter the same as in Figure 4a except K = 4,000 , α = 0.0000006 . Clearly, from Figure 5(a), we observe that if we choose r , q from red region (unstable region), then an increase in the parameter q eliminates the oscillations in populations and the endemic equilibrium point E 1 becomes stable. Therefore, if we enhance the non-cytolytic cure rate ( q ) of infected cells, it leads to better control of the viral load, reducing the oscillatory behaviour and stabilizing the infection dynamics around the endemic equilibrium point. This means that the infection is more predictable and does not fluctuate much. For better visualization, we have plotted time series solution for fixed r = 0.85 and three different values of q = 1 , 2.6 , and 3.6.

Figure 5

(a) Two-parameter stability region in r q -plane. The description of all different colour regions is same as shown in Figure 4(a) and (b) Time series solution for different values of q .

5 Parameter estimation

For parameter estimation, we have used MATLAB built-in function “fminsearch” for fitting the data with the numerical solution of the model. We have taken patient data from [6,31] and estimated the different parameters such as virus-to-cell infection rate ( k ) , cell-to-cell infection rate ( α ) , death rate of free virus ( δ ) , and non-cytolytic cure rate q , keeping other parameters fixed at Λ = 10 , μ = 0.04 , K = 5,000 , N = 500 , r = 0.03 , c = 2.4 , and E = 0.01 (Table 1) and simulate our model. Here, we are interested to know how the population of CD4+ T cells progresses as time passes, i.e., we have predicted CD4+ T cell population behaviour for 1300 days. The observed population of CD4+ T cells data is plotted with the blue circles and the solid red colour curve showing the model prediction (Figure 6 and Table 2). We can note that the predicted curve has approximately the same behaviour as the observed data of CD4+ T cell for the fitted set of parameters. Moreover, we found two-parameter values, i.e., ( k , α ) same as appears in Table 1, but we observe different values of δ and q in our parameter estimation, which is given in Table 3.

Figure 6

The solution with the estimated parameters alongside the patient data for the HIV model.

Table 2

Data for CD4+ T-cell [6,31]

Days	CD4+ T-cell count	Days	CD4+ T-cell count
168	629	618	526
179	547	678	513
199	954	735	459
241	578	804	539
261	886	865	462
305	712	935	479
366	438	1,046	492
432	420	1,059	420
464	511	1,183	479
495	401	1,191	347
561	360	1,246	383

Table 3

Estimated parameters

Parameter	Estimate
k	0.0000024
α	0.0000006
δ	2
q	0.001

6 Conclusion

We investigated a mathematical model of HIV infection that incorporates logistic growth of CD4+ T cell, cell-to-cell transmission, and non-cytolytic cure. We demonstrated the model’s positivity and boundedness and identified two steady states: a disease-free steady state and an endemic steady state. Furthermore, we calculate R 0 , i.e., the basic reproduction number, and perform the stability analysis of the model system, which shows that both steady states are locally asymptotically stable, provided Routh-Hurwitz criteria are stratified. We also show the endemic and disease-free steady state’s global stability by building a Lyapunov function and using geometric methods, respectively. We explore the bifurcation behaviour of the model system with respect to the proliferation rate of CD4+ T cells. Moreover, we observed that the model system goes through a forward transcritical bifurcation at the disease-free steady state and Hopf bifurcation and Hopf–Hopf bifurcation around the endemic steady state. We observed fascinating dynamics, i.e., the Paradox of enrichment: As carrying capacity grows, higher numbers of uninfected CD4+ T cells. Also, we have done data fitting and seen the same behaviour of observed data with predicted data.

Therefore, our findings provide valuable insights into the dynamics of HIV infection and may inform the development of strategies for controlling the spread of the virus.

Acknowledgment

The authors are thankful to the anonymous referees for their valuable suggestions which have significantly improved the manuscript.

Funding information: Surya Prakash acknowledges the financial support from the University Grant Commission (UGC) File No. 16-9 (June 2018)/2019(NET/CSIR). Anuj Kumar Umrao acknowledges the financial support from Council of Scientific and Industrial Research (CSIR) under Grant Number: 09/1023(0038)/2020-EMR-I.
Author contributions: The authors declare that they have contributed equally to the creation of this work.
Conflict of interest: The authors state that they have no competing interests to declare and have not engaged in any other activities that may possess a perceived conflict of interest.
Ethical approval: This research did not require ethical approval.

References

[1] Adak, D., Bairagi, N., & Hakl, R. (2020). Accounting for multi-delay effects in an HIV-1 infection model with saturated infection rate, recovery, and proliferation of host cells. BIOMATH, 9(2), 2012297. 10.11145/j.biomath.2020.12.297Search in Google Scholar

[2] Alvarez-Ramirez, J., Meraz, M., & Velasco-Hernandez, J. X. (2000). Feedback control of the chemotherapy of HIV. International Journal of Bifurcation and Chaos, 10(09), 2207–2219. 10.1142/S0218127400001377Search in Google Scholar

[3] Arenas, A. J., González-Parra, G., Naranjo, J. J., Cogollo, M., & DeLaEspriella, N. (2021). Mathematical analysis and numerical solution of a model of HIV with a discrete time delay. Mathematics, 9(3), 257. 10.3390/math9030257Search in Google Scholar

[4] Balasubramaniam, P., Tamilalagan, P., & Prakash, M. (2015). Bifurcation analysis of HIV infection model with antibody and cytotoxic T-lymphocyte immune responses and Beddington-De Angelis functional response. Mathematical Methods in the Applied Sciences, 38(7), 1330–1341. 10.1002/mma.3148Search in Google Scholar

[5] Castillo-Chavez, C., & Song, B. (2004). Dynamical models of tuberculosis and their applications. Mathematical Biosciences & Engineering, 1(2), 361. 10.3934/mbe.2004.1.361Search in Google Scholar PubMed

[6] Djomand, G., Duerr, A., Faulhaber, J. C., Struchiner, C. J., Pacheco, A.G., Barroso, P.F., …, Schechter, M. (2006). Viral load and CD4 count dynamics after HIV-1 seroconversion in homosexual and bisexual men in Rio de janeiro, Brazil. JAIDS Journal of Acquired Immune Deficiency Syndromes, 43(4), 401–404. 10.1097/01.qai.0000243117.21788.90Search in Google Scholar PubMed

[7] Dubey, B., Dubey, P., & Dubey, U. S. (2016). Modeling the intracellular pathogen-immune interaction with cure rate. Communications in Nonlinear Science and Numerical Simulation, 38, 72–90. 10.1016/j.cnsns.2016.02.007Search in Google Scholar

[8] Guidotti, L. G., & Chisari, F. V. (2001). Noncytolytic control of viral infections by the innate and adaptive immune response. Annual Review of Immunology, 19(1), 65–91. 10.1146/annurev.immunol.19.1.65Search in Google Scholar PubMed

[9] Hassard, B., Kazarinoff, N., & Wan, Y. (1981). Theory and applications of hopf bifurcation, Cambridge: Cambridge University Press. Search in Google Scholar

[10] Hübner, W., McNerney, G. P., Chen, P., Dale, B. M., Gordon, R. E., Chuang, F. Y., …, Chen, B. K. (2009). Quantitative 3d video microscopy of HIV transfer across T cell virological synapses. Science, 323(5922), 1743–1747. 10.1126/science.1167525Search in Google Scholar PubMed PubMed Central

[11] LaSalle, J.P. (1976). The stability of dynamical systems. Philadelphia, Pennsylvania: SIAM. 10.21236/ADA031020Search in Google Scholar

[12] Li, F., & Wang, J. (2015). Analysis of an HIV infection model with logistic target-cell growth and cell-to-cell transmission. Chaos, Solitons & Fractals, 81, 136–145. 10.1016/j.chaos.2015.09.003Search in Google Scholar

[13] Li, M. Y., & Muldowney, J. S. (1996). A geometric approach to global-stability problems. SIAM Journal on Mathematical Analysis, 27(4), 1070–1083. 10.1137/S0036141094266449Search in Google Scholar

[14] Lin, J., Xu, R., & Tian, X. (2017). Threshold dynamics of an HIV-1 virus model with both virus-to-cell and cell-to-cell transmissions, intracellular delay, and humoral immunity. Applied Mathematics and Computation, 315, 516–530. 10.1016/j.amc.2017.08.004Search in Google Scholar

[15] Liu, W.-M. (1994). Criterion of hopf bifurcations without using eigenvalues. Journal of Mathematical Analysis and Applications, 182(1), 250–256. 10.1006/jmaa.1994.1079Search in Google Scholar

[16] Martin, N., & Sattentau, Q. (2009). Cell-to-cell HIV-1 spread and its implications for immune evasion. Current Opinion in HIV and AIDS, 4(2), 143–149. 10.1097/COH.0b013e328322f94aSearch in Google Scholar PubMed

[17] Muldowney, J. S. (1990). Compound matrices and ordinary differential equations. The Rocky Mountain Journal of Mathematics, 20, 857–872. 10.1216/rmjm/1181073047Search in Google Scholar

[18] Nowak, M., & May, R. M. (2000). Virus dynamics: Mathematical principles of immunology and virology: Mathematical principles of immunology and virology. UK: Oxford University Press. 10.1093/oso/9780198504184.001.0001Search in Google Scholar

[19] Nowak, M. A., & Bangham, C. R. (1996). Population dynamics of immune responses to persistent viruses. Science, 272(5258), 74–79. 10.1126/science.272.5258.74Search in Google Scholar PubMed

[20] Pan, S., & Chakrabarty, S. P. (2018). Threshold dynamics of HCV model with cell-to-cell transmission and a non-cytolytic cure in the presence of humoral immunity. Communications in Nonlinear Science and Numerical Simulation, 61, 180–197. 10.1016/j.cnsns.2018.02.010Search in Google Scholar

[21] Perelson, A. S., Essunger, P., & Ho, D. D. (1997). Dynamics of HIV-1 and CD4+ lymphocytes in vivo. AIDS (London, England), 11, S17–S24. Search in Google Scholar

[22] Rebecca, V. C., & Ruan, S.. (2000). A delay-differential equation model of HIV infection of CD4+ T-cells. Mathematical Biosciences, 165(1), 27–39. 10.1016/S0025-5564(00)00006-7Search in Google Scholar

[23] Sigal, A., Kim, J. T., Balazs, A. B., Dekel, E., Mayo, A., Milo, R., & Baltimore, D. (2011). Cell-to-cell spread of HIV permits ongoing replication despite antiretroviral therapy. Nature, 477(7362), 95–98. 10.1038/nature10347Search in Google Scholar PubMed

[24] Srivastava, P. (2017). Deterministic and stochastic model for HTLV-I infection of CD4+ T cells. In: Mathematical Biology And Biological Physics (pp. 253–268). Singapore: World Scientific. 10.1142/9789813227880_0015Search in Google Scholar

[25] Srivastava, P., Banerjee, M., & Chandra, P. (2009). Modeling the drug therapy for HIV infection. Journal of Biological Systems, 17(2), 213–223. 10.1142/S0218339009002764Search in Google Scholar

[26] Srivastava, P., & Chandra, P. (2008). Hopf bifurcation and periodic solutions in a dynamical model for HIV and immune response. Differential Equations and Dynamical Systems, 16(1), 77–100. 10.1007/s12591-008-0006-2Search in Google Scholar

[27] Srivastava, P. K., & Chandra, P. (2010). Modeling the dynamics of HIV and CD4. T cells during primary infection. Nonlinear Analysis: Real World Applications, 11(2), 612–618. 10.1016/j.nonrwa.2008.10.037Search in Google Scholar

[28] Vanden Driessche, P., & Watmough, J. (2002). Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Mathematical Biosciences, 180(1–2), 29–48. 10.1016/S0025-5564(02)00108-6Search in Google Scholar

[29] Wang, J., Guo, M., Liu, X., & Zhao, Z. (2016). Threshold dynamics of HIV-1 virus model with cell-to-cell transmission, cell-mediated immune responses and distributed delay. Applied Mathematics and Computation, 291, 149–161. 10.1016/j.amc.2016.06.032Search in Google Scholar

[30] Wang, L., & Li, M.Y. (2006). Mathematical analysis of the global dynamics of a model for HIV infection of CD4+ T cells. Mathematical Biosciences, 200(1), 44–57. 10.1016/j.mbs.2005.12.026Search in Google Scholar PubMed

[31] Wang, S., Hottz, P., Schechter, M., & Rong, L. (2015). Modeling the slow CD4+ T cell decline in HIV-infected individuals. PLoS Computational Biology, 11(12), e1004665. 10.1371/journal.pcbi.1004665Search in Google Scholar PubMed PubMed Central

[32] Wang, X., Tang, S., Song, X., & Rong, L. (2017). Mathematical analysis of an HIV latent infection model including both virus-to-cell infection and cell-to-cell transmission. Journal of Biological Dynamics, 11(sup2), 455–483. 10.1080/17513758.2016.1242784Search in Google Scholar PubMed

[33] Weiss, R. A. (1993). How does HIV cause AIDS? Science, 260(5112), 1273–1279. 10.1126/science.8493571Search in Google Scholar PubMed

[34] Wodarz, D., Christensen, J. P., & Thomsen, A. R. (2002). The importance of lytic and nonlytic immune responses in viral infections. TRENDS in Immunology, 23(4), 194–200. 10.1016/S1471-4906(02)02189-0Search in Google Scholar

[35] Zanoni, M., Palesch, D., Pinacchio, C., Statzu, M., Tharp, G. K., Paiardini, M., …, Silvestri, G. (2020a). Innate, non-cytolytic CD8+ T cell-mediated suppression of HIV replication by MHC-independent inhibition of virus transcription. PLoS Pathogens, 16(9), e1008821. 10.1371/journal.ppat.1008821Search in Google Scholar PubMed PubMed Central

[36] Zhang, T., Meng, X., & Zhang, T. (2015). Global dynamics of a virus dynamical model with cell-to-cell transmission and cure rate. Computational and Mathematical Methods in Medicine, 2015, 758362. 10.1155/2015/758362Search in Google Scholar PubMed PubMed Central

Received: 2023-08-24

Revised: 2023-12-03

Accepted: 2023-12-03

Published Online: 2023-12-31

This work is licensed under the Creative Commons Attribution 4.0 International License.

Bifurcation analysis of HIV infection model with cell-to-cell transmission and non-cytolytic cure

Abstract

1 Introduction

2 Mathematical model

2.1 Positivity and boundedness

2.2 Basic reproduction number calculation

2.3 Feasibility and stability (local and global) of system’s steady states

Theorem 2.1

Proof

Theorem 2.2

Proof

Theorem 2.3

Proof

3 Bifurcation analysis

3.1 Transcritical bifurcation

Theorem 3.1

Proof

3.2 Existence of Hopf bifurcation at endemic steady state

Theorem 3.2

3.3 Direction and stability of Hopf bifurcation

4 Numerical simulation

Example 1

Example 2

4.1 Hopf–Hopf bifurcation and multiple stability switch

Example 3

4.2 Combined effect of cell to cell infection rate and proliferation rate of CD4+ T cell

Example 4

Example 5

4.3 Effect of non-cytolytic cure rate ( q )

5 Parameter estimation

6 Conclusion

Acknowledgment

References

Journal and Issue

Articles in the same Issue