Next Article in Journal
Light “You Only Look Once”: An Improved Lightweight Vehicle-Detection Model for Intelligent Vehicles under Dark Conditions
Next Article in Special Issue
A Computational Approach to Individual Cell-Based Decision Algorithms Involved in Bone Remodeling
Previous Article in Journal
A Task Orchestration Strategy in a Cloud-Edge Environment Based on Intuitionistic Fuzzy Sets
Previous Article in Special Issue
Effect of Impaired B-Cell and CTL Functions on HIV-1 Dynamics
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Periodic Behaviour of HIV Dynamics with Three Infection Routes

by
Miled El Hajji
* and
Rahmah Mohammed Alnjrani
Department of Mathematics and Statistics, Faculty of Science, University of Jeddah, P.O. Box 80327, Jeddah 21589, Saudi Arabia
*
Author to whom correspondence should be addressed.
Mathematics 2024, 12(1), 123; https://doi.org/10.3390/math12010123
Submission received: 30 November 2023 / Revised: 19 December 2023 / Accepted: 27 December 2023 / Published: 29 December 2023

Abstract

:
In this study, we consider a system of nonlinear differential equations modeling the human immunodeficiency virus type-1 (HIV-1) in a variable environment. Infected cells were subdivided into two compartments describing both latently and productively infected cells. Thus, three routes of infection were considered including the HIV-to-cell contact, latently infected cell-to-cell contact, and actively infected cell-to-cell contact. The nonnegativity and boundedness of the trajectories of the dynamics were proved. The basic reproduction number was determined through an integral operator. The global stability of steady states is then analyzed using the Lyapunov theory together with LaSalle’s invariance principle for the case of a fixed environment. Similarly, for the case of a variable environment, we showed that the virus-free periodic solution is globally asymptotically stable once R 0 1 , while the virus will persist once R 0 > 1 . Finally, some numerical examples are provided illustrating the theoretical investigations.

1. Introduction

Human immunodeficiency virus (HIV) infection is caused by two retroviruses (HIV-1 and HIV-2) that can destroy CD4 + lymphocytes and disrupt the functioning of cell-mediated immunity, leading to an increased risk of some infections and certain cancers. The initial infection may result in a nonspecific febrile illness. The risk of subsequent manifestations, linked to immunodeficiency, is related to the rate of CD4 + lymphocyte exhaustion. HIV destroys the brain, gonads, kidneys, and heart, provoking cognitive impairment, hypogonadism, kidney failure, and cardiomyopathy. Manifestations can range from asymptomatic forms to acquired immunodeficiency syndrome (AIDS).
HIV infection has attracted the attention of several researchers from various fields. In particular, mathematical modeling permits mathematicians to use a set of concepts, methods, amd mathematical theories that facilitate the description, understanding, and prediction of the evolution of a disease, making a link between reality and mathematics.
Mathematical models of HIV-immune dynamics are widely proposed in the literature [1,2,3,4,5,6]. The interest of a range of works has focused on the description of the immunological response to HIV infection, characterizing interactions of HIV with CD 4 + expressing non-conventional T cells. Most of these references proposed mathematical models describing the host–pathogen dynamics for the principal virus strain, HIV-1, based on biological and clinical observations. All of the proposed models and studies are based on systems of ordinary-differential equations. In particular, the mathematical modeling of the behavior of HIV dynamics was studied in several recent works [7,8,9]. In [10], an HIV mutant infection model was proposed, and the feasibility and effectiveness of the optimal therapies were studied for this model in which virtual therapies for wild-type infections are incorporated. Later, in [11], the possibilities of eradication at an early stage of infection were studied on the model proposed in [10]. In [12], the authors studied an HIV infection model with fixed parameters taking into account three routes of infection, namely virus-to-cell, latently infected cell-to-cell, and actively infected cell-to-cell.
A person is more likely to have a heart attack in winter or catch a sexually transmitted infection in summer. Indeed, the frequency of many diseases varies with the seasons. Several published studies explain that the prevalence of several infectious diseases is linked to the season [13]. To establish a timeline, researchers analyzed epidemiological data from the Centers for Disease Control and Prevention (the main federal agency in the United States for protecting public health), the World Health Organization (WHO), and the European Centre for Disease Prevention and Control (ECDC). Generally, there are, on the one hand, long-standing epidemics with low morbidity rates and, on the other hand, seasonal epidemics with high morbidity rates. Pathogens in seasonal circulation are of more interest to the population since everyone, or almost everyone, has been or will be infected one day or another. For example, seasonal flu occurs during the cold season. In spring and summer, epidemics of borreliosis (Lyme disease) and tick-borne encephalitis (summer meningoencephalitis) occur. The number of campylobacter food poisoning cases increases in winter during holidays as well as in summer during grilling time. Measles cases also follow a seasonal pattern; in this case, the elimination strategy helped reduce the chains of transmission. Thus, the season plays a major role in the occurrence of certain infectious diseases, which is evoked in several studies [14,15,16,17]. Mechanisms linked to seasonality in the occurrence of infectious diseases should be considered by the centers for disease control and prevention. This would strengthen public health policies aimed at controlling them. Several studies considered mathematical models that describe the influence of seasonality on the dynamics of several diseases. In [18], the seasonal dynamics of an SVEIR model was considered and analyzed. [19,20] consider the seasonal dynamics of HIV and chikungunya virus.
In our case, the dynamics of HIV in relation with cells will be studied for both fixed and variable environments. We provide a thorough exploration of HIV-1 dynamics through a nonlinear differential equation model. The subdivision of infected cells into two compartments, the consideration of three infection routes, and the variable environment add realism to the proposed model.
Thus, the present paper is organized as follows: in Section 2, we consider a simple mathematical model for HIV spread by considering both the variable environment and three routes of infection. In Section 3, we consider firstly the case of the fixed environment, and we express R 0 and investigate the global analysis of both the virus-free and the virus-infected steady states. Stability analyses were conducted using Lyapunov theory and LaSalle’s invariance principle. In Section 4, we investigate the stability of virus-free and virus-present periodic trajectories for the case of the variable environment. Some numerical tests are provided in Section 5, validating the given theoretical results. Finally, in Section 6, some conclusions are given.

2. Mathematical Model for HIV Dynamics

We aim in this paper to propose a dynamical system generalizing the one given in [12] to a variable environment. Let Y s , Y l , and Y i be the concentration of susceptible, latently infected, and productively infected cells, respectively. Let Y v and Y c be the free virions (HIV-1) and B cells and cytotoxic T lymphocytes (CTLs), respectively. The proposed mathematical model is an association of five ordinary differential equations and takes the following form.
Y s ˙ ( t ) = d s ( t ) Λ i n ( t ) d s ( t ) Y s ( t ) [ β 1 ( t ) Y v ( t ) + β 2 ( t ) Y l ( t ) + β 3 ( t ) Y i ( t ) ] Y s ( t ) , Y l ˙ ( t ) = [ β 1 ( t ) Y v ( t ) + β 2 ( t ) Y l ( t ) + β 3 ( t ) Y i ( t ) ] Y s ( t ) ( η ( t ) + d l ( t ) ) Y l ( t ) , Y i ˙ ( t ) = η ( t ) Y l ( t ) d i ( t ) Y i ( t ) γ ( t ) Y c ( t ) Y i ( t ) , Y v ˙ ( t ) = ε ( t ) Y i ( t ) d v ( t ) Y v ( t ) , Y c ˙ ( t ) = λ ( t ) Y i ( t ) d c ( t ) Y c ( t ) δ ( t ) Y c ( t ) Y i ( t ) ,
with positive initial condition ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) R + 5 .
The model’s variables and parameters are further described hereafter in Table 1.
Note that the incidence rate Y s increases once the number of susceptible cells, Y s , increases and it is zero in the absence of the susceptible cells ( Y s = 0 ). Thus, we impose some assumptions on the incidence rate f in the following assumption. We suppose also that the death rates of cells depend on the compartment to which the cell belongs. Thus, we assume that the death rate of susceptible cells is smaller than the death rate of latently infected cells, and that the death rate of latently infected cells is itself smaller than the death rate of productively infected cells. Therefore, we assume that d s ( t ) d l ( t ) d i ( t ) , t R + .

3. Autonomous System

In this section, we suppose that the parameters are not influenced by the environment and thus are all constants. The dynamics of the model (1) takes the following more simple form.
Y s ˙ = d s Λ i n d s Y s [ β 1 Y v + β 2 Y l + β 3 Y i ] Y s Y l ˙ = [ β 1 Y v + β 2 Y l + β 3 Y i ] Y s ( η + d l ) Y l , Y i ˙ = η Y l d i Y i γ Y c Y i , Y v ˙ = ε Y i d v Y v , Y c ˙ = λ Y i d c Y c δ Y c Y i ,
with positive initial condition ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) R + 5 . In this section, we make the following assumption.
Assumption 1.
β 2 Λ i n η + d l .

3.1. Basic Results

