Periodic Trajectories for HIV Dynamics in a Seasonal Environment With a General Incidence Rate

. In this paper, we propose a ﬁve-dimensional nonlinear system of diﬀerential equations for the human immunodeﬁciency virus (HIV) including the B-cell functions with a general nonlinear incidence rate. The compartment of infected cells was subdivided into three classes representing the latently infected cells, the short-lived productively infected cells, and the long-lived productively infected cells. The basic reproduction number was established, and the local and global stability of the equilibria of the model were studied. A sensitivity analysis with respect to the model parameters was undertaken. Finally, some numerical simulations are presented to illustrate the theoretical ﬁndings.


Introduction
Human immunodeficiency virus (HIV) is a type of virus that attacks the body's immune system [1].
Although HIV infection is a manageable chronic condition, if left untreated, it can weaken the immune system or progress to acquired immune deficiency syndrome (AIDS). Some people may have no symptoms after contracting HIV so the infection is not diagnosed until symptoms of AIDS appear. It can take up to 10 years. Symptoms of an HIV infection can last for a few days or weeks. They can go away on their own. It is common for HIV infection to be misdiagnosed as another illness first.
Historically, mathematics (in this context, we refer in particular to mathematical modeling and analysis) has been used to better understand the dynamics of the transmission of infectious diseases and to learn how to control them. This application of mathematics dates back to the work of Daniel Bernoulli, who used mathematical and statistical methods to study the potential impact of the smallpox vaccine in 1760 [2]. In the 1920s, Sir Ronald Ross, a physician by training, used mathematical modeling to propose effective methods of malaria control. In particular, he showed that the disease can be eradicated if the mosquito population is kept below a certain threshold, a discovery that won him the Nobel Prize in medicine. More recently, mathematics has helped shape effective public health policies against the spread of emerging and re-emerging diseases that pose a significant threat to public health, such as HIV/AIDS, influenza (e.g., the recent pandemics of bird flu and swine flu), malaria, severe acute respiratory syndrome (SARS) and tuberculosis. The mathematical modelling in epidemiology is a way to study how a disease is spread, predict the future behaviour, and propose control strategies.
However, seasonality is very repetitive in each of the ecological, biological, and human systems [14].
In particular, in the climate variation patterns repeated every year by the same way, bird migration is repeated according to the repeated season variation, schools open and close almost periodically each year, etc. Among other things, these seasonal factors affect the pathogens' survival in the environment, host behaviour, and the abundance of vectors and non-human hosts. Therefore, several diseases show seasonal behaviours. Taking into account the seasonality in mathematical modelling becomes very important. Note that even the simplest mathematical models that take into account seasonality present many difficulties to study [5]. In [15], Bacaër and Gomes discussed the periodic S-I-R model, a simple generalization of the classical model of Kermack and McKendrick [16]. In [17], the authors studied a SEIRS epidemic model with periodic fluctuations. They calculated the basic reproduction number R 0 of the time-averaged system (autonomous). Then, they proved a sufficient but not necessary condition (R 0 < 1) such that the disease could not persist in the population in a seasonal environment. In [18], Guerrero-Flores et al. considered a class of SIQRS models with periodic variations in the contact rate. They proved the existence of periodic orbits by using Leray-Schauder degree theory. Zhang and Teng [19] studied an alternative SEIRS epidemic model in a seasonal environment and established some sufficient equivalent conditions for the persistence and the extinction of the disease. These results were improved by Nakata and Kuniya in [3] by giving a threshold value between the uniform persistence and the extinction of the disease. In [6], Bacaër and Guernaoui gave the definition of the basic reproduction number in seasonal environments. In 2008, Wang and Zhao [20] defined R 0 for several compartmental epidemic models in seasonal environments. All these definitions were different, in several cases, from the basic reproduction number defined for the timeaveraged system. By considering general compartmental epidemic models in seasonal environments, Wang and Zhao [20] showed that R 0 was the threshold value for proving or not the local stability of the disease-free periodic trajectory. In [21], the authors studied the periodic behaviour of an "SVEIR" epidemic in a seasonal environment with vaccination.
As seasonality is very repetitive in the environment, which affects several diseases that show seasonal behaviours, taking seasonality into account in mathematical modelling becomes a necessity. In this paper, we proposed an extension of the HIV model proposed in [10,11] by reducing the system from six-dimensional to five-dimensional system and by taking into account the seasonal environment. The mathematical model includes the B-cell functions for HIV dynamics with a general nonlinear incidence rate. The infected compartment was subdivided into two classes, namely the latently infected class and the productively infected class. We studied, in a first step, the autonomous system by investigating the global stability of the steady states. In a second step, we showed that the disease-free periodic solution is globally asymptotically stable if R 0 is less than 1, and, if R 0 is greater than 1, the disease persists. The rest of the paper is organized as follows. In Section 2, we introduced the mathematical model. In Section 3, we studied the case of an autonomous system, where all parameters are supposed to be constants. In Section 4, we considered the non-autonomous system, gave some basic results, and gave the definition of R 0 . We showed that the value of R 0 around one was a threshold value between the disease's extinction and the disease's uniform persistence. We gave numerical examples that supported the theoretical findings in Section 5. Section 6 provided a brief conclusions of our obtained results.

