A Time-Continuous Model for an Untreated HIV-Infection and a Novel Non-Standard Finite-Difference-Method for Its Discretization

1  Introduction

For more than four decades, infections with human immunodeficiency virus (HIV) have become a global epidemic, with approximately 40 million people currently living with this infection [1–3]. However, nearly 5.5 million people do not know about their HIV infection [3]. In the absence of medical treatment of such infections, cellular mechanisms and their time development within human individuals have been well studied and understood [4–7]. Over time, many scientists have developed mathematical models for the evolution of an infection with HIV which are summarized in different reviews [8–10].

1.1 Motivation on Mathematical Modeling

Mathematical models and, in particular, differential equations play an important role in many branches of sciences. A growing field is a mathematical or theoretical biology where, for example, models for brain cancer [11] or models for competing species [12] are considered. Ethanol metabolism in the human body can be described by differential equations as well [13,14]. Temporal development of calcium-ion concentrations in liver cells can be modeled in a similar way as well [15]. In recent years, epidemiology has seen a large growth due to modeling approaches [16–18]. Further fields are chemistry [19–21] or computer science [22,23]. Mathematical modeling has also been an invaluable tool for understanding HIV-1 dynamics within human individuals. Let us first consider one basic mathematical model of untreated HIV dynamics within human hosts described by Stafford and co-workers [24], Perelson [25] or presented in Alizon’s and Magnus’s review article [8]. Here, we assume human blood as our compartment for our modeling idea as a non-linear system of ordinary differential equations. Since CD4+ T-cells are mainly affected by an HIV infection [4], we denote the density of susceptible, uninfected CD4+ T-cells by T(t), the density of infected CD4+ T-cells by I(t) and the density of virus particles by V(t). Uninfected target CD4+ T-cells are produced by a constant rate λ in the bone marrow, transfer to infected CD4+ T-cells by a constant rate β and die at a constant natural death rate dT. These infected CD4+ T-cells are also destroyed or die by a constant rate dI. Those infected cells can burst and expose new virus particles by a constant rate p. Finally, virus particles are cleared or die naturally from the human system by a constant rate dV. Hence, this fundamental model for the description of untreated HIV dynamics within human individuals reads