We start by proving that the state variables Y s ( t ) , Y l ( t ) , Y i ( t ) , Y v ( t ) , and Y c ( t ) remain non-negative all the time. Let us denote d = min ( d v , d c ) , and then we obtain the following lemma.
Lemma 1.
The bounded positive set
Σ = ( Y s , Y l , Y i , Y v , Y c ) R + 5 ; Y s + Y l + Y i Λ i n , Y v + Y c ε + λ d Λ i n
is invariant and an attractor of all solutions of (2).
Proof. 
R + 5 is invariant for (2) because Y s ˙ Y s = 0 = d s Λ i n > 0 , Y l ˙ Y l = 0 = ( β 1 Y v + β 3 Y i ) Y s 0 , Y i ˙ Y i = 0 = η Y l 0 , Y v ˙ Y v = 0 = ε Y i 0 , and Y c ˙ Y c = 0 = λ Y i 0 .
Consider the two variables S 1 ( t ) = Y s ( t ) + Y l ( t ) + Y i ( t ) Λ i n and S 2 ( t ) = Y v ( t ) + Y c ( t ) Λ i n d ( ε + λ ) . From Equation (2), we have S ˙ 1 ( t ) = d s Λ i n ( d s Y s ( t ) + d l Y l ( t ) + d i Y i ( t ) ) γ Y c ( t ) Y i ( t ) d s S 1 ( t ) since d s d l d i . Therefore, S 1 ( t ) S 1 ( 0 ) e d s t , and then Y s ( t ) + Y l ( t ) + Y i ( t ) Λ i n + ( Y s ( 0 ) + Y l ( 0 ) + Y i ( 0 ) Λ i n ) e d s t . Similarly, we have
S ˙ 2 ( t ) = Y ˙ v ( t ) + Y ˙ c ( t ) = ε Y i ( t ) d v Y v ( t ) + λ Y i ( t ) d c Y c ( t ) δ Y c ( t ) Y i ( t ) ( ε + λ ) Y i ( t ) ( d v Y v ( t ) + d c Y c ( t ) ) ( ε + λ ) Λ i n d ( Y v ( t ) + Y c ( t ) ) d S 2 ( t ) .
Therefore, S 2 ( t ) S 2 ( 0 ) e d t , and thus
Y v ( t ) + Y c ( t ) ε + λ d Λ i n + Y v ( 0 ) + Y c ( 0 ) ε + λ d Λ i n e d t .
One can see that if Y s ( 0 ) + Y l ( 0 ) + Y i ( 0 ) Λ i n , then Y s ( t ) + Y l ( t ) + Y i ( t ) Λ i n , and if Y v ( 0 ) + Y c ( 0 ) ε + λ d Λ i n , then Y v ( t ) + Y c ( t ) Λ i n d ( ε + λ ) , which confirms that Σ is a positively invariant attractor for system (2). □

3.2. Basic Reproduction Number and Equilibria of the Dynamics

We denote by R 0 , the basic reproduction number, that we used the next-generation matrix method [21,22] to calculate it. Let F = ( F i j ) 1 i , j 3 and V = ( V i j ) 1 i , j 3 with F 11 = β 2 Λ i n , F 12 = β 3 Λ i n , F 13 = β 1 Λ i n , F 21 = F 22 = F 23 = F 31 = F 32 = F 33 = 0 , V 11 = η + d l , V 21 = η , V 22 = d i , V 32 = ε , V 33 = d v , and V 12 = V 13 = V 23 = V 31 = 0 . Then, the next-generation matrix can be calculated as follows: F V 1 = ( G i j ) 1 i , j 3 with G 11 = ( η ε β 1 + d i d v β 2 + η d v β 3 ) Λ i n d i d v ( η + d l ) , G 12 = ( ε β 1 + d v β 3 ) Λ i n d i d v , G 13 = β 1 Λ i n d v , and G 21 = G 22 = G 23 = G 31 = G 32 = G 33 = 0 . Therefore, R 0 is given by
R 0 = ( η ε β 1 + d i d v β 2 + η d v β 3 ) Λ i n d i d v ( η + d l ) = η ε β 1 Λ i n d i d v ( η + d l ) + β 2 Λ i n ( η + d l ) + η β 3 Λ i n d i ( η + d l ) = R 1 + R 2 + R 3
where R 1 = η ε β 1 Λ i n d i d v ( η + d l ) , R 2 = β 2 Λ i n ( η + d l ) , and R 3 = η β 3 Λ i n d i ( η + d l ) .
Lemma 2.
  • If R 0 1 , then the model (2) admits a unique steady state represented by E 0 = ( Λ i n , 0 , 0 , 0 , 0 ) .
  • If R 0 > 1 , then the dynamics (2) admits two equilibria E 0 and E = ( Y s , Y l , Y i , Y v , Y c ) .
Proof. 
Let ( Y s , Y l , Y i , Y v , Y c ) be any steady state of dynamics (2) satisfying:
0 = d s Λ i n d s Y s [ β 1 Y v + β 2 Y l + β 3 Y i ] Y s , 0 = [ β 1 Y v + β 2 Y l + β 3 Y i ] Y s ( η + d l ) Y l , 0 = η Y l d i Y i γ Y c Y i , 0 = ε Y i d v Y v , 0 = λ Y i d c Y c δ Y c Y i .
By solving Equation (4), we obtain Y s = Λ i n ( η + d l ) ( d i d c Y i + ( d i δ + γ λ ) Y i 2 ) d s η ( d c + δ Y i ) , Y l = d i d c Y i + ( d i δ + γ λ ) Y i 2 η ( d c + δ Y i ) , Y v = ε Y i d v , and Y c = λ Y i ( d c + δ Y i ) .
From the second equation of (4), we have
( β 1 ε Y i d v + β 2 d i d c Y i + ( d i δ + γ λ ) Y i 2 η ( d c + δ Y i ) + β 3 Y i ) ( Λ i n ( η + d l ) ( d i d c Y i + ( d i δ + γ λ ) Y i 2 ) d s η ( d c + δ Y i ) ) = ( η + d l ) d i d c Y i + ( d i δ + γ λ ) Y i 2 η ( d c + δ Y i ) .
Therefore, if Y i = 0 then Y s = Λ i n and Y l = Y i = Y v = Y c = 0 , we obtain the HIV-free equilibrium point E 0 = ( Λ i n , 0 , 0 , 0 , 0 ) . However, if Y i 0 , we divide the previous equation by Y i , and we obtain
0 = ( β 1 ε d v + β 2 d i d c + ( d i δ + γ λ ) Y i η ( d c + δ Y i ) + β 3 ) ( Λ i n ( η + d l ) ( d i d c Y i + ( d i δ + γ λ ) Y i 2 ) d s η ( d c + δ Y i ) ) ( η + d l ) d i d c + ( d i δ + γ λ ) Y i η ( d c + δ Y i ) .
Consider the function g defined by
g ( Y i ) = ( β 1 ε d v + β 2 d i d c + ( d i δ + γ λ ) Y i η ( d c + δ Y i ) + β 3 ) ( Λ i n ( η + d l ) ( d i d c Y i + ( d i δ + γ λ ) Y i 2 ) d s η ( d c + δ Y i ) ) ( η + d l ) d i d c + ( d i δ + γ λ ) Y i η ( d c + δ Y i ) .
One can easily obtain
g ( 0 ) = β 1 ε d v + β 2 d i η + β 3 Λ i n ( η + d l ) d i η = ( η + d l ) d i η ( R 0 1 )
Let Y i 0 be the solution of the equation:
( η + d l ) ( d i d c Y i + ( d i δ + γ λ ) Y i 2 ) d s η ( d c + δ Y i ) = Λ i n ,
which can be written in the following form: d s η ( d c + δ Y i ) Λ i n = ( η + d l ) ( d i d c Y i + ( d i δ + γ λ ) Y i 2 ) . We then obtain:
a Y i 2 + b Y i + c = 0
where a = ( η + d l ) ( d i δ + γ λ ) , b = ( η + d l ) d i d c η δ d s Λ i n , and c = η d c d s Λ i n . Then, the discriminant is given by: Δ = b 2 4 a c = ( η + d l ) d i d c + η δ d s Λ i n 2 + 4 η γ λ d s d c ( η + d l ) Λ i n > 0 . Therefore, Equation (5) admits two real solutions. One of them is nonpositive and then excluded and we keep only the positive solution given by
Y i 0 = ( η + d l ) d i d c + η δ d s Λ i n + ( η + d l ) d i d c + η δ d s Λ i n 2 + 4 η γ λ d s d c ( η + d l ) Λ i n 2 ( η + d l ) ( d i δ + γ λ ) > 0 .
One can see easily that g ( Y i 0 ) = ( η + d l ) d i d c + ( d i δ + γ λ ) Y i 0 η ( d c + δ Y i 0 ) < 0 . The derivative of g is calculated as follows:
g ( Y i ) = β 1 ε d v + β 2 d i d c + ( d i δ + γ λ ) Y i η ( d c + δ Y i ) + β 3 × ( η + d l ) ( δ 2 d i Y i 2 + λ γ δ Y i 2 + 2 δ d i d c Y i + 2 λ γ d c Y i + d i d c 2 ) η d s ( d c + δ Y i ) 2 + β 2 d c γ λ η ( d c + δ Y i ) 2 Λ i n ( η + d l ) ( d i d c Y i + ( d i δ + γ λ ) Y i 2 ) d s η ( d c + δ Y i ) ( η + d l ) β 2 .
Therefore, it is clear that the first term of g ( Y i ) is negative. According to Assumption 1, we have Λ i n ( η + d l ) β 2 < 0 , and then the second term of g ( Y i ) is negative. Thus, g ( Y i ) < 0 which permits us to conclude the existence and uniqueness of Y i ( 0 , Y i 0 ) satisfying g ( Y i ) = 0 . Then, we obtain
Y s = Λ i n ( η + d l ) ( d i d c Y i + ( d i δ + γ λ ) ( Y i ) 2 ) η d s ( d c + δ Y i ) , Y l = d i d c Y i + ( d i δ + γ λ ) ( Y i ) 2 η ( d c + δ Y i ) , Y v = ε Y i d v , Y c = λ Y i ( d c + δ Y i ) .
Then, the persistence equilibrium point E = ( Y s , Y l , Y i , Y v , Y c ) exists and is unique once R 0 > 1 . □
Next, we investigate the local stability of the virus-free steady state E 0 using the Jacobian matrix. Note that, generally, the value of R 0 with respect to the unity will be decisive in whether the virus persists or goes to extinction. We give the stability result concerning the virus-free steady state E 0 .
Theorem 1.
If R 0 < 1 , then the virus-free steady state E 0 is locally asymptotically stable.
Proof. 
The Jacobian matrix of (2) at the virus-free steady state E 0 is:
J 0 = d s β 2 Λ i n β 3 Λ i n β 1 Λ i n 0 0 β 2 Λ i n ( η + d l ) β 3 Λ i n β 1 Λ i n 0 0 η d i 0 0 0 0 ε d v 0 0 0 λ 0 d c
J 0 admits λ 1 = d s < 0 and λ 2 = d c < 0 as eigenvalues. The other three eigenvalues are the roots of:
P 0 ( λ ) = β 2 Λ i n ( η + d l ) λ ) β 3 Λ i n β 1 Λ i n η ( d i + λ ) 0 0 ε ( d v + λ ) = ( β 2 Λ i n ( η + d l ) λ ) ( d i + λ ) 0 ε ( d v + λ ) η β 3 Λ i n β 1 Λ i n ε ( d v + λ ) = ( β 2 Λ i n ( η + d l ) λ ) ( d i + λ ) ( d v + λ ) η ( d v + λ ) β 3 Λ i n β 1 Λ i n ε = ( β 2 Λ i n ( η + d l ) λ ) ( d i d v + d i λ + d v λ + λ 2 ) + η ( d v + λ ) β 3 Λ i n + η ε β 1 Λ i n = ( λ 3 + a 2 λ 2 + a 1 λ + a 0 )
where
a 2 = d i + d v + ( η + d l ) β 2 Λ i n = d i + d v + ( η + d l ) ( 1 R 2 ) > 0 , for   R 0 < 1 . a 1 = d i d v + [ ( η + d l ) β 2 Λ i n ] ( d i + d v ) η β 3 Λ i n = d i d v + ( η + d l ) d v ( 1 R 2 ) + d i ( η + d l ) ( 1 ( R 2 + R 3 ) ) > 0 , for   R 0 < 1 . a 0 = [ ( η + d l ) β 2 Λ i n ] d i d v η β 3 Λ i n d v η ε β 1 Λ i n = ( η + d l ) d i d v + ( 1 R 0 ) . a 2 a 1 a 0 = d i + d v + ( η + d l ) β 2 Λ i n d i d v + [ ( η + d l ) β 2 Λ i n ] ( d i + d v ) η β 3 Λ i n [ ( η + d l ) β 2 Λ i n ] d i d v + η β 3 Λ i n d v + η ε β 1 Λ i n = d i + d v + ( η + d l ) β 2 Λ i n [ ( η + d l ) β 2 Λ i n ] ( d i + d v ) η β 3 Λ i n + ( d i + d v ) d i d v + η β 3 Λ i n d v + η ε β 1 Λ i n = d i + d v + ( η + d l ) β 2 Λ i n ( η + d l ) d v ( 1 R 2 ) + d i ( η + d l ) ( 1 ( R 2 + R 3 ) ) + ( d i + d v ) d i d v + η β 3 Λ i n d v + η ε β 1 Λ i n > 0 , f o r   R 0 < 1 .
It is easy to see that for R 0 < 1 , we have a 2 > 0 , a 1 > 0 , a 0 > 0 , a 2 a 1 > a 0 . By applying Routh-Hurwitz criteria [23,24], all eigenvalues admit negative real parts, and E 0 is locally asymptotically stable if R 0 < 1 , and it is unstable if R 0 > 1 (See [25] for another application). □