Mathematical Model for HIV Dynamics
In this paper, we generalized the mathematical model studied in [10,11]. The mathematical model is compartmental model since it describes the transfer of molecules through different compartments of the body. Let X u , X l , X i , X p , and X c to be the number of uninfected cells, latently infected cells, productively infected cells, free virions, and B cells, respectively. The change in the number or concentration of CD4 + T lymphocytes (dX u ) over a small time interval (dt) is a function of a constant (rate of de novo cell production ), the rate of cell death (d u ) proportional to the number of cells present (X u ), and the number of infected cells that leave this "compartment" to join the compartment of infected cells (X l , X i ). The number of infected cells is distributed into two compartments depending on the type of infected cell (X l , X i ).
The number of infected cells is, therefore, a function of the number of target cells, the number of circulating virions (X p ), and the incidence rate of infection (f 1 (X p )X u ), also called infectivity. The incidence rate of infection is very important for understanding the dynamics of the system. The variation in the number of virions (dX p ) per unit of time depends on the number of virions produced (d i r 1 X i ), the clearance d p of the virus, and the neutralized part of the HIV particles, f 2 (X p )X c .
The B-cell immune response is assumed to be proportional to the free virions' population (εX p ). The B cell impairment rate is assumed to be proportional to the contact with the free virions' population (υf 2 (X p )X c , where υ is a positive constant).
(p 1 + p 2 )f 1 (X p )X u is the incidence rate of infection and f 2 (X p )X c is the neutralization rate of HIV particles. Note that the incidence rate (f 1 ) and the neutralization rate (f 2 ), increase once the free viruses increase and are neutral in the absence of the virus. Thus, the functions f 1 and f 2 satisfy the following assumption.
The model's parameters are positive and are given hereafter as in Table 1.