{T′(t)=λ−β⋅T(t)⋅V(t)−dT⋅T(t),I′(t)=β⋅T(t)⋅V(t)−dI⋅I(t),V′(t)=p⋅I(t)−dV⋅V(t),T(0)=T0>0,I(0)=I0≥0,V(0)=V0>0,(1)

where T0, I0 and V0 denote non-negative initial conditions. Here, we assume T0 and V0 to be positive.

From a mathematical modeling point of view, different adaptations of (1) exist in the literature [8–10]. Let us chronologically summarize some important studies in this regard. In 1991, Nowak and May mainly investigated numerically a simple mathematical model for antigenic variation with multiple virus strains in HIV infections [26]. Based on (1), Perelson, Kirschner and De Boer focused on stability and simulation of an extension of (1) by latently and actively infected CD4+ T-cells and a logistic depletion of susceptible, uninfected CD4+ T-cells [27] two years later. In 1994, Bonhoeffer and Nowak concentrated on the simulation and stability of intra-host time evolution of two HIV strains [28].

Before we continue our chronological summary of different contributions to the field of mathematical modeling of primary HIV infection, we want to introduce our model under consideration which is excellently described by Kirschner [29]. Again, we denote the density of susceptible, uninfected CD4+ T-cells by T(t), the density of infected CD4+ T-cells by I(t) and the density of virus particles by V(t). In contrast to (1), Kirschner modified (1) by a term which represents stimulation of CD4+ T-cells due to presence of virus particle. This adjusted model reads

{T′(t)=k1(t)−k2⋅T(t)+k3⋅T(t)⋅V(t)k4+V(t)−k5⋅T(t)⋅V(t),I′(t)=k5⋅T(t)⋅V(t)−k6⋅I(t)−k3⋅I(t)⋅V(t)k4+V(t),V′(t)=k7⋅k3⋅I(t)⋅V(t)k4+V(t)−k8⋅T(t)⋅V(t)+k9⋅V(t)k10+V(t),T(0)=T0>0,I(0)=I0≥0,V(0)=V0>0,(2)

where T0>0, I0≥0 and V0>0 stand for initial conditions. All meanings of time-dependent functions are presented in Table 1.

In the first equation of (2), k1(t) denotes a time-dependent production rate of susceptible CD4+ T-cells while k2 is a constant natural death rate of these cells. In particular, k1(t) should be a decreasing function if the virus load becomes larger which means that k1(t) should be bounded below by zero, bounded above by a positive constant K1 for all t≥0 and it converges to a non-negative k1~ as t→∞. In addition to that, k3 symbolizes a constant stimulation rate of the susceptible cells by the presence of virus particles. Here, k2>k3 is assumed. k5 is a constant transfer rate of susceptible cells to infected CD4+ T-cells by infiltration of susceptible CD4+ T-cells from virus particles. k4 is an adjustable parameter. In the second equation of (2), k6 symbolizes the clearance and death rate of infected CD4+ T-cells, while the last term −k3⋅I(t)⋅V(t)k4+V(t) represents a bursting term of infected CD4+ T-cells after their infiltration by virus particles. In the last equation of (2), k7 stands for the replication rate of virus particles after division of infected CD4+ T-cells while k8 represents a loss rate of virus particles by immune reactions. At last, k9 is the growth rate of virus particles, whereas k10 stands for a saturation rate of this proliferation process. All meanings of constant, positive problem parameters are summarized in Table 2.

Now, we can complete our brief overview on mathematical modeling of primary HIV infections. In 1996, Perelson and co-authors suggested an adjustment of (1) by an infectious pool and a non-infectious pool of virions after introducing a pharmacological therapy [30]. One year later, Perelson and co-workers proposed a changed dynamical model of HIV infection under combination therapy [31]. In 1999, Perelson and Nelson reviewed some fundamental mathematical models, especially regarding the stability of systems of ordinary differential equations concerning primary HIV-infections [32]. In 2000, Stafford and co-authors specifically investigated parameter estimation of model (1) concerning data from different patients [24]. Two years later, Perelson reviewed some basic models regarding primary HIV infection with and without medical treatment [25].

A delay-based system of integro-differential equations of cell-to-cell spread of primary HIV infections was introduced by Culshaw et al. in 2003 [33]. Based on (1), Stilianakis and Schenzle recommended a time-dependent infection rate [34]. In 2009, Kamp examined different strains under HIV coreceptor switch from a dynamical perspective, where investigation of stability was mainly focused [35]. One year later, Elaiw proposed an alternation of (1) by latently infected and actively infected CD4+ T-cells. The main focus of this research work was the proof of globally asymptotic stability [36]. Lythgoe and Fraser came up with a model of three host-cell types, namely susceptible CD4+ T-cells, latently infected memory CD4+ T-cells and macrophages [37]. In the same year of 2012, Alizon and Magnus reviewed the target-cell limitation model (1) and general schemes of the HIV within-host model [8]. Perelson and Ribeiro reconsidered (1) and modifications of within-host HIV models [9]. In 2014, Pourbashash and co-authors examined a global analysis of within-host virus models with cell-to-cell viral transmission [38]. Shen et al. analyzed the global stability of an infection-age structured HIV-1 model linking within-host and between-host dynamics in 2015 [39]. In the same year, Wang and co-authors developed a mathematical model for slow CD4+ T-cell decline in HIV-infected individuals [40]. Based on model (1), Pankavich and Parkinson analyzed a modification with Laplacian operators to include spatial heterogeneity in within-host viral dynamics models [41]. In 2017, Wang and co-authors presented an age-structured within-host HIV model with T-cell competition [42]. One year later, Almocera and co-authors established one multiscale within-host and between-host model for viral infectious diseases [43]. In 2023, D’Orso and Forst reviewed different models for HIV infections [10]. In 2024, Clarke and Pankavich proposed a three-stage model for HIV infections and its implications for antiretroviral therapy [44]. Wacker further investigated a simple mathematical model for primary HIV infection concerning classical mathematical questions [45] in 2024. In 2023, Li and co-authors investigated the impact of standard and non-standard finite-difference methods for between-host dynamics of HIV infections within a population [46], which motivates a quick overview of numerical discretizations of different mathematical models of HIV infection. In 2024, Meetei and co-authors investigated one HIV model against data from real cases [47]. In the same year, Ayele and co-authors examined a population-based model for the co-infection of HIV with tuberculosis [48]. In 2025, Hao and co-workers proposed a mathematical model in order to evaluate an international target strategy against HIV [49]. Recently, Teklu and co-workers have suggested a population-based and age-structured HIV model with optimal control [50].

We also like to briefly highlight stochastic approaches for the solution of within-host HIV models, such as those presented by Mbogo [51] and co-authors in 2013 or by Umar and co-workers [52] in 2022.

1.2 Motivation on Numerical Discretizations

While many works on mathematical models of HIV infections exist, numerical solution of these models, such as time-discrete variants of time-continuous models, is an important research topic as well. We want to summarize some studies with regard to this matter.

Many classical time integration methods such as Runge-Kutta schemes can be applied to such systems of ordinary differential equations [53,54]. As an alternative approach, Mickens popularized the concept of non-standard finite-difference methods for solving differential equations [55,56]. Let

y′(t)=f(t,y(t))

be a system of ordinary differential equations. Often, we want to find preferably a finite-difference-scheme of the form

Dh(yn)=Fh(f,yn)

where Dh(yn)≈y′(tn), Fh(f,yn)≈f(tn,yn) and yn≈y(tn) are appropriate approximations at time grid points tn for numerical simulations. A discretization scheme is called a non-standard finite-difference-method if at least one of the conditions

1.   Dh(yn)=yn+1−ynφ(hn) is an approximation of the first derivatives of our dynamical system where φ(h)=h+𝒪(h2) holds and hn=tn+1−tn is a possibly non-equidistant time-step size;

2.   Fh(f,yn)=g(yn+1,yn,…,h) is a non-local approximation of the right-hand-side vectorial function hold.

In 2013, Obaid et al. established a non-standard finite-difference method of development of HIV infections between hosts in a population [57]. In 2016, Yang and co-authors presented one non-standard finite-difference method for a diffusive within-host dynamics model with both virus-to-cell and cell-to-cell transmissions [58]. These authors mainly focused on analyzing stability properties for their suggested discretization method. Elaiw and Alshaikh examined the global stability of different time-discrete HIV dynamical models with three categories of infected CD4+ T-cells [59]. Salman introduced a non-standard finite-difference method and an optimal control problem for one HIV model with Beddington-DeAngelis incidence and cure rates [60]. Attaullah and Sohaib suggested one continuous Petrov-Galerkin method and one Legendre Wavelet Collocation method for the discretization of the within-host HIV model [61]. In 2023, Elaiw and co-authors developed a non-standard finite-difference-method for a co-infection model of HIV and another illness where the vectorial right-hand-side function of the dynamical system is a non-linear polynomial [62]. In that work, the authors established the global asymptotic stability of their four different equilibria by providing different Lyapunov functions. In 2025, Namjoo and co-workers proposed a non-standard finite-difference-method for a bounded, time-continuous HIV model of CD+ T-cells [63].

1.3 Contributions

All in all, many studies focused mainly on stability properties of (2) and only a few non-standard finite-difference methods exist in the literature. Henceforth, our contributions can be stated as follows:

1.   First, we extend some analytical results of (2) by proofs of non-negativity, boundedness of T(t) and I(t), at most linear growth of V(t), global existence in time, global uniqueness in time on a time interval [0,T] for an arbitrary T>0 and discuss further analytical and numerical investigations on stability properties of model (2). Here, we especially mention all results on the boundedness of solutions, global existence in time, and global uniqueness in time, which were not proved in Kirschner’s original work [29]. In particular, we would like to highlight that we, in contrast to Kirschner, establish local stability of equilibria analytically or numerically;

2.   As our main contribution, we suggest a novel non-standard finite-difference-method

{Tj+1−Tjhj=k1,j+1−k2⋅Tj+1+k3⋅Tj⋅Vjk4+Vj−k5⋅Tj+1⋅Vj,Ij+1−Ijhj=k5⋅Tj+1⋅Vj−k6⋅Ij+1−k3⋅Ij+1⋅Vjk4+Vj,Vj+1−Vjhj=k7⋅k3⋅Ij+1⋅Vjk4+Vj−k8⋅Tj+1⋅Vj+1+k9⋅Vjk10+Vj,T1=T0,I1=I0,V1=V0,(3)

with non-equidistant time-step sizes hj=tj+1−tj and, for example, Tj denotes time-discrete solutions at grid points tj on a time mesh and k1,j+1 is an appropriate time-discretization of k1(t). Additionally, we show that many properties of the time-continuous case such as non-negativity, boundedness or at most linear growth transfer to the time-discrete setting and we demonstrate linear convergence of the time-discrete solution components towards the time-continuous ones;

3.   Finally, we illustrate our theoretical findings by numerical examples.

Our work is organized as follows. After stating our main goals in Section 1, we provide analytical results of the time-continuous model (2) in Section 2. In Section 3, we examine our proposed non-standard finite-difference-method (3) thoroughly and show evidence of our theoretical results by numerical experiments in Section 4. Conclusively, we summarize our findings and possible future research directions in Section 5.

2  Analysis of the Time-Continuous Problem Formulation

Since the vectorial right-hand-side function of (2) is at least continuous, we can conclude that all solution components are at least continuously differentiable.

2.1 Preliminaries

Before providing our analytical results, we want to summarize some important results that we need for our investigation. Let us consider the general initial-value problem

{z′(t)=G(t,z(t)),z(0)=z0,(4)

where all state variables are included in z(t)=(z1(t),…,zn(t)) and all initial values are given by the initial vector z0∈Rn. We define a function

G:[0,∞)×(𝒞([0,∞)×Rn;R))n⟶(𝒞([0,∞)×Rn;R))n

where X:=(𝒞([0,∞)×Rn;R))n denotes the metric space of continuous functions on [0,∞) equipped with the classical supremum norm ‖⋅‖∞. Now, we want to summarize important definitions such as local and global Lipschitz-continuity which can be found in [[64], Subsection 3.2.1]. Let (X,‖⋅‖X) and (Y,‖⋅‖Y) be normed spaces. A function f:X⟶Y is called globally Lipschitz-continuous on X if there exists a constant K≥0 such that the inequality

‖f(x1)−f(x2)‖Y≤K⋅‖x1−x2‖X

holds for all x1,x2∈X. A function f:X⟶Y is called locally Lipschitz-continuous if, for every x∈X, there exists a neighborhood U of x such that f restricted to U globally Lipschitz-continuous on U. Let us recapitulate one important theorem from [[64], Theorem 4.7.1].

Theorem 1. If G:[0,∞)×Rn⟶Rn is locally Lipschitz-continuous in both time and state variables and if there exist non-negative real functions D:[0,∞)⟶[0,∞) and K:[0,∞)⟶[0,∞) such that

‖G(t,z(t))‖Rn≤K(t)⋅‖z(t)‖Rn+D(t)(5)

holds for all z(t)∈Rn where ‖⋅‖ denotes an appropriate norm on the corresponding space. Then it follows a solution of the initial-value problem (4) exists for all t≥0 and moreover, for all finite T≥0, it holds

‖z(t)‖Rn≤‖z0‖Rn⋅exp(Kmax⋅t)+DmaxKmax⋅(exp(Kmax⋅t)−1)(6)

for all t∈[0,T] where

Dmax:=max0≤s≤T|D(s)|andKmax:=max0≤s≤T|K(s)|.

2.2 Non-Negativity, Boundedness of T(t) and I(t) and Linear Growth of V(t) in Time

First, we want to show that possible solutions of (2) remain non-negative for all t≥0.

Theorem 2. Solutions of (2) stay non-negative for all t≥0.

Proof. Since all solution components of (2) are continuous, there exists a time point t⋆>0 such that

t⋆=supt>0{t>0:T(t)≥0∧I(t)≥0∧V(t)≥0}

holds. Let us assume that t⋆>0 is finite. We can conclude that

T′(t)=k1(t)−k2⋅T(t)+k3⋅T(t)⋅V(t)k4+V(t)−k5⋅T(t)⋅V(t)≥−k2⋅T(t)−k5⋅T(t)⋅V(t)=−T(t)⋅(k2+k5⋅V(t))

is valid for all t∈[0,t⋆] and by the comparison principle, we obtain

T(t)≥T0⋅exp(∫0t(−k2−k5⋅V(τ))dτ).

Additionally, we notice that

I′(t)≥−k6⋅I(t)−k3k4⋅I(t)=−(k6+k3k4)⋅I(t)

holds for all t∈[0,t⋆] and this implies

I(t)≥I0⋅exp(∫0t−(k6+k3k4)dτ)=I0⋅exp(−(k6+k3k4)⋅t)

for all t∈[0,t⋆]. Furthermore, we see that

V′(t)≥−k8⋅T(t)⋅V(t)

holds for all t∈[0,t⋆] and the comparison principle yields

V(t)≥V0⋅exp(−∫0tk8⋅T(τ)dτ)

for all t∈[0,t⋆]. All these inequalities contradict the maximality of t⋆ and consequently, all solution components of (2) remain non-negative for all t≥0. □

Next, we want to prove the boundedness of T(t) globally in time.

Lemma 1. The solution components T(t) and I(t) are bounded for all t≥0.

Proof. We separate this proof into two parts.

1.   It has already been proven that the inequality T(t)≥0 is valid for all t≥0. For that reason, we only need to demonstrate that there exists a constant A>0 such that T(t)≤A for all t≥0. Since k1(t)≤K1 for all t≥0 and k2>k3 hold, we obtain

T′(t)=k1(t)−k2⋅T(t)+k3⋅T(t)⋅V(t)k4+V(t)−k5⋅T(t)⋅V(t)≤K1−k2⋅T(t)+k3⋅T(t)≤K1−(k2−k3)⋅T(t)

for all t≥0 due to the non-negativity of all solution components. By the comparison principle, we get

T(t)≤K1k2−k3+(T0−K1k2−k3)⋅exp(−(k2−k3)⋅t)≤max{T0;K1k2−k3}=:Tmax

which proves our first assertion.

2.   It has also been shown that the inequality I(t)≥0 holds for all t≥0. Consequently, we solely need to find a constant B>0 such that I(t)≤B for all t≥0. First, by adding T′(t) and I′(t) and applying k1(t)≤K1 and k2>k3 for all t≥0, we conclude that

N′(t):=T′(t)+I′(t)={k1(t)−k2⋅T(t)+k3⋅T(t)⋅V(t)k4+V(t)−k5⋅T(t)⋅V(t)}+{k5⋅T(t)⋅V(t)−k6⋅I(t)−k3⋅I(t)⋅V(t)k4+V(t)}≤K1−(k2−k3)⋅T(t)−k6⋅I(t)≤K1−min{k2−k3;k6}⋅(T(t)+I(t))

holds because of the non-negativity of all solution components. Again, applying the comparison principle, we obtain

T(t)+I(t)≤K1min{k2−k3;k6}+(T0+I0−K1min{k2−k3;k6})⋅exp(−min{k2−k3;k6}⋅t)≤max{T0+I0;K1min{k2−k3;k6}}=:Nmax

for all t≥0 which shows our second assertion due to the non-negativity of both solution components T(t) and I(t).

This finishes our proof. □

As our last result of this subsection, we prove that V(t) grows at most linearly in time.

Lemma 2. The solution component V(t) grows at most linearly in time.

Proof. The differential equation for V(t) yields the estimate

V′(t)=k7⋅k3⋅I(t)⋅V(t)k4+V(t)−k8⋅T(t)⋅V(t)+k9⋅V(t)k10+V(t)≤k3⋅k7⋅Nmax+k9

due to the non-negativity of all solution components. Consequently, it holds

V(t)≤V0+(k3⋅k7⋅Nmax+k9)⋅t

by the comparison principle for all t≥0 which proves our claim. □

2.3 Existence Globally in Time

Let us define

F(t,(T(t),I(t),V(t))):=(F1(t,(T(t),I(t),V(t)))F2(t,(T(t),I(t),V(t)))F3(t,(T(t),I(t),V(t))))=(k1(t)−k2⋅T(t)+k3⋅T(t)⋅V(t)k4+V(t)−k5⋅T(t)⋅V(t)k5⋅T(t)⋅V(t)−k6⋅I(t)−k3⋅I(t)⋅V(t)k4+V(t)k7⋅k3⋅I(t)⋅V(t)k4+V(t)−k8⋅T(t)⋅V(t)+k9⋅V(t)k10+V(t))

as the right-hand-side function of (2). We notice that this function is locally Lipschitz-continuous on [0,∞)×(𝒞([0,∞);R))3 because it is continuously differentiable on this set concerning t, T, I and V by [[64], Proposition 3.2.2]. Consequently, we can state the following result.

Theorem 3. Solutions of (2) exist for all t≥0.

Proof. We need to estimate all components of F(t,(T(t),I(t),V(t))) separately.

1. We see that

‖F1(t,(T(t),I(t),V(t)))‖∞=‖k1(t)−k2⋅T(t)+k3⋅T(t)⋅V(t)k4+V(t)−k5⋅T(t)⋅V(t)‖∞≤K1+k2⋅Tmax+k3⋅Tmax+k5⋅Tmax⋅‖V(t)‖∞≤(K1+(k2+k3)⋅Tmax)+k5⋅Tmax⋅‖(T(t),I(t),V(t))‖∞

holds.

2. Additionally, we notice that

‖F2(t,(T(t),I(t),V(t)))‖∞=‖k5⋅T(t)⋅V(t)−k6⋅I(t)−k3⋅I(t)⋅V(t)k4+V(t)‖∞≤k5⋅Tmax⋅‖V(t)‖∞+k6⋅Nmax+k3⋅Nmax≤(k3+k6)⋅Nmax+k5⋅Tmax⋅‖(T(t),I(t),V(t))‖∞

is valid.

3. Furthermore, we obtain

‖F3(t,(T(t),I(t),V(t)))‖∞=‖k7⋅k3⋅I(t)⋅V(t)k4+V(t)−k8⋅T(t)⋅V(t)+k9⋅V(t)k10+V(t)‖∞≤k3⋅k7⋅Nmax+k8⋅Tmax⋅‖V(t)‖∞+k9≤(k3⋅k7⋅Nmax+k9)+k8⋅Tmax⋅‖(T(t),I(t),V(t))‖∞.

4. Conclusively, these three inequalities yield

‖F(t,(T(t),I(t),V(t)))‖∞≤(K1+k9+(k2+2⋅k3+k6+k3⋅k7)⋅Nmax)+(2⋅k5+k8)⋅Tmax⋅‖(T(t),I(t),V(t))‖∞.

Finally, this proves the existence of solutions of (2) globally in time since all assumptions of Theorem 1 are fulfilled. □

2.4 Uniqueness in Time

By applying [[65], Corollary 2, Subchapter 2.4] and [[65], Theorem 4, Subchapter 2.4], we obtain a uniqueness result globally in time on [0,T] for an arbitrary T>0.

Theorem 4. Solutions of (2) exist uniquely for all t∈[0,T].

Proof. Let us construct the open set E:=(−12⋅min{k4;k10},∞)3 for the solution variables (T,I,V). Our initial conditions (T0,I0,V0)∈[0,∞)3⊂E are assumed to be non-negative. Furthermore, we have already observed that our right-hand-side function is continuously differentiable on E. By Lemmas 1 and 2, we see that

(T(t),I(t),V(t))⊂K:=[0,Nmax]×[0,Nmax]×[0,V0+(k3⋅k7⋅Nmax+k9)⋅T]⊂E

holds for all t∈[0,T] for an arbitrary T>0. If we assume two different solution vectors starting at the same initial condition, we conclude from [[65], Theorem 4, Subchapter 2.4] that all solution components exist uniquely for all t∈[0,T] for an arbitrary T>0. As a consequence, uniqueness on [0,T] is established. □

2.5 Equilibrium States

For searching equilibrium states of system (2), we need to solve the system of non-linear equations

0=k1~−k2⋅T⋆+k3⋅T⋆⋅V⋆k4+V⋆−k5⋅T⋆⋅V⋆,0=k5⋅T⋆⋅V⋆−k6⋅I⋆−k3⋅I⋆⋅V⋆k4+V⋆,0=k3⋅k7⋅I⋆⋅V⋆k4+V⋆−k8⋅T⋆⋅V⋆+k9⋅V⋆k10+V⋆,(7)

where (T⋆,I⋆,V⋆) denotes a possible equilibrium state. Here, we assume that k1(t)=k1~ is a constant for all t≥0 for simplicity.

2.5.1 Virus-Free Equilibrium State

It can be seen from (7) that the virus-free equilibrium state reads

(T1⋆,I1⋆,V1⋆)=(k1~k2,0,0).(8)

Let us define the basic reproduction number R0 by

R0:=k2⋅k9k1~⋅k8⋅k10(9)

which is going to be motivated in the proof of our next statement.

Theorem 5. If R0<1 of (9) holds, our virus-free equilibrium state (8) of (2) is locally asymptotically stable.

Proof. We divide our proof into different parts.

1. We first calculate the Jacobian matrix, denoted by DF(T⋆,I⋆,V⋆), of the right-hand-side function of our dynamical system at possible equilibria. We notice that

DF(T⋆,I⋆,V⋆)=(AF,1,1AF,1,2AF,1,3AF,2,1AF,2,2AF,2,3AF,3,1AF,3,2AF,3,3)(10)

with entries

AF,1,1:=−k2+k3⋅V⋆k4+V⋆−k5⋅V⋆,AF,1,2:=0,AF,1,3:=k3⋅T⋆k4+V⋆−k3⋅T⋆⋅V⋆(k4+V⋆)2−k5⋅T⋆,AF,2,1:=k5⋅V⋆,AF,2,2:=−k6−k3⋅V⋆k4+V⋆,AF,2,3:=k5⋅T⋆−k3⋅V⋆k4+V⋆+k3⋅I⋆⋅V⋆(k4+V⋆)2,AF,3,1:=−k8⋅V⋆,AF,3,2:=k3⋅k7⋅V⋆k4+V⋆,AF,3,3:=k3⋅k7⋅I⋆k4+V⋆−k3⋅k7⋅I⋆⋅V⋆(k4+V⋆)2−k8⋅T⋆+k9⋅1k10+V⋆−k9⋅V⋆(k10+V⋆)2,

holds.

2. If we plug our virus-free equilibrium state (8), we obtain

DF(T1⋆,0,0)=(−k20k3k4⋅T1⋆−k5⋅T1⋆0−k6k5⋅T1⋆00−k8⋅T1⋆+k9k10).

Calculating the characteristic polynomial χ(λ) of this Jacobian at the virus-free equilibrium state, we see that

χ(λ)=det(DF(T1⋆,0,0)−λ⋅(100010001))=(−k2−λ)⋅(−k6−λ)⋅(−k8⋅T1⋆+k9k10−λ)=0

is valid. Both eigenvalues λ1=−k2<0 and λ2=−k6<0 are negative. The last eigenvalue λ3=k9k10−k8⋅T⋆=k9k10−k1~⋅k8k2 is negative if R0:=k2⋅k9k1~⋅k8⋅k10<1 holds.

This finishes our argument. □

2.5.2 Virus-Endemic Equilibrium State

Here, we assume that

R0:=k2⋅k9k1~⋅k8⋅k10>1

holds. From (7), we obtain

T2⋆=k1~k2+k5⋅V2⋆−k3⋅V2⋆k4+V2⋆≥0,I2⋆=k5⋅T2⋆⋅V2⋆k3⋅V2⋆k4+V2⋆+k6≥0,V2⋆⋅(k3⋅k7⋅I2⋆k4+V2⋆−k8⋅T2⋆+k9k10+V2⋆)=0.(11)

Since V2⋆>0 should hold in a virus-endemic equilibrium state (T2⋆,I2⋆,V2⋆), it must hold

k3⋅k7⋅I2⋆k4+V2⋆−k8⋅T2⋆+k9k10+V2⋆=0.

Plugging all results of (11) into our previous equation, we get

0=k1~⋅k3⋅k5⋅k7⋅V2⋆(k4+V2⋆)⋅(k3⋅V2⋆k4+V2⋆+k6)⋅(k2+k5⋅V2⋆−k3⋅V2⋆k4+V2⋆)−k1~⋅k8k2+k5⋅V2⋆−k3⋅V2⋆k4+V2⋆+k9k10+V2⋆(12)

and we can re-establish zeros of (12) by finding zeros of

0=k1~⋅k3⋅k5⋅k7⋅V2⋆(k3+k6)⋅V2⋆+k4⋅k6−k1~⋅k8+k9⋅(k5⋅V2⋆−k3⋅V2⋆V2⋆+k4+k2)V2⋆+k10=:G(V2⋆).(13)

Let us define the coefficients

a3:=k1~⋅k3⋅k5⋅k7−k1~⋅k3⋅k8−k1~⋅k6⋅k8+k3⋅k5⋅k9+k5⋅k6⋅k9,a2:=k1~⋅k3⋅k4⋅k5⋅k7−k1~⋅k3⋅k4⋅k8−2⋅k1~⋅k4⋅k6⋅k8−k32⋅k9+k2⋅k3⋅k9+k3⋅k4⋅k5⋅k9+k2⋅k6⋅k9−k3⋅k6⋅k9+2⋅k4⋅k5⋅k6⋅k9+k1~⋅k3⋅k5⋅k7⋅k10−k1~⋅k3⋅k8⋅k10−k1~⋅k6⋅k8⋅k10,a1:=−k1~⋅k42⋅k6⋅k8+k2⋅k3⋅k4⋅k9+2⋅k2⋅k4⋅k6⋅k9−k3⋅k4⋅k6⋅k9+k42⋅k5⋅k6⋅k9+k1~⋅k3⋅k4⋅k5⋅k7⋅k10−k1~⋅k3⋅k4⋅k8⋅k10−2⋅k1~⋅k4⋅k6⋅k8⋅k10,a0:=k2⋅k42⋅k6⋅k9−k1~⋅k42⋅k6⋅k8⋅k10=k1~⋅k4⋅k6⋅k8⋅k10⋅(R0−1).(14)