3.3. Global Analysis

Consider the function G defined by G ( x ) = x 1 ln ( x )
Theorem 2.
The virus-free steady state E 0 is globally asymptotically stable if R 0 1 .
Proof. 
Consider the case where R 0 1 and consider the Lyapunov function L 0 given as follows:
L 0 ( Y s , Y l , Y i , Y v , Y c ) = Λ i n G Y s Λ i n + Y l + ( η + d l ) η ( 1 R 2 ) Y i + β 1 Λ i n d v Y v + d i ( η + d l ) λ η ( 1 R 0 ) Y c .
Clearly, L 0 ( Y s , Y l , Y i , Y v , Y c ) > 0 for all Y s , Y l , Y i , Y v , Y c > 0 , and L 0 ( Λ i n , 0 , 0 , 0 , 0 ) = 0 . Furthermore, L 0 admits the following derivative along (2).
d L 0 d t = 1 Λ i n Y s d s Λ i n d s Y s β 1 Y s Y v β 2 Y s Y l β 3 Y s Y i + β 1 Y s Y v + β 2 Y s Y l + β 3 Y s Y i η + d l Y l + η + d l 1 R 2 η η Y l d i Y i γ Y c Y i + β 1 Λ i n d v ε Y i d v Y v + d i η + d l 1 R 0 λ η λ Y i d c Y c δ Y c Y i = d s Y s Λ i n 2 Y s + β 2 Λ i n η + d l + η + d l 1 R 2 Y l + β 3 Λ i n d i η + d l 1 R 2 η + Λ i n β 1 ε d v + d i η + d l 1 R 0 η Y i d i d c η + d l 1 R 0 λ η Y c γ η + d l 1 R 2 η + d i δ η + d l 1 R 0 λ η Y i Y c = d s Y s Λ i n 2 Y s d i d c η + d l 1 R 0 λ η Y c γ η + d l 1 R 2 η + d i δ η + d l 1 R 0 λ η Y i Y c .
If R 0 1 , then d L 0 d t 0 for all Y s , Y l , Y i , Y v , Y c > 0 . Let W 0 = { ( Y s , Y l , Y i , Y v , Y c ) : d L 0 d t = 0 } . It can be easily shown that W 0 = { E 0 } . By applying LaSalle’s invariance principle [26], one obtain that E 0 is globally asymptotically stable if R 0 1 (see [27] for an other application). □
We give now the following result regarding the global stability of the virus-present steady state, E when R 0 > 1 .
Theorem 3.
By considering the dynamics (2), if R 0 > 1 , then E is globally asymptotically stable.
Proof. 
Consider the Lyapunov function L ( Y s , Y l , Y i , Y v , Y c ) defined as:
L ( Y s , Y l , Y i , Y v , Y c ) = Y s G Y s Y s + Y l G Y l Y l + Y s ( β 1 ε + β 3 d v ) Y i d v ( d i + γ Y c ) G Y i Y i + β 1 Y s Y v d v G Y v Y v + γ Y s ( β 1 ε + β 3 d v ) 2 d v ( d i + γ Y c ) ( λ δ Y c ) ( Y c Y c ) 2 .
Clearly, L ( Y s , Y l , Y i , Y v , Y c ) > 0 for all Y s , Y l , Y i , Y v , Y c > 0 and L ( Y s , Y l , Y i , Y v , Y c ) = 0 . By calculating d L d t along the dynamics (2), we obtain:
d L d t = 1 Y S Y s d S Λ i n d S Y S β 1 Y S Y v β 2 Y S Y l β 3 Y S Y i + β 1 Y S d v 1 Y v Y v ε Y i d v Y v + 1 Y l Y l β 1 Y S Y v + β 2 Y S Y l + β 3 Y S Y i η + d l Y l + Y S β 1 ε + β 3 d v d v d i + γ Y c 1 Y i Y i η Y l d i Y i γ Y c Y i + γ Y S β 1 ε + β 3 d v d v d i + γ Y c λ δ Y c Y c Y c λ Y i d c Y c δ Y c Y i = d S 1 Y S Y s Λ i n Y S + β 2 Y S η + d l + η Y S β 1 ε + β 3 d v d v d i + γ Y c Y l + β 3 Y S Y i β 1 Y S Y v + β 2 Y S Y l + β 3 Y S Y i Y l Y l + η + d l Y l η β 1 ε + β 3 d v d v d i + γ Y c Y S Y l Y i Y i d i β 1 ε + β 3 d v d v d i + γ Y c Y S Y i Y i γ β 1 ε + β 3 d v d v d i + γ Y c Y S Y c Y i Y i β 1 ε d v Y S Y i Y v Y v + β 1 ε d v Y S Y i + β 1 Y S Y v + γ Y S β 1 ε + β 3 d v d v d i + γ Y c λ δ Y c Y c Y c λ Y i d c Y c δ Y c Y i .
Now, since we have d s Λ i n = d s Y s + β 1 Y s Y v + β 2 Y s Y l + β 3 Y s Y i , β 1 Y s Y v + β 2 Y s Y l + β 3 Y s Y i = ( η + d l ) Y l , Y i = η Y l d i + γ Y c , Y v = ε Y i d v and λ Y i d c Y c δ Y c Y i , we obtain
β 1 Y s Y v + β 3 Y s Y i = Y s ( β 1 ε + β 3 d v ) Y i d v = η Y s ( β 1 ε + β 3 d v ) d v ( d i + γ Y c )
and
β 2 Y s ( η + d l ) + η Y s ( β 1 ε + β 3 d v ) d v ( d i + γ Y c ) Y l = 0 .
Therefore, we obtain
d L d t = d s 1 Y s Y s Y s Y s + β 1 Y s Y v + β 2 Y s Y l + β 3 Y s Y i 1 Y s Y s + β 3 Y s Y i β 1 Y s Y v Y l Y l β 2 Y s Y l β 3 Y s Y i Y l Y l + 2 β 1 Y s Y v + β 2 Y s Y l + β 3 Y s Y i + β 1 ε d v Y s Y i η β 1 ε + β 3 d v d v d i + γ Y c Y s Y l Y i Y i d i β 1 ε + β 3 d v d v d i + γ Y c Y s Y i Y i β 1 ε d v Y s Y i Y v Y v γ β 1 ε + β 3 d v d v d i + γ Y c Y s Y c Y i Y i + γ β 1 ε + β 3 d v d v d i + γ Y c λ δ Y c Y s Y c Y c ×          λ Y i d c Y c δ Y c Y i λ Y i + d c Y c + δ Y c Y i = d s Y s Y s 2 Y s + β 1 Y s Y v + β 2 Y s Y l + β 3 Y s Y i 2 Y s Y s + β 3 Y s Y i β 2 Y s Y l β 1 Y s Y v Y l Y l β 3 Y s Y i Y l Y l β 1 Y s Y v + β 3 Y s Y i Y l Y i Y l Y i + β 1 Y s Y v β 1 ε + β 3 d v d v d i + γ Y c Y s d i + γ Y c Y i Y i + β 1 Y s Y v Y i Y i β 1 Y s Y i Y v 2 Y i Y v γ β 1 ε + β 3 d v d c + δ Y i d v d i + γ Y c λ δ Y c Y s Y c Y c 2 = d s Y s Y s 2 Y s + β 1 Y s Y v + β 2 Y s Y l + β 3 Y s Y i 2 Y s Y s + β 3 Y s Y i β 1 Y s Y v Y l Y l β 2 Y s Y l β 3 Y s Y i Y l Y l β 1 Y s Y v + β 3 Y s Y i Y l Y i Y l Y i β 1 Y s Y v + β 3 Y s Y i Y i Y i 1 + β 1 Y s Y i Y v Y i β 1 Y s Y i Y v 2 Y i Y v + β 1 Y s Y v γ β 1 ε + β 3 d v d c + δ Y i d v d i + γ Y c λ δ Y c Y s Y c Y c 2 = d s + β 2 Y l Y s Y s 2 Y s γ β 1 ε + β 3 d v d c + δ Y i d v d i + γ Y c λ δ Y c Y s Y c Y c 2 + β 1 Y s Y v 4 Y s Y s Y l Y i Y l Y i Y i Y v Y i Y v Y s Y v Y l Y s Y v Y l + β 3 Y s Y i 3 Y s Y s Y l Y i Y l Y i Y s Y i Y l Y s Y i Y l .
Since we have z 1 + z 2 + z 3 + + z n n z 1 · z 2 · z 3 z n n for all z 1 , z 2 , z 3 , , z n 0 with positive integer n. Then we deduce that
3 Y s Y s Y l Y i Y l Y i Y s Y i Y l Y s Y i Y l 0
and
4 Y s Y s Y l Y i Y l Y i Y i Y v Y i Y v Y s Y v Y l Y s Y v Y l 0 .
Finally, we obtain d L d t ( Y s , Y l , Y i , Y v , Y c ) 0 for all Y s , Y l , Y i , Y v , Y c > 0 and d L d t = 0 if and only if ( Y s , Y l , Y i , Y v , Y c ) = ( Y s , Y l , Y i , Y v , Y c ) . By applying the LaSalle’s invariance principle [26], we conclude that E is globally asymptotically stable (see [27] for other application). □