Parameter
Description Incidence rate between X p and X l p 2 Incidence rate between X p and X i   ( (2) Proof.

Autonomous system
Consider the case where all parameters of dynamics (2.1) are constant and the system becomes autonomous and takes the following form 3.1. Basic Results. It is necessary that the state variables X u (t), X l (t), X i (t), X p (t) and X c (t) remain non-negative for all t ≥ 0.
is a positively invariant attractor set for the dynamics (3.1).
Proof. R 5 + is positively invariant set for (3.1), since we havė In order to prove that the solution is bounded, let us define T 1 (t) = X u (t) + X l (t) + X i (t) − d 1 and Similarly, we havė then, and if then, therefore, Γ is positively invariant for system (3.1).

Basic reproduction number and steady states.
Let R 0 to be the basic reproduction number, describing the average number of new cases of a disease that a single infected and contagious person will generate on average in a susceptible population. Diekmann et al. [22] were the first to propose the next-generation matrix method to calculate R 0 . Late, Van den Driessche and Watmough [23] elaborated this method.
Then, the next-generation matrix is given by and its characteristic polynomial is given by: Therefore, the spectral radius representing the basic reproduction number is: Lemma 3.2.
• If R 0 ≤ 1, then the model (3.1) admits a unique equilibrium point given by Proof. Let (X u , X l , X i , X p , X c ) be any equilibrium point of the model (3.1) satisfying: By solving equation (3.3), we obtain a steady state given by the HIV-free steady state E 0 = ( d u , 0, 0, 0, 0). Moreover, we have: .
Using the third equation of the dynamics (3.3), we obtain Therefore, X p = 0 is a solution and then X u = d u and X l = X i = X p = X c = 0, called the HIV-free . Furthermore, assume that X p = 0, and divide the equation (3.5) by X p , we obtain Consider the function g defined by One can easy obtain Since Furthermore, the derivative of g is given by By Lemma 2.1, the first term of g (X p ) is negative and since the second and the third terms are negative then g (X p ) < 0, and thus, g is a decreasing function. Then, the existence and uniqueness Then, the persistence equilibrium point

Local Stability.
In this section, we aim to investigate the local stability of the steady states of system (3.1) using the linearization approach. Recall that the value of R 0 with respect to the unit is important in concluding if the disease persists or not. Hereafter, we give the results concerning the local stability of the equilibrium points E 0 and E * .
Theorem 3.1. If R 0 < 1, then the steady state E 0 is locally asymptotically stable.
Proof. The Jacobian matrix of the linear approximation of system (3.1) at the trivial steady state E 0 is: J 0 admits five eigenvalues. The first two eigenvalues are given by The other three eigenvalues are the roots of: It is easy to see that for R 0 < 1, we have Then, by using the Routh-Hurwitz criteria [24,25], the eigenvalues have negative real parts. Then, the trivial steady state E 0 is locally asymptotically stable once R 0 < 1; however, it is a saddle point Theorem 3.2. The infected steady state E * is locally asymptotically stable once R 0 > 1.
Proof. The value of the Jacobian matrix at the infected equilibrium point E * is: The eigenvalues of J * are the roots of: The characteristic polynomial X * p (λ) can be written in the form X * p (λ) = λ 5 + a 4 λ 4 + a 3 λ 3 + a 2 λ 2 + a 1 λ + a 0 . Then, by using the Maple software, we can prove that for R 0 < 1, we have Therefore, the roots of X * p (eigenvalues) have negative real parts by the Routh-Hurwitz criteria [24,25]. The steady state, E * , is then locally asymptotically stable if R 0 > 1.

Global Stability.
Define G to be the function G(z) = z − 1 − ln z and the constant: Theorem 3.3. The trivial equilibrium point E 0 is globally asymptotically stable once R 0 ≤ 1.
Proof. Assume that R 0 ≤ 1, and let define the Lyapunov function L 0 (X u , X l , X i , X p , X c ) by The derivative of L 0 along Model (3.1) is: It can be easily shown that W 0 = {E 0 }. Using LaSalle's invariance principle [26] (see [8,11,12,27] for some examples), one deduces that E 0 is globally asymptotically stable once R 0 ≤ 1.
Proof. Let a function L * (X u , X l , X i , X p , X c ) be defined as: Clearly, L * (X u , X l , X i , X p , X c ) > 0 for all X u , X l , X i , X p , X c > 0 and L * (X * u , X * l , X * i , X * p , X * c ) = 0. Calculating dL * dt along the trajectories of (3.1) and using the fact that = d u X * u + (p 1 +p 2 )f (X * p )X * u , we obtain: Now, since: .
Based on the rule: 10) and the Lemma 2.1, we obtain dL * dt (X u , X l , X i , X p , X c ) ≤ 0 for all X u , X l , X i , X p , X c > 0 and dL * dt (X u , X l , X i , X p , X c ) = 0 if and only if (X u , X l , X i , X p , X c ) = (X * u , X * l , X * i , X * p , X * c ). From LaSalle's invariance principle [26], we deduce the global stability of E * (see [11,28,29] for other applications).

Periodic System
Return now to the main model (2.1) where all parameters are continuous and T -periodic positive functions: to be the ordered m-dimensional Euclidean space associated with norm · . For X 1 , X 2 ∈ R m , we denote by . Consider a T -periodic m × m matrix function denoted by C(t) which is continuous, irreducible and cooperative. Let us denote by φ C (t) the fundamental matrix, solution of the following systemẋ (t) = C(t)x(t). (4.2) Let denote the spectral radius of the matrix φ C (T ) by r (φ C (T )). Therefore, all entries of φ C (t) are positive for each t > 0. Let apply the theorem of Perron-Frobenius to deduce that r (φ C (T )) is the principal eigenvalue of φ C (T ) (simple and admits an eigenvector y * 0). For the rest of the paper, the following lemma will be useful.
Lemma 4.1. [30]. There exists a positive T -periodic function y (t) such that x(t) = y (t)e kt will be a solution of system (4.2) where k = 1 T ln(r (φ C (T ))).
Let start by proving the existence (and uniqueness) of the disease free periodic trajectory of model (4.1). Let consider the following equatioṅ with initial condition X 0 u ∈ R+. (4.3) admits a unique T -periodic solution X * u (t) with X * u (t) > 0 which is globally attractive in R+ and hence, system (4.1) has a unique disease free periodic solution (0, 0, 0, 0, X * u (t)). For a continuous, positive T -periodic function g(t), we set g u = max t∈[0,T ) g(t), 2ε u X c ≤S 2 } is a positively invariant attractor set for system (4.1). Furthermore, we have Proof. From system (4.1), we havė This means that Ω u is a forward invariant compact absorbing set of all solutions of system (4.1). Let , and this means that lim t→∞ y (t) = lim t→∞ (X u (t)+X l (t)+X i (t)−X * u (t)) = 0.
Next, in subsection 4.1, we define R 0 , the basic reproduction number and we will prove that the disease free periodic trajectory (0, 0, 0, 0, X * u (t)) is globally asymptotically stable (and therefore, the disease dies out) once R 0 < 1. Then, in subsection 4.2, we will prove that X i (t) is uniform persistence (and then the disease persists) once R 0 > 1. Therefore, we deduce that R 0 is the threshold parameter between the uniform persistence and the extinction of the disease.

