Revisiting the classical target cell limited dynamical within-host HIV model - Basic mathematical properties and stability analysis

1.   Introduction

1.1.   Motivation

Today, approximately 40 million people are infected with the human immunodeficiency virus (HIV) and nearly 5.5 million people are unaware of it ^[1,2]. For that reason, research on this infectious disease without treatment still can be regarded as an important topic from a biological and clinical point of view.

Since HIV was found to be the main reason for the acquired immune deficiency syndrome (AIDS), many modeling approaches have been explored over the course of the last decades to simulate its time development. At the end of the twentieth century, different approaches such as $\text{CD}4^{+}$ T cell subpopulations ^[3], experimental data ^[4] or simpler fundamental models ^[5,6,7,8,9] were applied to better understand the dynamics of primary HIV infections. Reviews and fundamental models were proposed at the beginning of the twenty-first century ^[10,11,12]. Afterward, some works on global stability of fundamental models on viral dynamics were published ^[13,14,15]. Later, main fundamental models were reviewed in ^[16,17,18]. Furthermore, different aspects such as drug therapy or treatment can be implemented to obtain realistically dynamical models ^{[19,20,21,22,23]}. Furthermore, agent-based models can be applied as an alternative ^[24,25]. In this work, we specifically consider the well-known, classical target cell within-host HIV model

where all model parameters, also called constant problem parameters, and all variables, also known as solution components, are described in Table 1; quantities with index $0$ represent initial conditions. We briefly want to remark that the state variable $T_{i} \left(t \right)$ indicates infected ${\rm{CD}}4^{+}$ T cells and the index $i$ reflects this infected state.

As later explained, this model accurately describes viral load during primary HIV infection in the acute phase. Additionally, our ideas for proofs of basic mathematical properties, given in the later part of this work, can be transferred to more complex models. For these reasons, we choose to mainly concentrate our mathematical examination on model (1.1), although we are aware of dynamical models including more tissue and mechanisms, as presented in our previous discussion.

Throughout this article, all model parameters are assumed to be positive. Units are taken from the work by Alizon and Magnus ^[17]. Further, we want to shed some light on the dynamical system's structure. $r$ represents a constant production rate of target ${\rm{CD}}4^{+}$ T cells while the linear term $-d \cdot T \left(t \right)$ stands for the constant elimination of ${\rm{CD}}4^{+}$ T cells due to non-disease natural reasons. The main non-linear term of this first-order, non-linear dynamical system is given by $- \beta \cdot T \left(t \right) \cdot V \left(t \right)$, which models the reduction or target ${\rm{CD}}4^{+}$ T cells by being infiltrated by virus particles. Consequently, $- \delta \cdot T_{i} \left(t \right)$ and $- c \cdot V \left(t \right)$ model death or elimination processes of infected ${\rm{CD}}4^{+}$ T cells or of virions. Hence, system (1.1) describes a non-linear dynamical system of first order. Average lifetimes of infected ${\rm{CD}}4^{+}$ T cells and of virus particles are represented by $\frac{1}{d}$ and $\frac{1}{c}$, as described by Nowak and Bangham ^[6]. For further details regarding biological interpretations, we refer interested readers to ^[5,6,7].

1.2.   Literature review