4. Variable Environment and Periodic Solution

Let us go back to the main dynamics (1) such that the positive parameters of the model are continuous and T-periodic functions. Let us denote by A ( t ) a continuous T-periodic m × m irreducible and cooperative matrix function and let κ A ( t ) be the matrix solution of the equation hereafter
z ˙ ( t ) = A ( t ) z ( t ) .
We denote by r ( κ A ( T ) ) the spectral radius of κ A ( T ) , which has positive entries for all t > 0 . By using the theorem of Perron–Frobenius [28], we can obtain that κ A ( T ) admits a principal eigenvalue r ( κ A ( T ) ) , and we consider the lemma hereafter that will be used later frequently.
Lemma 3
([29]). Equation (7) admits a solution of the form z ( t ) = y ( t ) e t where y ( t ) is a positive T-periodic function and = 1 T ln ( r ( κ A ( T ) ) ) .
With the aim to prove that dynamics (1) admits a unique virus-free periodic solution, let the equation
Y s ˙ ( t ) = d s ( t ) Λ i n ( t ) d s ( t ) Y s ( t ) ,
where the initial condition Y s 0 R + . Equation (8) admits a unique T-periodic solution Y s ( t ) globally attractive in R + , such that Y s ( t ) > 0 . Thus, the dynamics (1) admits a unique virus-free periodic solution E 0 ( t ) = ( Y s ( t ) , 0 , 0 , 0 , 0 ) . Consider a continuous, positive T-periodic function σ ( t ) , let σ u = max t [ 0 , T ) σ ( t ) and σ l = min t [ 0 , T ) σ ( t ) . Let us define d ( t ) = min t 0 ( d v ( t ) , d c ( t ) ) , and then we obtain the following Proposition.
Proposition 1.
The compact positive set
Σ u = ( Y s , Y l , Y i , Y v , Y c ) R + 5 / Y s + Y l + Y i Λ i n u ; Y v + Y c ( ε u + λ u ) Λ i n u d l
is invariant and an attractor of all trajectories of model (1) with
lim t Y s ( t ) + Y l ( t ) + Y i ( t ) Y s ( t ) = 0 .
Proof. 
By considering the model (1), we obtain
Y ˙ s ( t ) + Y ˙ l ( t ) + Y ˙ i ( t ) = d s ( t ) Λ i n ( t ) d s ( t ) Y s ( t ) d l ( t ) Y l ( t ) d i ( t ) Y i ( t ) η ( t ) Y l ( t ) γ ( t ) Y c ( t ) Y i ( t ) d s ( t ) Λ i n ( t ) ( Y s ( t ) + Y l ( t ) + Y i ( t ) ) 0 , if ( Y s ( t ) + Y l ( t ) + Y i ( t ) ) Λ i n u ,
and
Y ˙ v ( t ) + Y ˙ c ( t ) = ( ε ( t ) + λ ( t ) ) Y i ( t ) d v ( t ) Y v ( t ) d c ( t ) Y c ( t ) δ ( t ) Y c ( t ) Y i ( t ) ( ε ( t ) + λ ( t ) ) Y i ( t ) d v ( t ) Y v ( t ) d c ( t ) Y c ( t ) ( ε u + λ u ) Λ i n u d ( t ) ( Y v ( t ) + Y c ( t ) ) 0 , if d ( t ) ( Y v ( t ) + Y c ( t ) ) ( ε u + λ u ) Λ i n u .
In Section 4.1, we define the basic reproduction number R 0 according to the theory in [30]. It will be proved that if R 0 < 1 , then the virus-free periodic solution E 0 ( t ) = ( Y s ( t ) , 0 , 0 , 0 , 0 ) is globally asymptotically stable and the virus goes to extinction. In Section 4.2, it will be proven that if R 0 > 1 , then the dynamics (1) will be uniformly persistent. This means that R 0 plays the role of the threshold parameter between uniform persistence and extinction of the virus.

4.1. Virus-Free Periodic Solution

According to the theory in [30] and to define the basic reproduction number R 0 , we rewrite the dynamics (1) in a more suitable form with the theory of [30]. Let Y ( t ) = Y l ( t ) , Y i ( t ) , Y v ( t ) , Y s ( t ) , Y c ( t ) T , F ( t , Y ( t ) ) = ( β 1 ( t ) Y v ( t ) + β 2 ( t ) Y l ( t ) + β 3 ( t ) Y i ( t ) ) Y s ( t ) , η ( t ) Y l ( t ) , ε ( t ) Y i ( t ) , 0 , 0 T , V ( t , Y ( t ) ) = ( η ( t ) + d l ( t ) ) Y l ( t ) , d i ( t ) Y i ( t ) + γ ( t ) Y c ( t ) Y i ( t ) , d v ( t ) Y v ( t ) , ( d s ( t ) + β 1 ( t ) Y v ( t ) + β 2 ( t ) Y l ( t ) + β 3 ( t ) Y i ( t ) ) Y s ( t ) , d c ( t ) Y c ( t ) + δ ( t ) Y c ( t ) Y i ( t ) T and V + ( t , Y ( t ) ) = 0 , 0 , 0 , d s ( t ) Λ i n ( t ) , λ ( t ) Y i ( t ) T . We aim to satisfy Assumptions (A1)–(A7) in [30]. The dynamics (1) can be written as follows
Y ˙ ( t ) = F ( t , Y ( t ) ) V ( t , Y ( t ) ) = F ( t , Y ( t ) ) V ( t , Y ( t ) ) + V + ( t , Y ( t ) ) .
One can easily see that Assumptions (A1)–(A5) are verified.
The dynamics (10) has a virus-free periodic solution, Y ( t ) = 0 , 0 , 0 , Y s ( t ) , 0 T .
Consider the functions f ( t , Y ( t ) ) = F ( t , Y ( t ) ) V ( t , Y ( t ) ) + V + ( t , Y ( t ) ) , and M ( t ) = f i ( t , Y ( t ) ) Y j 4 i , j 5 , such that f i ( t , Y ( t ) ) and Y i ( t ) are the i-th component of f ( t , Y ( t ) ) and Y ( t ) , respectively. A simple calculus allows us to obtain M ( t ) = d s ( t ) 0 0 d c ( t ) ; thus, r ( ϕ M ( T ) ) < 1 . Therefore, the virus-free solution Y ( t ) is asymptotically stable in Σ s ; given by
Σ s = ( 0 , 0 , 0 , Y s , 0 ) R + 5 ,
which means that condition (A6) of [30] is verified.
Let us consider F ( t ) and V ( t ) to be 3 by 3 matrices defined by F ( t ) = F i ( t , Y ( t ) ) Y j 1 i , j 3 , and V ( t ) = V i ( t , Y ( t ) ) Y j 1 i , j 3 , such that F j ( t , Y ( t ) ) and V j ( t , Y ( t ) ) are the j-th component of F ( t , Y ( t ) ) and V ( t , Y ( t ) ) , respectively. An easy calculus allows us to conclude from dynamics (10) that
F ( t ) = β 2 ( t ) Y s ( t ) β 3 ( t ) Y s ( t ) β 1 ( t ) Y s ( t ) η ( t ) 0 0 0 ε ( t ) 0 and   V ( t ) = η ( t ) + d l ( t ) 0 0 0 d i ( t ) 0 0 0 d v ( t ) .
Consider the equation d d t Z ( t 1 , t 2 ) = V ( t 1 ) Z ( t 1 , t 2 ) with t 1 t 2 and Z ( t 1 , t 1 ) = I , the three-by-three identity matrix, and let Z ( t 1 , t 2 ) be its three-by-three matrix solution. Then, condition (A7) of [30] is verified.
Consider K : C T C T to be the linear operator given hereafter
( K ϕ ) ( p ) = 0 Z ( p , p z ) F ( p z ) ϕ ( p z ) d z , p R , ϕ C T
where C T is the Banach space of functions that are T-periodic and defined on R R 3 , equipped with the norm . . Now, we can define the expression of the basic reproduction number, R 0 associated with the dynamics (1), as the spectral radius of K, defined previously,
R 0 = r ( K ) .
We give in the theorem hereafter, the local stability result for the virus-free periodic solution E 0 ( t ) = ( Y s ( t ) , 0 , 0 , 0 , 0 ) of the dynamics (1) according to the definition of the basic reproduction number [30].
Theorem 4
([30], Theorem 2.2).
  • R 0 < 1 r ( ϕ F V ( T ) ) < 1 .
  • R 0 = 1 r ( ϕ F V ( T ) ) = 1 .
  • R 0 > 1 r ( ϕ F V ( T ) ) > 1 .