Disease Free Periodic Solution .
We start by giving the definition of the basic reproduction number of model (4.1), by using the theory given in [20] where F(t, X) =  Our aim is to check the conditions (A1)-(A7) in [20, Section 1]. Note that system(4.1) can have the following formẊ = F(t, X) − V(t, X) = F(t, X) − V − (t, X) + V + (t, X). (4.7) The first five conditions (A1)-(A5) are fulfilled.
The system (4.7) admits a disease free periodic trajectory X * (t) = where f i (t, X(t)) and X i are the i -th component of f (t, X(t)) and X, respectively. By an easy calculus, we get M(t) = −(p 1 (t) + p 2 (t))f 1 (X p (t)) − d u (t) and then r (φ M (T )) < 1. Therefore X * (t) is linearly asymptotically stable in the subspace Γ s = (0, 0, 0, 0, X u ) ∈ R 5 + . Thus, the condition (A6) in [20, Section 1] is satisfied. Now, let us define F(t) and V(t) to be four by four matrices given by and where F i (t, X) and V i (t, X) are the i -th component of F(t, X) and V(t, X), respectively. By an easy calculus, we obtain from system (4.7) Consider Z(t 1 , t 2 ) to be the two by two matrix solution of the system d dt for any t 1 ≥ t 2 , with Z(t 1 , t 1 ) = I, the two by two identity matrix. Thus, condition (A7) was satisfied.
Let define C T to be the ordered Banach space of T -periodic functions defined on R → R 2 , associated to the maximum norm . ∞ and the positive cone C + T = {ψ ∈ C T : ψ(s) ≥ 0, for any s ∈ R}. Define the linear operator K : Let now define the basic reproduction number, R 0 , of model (4.1) by R 0 = r (K).
Consider the system hereafter system Then, we deduce that the disease free periodic solution E 0 (t) is globally attractive which complete the proof.
For the following subsection, we consider only the case where R 0 > 1.
Let us define the function P : R 5 + → R 5 + to be the Poincaré map associated to system (4.1) such that X 0 → u(T, X 0 ), where u(t, X 0 ) is the unique solution of the system (4.1) with the initial condition u(0, X 0 ) = X 0 ∈ R 5 + . Let us define (4.15) The fixed pointX uα0 of the function P associated to the perturbed system (4.13) is globally attractive such thatX uα (t) >X u (t) − η, then ∃ T 2 > 0 large enough and satisfying since the trajectories are bounded. Therefore, the inequality (4.14) is satisfied and P is weakly uniformly persistent respecting to (Γ 0 , ∂Γ 0 ). By applying Lemma 4.2, P has a global attractor. We We prove also by contradiction thatX 0 u > 0. Assume thatX 0 u = 0. Using the first equation of the reduced system (4.1),X u (t) verifieṡ withX 0 u =X u (mT ) = 0, m = 1, 2, 3, · · · . Applying Lemma 4.2, ∀ δ 3 > 0, there exists T 3 > 0 large enough and satisfyingX p (t) ≤S 2 + δ 3 , t > T 3 . Then, we obtaiṅ There existsm large enough and satisfying mT > T 3 for all m >m. Applying the comparison principle, we deducẽ for any m >m which is impossible. Therefore,X 0 u > 0 and (X 0 u ,X 0 l ,X 0 i ,X 0 p ,X 0 c ) is a positive T -periodic trajectory of the reduced system (4.1).
c , ν 1 , υ 1 , r 1 1 and ε 1 measure the amplitude(< 1) of the seasonal variation in each of the parameters. φ is the phase shift. Some fixed constants used for the numerical simulations are given in Table 2.
We will consider three cases. The first case is dedicated for the case of constant parameters (autonomous system) to validate the obtained theoretical results concerning the local and global stability of the equilibrium points E 0 and E * . The second case deals only with a seasonal contact (p 1 and p 2 are periodic functions) however the other parameters are constants (partially non-autonomous system).
The third case considers all parameters as periodic functions (non-autonomous system).