Let us present a concise history on mathematical modeling in primary HIV infection and related fields. In 1994, Essunger and Perelson proposed a model of HIV infection of ${\rm{CD}}4^{+}$ T cell subpopulations ^[3]. Their main mathematical interest was the possible stability of steady states, also known as equilibrium points. In 1996, Kirschner presented one simplified model for primary HIV infection similar to model (1.1), but modified by a Michaelis-Menten mechanism ^[5]. In her work, she mainly numerically investigated steady states or linear growth of virus particles in time. In the same year, Nowak and Bingham suggested different models, including (1.1) in ^[6] and they also mainly investigated plateau-phase equlibrium points. For example, they extended our basic, i.e., very fundamental, model by immune responses to virus particles. One year later, Bonhoeffer and co-authors discussed model (1.1) and modified versions in which they emphasized the discussion of equilibria and numerical simulations ^[7] (compare page 6971 and especially the section about "A Basic Model" in this reference). Finzi and Siliciano also discussed equilibria for a simplified primary HIV-infection model ^[8]. In 1998, De Boer and Perelson investigated different models similar to (1.1) in ^[9]. However, they were mainly concerned with steady states. One year later, Perelson and Nelson primarly summarized analytical investigations on (1.1) regarding stability of equilibria ^[10]. In 2000, Stafford and co-workers used model (1.1) for parameter estimation problems ^[11]. Two years later, Perelson reviewed some classical mathematical models of primary HIV infection ^[12]. In 2004, Korobeinikov introduced different Lyapunov functions for global stability of simple epidemiological and virus dynamical models ^[13,14]. Additionally, Wang and Li proved global stability with respect to a modified mathematical model of primary HIV infection ^[15]. One year later, Ribeiro investigated the basic mathematical model (1.1) numerically ^[16]. In 2012 and 2013, Alizon, Magnus, Perelson and Ribeiro reviewed basic and more sophisticated mathematical models on dynamics of primary HIV infections including, for example, immune responses as seen in ^[6] or different virus strains. In recent years, stability analysis of more sophisticated models on HIV infections was developed ^[20,21,22]. Recently, Xu modified the basic model (1.1) by a constant ${\rm{CD}}8^{+}$ T cell density in blood and examined stability properties of the suggested model ^[23].

The aforementioned articles mainly examined stability of the plateau-phase equilibrium states of mathematical models of HIV infections and the modeling of infection over time. Hence, we aim to provide proofs of existence and uniqueness of solutions of (1.1) globally in time as these are essential properties for biologically plausible models. Additionally, since this model (1.1) properly describes primary infection during the acute phase of HIV ^[11], it seems crucial to present thorough proofs of certain basic mathematical properties. For that reason, we need statements on boundedness and non-negativity of (1.1). Additionally, we want to derive all possible equilibria in a thorough manner. To the best of our knowledge, existence and uniqueness were not proven in the aforementioned articles. Furthermore, different properties such as boundedness, non-negativity, or derivation of the equilibrium points were only mentioned in the aforementioned works. Hence, our main goal is to collect these important properties with mathematical derivations.

1.3.   Contributions

Although system (1.1) is one of the simplest models for primary HIV infection, many studies have focused on applications and modelling with respect to the disease's time development ^{[4,5,10,12,18]}. For that reason, we want to thoroughly investigate and re-examine system (1.1), since a detailed analysis can be seen as a preparatory step for future research. In previous investigations ^{[26,27,28,29,30]}, different systems of differential equations were first analyzed and then numerical algorithms were developped and applied for its numerical solution.

Our two main contributions can be summarized as follows:

$1)$ In Section 2, we mainly prove analytical results with respect to the dynamical system (1.1). Here, we begin with the non-negativity of possible solutions for all $t \geq 0$. Afterward, we demonstrate that all solution components remain bounded for all $t \geq 0$. To establish these results, we need to examine subsystems of (1.1), which is done in the proof of Lemma 3. In addition, we provide results for existence and uniqueness globally in time for all $t \geq 0$ in Sections 2.3 and 2.4. We conclude this section with a stability result for the virus-free equilibrium point and the plateau-phase disease equilibrium point based on the Routh–Hurwitz criterion. Here, the basic reproduction number occurs as a by-product of this criterion. In addition, we obtain the basic reproduction number by considering the approach of van den Driessche ^[31,32,33].

$2)$ By applying typical Runge–Kutta time stepping schemes in Section 3, we investigate our theoretical findings using numerical simulations and thus underline our results.