Thus, E 0 ( t ) is locally asymptotically stable once R 0 < 1 , however it is unstable once R 0 > 1 .
Theorem 5.
The virus-free trajectory E 0 ( t ) is globally asymptotically stable if R 0 < 1 ; however, it is unstable once R 0 > 1 .
Proof. 
By applying Theorem 4, E 0 ( t ) is locally asymptotically stable once R 0 < 1 ; however, it is unstable once R 0 > 1 . Thus, we need to prove the global attractivity of the virus-free solution E 0 ( t ) once R 0 < 1 . Let us concentrate on the case where R 0 < 1 . By the limit (9) in Lemma 1, for any ς 1 > 0 , there exists a positive constant T 1 > 0 such that Y s ( t ) + Y l ( t ) + Y i ( t ) Y s ( t ) + ς 1 for t > T 1 . Therefore, Y s ( t ) Y s ( t ) + ς 1 , and we obtain the following subsystem
Y l ˙ ( t ) [ β 1 ( t ) Y v ( t ) + β 2 ( t ) Y l ( t ) + β 3 ( t ) Y i ( t ) ] ( Y s ( t ) + ς 1 ) ( η ( t ) + d l ( t ) ) Y l ( t ) , Y i ˙ ( t ) = η ( t ) Y l ( t ) d i ( t ) Y i ( t ) γ ( t ) Y c ( t ) Y i ( t ) , Y v ˙ ( t ) = ε ( t ) Y i ( t ) d v ( t ) Y v ( t ) ,
for t > T 1 . Let us define M 2 ( t ) to be the three-by-three matrix function given hereafter.
M 2 ( t ) = β 2 ( t ) ( Y s ( t ) + ς 1 ) β 3 ( t ) ( Y s ( t ) + ς 1 ) β 1 ( t ) ( Y s ( t ) + ς 1 ) η ( t ) 0 0 0 ε ( t ) 0 .
Similarly, we use Theorem 4, we obtain r ( φ F V ( T ) ) < 1 , and we choose ς 1 > 0 such that r ( φ F V + ς 1 M 2 ( T ) ) < 1 . Let the three-equations system hereafter be
Y ¯ ˙ l ( t ) = [ β 1 ( t ) Y ¯ v ( t ) + β 2 ( t ) Y ¯ l ( t ) + β 3 ( t ) Y ¯ i ( t ) ] ( Y s ( t ) + ς 1 ) ( η ( t ) + d l ( t ) ) Y ¯ l ( t ) , Y ¯ ˙ i ( t ) = η ( t ) Y l ( t ) d i ( t ) Y ¯ i ( t ) γ ( t ) Y ¯ c ( t ) Y ¯ i ( t ) , Y ¯ ˙ v ( t ) = ε ( t ) Y ¯ i ( t ) d v ( t ) Y ¯ v ( t ) .
By using Lemma 3 and by applying the comparison principle, it can easily deduced that there exists a T-periodic positive function y 1 ( t ) such that x ( t ) y 1 ( t ) e k 1 t with x ( t ) = Y l ( t ) , Y i ( t ) , Y v ( t ) and k 1 = 1 T ln r ( φ F V + ς 1 M 2 ( T ) < 0 . Therefore, lim t Y l ( t ) = lim t Y i ( t ) = lim t Y v ( t ) = 0 . Then, we deduce that lim t Y c ( t ) = 0 . Moreover, we obtain lim t ( Y s ( t ) Y s ( t ) ) = 0 . Therefore, it is deduced that the virus-free periodic solution E 0 ( t ) is globally attractive, which permits us to complete the proof. □
Let us now concentrate on the case for which the basic reproduction number satisfies R 0 > 1 .

4.2. Virus-Infected Periodic Solution

Let us define the Poincaré map Q : R + 5 R + 5 associated with the dynamics (1) with Y 0 w ( T , Y 0 ) and w ( t , Y 0 ) is the unique solution of (1) with the initial condition w ( 0 , Y 0 ) = Y 0 R + 4 . Let us consider the following sets Ω = ( Y s , Y l , Y i , Y v , Y c ) R + 5 , Ω 0 = I n t ( R + 5 ) and Ω 0 = Ω Ω 0 .
Note that both Ω and Ω 0 are positively invariant from Lemma 1. Q is point dissipative. Define
M = ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) Ω 0 : Q n ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) Ω 0 , n 0 .
To use the uniform persistence theory given by Zhao [29,31], we need to prove that
M = ( Y s , 0 , 0 , 0 , 0 ) , Y s 0 .
It easy to see that M ( Y s , 0 , 0 , 0 , 0 ) , Y s 0 . To show that M { ( Y s , 0 , 0 , 0 , 0 ) , Y s 0 } = , let ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) M ( Y s , 0 , 0 , 0 , 0 ) , Y s 0 .
If Y i 0 = 0 and 0 < Y l 0 , then Y l ( t ) > 0 for any t > 0 . Then, it holds that Y ˙ i ( t ) | t = 0 = η ( 0 ) Y l 0 > 0 . If Y i 0 > 0 and Y l 0 = 0 , then Y i ( t ) > 0 and Y s ( t ) > 0 for any t > 0 . Therefore, for any t > 0 , we have
Y l ( t ) = Y l 0 + 0 t [ β 1 ( ω ) Y v ( ω ) + β 2 ( ω ) Y l ( ω ) + β 3 ( ω ) Y i ( ω ) ] Y s ( ω ) e 0 ω ( η ( s ) + d l ( s ) ) d s d ω e 0 t ( η ( s ) + d l ( s ) ) d s > 0
for any t > 0 . Therefore, ( Y s ( t ) , Y l ( t ) , Y i ( t ) , Y v ( t ) , Y c ( t ) ) Ω 0 for 0 < t small enough. Since Ω 0 is positively invariant according to Proposition 1, then we deduced Equation (15). Therefore, there exists a unique fixed point ( Y s ( 0 ) , 0 , 0 , 0 , 0 ) of Q in M , which means that the disease will persist according to the following theorem:
Theorem 6.
Assume that R 0 > 1 . Dynamics (1) admits at least one positive periodic solution and there exists a positive constant ϖ > 0 , such that ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) R + × I n t ( R + 5 ) ,
lim inf t Y i ( t ) ϖ > 0 .
Proof. 
We aim to prove that Q is uniformly persistent, with respect to ( Ω 0 , Ω 0 ) . This permits us to prove the uniform persistence of the solution of model (1) with respect to ( Ω 0 , Ω 0 ) when applying ([31], Theorem 3.1.1). From Theorem 4, we have r ( φ F V ( T ) ) > 1 . Therefore, ζ > 0 sufficiently small such that r ( φ F V ζ M 2 ( T ) ) > 1 . Define the following perturbed system
Y ˙ s α ( t ) = d s ( t ) Λ i n ( t ) d s ( t ) Y s α ( t ) [ β 1 ( t ) + β 2 ( t ) + β 3 ( t ) ] Y s α ( t ) α ,
By associating model (16), function Q admits a unique fixed point Y ¯ s α 0 that it is positive and globally attractive. According to the implicit function theorem we can deduce that Y ¯ s α 0 is continuous with respect to α . Therefore, it is sufficient to choose α > 0 sufficiently small and satisfying Y ¯ s α ( t ) > Y ¯ s ( t ) ζ , t > 0 . Let M 1 = ( Y ¯ s 0 , 0 , 0 , 0 , 0 ) . As the solutions are continuous with respect to the initial condition, α satisfies ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) Ω 0 with ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) M 1 α , and we deduce that
w ( t , ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) ) w ( t , M 1 ) < α   for   0 t T .
Let us prove by contradiction that
lim sup n d ( Q n ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) , M 1 ) α   for   any   ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) Ω 0 .
Assume that lim sup n d ( Q n ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) , M 1 ) < α for some ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) Ω 0 . In particular, assume that d ( Q n ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) , M 1 ) < α , n > 0 . Then, we deduce that w ( t , Q n ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) ) w ( t , M 1 ) < α for all n > 0 and 0 t T . For any t 0 , assume that t = n T + t 1 , with t 1 [ 0 , T ) and n t T is the greatest integer of t T . Then, we obtain w ( t , ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) ) w ( t , M 1 ) = w ( t 1 , Q n ( S 0 , V 0 , E 0 , I 0 ) ) w ( t 1 , M 1 ) < α , t 0 . Let ( Y s ( t ) , Y l ( t ) , Y i ( t ) , Y v ( t ) , Y c ( t ) ) = w ( t , ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) ) . Therefore, 0 Y l ( t ) , Y i ( t ) , Y v ( t ) α for all t 0 and
Y ˙ s ( t ) d s ( t ) Λ i n ( t ) d s ( t ) Y s ( t ) α ( β 1 ( t ) + β 2 ( t ) + β 3 ( t ) ) Y s ( t ) .
The function Q admits the fixed point Y ¯ s α 0 when associated with model (16), such that it is globally attractive with Y ¯ s α ( t ) > Y ¯ s ( t ) ζ , then T 2 > 0 large enough such that
Y ¯ s ( t ) > Y ¯ s ( t ) ζ , t > T 2 .
Thus, for all t > T 2
Y l ˙ ( t ) [ β 1 ( t ) Y v ( t ) + β 2 ( t ) Y l ( t ) + β 3 ( t ) Y i ( t ) ] ( Y ¯ s ( t ) ζ ) ( η ( t ) + d l ( t ) ) Y l ( t ) , Y i ˙ ( t ) = η ( t ) Y l ( t ) d i ( t ) Y i ( t ) γ ( t ) Y c ( t ) Y i ( t ) , Y v ˙ ( t ) = ε ( t ) Y i ( t ) d v ( t ) Y v ( t ) .
Since we have r ( φ F V η M 2 ( T ) ) > 1 , then by applying again Lemma 3 and the comparison principle, there exists a T-periodic positive solution y 2 ( t ) that satisfies J ( t ) e k 2 t y 2 ( t ) and k 2 = 1 T ln r φ F V η M 2 ( T ) > 0 . Thus, it can be deduced that lim t Y i ( t ) = , which contradicts the fact that the solution is bounded. Then, the inequality (17) is verified, and then the function Q is weakly uniformly persistent with respect to ( Ω 0 , Ω 0 ) . By using Lemma 1, the function Q admits a global attractor. Therefore, we obtain that the set M 1 = ( Y ¯ s 0 , 0 , 0 , 0 , 0 ) is invariant in Ω and W s ( M 1 ) Ω 0 = . All trajectories in M converge to M 1 , such that it is acyclic in M . According to ([31], Theorem 1.3.1 and Remark 1.3.1), we obtain the uniform persistence of the function Q with respect to ( Ω 0 , Ω 0 ) . Moreover, by applying ([31], Theorem 1.3.6), the function Q has a fixed point ( Y ˜ s 0 , Y ˜ l 0 , Y ˜ i 0 , Y ˜ v 0 , Y ˜ c 0 ) Ω 0 . Note that ( Y ˜ s 0 , Y ˜ l 0 , Y ˜ i 0 , Y ˜ v 0 , Y ˜ c 0 ) R + × I n t ( R + 4 ) . We aim to prove also by contradiction that Y ˜ s 0 > 0 . Suppose that Y ˜ s 0 = 0 . From the first equation of dynamics (1), Y ˜ s ( t ) satisfies
Y ˜ ˙ s ( t ) d s ( t ) Λ i n ( t ) d s ( t ) Y ˜ s ( t ) ( β 1 ( t ) Y ˜ v ( t ) + β 2 ( t ) Y ˜ l ( t ) + β 3 ( t ) Y ˜ i ( t ) ) Y ˜ s ( t ) ,
with Y ˜ s 0 = Y ˜ s ( m T ) = 0 , m = 1 , 2 , 3 , . By using Lemma 1, ς 3 > 0 , there exists T 3 > 0 large enough that satisfies Y ˜ l ( t ) , Y ˜ i ( t ) , Y ˜ v ( t ) N ¯ + ς 3 , t > T 3 . Therefore, we obtain
Y ˜ ˙ s ( t ) d s ( t ) Λ i n ( t ) d s ( t ) Y ˜ s ( t ) ( β 1 ( t ) + β 2 ( t ) + β 3 ( t ) ) Y ˜ s ( t ) ( N ¯ + ς 3 ) , for   t T 3 .
There exists a large enough m ¯ that satisfies m T > T 3 for all m > m ¯ . By using the comparison principle, it can be deduced that
Y ˜ s ( m T ) = Y ˜ s 0 + 0 m T d s ( ω ) Λ i n ( ω ) e 0 ω [ β 1 ( u ) + β 2 ( u ) + β 3 ( u ) ] ( N ¯ + ς 3 ) + d s ( u ) d u d ω e 0 m T ( [ β 1 ( u ) + β 2 ( u ) + β 3 ( u ) ] ( N ¯ + ς 3 ) + d s ( u ) ) d u > 0
for all m > m ¯ , and this is impossible. Thus, Y ˜ s 0 > 0 and the T-periodic trajectory ( Y ˜ s 0 , Y ˜ l 0 , Y ˜ i 0 , Y ˜ v 0 , Y ˜ c 0 ) is a positive solution of the dynamics (1). □