5.1.
The case of the autonomous system. In a first step, we consider that all parameters of the system (2.1) are constants (ρ 1 = p 1 Thus, the model is given by with positive initial condition (X 0 u , X 0 l , X 0 i , X 0 p , X 0 c ) ∈ R 5 + . We give some numerical results that confirm the stability of the equilibrium points of (5.2). In Fig. 1, we give the results for the case where R 0 > 1. The approximated solution of the given model (5.2) approaches asymptotically to E * , which confirms that E * is globally asymptotically stable once R 0 > 1. In Fig. 2, we give the results for the case where R 0 < 1. The approximated solution of the given model (5.2) approaches the equilibrium E 0 , which confirms that E 0 is globally asymptotically stable once R 0 ≤ 1.

5.2.
The case of the partially non-autonomous system. In a second step, we performed numerical simulations on the system (2.1) using linear function to express the transmission rate where only the seasonally forced T -periodic function p 1 (t) and p 2 (t) are depending on time, t. The other parameters are constants.
Thus the model is given by with positive initial condition (X 0 u , X 0 l , X 0 i , X 0 p , X 0 c ) ∈ R 5 + where β 1 = 0.8. The basic reproduction number, R 0 , was approximated using the time-averaged system. We give some numerical results that confirm the asymptotic behaviour of the solution of (5.3). In Fig. 3, we give the results for the case where R 0 > 1. The approximated solution of the given model (5.3) approaches asymptotically to a periodic solution with persistence of the disease. In Fig. 4, we give the results for the case where

5.3.
The case of totally non-autonomous system. In a third step, we performed numerical simulations on the system (2.1) using classical Monod function to express the transmission rate where all parameters are T -periodic functions. Thus the model is given by X l (t) = p 1 (t)f 1 (X p (t))X u (t) − (d l (t) + ν(t))X l (t), X i (t) = p 2 (t)f 1 (X p (t))X u (t) + ν(t)X l (t) − d i (t)X i (t), X p (t) = d i (t)r 1 (t)X i (t) − d p (t)X p (t) − f 2 (X p (t))X c (t), X c (t) = ε(t)X p (t) − d c (t)X c (t) − υ(t)f 2 (X p (t))X c (t), , with positive initial condition (X 0 u , X 0 l , X 0 i , X 0 p , X 0 c ) ∈ R 5 + . Additional constants used for the numerical simulations in this step are given in Table 3. The basic reproduction number, R 0 , was approximated using the time-averaged system.  Table 3. Additional parameters for numerical simulations of the totally nonautonomous system.
We give some numerical results that confirm the asymptotic behaviour of the solution of (5.4). In In Fig. 6, we give the results for the case where R 0 < 1. The approximated solution of the given model (5.4) approaches the disease-free periodic trajectory E 0 (t) = (X * u (t), 0, 0, 0, 0) once R 0 ≤ 1.

Conclusions
The human immunodeficiency virus (HIV) is the pathogen responsible for acquired immunodeficiency syndrome (AIDS). In this work, we proposed an extension of the HIV epidemic model given in [10,11] in a seasonal environment. In the first step we studied the case of autonomous system where all parameters are supposed to be constants. In the second step, we considered the non-autonomous system and we give some basic results and we defined the basic reproduction number, R 0 . We show that R 0 -value compared to the unit is the threshold value between uniform persistence and extinction of the considered disease. More precisely, we showed that if R 0 is less than 1, then the disease free periodic solution is globally asymptotically stable and if R 0 is greater than 1, then the disease persists. Finally, we gave some numerical examples that supports the theoretical findings, including the autonomous system, the partially non-autonomous system and the full non-autonomous system.
It is deduced that if the system is autonomous, the trajectories converge to one of the equilibrium of the system (2.1) according to theorems 3.3 and 3.4. However, if at least one of the model parameters is periodic, the trajectories converge to a limit cycle according to theorems 4.2 and 4.3.

Conflicts of Interest:
The authors declare that there are no conflicts of interest regarding the publication of this paper.