To summarize, we aim to provide thorough proofs of analytical results in order to stress the biological usefulness of model (1.1).

2.   Analytical results of (1.1)

In this section, we prove some analytical results with respect to (1.1). Time-continuous solution components are assumed because of the dynamical system's structure, since all problem constants are positive and, as a consequence, no jumps in problem parameters exist. For additional interpretation, we refer our readers back to the introduction.

We want to note that all right-hand side functions $f_{j} \left(t, T \left(t \right), T_{i} \left(t \right), V \left(t \right) \right)$ of (1.1) for each index $j \in \left\{ 1, 2, 3 \right\}$ are continuous with respect to all state variables. Since (1.1) can be equivalently written as an integral equation, all state variables are continuously differentiable functions with respect to time.

2.1.   Non-negativity

Here, we discuss the non-negativity of possible solutions for system (1.1). This is of importance since only non-negative solution components of (1.1) have biological relevance.

Lemma 2.1. All solution components of (1.1) remain non-negative for all $t \geq 0$.

Proof. Let us assume that there might be a time where at least one solution component becomes negative. Due to the continuity of all solution components, there exists a time point $t_{0} \geq 0$ where $T \left(t_{0} \right) = 0$, $T_{i} \left(t_{0} \right) = 0$ or $V \left(t_{0} \right) = 0$ hold. Here, it is important to keep in mind that our initial conditions need to be non-negative as stated in (1.1).

Let us first assume that $T \left(t_{0} \right) = 0$ holds while all other solution components are non-negative due to continuity. Then

holds and this implies $T^{\prime} \left(t_{0} \right) > 0$.

Now, let us assume that $T_{i} \left(t_{0} \right) = 0$ holds while all other solution components are non-negative due to continuity. Then

holds and this implies $T_{i}^{\prime} \left(t_{0} \right) \geq 0$.

Let us assume that $V \left(t_{0} \right) = 0$ holds while all other solution components are non-negative due to continuity. Then

holds and this implies $V^{\prime} \left(t_{0} \right) \geq 0$.

Inductively, for later time points where at least one solution component is zero, we can apply the same argument. This means that no state variable can become negative. In conclusion, all solution components remain non-negative for all $t \geq 0$ due to continuity and non-negativity of initial conditions, which finishes our proof.

2.2.   Boundedness of all solution components

In this subsection, we investigate the boundedness of all solution components of system (1.1). For this proof, we need one comparison principle from differential equations, stated in [34, Lemma 3.4].

Lemma 2.2. Consider the scalar differential equation

where $f \left(t, u \right)$ is continuous in $t$ and locally Lipschitz in $u$ for all $t \geq t_{0}$ and all $u \in J \subset \mathbb{R}$. Let $\left[t_{0}, T \right)$ (T could be infinity) be the maximal interval of existence of the solution $u \left(t \right)$. Suppose $u \left(t \right) \in J$ for all $t \in \left[t_{0}, T \right)$. Let $v \left(t \right)$ be a continuous functions whose upper right-hand derivative $D^{+} v \left(t \right)$ satisfies the differential inequality

with $v \left(t \right) \in J$ for all $t \in \left[t_{0}, T \right)$. Then, it holds $v \left(t \right) \leq u \left(t \right)$ for all $t \in \left[t_{0}, T \right)$.

Now, we are able to prove the boundedness of all states variables of (1.1).

Lemma 2.3. All solution components of (1.1) remain bounded for all $t \geq 0$.

Proof. 1) Let us first consider $T^{\prime} \left(t \right) = r - \beta \cdot V \left(t \right) \cdot T \left(t \right) - d \cdot T \left(t \right)$. Due to Lemma 2.1, we conclude that

holds. Consequently, by application of the comparison principle from Lemma 2.2 and by non-negativity, we notice that

is valid for all $t \geq 0$ because all assumptions of Lemma 2.2 are fulfilled in this situation. Finally, this implies