5. Numerical Examples

We illustrate the obtained theoretical results on the dynamics (1) by several numerical examples. The T-periodic functions that we used as model parameters are given hereafter.
Λ i n ( t ) = Λ i n 0 ( 1 + Λ i n 1 cos ( 2 π ( t + ψ ) ) ) , d s ( t ) = d s 0 ( 1 + d s 1 cos ( 2 π ( t + ψ ) ) ) , β 1 ( t ) = β 10 ( 1 + β 11 cos ( 2 π ( t + ψ ) ) ) , d l ( t ) = d l 0 ( 1 + d l 1 cos ( 2 π ( t + ψ ) ) ) , β 2 ( t ) = β 20 ( 1 + β 21 cos ( 2 π ( t + ψ ) ) ) , d i ( t ) = d i 0 ( 1 + d i 1 cos ( 2 π ( t + ψ ) ) ) , β 3 ( t ) = β 30 ( 1 + β 31 cos ( 2 π ( t + ψ ) ) ) , d v ( t ) = d v 0 ( 1 + d v 1 cos ( 2 π ( t + ψ ) ) ) , η ( t ) = η 0 ( 1 + η 1 cos ( 2 π ( t + ψ ) ) ) , d c ( t ) = d c 0 ( 1 + d c 1 cos ( 2 π ( t + ψ ) ) ) , γ ( t ) = γ 0 ( 1 + γ 1 cos ( 2 π ( t + ψ ) ) ) , δ ( t ) = δ 0 ( 1 + δ 1 cos ( 2 π ( t + ψ ) ) ) , ε ( t ) = ε 0 ( 1 + ε 1 cos ( 2 π ( t + ψ ) ) ) , λ ( t ) = λ 0 ( 1 + λ 1 cos ( 2 π ( t + ψ ) ) ) .
Λ i n 1 , d s 1 , d l 1 , d i 1 , d v 1 , d c 1 , β 11 , β 21 , β 31 , η 1 , γ 1 , δ 1 , ε 1 and λ 1 are the amplitudes reflecting the influence of the variations on all of the parameters where | Λ i n 1 |   < 1 , | d s 1 |   < 1 , | d l 1 |   < 1 , | d i 1 |   < 1 , | d v 1 |   < 1 , | d c 1 |   < 1 , | β 11 |   < 1 , | β 21 |   < 1 , | β 31 |   < 1 , | η 1 |   < 1 , | γ 1 |   < 1 , | δ 1 |   < 1 , | ε 1 |   < 1 and | λ 1 |   < 1 . The phase shift is given by ψ . All fixed constants used in the numerical examples are summarized in Table 2.
Three environmental situations were considered. The first situation was dedicated to a totally fixed environment and all parameters are constants, and thus we obtain an autonomous system. We aim then to confirm the theoretical findings regarding the stability of the steady states E 0 and E . The second situation was dedicated to a variable contact, which means that only the contact rates β 1 ( t ) , β 2 ( t ) and β 3 ( t ) are T-periodic, and the rest of the parameters are fixed. The third situation deals with a total variable environment such that all the parameters are T-periodic functions.

5.1. Totally Fixed Environment

We first consider the case of a totally fixed environment, and thus all the parameters of system (1) are constants as follows:
Y s ˙ ( t ) = d s 0 Λ i n 0 d s 0 Y s ( t ) [ β 10 Y v ( t ) + β 20 Y l ( t ) + β 30 Y i ( t ) ] Y s ( t ) , Y l ˙ ( t ) = [ β 10 Y v ( t ) + β 20 Y l ( t ) + β 30 Y i ( t ) ] Y s ( t ) ( η 0 + d l 0 ) Y l ( t ) , Y i ˙ ( t ) = η 0 Y l ( t ) d i 0 ( t ) Y i ( t ) γ 0 Y c ( t ) Y i ( t ) , Y v ˙ ( t ) = ε 0 Y i ( t ) d v 0 Y v ( t ) , Y c ˙ ( t ) = λ 0 Y i ( t ) d c 0 Y c ( t ) δ 0 Y c ( t ) Y i ( t ) ,
with initial condition ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) = ( 10 , 10 , 6 , 8 , 10 ) R + 5 . We provide several numerical examples that confirm the stability of the steady states of (23). In Figure 1, we display the variables’ behavior when R 0 < 1 . The trajectory of the system (23) converges to the steady state E 0 , and this confirms the global stability of E 0 when R 0 1 . In Figure 2, we display the behavior of all components of the solution for several initial conditions, confirming the global stability of E 0 if R 0 < 1 . However, in Figure 3, we display the variables’ behavior when R 0 > 1 . The trajectory of the system (23) converges to the steady state E , confirming the global stability of E when R 0 > 1 . In Figure 4, we display the behavior of all components of the solution for different initial conditions, confirming the global stability of E if R 0 > 1 .

5.2. Variable Contact Rates

Let us now perform numerical simulations on System (1) where only the variable forced T-periodic function β depends on time, t. The other parameters are constant. Thus, the model is given by
Y s ˙ ( t ) = d s 0 Λ i n 0 d s 0 Y s ( t ) [ β 1 ( t ) Y v ( t ) + β 2 ( t ) Y l ( t ) + β 3 ( t ) Y i ( t ) ] Y s ( t ) , Y l ˙ ( t ) = [ β 1 ( t ) Y v ( t ) + β 2 ( t ) Y l ( t ) + β 3 ( t ) Y i ( t ) ] Y s ( t ) ( η 0 + d l 0 ) Y l ( t ) , Y i ˙ ( t ) = η 0 Y l ( t ) d i 0 ( t ) Y i ( t ) γ 0 Y c ( t ) Y i ( t ) , Y v ˙ ( t ) = ε 0 Y i ( t ) d v 0 Y v ( t ) , Y c ˙ ( t ) = λ 0 Y i ( t ) d c 0 Y c ( t ) δ 0 Y c ( t ) Y i ( t ) ,
with positive initial condition ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) = ( 10 , 10 , 6 , 8 , 10 ) R + 5 . Note that R 0 was calculated by using the time-averaged system as in [32,33]. Other definitions of R 0 can be found in [34,35]. In Figure 5, we display the variables’ behavior when R 0 < 1 . The trajectory of the system (24) converges to the steady state E 0 = Λ i n 0 , 0 , 0 , 0 , 0 , 0 , and this confirms the global stability of E 0 when R 0 1 . In Figure 6, we display the behavior of all components of the solution for different initial conditions, confirming the global stability of E 0 if R 0 < 1 . However, in Figure 7, we display the variables’ behavior when R 0 > 1 . The trajectory of the system (24) converges to the steady state E , and this confirms the global stability of E when R 0 > 1 . In Figure 8, we display the behavior of all components of the solution for different initial conditions, confirming the global stability of E if R 0 > 1 .

5.3. Totally Variable Environment

In this subsection, we consider the dynamics (1) such that all parameters are T-periodic functions as follows:
Y s ˙ ( t ) = d s ( t ) Λ i n ( t ) d s ( t ) Y s ( t ) [ β 1 ( t ) Y v ( t ) + β 2 ( t ) Y l ( t ) + β 3 ( t ) Y i ( t ) ] Y s ( t ) , Y l ˙ ( t ) = [ β 1 ( t ) Y v ( t ) + β 2 ( t ) Y l ( t ) + β 3 ( t ) Y i ( t ) ] Y s ( t ) ( η ( t ) + d l ( t ) ) Y l ( t ) , Y i ˙ ( t ) = η ( t ) Y l ( t ) d i ( t ) Y i ( t ) γ ( t ) Y c ( t ) Y i ( t ) , Y v ˙ ( t ) = ε ( t ) Y i ( t ) d v ( t ) Y v ( t ) , Y c ˙ ( t ) = λ ( t ) Y i ( t ) d c ( t ) Y c ( t ) δ ( t ) Y c ( t ) Y i ( t ) ,
with positive initial condition ( Y s 0 , Y l 0 , Y i 0 , Y v 0 , Y c 0 ) = ( 10 , 10 , 6 , 8 , 10 ) R + 5 . Similarly to the case of system (24), R 0 was calculated using a time-averaged system as in [32,33]. We presented several numerical examples confirming the asymptotic stability of dynamics (25). In Figure 9, we display the variables’ behavior when R 0 < 1 . The trajectory of the system (25) converges to the free-virus periodic state E 0 ( t ) = ( Y s ( t ) , 0 , 0 , 0 , 0 , 0 ) , confirming the global stability of E 0 ( t ) when R 0 1 . In Figure 10, we display the behavior of all components of the solution with several initial conditions, confirming the global stability of E 0 ( t ) if R 0 < 1 . However, in Figure 11, we display the variables’ behavior when R 0 > 1 . The trajectory of the system (25) converges asymptotically to a periodic solution reflecting the virus’ persistence when R 0 > 1 . In Figure 12, we display the behavior of all components of the solution with several initial conditions. The solution converges to a periodic trajectory with persistence of the virus. when R 0 > 1 .

6. Conclusions

In this study, we extended the HIV dynamics model already studied in [12] by considering a variable environment. The model explored three routes of infection including the HIV-to-cell contact, latently infected cell-to-cell contact, and actively infected cell-to-cell contact. In the first stage, we analyzed the case of a fixed environment, such that all parameters are fixed positive constants. In the second stage, we considered the case of a variable environment, where we defined the basic reproduction number R 0 . We proved the global stability of the virus-free periodic trajectory, E 0 ( t ) if R 0 < 1 . However, we proved that the virus will persist if R 0 > 1 and the solution converges to a periodic trajectory. We confirmed the obtained results by several numerical examples, including the fixed environment, the variable contact rates, and the fully variable environment. The trajectories converged to one of the steady states of (2), confirming Theorems 2 and 3, in the case of the fixed environment. Since at least the contact rates are periodic, the solution of (1) will be periodic and converges to a limit cycle confirming Theorems 5 and 6.

Author Contributions

Conceptualization, M.E.H. and R.M.A.; methodology, M.E.H. and R.M.A.; writing—original draft, M.E.H. and R.M.A.; writing—review and editing, M.E.H. and R.M.A. All authors have read and agreed to the published version of the manuscript.

Funding

This work was funded by the University of Jeddah, Jeddah, Saudi Arabia under Grant No. (UJ-23-FR-59).

Data Availability Statement

Data are contained within the article.

Acknowledgments

This work was funded by the University of Jeddah, Jeddah, Saudi Arabia under Grant No. (UJ-23-FR-59). Therefore, the authors thank the University of Jeddah for its technical and financial support. The authors are also grateful to the unknown referees for the many constructive suggestions, which helped to improve the presentation of the paper.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Perelson, A.S.; Neumann, A.U.; Markowitz, M.; Leonard, J.M.; Ho, D.D. HIV-1 dynamics in vivo: Virion clearance rate, infected cell life-span, and viral generation time. Science 1996, 271, 1582–1586. [Google Scholar] [CrossRef] [PubMed]
  2. Nowak, M.A.; Bonhoeffer, S.; Shaw, G.M.; May, R.M. Anti-viral Drug Treatment: Dynamics of Resistance in Free Virus and Infected Cell Populations. J. Theor. Biol. 1997, 184, 203–217. [Google Scholar] [CrossRef] [PubMed]
  3. Perelson, A.S.; Essunger, P.; Cao, Y.; Vesanen, M.; Hurley, A.; Saksela, K.; Markowitz, M. Decay characteristics of HIV-1-infected compartments during combination therapy. Nature 1997, 387, 188–191. [Google Scholar] [CrossRef] [PubMed]
  4. De Boer, R.; Perelson, A. Target Cell Limited and Immune Control Models of HIV Infection: A Comparison. J. Theor. Biol. 1998, 190, 201–214. [Google Scholar] [CrossRef] [PubMed]
  5. Wodarz, D.; Lloyd, A.L.; Jansen, V.A.; Nowak, M.A. Dynamics of Macrophage and T Cell Infection by HIV. J. Theor. Biol. 1999, 196, 101–113. [Google Scholar] [CrossRef] [PubMed]
  6. Bajaria, S.H.; Webb, G.; Kirschner, D.E. Predicting differential responses to structured treatment interruptions during HAART. Bull. Math. Biol. 2004, 66, 1093–1118. [Google Scholar] [CrossRef] [PubMed]
  7. Elaiw, A.; AlShamrani, N. Stability of a general CTL-mediated immunity HIV infection model with silent infected cell-to-cell spread. Adv. Differ. Equ. 2020, 2020, 335. [Google Scholar] [CrossRef]
  8. AlShamrani, N. Stability of a general adaptive immunity HIV infection model with silent infected cell-to-cell spread. Chaos Solitons Fractals 2021, 150, 110422. [Google Scholar] [CrossRef]
  9. Alsahafi, S.; Woodcock, S. Exploring HIV Dynamics and an Optimal Control Strategy. Mathematics 2022, 10, 749. [Google Scholar] [CrossRef]
  10. Stengel, R.F. Mutation and control of the human immunodeficiency virus. Math. Biosci. 2008, 213, 93–102. [Google Scholar] [CrossRef]
  11. Starkov, K.E.; Kanatnikov, A.N. Eradication Conditions of Infected Cell Populations in the 7-Order HIV Model with Viral Mutations and Related Results. Mathematics 2021, 9, 1862. [Google Scholar] [CrossRef]
  12. AlShamrani, N.H.; Halawani, R.H.; Shammakh, W.; Elaiw, A.M. Global Properties of HIV-1 Dynamics Models with CTL Immune Impairment and Latent Cell-to-Cell Spread. Mathematics 2023, 11, 3743. [Google Scholar] [CrossRef]
  13. Xiao, D. Dynamics and bifurcations on a class of population model with seasonal constant-yield harvesting. Discret. Contin. Dyn. Syst. -B 2016, 21, 699–719. [Google Scholar] [CrossRef]
  14. Ibrahim, M.A.; Dénes, A. Stability and Threshold Dynamics in a Seasonal Mathematical Model for Measles Outbreaks with Double-Dose Vaccination. Mathematics 2023, 11, 1791. [Google Scholar] [CrossRef]
  15. Bacaër, N.; Guernaoui, S. The epidemic threshold of vector-borne diseases with seasonality. J. Math. Biol. 2006, 53, 421–436. [Google Scholar] [CrossRef]
  16. Bacaër, N. Approximation of the basic reproduction number R0 for vector-borne diseases with a periodic vector population. Bull. Math. Biol. 2007, 69, 1067–1091. [Google Scholar] [CrossRef]
  17. Ibrahim, M.A.; Dénes, A. A mathematical model for Lassa fever transmission dynamics in a seasonal environment with a view to the 2017–20 epidemic in Nigeria. Nonlinear Anal. Real World Appl. 2021, 60, 103310. [Google Scholar] [CrossRef]
  18. El Hajji, M.; Alshaikh, D.M.; Almuallem, N.A. Periodic behaviour of an epidemic in a seasonal environment with vaccination. Mathematics 2023, 11, 2350. [Google Scholar] [CrossRef]
  19. El Hajji, M.; Alnjrani, R.M. Periodic Trajectories for HIV Dynamics in a Seasonal Environment With a General Incidence Rate. Int. J. Anal. Appl. 2023, 21, 96. [Google Scholar] [CrossRef]
  20. El Hajji, M. Periodic solutions for chikungunya virus dynamics in a seasonal environment with a general incidence rate. AIMS Math. 2023, 8, 24888–24913. [Google Scholar] [CrossRef]
  21. Diekmann, O.; Heesterbeek, J. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J. Math. Biol. 1990, 28, 365–382. [Google Scholar] [CrossRef]
  22. Den Driessche, P.V.; Watmough, J. Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Math. Biosci. 2002, 180, 29–48. [Google Scholar] [CrossRef]
  23. Hurwitz, A. Ueber die Bedingungen, unter welchen eine Gleichung nur Wurzeln mit negativen reellen Theilen besitzt. (English translation “On the conditions under which an equation has only roots with negative real parts” by H. G. Bergmann in Selected Papers on Mathematical Trends in Control Theory R. Bellman and R. Kalaba Eds. New York: Dover, 1964; pp. 70–82). Math. Ann. 1895, 46, 273–284. [Google Scholar]
  24. Routh, E.J. A Treatise on the Stability of a Given State of Motion: Particularly Steady Motion; Macmillan: New York, NY, USA, 1877. [Google Scholar]
  25. El Hajji, M. Mathematical modeling for anaerobic digestion under the influence of leachate recirculation. AIMS Math. 2023, 8, 30287–30312. [Google Scholar] [CrossRef]
  26. LaSalle, J. The Stability of Dynamical Systems; SIAM: Philadelphia, PA, USA, 1976. [Google Scholar]
  27. Alshehri, A.; El Hajji, M. Mathematical study for Zika virus transmission with general incidence rate. AIMS Math. 2022, 7, 7117–7142. [Google Scholar] [CrossRef]
  28. Frobenius, G. Uber Matrizen aus Nicht Negativen Elementen; Sitzungsberichte Preussische Akademie der Wissenschaft: Berlin, Germany, 1912; pp. 456–477. [Google Scholar]
  29. Zhang, F.; Zhao, X. A periodic epidemic model in a patchy environment. J. Math. Anal. Appl. 2007, 325, 496–516. [Google Scholar] [CrossRef]
  30. Wang, W.; Zhao, X. Threshold dynamics for compartmental epidemic models in periodic environments. Dynam. Differ. Equ. 2008, 20, 699–717. [Google Scholar] [CrossRef]
  31. Zhao, X. Dynamical Systems in Population Biology; CMS Books in Mathematics (CMSBM); Springer: New York, NY, USA, 2003; Volume 16. [Google Scholar]
  32. Ma, J.; Ma, Z. Epidemic threshold conditions for seasonally forced SEIR models. Math. Biosci. Eng. 2006, 3, 161–172. [Google Scholar] [CrossRef]
  33. Zhang, T.; Teng, Z. On a nonautonomous SEIRS model in epidemiology. Bull. Math. Biol. 2007, 69, 2537–2559. [Google Scholar] [CrossRef]
  34. Guerrero-Flores, S.; Osuna, O.; de Leon, C.V. Periodic solutions for seasonal SIQRS models with nonlinear infection terms. Electron. J. Differ. Equations 2019, 2019, 1–13. [Google Scholar]
  35. Nakata, Y.; Kuniya, T. Global dynamics of a class of SEIRS epidemic models in a periodic environment. J. Math. Anal. Appl. 2010, 363, 230–237. [Google Scholar] [CrossRef]