for all $t \geq 0$.

2) By investigating the subsystem

of (1.1) and defining

as the complete number of target ${\rm{CD}}4^{+}$ T cells, we obtain the following differential equation

with initial condition $M_{0} : = T_{0} + T_{i, 0}$. Hence, it follows

for all $t \geq 0$, which implies the validity of

for all $t \geq 0$.

3) Due to the boundedness of $T \left(t \right)$ and $T_{i} \left(t \right)$ for all $t \geq 0$, we obtain the following differential inequality

and we can conclude

for all $t \geq 0$, which has

for all $t \geq 0$ as a consequence.

This proves our assertion that all solution components of (1.1) remain bounded for all $t \geq 0$. This finishes our proof.

2.3.   Existence of all solution components of (1.1) for all $t \geq 0$

Here, we want to give one statement from [35, Theorem 4.7.1] or [36, Theorem 2.2]. We consider a general initial-value problem

where ${\bf z} \left(t \right) = \left(z_{1} \left(t \right), \ldots, z_{n} \left(t \right) \right)^{T}$ denotes our solution vector, the vectorial function is given by ${\bf G} \left(t, {\bf z} \left(t \right) \right) = \left(g_{1} \left(t, {\bf z} \left(t \right) \right), \ldots, g_{n} \left(t, {\bf z} \left(t \right) \right) \right)^{T}$, and initial conditions are given by ${\bf z}_{0} \in \mathbb{R}^{n}$. By $\lVert \cdot \rVert_{\mathbb{R}^{n}}$, we denote a suitable vector norm on $\mathbb{R}^{n}$. Here, we apply the supremum norm

on the space of bounded, continuous functions on $\mathbb{R}^{n}$ since this leads to a Banach space.

Theorem 2.1. If ${\bf G} \, \colon \, \left[0, \infty \right) \times \mathbb{R}^{n} \longrightarrow \mathbb{R}^{n}$ is locally Lipschitz-continuous in both time and state variables and if there are non-negative real functions $D \, \colon \, \left[0, \infty \right) \longrightarrow \left[0, \infty \right)$ and $K \, \colon \, \left[0, \infty \right) \longrightarrow \left[0, \infty \right)$ such that

holds for all ${\bf z} \left(t \right) \in \mathbb{R}^{n}$, then the solution of the initial-value problem (2.7) exists for all $t \geq 0$. Moreover, for every finite $T \geq 0$, we have

for all $t \in \left[0, T \right]$ where

are described.

Now, we are able to show the global existence of all solution components for our dynamical system in (1.1) for all $t \geq 0$. Existence is one main property that reliable models in natural sciences should fulfill. This also holds for uniqueness, later analyzed in this work.

Theorem 2.2. There exists a solution of (1.1) globally in time for all $t \geq 0$.

Proof. We define

for our vectorial function of (1.1).

1) At first, we prove Lipschitz-continuity locally for $g_{1} \left(t, T \left(t \right), T_{i} \left(t \right), V \left(t \right) \right)$ because the other component functions are estimated in a similar manner. By the boundedness of all solution components by Lemma 2.3, we can assume that

exist. We obtain

This result implies

As a consequence, we conclude that our defined vectorial function of (1.1) is locally Lipschitz-continuous.

2) By our assumptions, we obtain

for the first vectorial component,

for the second vectorial component and

for the final vectorial component. Set

Furthermore, we see that

is valid for all $t \geq 0$. This inequality implies the existence of all solution components globally in time for all $t \geq 0$ due to the application of Theorem 2.1.

Hence, our proof is complete.

2.4.   Uniqueness of all solution components of (1.1) for all $t \geq 0$

To show the uniqueness of system (1.1), we need Banach's fixed point theorem; compare [37, Theorem Ⅴ.18].

Theorem 2.3. Let $\left(X, \varrho \right)$ be a complete metric space with metric $\varrho \, \colon \, X \times X \longrightarrow \left[0, \infty \right)$. Let $T \, \colon \, X \longrightarrow X$ be a strict contraction, i.e., there exists a constant $K \in \left[0, 1 \right)$ such that $\varrho \left(Tx, Ty \right) \leq K \cdot \varrho \left(x, y \right)$ holds for all $x, y \in X$. Then, the mapping $T$ has a unique fixed point in $X$.

Now, we are able to formulate our statement on uniqueness of system (1.1) globally in time.

Theorem 2.4. System (1.1) possesses a unique solution globally in time for all $t \geq 0$.

Proof. We first define the equivalent system of integral equations

to system (1.1) for the application of Banach's fixed point theorem.

1) The function space of bounded, continuous functions on the interval $\left[0, \infty \right)$ is a complete metric space with the supremum norm.

2) First, we estimate

by the boundedness of all solution components. If we choose $\tau_{1} \leq \dfrac{1}{2 \cdot \left(\beta \cdot T_{sup} + \beta \cdot V_{sup} + d \right)}$, we obtain

3) Second, we see that

holds by the boundedness of all solution components. If we choose $\tau_{2} \leq \dfrac{1}{2 \cdot \left(\beta \cdot T_{sup} + \beta \cdot V_{sup} + \delta \right)}$, we obtain

4) Third, we notice that

is valid by the boundedness of all solution components. If we choose $\tau_{3} \leq \dfrac{1}{2 \cdot \left(\pi + c \right)}$, we obtain

$5)$ If we choose $\tau \leq \min \left\{ \tau_{1}, \tau_{2}, \tau_{3} \right\}$, we get

as an estimate. Hence, system (1.1) has a unique fixed point on $\left[0, \tau \right]$.

Inductively, we can conclude that this fixed point is unique on every interval $\left[k \cdot \tau, \left(k + 1 \right) \cdot \tau \right]$ for all $k \in \left\{ 0 \right\} \cup \mathbb{N}$. This implies the uniqueness of a solution of (1.1) globally in time for all $t \geq 0$.

2.5.   Stability analysis of equilibrium states

Denote possible equilibrium states by $\left(T^{\star}, T_{i}^{\star}, V^{\star} \right)$. From (1.1), we obtain the system of equations

for possible equilibrium states.

Lemma 2.4. The two possible equilibrium points read

Proof. We want to split our proof into two parts.

1) We can easily check by plugging

into (2.9) that this point is definitely one possible equilibrium state.

2) The third equation of (2.9) yields

Considering the second equation of (2.9), we obtain

and this implies

Looking at the first equation of (2.9), we get

and it follows

As a consequence, it holds

and this results in

Hence, the second possible equilibrium state reads

This proves our proposition.

We note that these two equilibrium points are the same as solely mentioned in ^[11]. Here, we give one statement for locally asymptotic stable equilibria of an autonomous dynamical system

from [35, Theorem 6.1.1].

Theorem 2.5. Suppose that ${\bf b_{\star}}$ is an equilibrium point for ${\bf x}^{\prime} \left(t \right) = {\bf G} \left({\bf x} \left(t \right) \right)$ where ${\bf G} \in \mathcal{C}^{1} \left(\mathcal{U} \right)$ with a domain $\mathcal{U} \subset \mathbb{R}^{d}$. Furthermore, we assume that

holds for all $j \in \left\{ 1, 2, \ldots, d \right\}$ where ${\bf D} {\bf G_{\star}}$ is the Jacobian of ${\bf G}$. Then, there is a neighborhood $\mathcal{V}$ of ${\bf b_{\star}}$ in $\mathbb{R}^{d}$ such that, for any initial data ${\bf b} \in \mathcal{V}$, the initial value problem

has a solution for all $t \geq 0$ and $\lim\limits_{t \to \infty} {\bf x} \left(t \right) = {\bf b_{\star}}$.

We consider a special case of the Routh–Hurwitz criterion; compare ^[38].

Lemma 2.5. Consider the characteristic equation

of a corresponding Jacobian of the linearization of a dynamical system. Its eigenvalues all have negative real parts and it is locally asymptotically stable if and only if

Hence, we can prove that the virus-free equilibrium state $\left(T_{1}^{\star}, T_{1, i}^{\star}, V_{1}^{\star} \right)$ and the plateau-phase equilibrium state $\left(T_{2}^{\star}, T_{2, i}^{\star}, V_{2}^{\star} \right)$ from Lemma 2.4 are locally asymptotically stable. For that reason, we define the basic reproduction number $R_{0} : = \dfrac{\beta \cdot \pi \cdot r}{c \cdot d \cdot \delta}$ of primary HIV infections, which can be regarded as a transition point from the virus-free to the plateau-phase equilibrium state. In ^[6], this threshold was only defined. Here, we want to give a derivation based on an approach by van den Driessche ^[32]. Fore more details, we refer interested readers to that article. We mainly follow concise ideas from ^[33] for a within-host model of COVID-19.

We reorganize (1.1) as follows

where we consider the subsystem

of infected and viral particles. We define two vectors

such that

holds. For the approach by van den Driessche, one needs the virus-free equilibrium state

as later given in Theorem 2.6. By computing both Jacobians

of these vectorial functions at this virus-free equilibrium state, the basic reproduction number is given by

Here, $\varrho$ defines the spectral radius of the considered matrix. Since we obtain

this yields

This basic reproduction number helps us to determine and distinguish the stability of equilibrium states. More specifically, if $R_{0} < 1$ holds, the disease's progress settles to the virus-free equilibrium state while it settles in the disease-plateau-phase equilibrium state if $R_{0} > 1$ holds. This threshold is also an important number in mathematical epidemiology for a disease's spread ^[13,14]. Additionally, although Korobeinikov proved the global stability ^[13,14], we want to show locally asymptotic stability by using the Routh–Hurwitz criterion, which, to the best of our knowledge, cannot be found in the aforementioned articles.

Theorem 2.6. (i) Suppose that

holds. The virus-free equilibrium state

of our dynamical system

from (1.1) is locally asymptotically stable.

(ii) Suppose that

holds. The plateau-phase equlibrium state

of our dynamical system

from (1.1) is locally asymptotically stable.

Proof. We divide our proof into multiple steps.

1)$ Let

Its Jacobian reads

2) We compute the characteristic equation of the previous Jacobian. At first, we obtain

This yields

3)ⅰ. Let us first consider the virus-free equilibrium state $\left(T_{1}^{\star}, T_{1, i}^{\star}, V_{1}^{\star} \right)$ and let $R_{0} < 1$. We can define

Obviously, $a_{0}$ and $a_{1}$ are positive. By plugging in the definition of the virus-free equilibrium point, we conclude that

and

hold. Hence, if $R_{0} < 1$ is valid, we have locally asymptotic stability by the Routh–Hurwitz criterion.

3)ⅱ. Here, we consider the plateau-phase equilibrium state $\left(T_{2}^{\star}, T_{2, i}^{\star}, V_{2}^{\star} \right)$ and let $R_{0} > 1$. Now, we can again define

Obviously, $a_{0}$ and $a_{1}$ are positive. By plugging in the definition of the plateau-phase equilibrium point, we conclude that

and

hold. Hence, all coefficients $a_{j}$ are positive for all $j \in \left\{ 0, 1, 2, 3 \right\}$. Since we want to apply the Routh–Hurwitz criterion from Lemma 2.5, we still need to show

We obtain

because

holds, which shows this assertion.

Since all assumptions of the Routh–Hurwitz criterion from Lemma 2.5 are fulfilled, its application for both cases yields the desired stability results and finishes our proof.

3.   Numerical simulations

We give two numerical examples to provide evidence for our theoretical findings. First, we consider the case $R_{0} > 1$ with all model parameters taken from ^[11,17], which means that we obtain the plateau-phase equilibrium point. Additionally, we also provide the case $R_{0} < 1$, which corresponds to the virus-free equilibrium state.

3.1.   Case 1: $R_{0} > 1$

In this subsection, we apply the following initial conditions and constant problem parameters, taken from ^[11,17] and presented in Table 2. We want to remark again that all constant problem parameters are assumed to be positive. In addition, we note that the initial conditions of (1.1) are non-negative for numerical simulations. Taking all values of Table 2, we can compute the basic reproduction number $R_{0}$, which reads

in this case. Hence, $R_{0} > 1$ holds and we expect an plateau-phase equilibrium state.

Here, we use the standard function ode45 of GNU Octave ^[39]. For further information regarding these modified Runge-Kutta time stepping methods, we refer interested readers to ^[40,41]. Our GNU Octave code can be found in the supplementary file for reproducibility. We must note that we shortened all time vectors and all solution components vectors for the plotting of the figure due to representation problems on the author's computer.

Using the given initial conditions and problem parameters from Table 2, we obtain

for the coordinates of the plateau-phase equilibrium point. Simulation results of (1.1) with parameters from Table 2 can be seen in Figure 1. Here, we can see that the model of primary HIV infection converges towards the plateau-phase disease equilibrium point after a certain amount of time. This shows that system (1.1) is especially appropriate for the acute and asymptotic phase of HIV infections ^[5,17] and it settles into the correct equilibrium state as we can see in Figure 1. Additionally, we also notice that our theoretical results of boundedness and non-negativity for all solution components of system (1.1) hold in numerical simulations. However, we address that preservation of boundedness or non-negativity is not intrinsically fulfilled by explicit time integration methods; compare with ^[30].

3.2.   Case 2: $R_{0} < 1$

Here, we present a case for $R_{0} < 1$. Let us take the following model parameters as presented in Table 3. It holds that

and we notice that the graphs of all solution components converge to the correct equilibrium point. Using all model parameters from Table 3, we obtain

as the correct virus-free equilibrium state, as seen in Figure 2.

4.   Conclusions and outlook

In this work, we examined and re-investigated analytical properties of the classical target cell limited dynamical within-host HIV model (1.1). Here, we focused on important facts such as non-negativity, boundedness, global existence and global uniqueness in time for all solution components. These are all important properties regarding biological relevance of model (1.1). Furthermore, we showed that the virus-free equilibrium point and the plateau-phase disease equilibrium point of (1.1) are locally asymptotically stable. We proved that they can be distinguished by the basic reproduction number $R_{0} : = \dfrac{\beta \cdot \pi \cdot r}{c \cdot d \cdot \delta}$. Finally, we highlighted our thereotical findings by numerical simulations, based on Runge–Kutta time stepping methods, and demonstrated that our results hold true in both cases $R_{0} < 1$ and $R_{0} > 1$.

We want to remark that our basic model (1.1) does not include more complex factors of HIV infections such as different virus strains, immune responses, or other cells such as ${\rm{CD}}8^{+}$ T cells. Hence, it could be of interest to examine the analytical properties of more complex models ^[5,17] that include additional aspects of HIV infections. The model in the work by Kirschner ^[5] is better suited for long-time modeling of HIV infection, since it also includes the rapid decline of ${\rm{CD}}4^{+}$ T cells approximately ten years after infection. It is also worth noting that classical explicit time stepping schemes do not preserve boundedness or non-negativity ^[30,42]. These methods are based on a methodology proposed by Mickens ^[43]. Hence, the development of non-negativity-preserving time integration methods for model (1.1) might be a future research direction, adding to this research article, where we mainly focused on analytical aspects.

Use of AI tools declaration

The author declares that he has not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The author declares there is no conflict of interest.