Figure 1. Behavior of dynamics (23) for Λ i n 0 = 0.2 then R 0 0.6364 < 1 . The red star describe the limit (final state) of the variable ( Y l , Y i , Y v ) .
Figure 1. Behavior of dynamics (23) for Λ i n 0 = 0.2 then R 0 0.6364 < 1 . The red star describe the limit (final state) of the variable ( Y l , Y i , Y v ) .
Mathematics 12 00123 g001
Figure 2. Behavior of dynamics (23) for several initial conditions where Λ i n 0 = 0.2 ( R 0 0.6364 < 1 ). The different graphs (colors) reflect different initial conditions.
Figure 2. Behavior of dynamics (23) for several initial conditions where Λ i n 0 = 0.2 ( R 0 0.6364 < 1 ). The different graphs (colors) reflect different initial conditions.
Mathematics 12 00123 g002
Figure 3. Behavior of dynamics (23) for Λ i n 0 = 2 then R 0 6.3642 > 1 . The red star describe the limit (final state) of the variable ( Y l , Y i , Y v ) .
Figure 3. Behavior of dynamics (23) for Λ i n 0 = 2 then R 0 6.3642 > 1 . The red star describe the limit (final state) of the variable ( Y l , Y i , Y v ) .
Mathematics 12 00123 g003
Figure 4. Behavior of dynamics (23) for several initial conditions where Λ i n 0 = 2 ( R 0 6.3642 > 1 ). The different graphs (colors) reflect different initial conditions.
Figure 4. Behavior of dynamics (23) for several initial conditions where Λ i n 0 = 2 ( R 0 6.3642 > 1 ). The different graphs (colors) reflect different initial conditions.
Mathematics 12 00123 g004
Figure 5. Behavior of dynamics (24) for Λ i n 0 = 0.2 then R 0 0.6364 < 1 . The red star describe the limit (final state) of the variable ( Y l , Y i , Y v ) .
Figure 5. Behavior of dynamics (24) for Λ i n 0 = 0.2 then R 0 0.6364 < 1 . The red star describe the limit (final state) of the variable ( Y l , Y i , Y v ) .
Mathematics 12 00123 g005
Figure 6. Behavior of dynamics (24) for several initial conditions where Λ i n 0 = 0.2 ( R 0 0.6364 < 1 ). The different graphs (colors) reflect different initial conditions.
Figure 6. Behavior of dynamics (24) for several initial conditions where Λ i n 0 = 0.2 ( R 0 0.6364 < 1 ). The different graphs (colors) reflect different initial conditions.
Mathematics 12 00123 g006
Figure 7. Behavior of dynamics (24) for Λ i n 0 = 2 then R 0 6.3642 > 1 . The red curve describe the limit-cycle of the variable ( Y l , Y i , Y v ) .
Figure 7. Behavior of dynamics (24) for Λ i n 0 = 2 then R 0 6.3642 > 1 . The red curve describe the limit-cycle of the variable ( Y l , Y i , Y v ) .
Mathematics 12 00123 g007
Figure 8. Behavior of dynamics (24) for several initial conditions where Λ i n 0 = 2 ( R 0 6.3642 > 1 ). The different graphs (colors) reflect different initial conditions.
Figure 8. Behavior of dynamics (24) for several initial conditions where Λ i n 0 = 2 ( R 0 6.3642 > 1 ). The different graphs (colors) reflect different initial conditions.
Mathematics 12 00123 g008
Figure 9. Behavior of dynamics (25) for Λ i n 0 = 0.2 then R 0 0.6364 < 1 . The red star describe the limit (final state) of the variable ( Y l , Y i , Y v ) .
Figure 9. Behavior of dynamics (25) for Λ i n 0 = 0.2 then R 0 0.6364 < 1 . The red star describe the limit (final state) of the variable ( Y l , Y i , Y v ) .
Mathematics 12 00123 g009
Figure 10. Behavior of dynamics (25) for several initial conditions where Λ i n 0 = 0.2 ( R 0 0.6364 < 1 ). The different graphs (colors) reflect different initial conditions.
Figure 10. Behavior of dynamics (25) for several initial conditions where Λ i n 0 = 0.2 ( R 0 0.6364 < 1 ). The different graphs (colors) reflect different initial conditions.
Mathematics 12 00123 g010
Figure 11. Behavior of dynamics (25) for Λ i n 0 = 2 ( R 0 6.3642 > 1 ). The red curve describe the limit-cycle of the variable ( Y l , Y i , Y v ) .
Figure 11. Behavior of dynamics (25) for Λ i n 0 = 2 ( R 0 6.3642 > 1 ). The red curve describe the limit-cycle of the variable ( Y l , Y i , Y v ) .
Mathematics 12 00123 g011
Figure 12. Behavior of dynamics (25) for several initial conditions where Λ i n 0 = 2 ( R 0 6.3642 > 1 ). The different graphs (colors) reflect different initial conditions.
Figure 12. Behavior of dynamics (25) for several initial conditions where Λ i n 0 = 2 ( R 0 6.3642 > 1 ). The different graphs (colors) reflect different initial conditions.
Mathematics 12 00123 g012
Table 1. Variables and parameters.
Table 1. Variables and parameters.
VariableDescription
Y s Susceptible cells
Y l Latently infected cells
Y i Productively infected cells
Y v Free virions
Y c B cells and cytotoxic T lymphocytes
ParameterDescription
Λ i n Generation rate of Healthy (susceptible) cells Y s
β 1 Periodic contact rate between Y s and Y v
β 2 Periodic contact rate between Y s and Y l
β 3 Periodic contact rate between Y s and Y i
d s Mortality rate of Healthy (susceptible) cells
d l Mortality rate of latently infected cells
d i Mortality rate of productively infected cells
d v Mortality rate of free virions (HIV-1 particles)
d c Mortality rate of B cells and cytotoxic T lymphocytes (CTLs)
η Conversion rate of Y l cells into Y i cells
ε Generated rate of HIV from Y i cells
λ B-cell immune rate produced by Y i cells
γ Y c Y i Neutralization rate
δ Y c Y i Impairment rate
Table 2. Constants used in the numerical examples.
Table 2. Constants used in the numerical examples.
Λ i n 0 d s 0 d l 0 d i 0 d v 0 d c 0 β 10 β 20 β 30 η 0 γ 0 δ 0 ε 0 λ 0
10 0.8 0.7 2 0.5 1 0.2 0.8 42 0.5 1 0.2 0.8
Λ i n 1 d s 1 d l 1 d i 1 d v 1 d c 1 β 11 β 21 β 31 η 1 γ 1 δ 1 ε 1 λ 1 ψ
10 0.8 0.7 2 0.5 1 0.2 0.8 42 0.5 1 0.2 0.8 4
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

El Hajji, M.; Alnjrani, R.M. Periodic Behaviour of HIV Dynamics with Three Infection Routes. Mathematics 2024, 12, 123. https://doi.org/10.3390/math12010123

AMA Style

El Hajji M, Alnjrani RM. Periodic Behaviour of HIV Dynamics with Three Infection Routes. Mathematics. 2024; 12(1):123. https://doi.org/10.3390/math12010123

Chicago/Turabian Style

El Hajji, Miled, and Rahmah Mohammed Alnjrani. 2024. "Periodic Behaviour of HIV Dynamics with Three Infection Routes" Mathematics 12, no. 1: 123. https://doi.org/10.3390/math12010123

APA Style

El Hajji, M., & Alnjrani, R. M. (2024). Periodic Behaviour of HIV Dynamics with Three Infection Routes. Mathematics, 12(1), 123. https://doi.org/10.3390/math12010123

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop