Analysis of a time-delayed HIV/AIDS epidemic model with education campaigns
AAnalysis of a time-delayed HIV/AIDS epidemic model witheducation campaigns.
Dawit Denu ∗ Sedar Ngoma † Rachidi B. Salako ‡ Abstract
We consider a time-delayed HIV/AIDS epidemic model with education dissemination andstudy the asymptotic dynamics of solutions as well as the asymptotic behavior of the endemicequilibrium with respect to the amount of information disseminated about the disease. Underappropriate assumptions on the infection rates, we show that if the basic reproduction numberis less than or equal to one, then the disease will be eradicated in the long run and any solutionto the Cauchy problem converges to the unique disease-free equilibrium of the model. On theother hand, when the basic reproduction number is greater than one, we prove that the diseasewill be permanent but its impact on the population can be significantly minimized as theamount of education dissemination increases. In particular, under appropriate hypothesis onthe model parameters, we establish that the size of the component of the infected populationof the endemic equilibrium decreases linearly as a function of the amount of informationdisseminated. We also fit our model to a set of data on HIV/AIDS from Uganda within theperiod 1992-2005 in order to estimate the infection, effective response, and information ratesof the disease. We then use these estimates to present numerical simulations to illustrate ourtheoretical findings.
Key words.
Delay differential equations, Epidemic model, Asymptotic behavior, Bifurcation,Numerical approximations.
The human immunodeficiency virus (HIV) is the virus which is responsible for causing acquiredimmunodeficiency syndrome (AIDS). The virus cannot replicate by its own and thus in orderto reproduce it will hijack cells (CD4) that play a central role in responding to the invasion ofinfections in the body. The virus can spread through some body fluids and once it is inside thebody it can destroy the immune system, making the person too weak to fight other infections.The elevated viral load of someone recently infected with HIV is the main biological reason thatthey are more likely to transmit HIV to others. The higher the viral load, the greater the riskis of transmitting HIV to others. When the virus has destroyed a certain number of CD4 cellsand the CD4 count drops below 200 cells per mm , a person will have progressed to AIDS. Untilnow there is no vaccine to prevent HIV infection nor a cure for AIDS and if HIV is untreated ∗ Department of Mathematics, Georgia Southern University, Savannah, GA 31419 ([email protected]) † Department of Mathematics, State University of New York at Geneseo, NY 14454 ([email protected]) ‡ Department of Mathematics, Ohio State University, Columbus, OH 43210 ([email protected]) a r X i v : . [ m a t h . D S ] J a n arlier it will lead to AIDS in about 8 to 10 years. However, there are some medications andantiviral treatments to slow down the progression of the disease [1, 2, 21]. According to UNAIDS,in 2018 there were approximately 38 million people with HIV/AIDS worldwide. Out of thesemany people the majority (95 . S AB group. Similarly, the interaction of thesusceptible group with the information will result in changing the behavior of some of the peopleto use condom and as a result they will join the S C group.2ne way to better capture the dynamics of HIV starting from the time of cell’s infection tothe time it produces new virions is to introduce a latency (or delay) time u >
0. Time delay canarise for various practical reasons in epidemiology. For example, pharmacological delay occursbetween the ingestion of drug and its appearance within cells, while intracellular delay is the timebetween initial infection of a cell by HIV and the release of new virions. In this perspective, thereare several HIV epidemic models that have been developed to study the effect of the delays on theoverall dynamics of the infection [4,7,14,16,18,22]. In particular, Z. Mukandavire et al. presenteda mathematical model for HIV/AIDS with explicit incubation period as a system of discrete timedelay differential equations [16]. They analyzed the stability of equilibrium points using theLyapunov functional approach. Moreover, they investigated the positivity and boundedness ofthe solutions, and provided conditions for the permanence of the infection.Thus, motivated by the above discussion, due to the advantage of the positive effect of infor-mation and education campaigns to control the transmission of HIV/AIDS and the fact that delaymodels provide a more realistic dynamic of the infection, we propose and analyze a time-delayedHIV/AIDS model with the presence of information and education campaigns.
In this paper we consider a population of total size N ( t ) at time t , with a constant recruitmentrate of B . The population size N ( t ) is divided into five compartments S ( t ), S ( t ), S ( t ), I ( t ),and R ( t ). The S -class, also called the general susceptible class, denotes the sub-group of thepopulation which are not yet infected but are susceptible to contract the disease upon interactionwith infected people, and they do not have any immunity with respect to the virus. Furthermore,the S -class individuals have direct access to information and education campaigns about thedisease, which includes practicing abstinence and faithfulness and the use of condom. We denoteby Z ( t ) the amount of the educational information from the information and education campaigns(IEC). Thanks to the access to education campaigns, part of the S population will change theirbehavior and join other compartments according to the type of information they receive andchoose to practice. The group S ( t ) denotes those individuals who change their behavior andpractice abstinence and faithfulness as a result of the information and education campaigns.Similarly, S ( t ) represents the group that changes their behavior and started practicing safe sexwith the use of condoms. There are some individuals who are not only exposed to the informationand education campaigns, but also used the provided information to the fullest that they willnever get infected by HIV throughout their life. These individuals are grouped under the removedcompartment R ( t ). It is important to note that this change in behavior for the R ( t ) class isindependent of what the susceptible people in groups S and S practice. For example, it couldbe having fewer sexual partners, and never sharing needles, etc. We assume that the interactionbetween S and Z will have rates of γ to go to the R group, γ to go to the S group, and γ togo to the S class.The density of information Z ( t ) is assumed to be proportional to the number of regional andglobal organizations providing information and educational campaigns on HIV prevention, andthe number of such organizations is in general proportional to the number of infections in thepopulation. We suppose that there is intracellular constant delay time u ≥ I ( t ) are denoted by β ,u , β ,u , and β ,u and cor-respond respectively, to the groups S , S and S . In addition to the natural death rates µ i ,3 = S , S , S , I, R of each group, we denote by d > S (t) S (t)S (t) I(t)Z(t)R B Figure 1: A flow diagram describing the dynamics of HIV/AIDS epidemic model where S ( t ) , S ( t ) , and S ( t ) are the different susceptible groups of the human population, I ( t ) and R ( t )are the number of infected and removed groups, respectively. Also, Z ( t ) is the information den-sity at time t , and the broken blue arrow indicates the interaction of the S ( t ) group with theinformation Z ( t ) to produce S ( t ) and S ( t ). The solid arrows represent transitions between thedifferent epidemiological classes, whereas the red broken arrows represent the death rates of thevarious groups in the model.Table 1 provides a summary description of the model parameters. The values of the parametersin Table 1 together with the parameters to be estimated are provided in Section 5.Table 1: A list of model parameters and their interpretationsSymbols Description u Intracellular time delay γ , γ , γ Effective response rate of S ( t ) , S ( t ) , S ( t ) β ,u , β ,u , β ,u Infection transmission rates of S ( t ) , S ( t ) , S ( t ) µ S , µ S , µ S , µ I , µ R Natural death rates of S ( t ) , S ( t ) , S ( t ) , I ( t ) , R ( t ) d Death rate due to the disease q Rate of increase of information w.r.t. I ( t ) η Growth rate of information B Recruitment rate or rate of entry per person into S ( t )Taking into account the above considerations, the model dynamics is assumed to be governed4y the following system of nonlinear delay differential equations ddt S ( t ) = B − (cid:80) j =0 γ j S ( t ) Z ( t ) − β ,u S ( t ) I ( t − u ) − µ S S ( t ) , ddt S j ( t ) = γ j S ( t ) Z ( t ) − β j,u S j ( t ) I ( t − u ) − µ S j S j ( t ) , j = 1 , , ddt Z ( t ) = qI ( t ) − ηZ ( t ) , ddt I ( t ) = (cid:80) j =0 β j,u S j ( t ) I ( t − u ) − µ I I ( t ) − dI ( t ) , ddt R ( t ) = γ S ( t ) Z ( t ) − µ R R ( t ) . (1.1)In the remaining part of this paper, the natural death rate of all the susceptible, infected,and removed groups is assumed to be equal and denoted by µ . That is, µ S j = µ I = µ R = µ , j = 0 , ,
2. It is the aim of the current work to study the dynamics of solutions of (1.1). We alsopresent numerical simulations to support of findings.
In this subsection we state our main results. We observe that the first five equations of (1.1)decouple from the last equation. Furthermore, we note that the solution to the last equation in(1.1) is completely determined by the solution of the decouple sub-system formed by the first fiveequations. For this reason, we shall focus on the study of the dynamics of the sub-system ddt S ( t ) = B − (cid:80) j =0 γ j S ( t ) Z ( t ) − β ,u S ( t ) I ( t − u ) − µS ( t ) , ddt dS j ( t ) = γ j S ( t ) Z ( t ) − β j,u S j ( t ) I ( t − u ) − µS j ( t ) , j = 1 , , ddt Z ( t ) = qI ( t ) − ηZ ( t ) , ddt I ( t ) = (cid:80) j =0 β j,u S j ( t ) I ( t − u ) − µI ( t ) − dI ( t ) , (1.2)where q, η, β , µ, d, B > β j,u > j ∈ { , } and γ j > j ∈ { , , } , andsubject to the initial conditions( S (0) , S (0) , S (0) , Z (0)) ∈ R , S (0) > , I ( · ) ∈ C ([ − u,
0] : R + ) , (1.3)where R + = { x ∈ R : x ≥ } . Recall that u ≥ C ([ − u,
0] : R + ) withthe sup-norm (cid:107) φ (cid:107) ∞ = sup − u ≤ s ≤ | φ ( s ) | . Due to biological interpretation, we suppose that theinfections rates β j,u is a non-increasing function of u for every j = 0 , ,
2. It follows from standardtheory on Delay-Differential Equations, that (1.2) subject to an initial conditions satisfying (1.3)has a unique non-negative solution, which is defined for all t ≥
0. Moreover, with the notations S ( t ) := (cid:80) j =0 S j ( t ) and N ( t ) := S ( t ) + I ( t ), it holds that ddt N ( t ) = B − µN ( t ) − dI ( t ) − γ Z ( t ) S ( t ) ≤ B − µN ( t ) , t > . Hence, it follows from comparison principle for ODEs thatlim sup t →∞ N ( t ) ≤ Bµ and N ( t ) ≤ max (cid:26) N (0) , Bµ (cid:27) , ∀ t ≥ . (1.4)Thus, using again the comparison principle for ODEs, we getlim sup t →∞ Z ( t ) ≤ qη lim sup t →∞ I ( t ) ≤ qBηµ and Z ( t ) ≤ max (cid:26) Z (0) , qη max (cid:26) N (0) , Bµ (cid:27)(cid:27) . X u = ( S , S , S , Z, I ( · )) ∈ R × C ([ − u,
0] : R + ) : (cid:88) j =0 S j + I (0) ≤ Bµ and Z ≤ qBµη is forward invariant for the flow generated by solutions of (1.2). Observe that when I ( s ) ≡ − u,
0] then I ( t ) ≡ t ≥ Z ( t ) = Z (0) e − ηt → t → ∞ . In which case it is also easilyseen that ( S ( t ) , S ( t ) , S ( t )) → ( Bµ , ,
0) as t → ∞ . Hence, due to biological interpretations, weshall always suppose that our initial conditions belong to X u with the additional properties that I ( s ) > s ∈ [ − u,
0] and S (0) > E = (cid:16) Bµ , , , , (cid:17) T is always an equilibrium solution of (1.2). Theequilibrium state E corresponds to the disease-free equilibrium. We note that any equilibriumsolution ˜E = ( ˜ S , ˜ S , ˜ S , ˜ Z, ˜ I ) T of (1.2) is uniquely determined by its ˜ I coordinate as the remainingcomponents are given by the formulas˜ S = Bµ + b ,u ˜ I , ˜ Z = q ˜ Iη , and ˜ S j = qγ j ˜ IBη ( µ + β j,u ˜ I )( µ + b ,u ˜ I ) , j = 1 , , (1.5)where γ := (cid:80) j =0 γ j and b ,u = β ,u + qη γ . Any equilibrium solution E = ( ˜ S , ˜ S , ˜ S , ˜ Z, ˜ I ) T of(1.2) for which ˜ I > R ,u defined as R ,u = Bβ ,u µ ( µ + d ) . (1.6)To simplify notations in our proofs, we introduce the following new quantity D = µ + dB/µ , (1.7)which is the ratio of the death rate µ + d of the infected class (i.e., the rate at which an individualleaves the infective class) by the total size of the population Bµ when at the disease-free equilibrium.In other words, D represents the per capita death rate of an infective individual in an entirelysusceptible population. We will refer to D as the critical death rate of the model (1.1). Observethat R ,u ≤ β ,u ≤ D and R ,u > β ,u > D .For this reason, we will use D in the remainder of this paper. Theorem 1.1 and Theorem 1.2below are good illustrations that justify the introduction of the critical death rate D . Theorem 1.1 (Stability of the disease-free equilibrium E ) . For every time delay u ≥ ,the disease-free equilibrium E is linearly stable if and only if β ,u < D . Furthermore if max { β ,u , β ,u , β ,u } ≤ D then the disease-free equilibrium is globally stable. By Theorem 1.1, there is no endemic equilibrium whenever β ,u = max { β ,u , β ,u , β ,u } ≤ D .In reality, it is expected that when people start adopting some sort of protection behaviors to avoidcontracting the disease, then the infection rate decreases as well. Then it is natural to supposethat β ,u > max { β ,u , β ,u } . However by curiosity, one might want to know from a mathematical6oint of view what happens if β ,u < max { β ,u , β ,u } . As we shall see from Theorem 3.2 ( ii )later, there is a range of parameters satisfying β ,u ≤ D < max { β ,u , β ,u } for which E is notasymptotically stable. In such case there are two endemic equilibrium solutions of (1.2).Our next result is about the existence and uniqueness of the endemic equilibrium when β ,u > D . Theorem 1.2 (Existence and uniqueness of endemic equilibrium point) . System (1.2) has aunique endemic equilibrium solution whenever β ,u > D . We first remark that Theorem 1.2 follows from Theorem 3.2 below. Theorem 3.2, which is ofindependent interest, provides a complete picture of the existence of endemic solutions of (1.2).Suppose β ,u > D and let E ∗ = ( S ∗ , S ∗ , S ∗ , Z ∗ , I ∗ ) T denotes the unique endemic equilibriumsolution of (1.2) where S ∗ , S ∗ , S ∗ and Z ∗ are given by (1.5) with ˜ I = I ∗ . The existence of I ∗ isuniquely determined by the unique positive solution of the algebraic equation G ( I ) = µ + d = B D µ , (1.8)where G ( I ) := β ,u Bµ + b ,u I + Bqη (cid:88) j =1 β j,u γ j I ( µ + β j,u I )( µ + b ,u I ) . (1.9)It is clear that an explicit expression is very difficult to derive for I ∗ . This makes the analysis ofthe endemic equilibrium solution more difficult. In fact, the study of the sign of the real parts ofthe roots of the characteristic equation of the linearized system at E ∗ , see (2.2), is more complexdue to its form. A plot of the nonlinear function G ( I ) and the line y = µ + d together with thevalue of I ∗ , their intersection point, is given in Figure 2. Note that this graph has been producedwith parameter values estimated in Section 5. I G ( I ) and ( + d ) I* = 1.2553
G(I) + d
Figure 2: Graphs of G ( I ) and µ + d together with their intersection point, I ∗ .We note that while the expression of the disease-free equilibrium depends only on the recruit-ment rate B and the natural death rate µ , the endemic equilibrium however is a function of all theparameters of the model. In reality, health officials are concerned with finding ways to control theoutbreak of a disease so as to minimize its impact on the population, and if possible completely7radicate it. In most cases, the concerned population is provided with more information on thecauses and ways to limit infections. Whence it becomes relevant to know how the parameters ofthe model affect the size of the infected class. In particular, studying the effect of the informationor education dissemination on the control of the disease is of interest. Theorem 1.1 shows thatthe (local) stability of the disease-free equilibrium is independent of the amount of educationavailable about the disease, and the different infections ( S and S classes) and education rates.In order to provide some information on the effect of education on the control of the disease,we decide to study how the magnitude of the component I ∗ of the endemic equilibrium solution,when exists, is affected by the change in populations behavior due to education dissemination.To this end we introduce the quantity τ = qη which measures the “ effective rate of education dissemination ” in the sense that τ > τ < I ∗ with respect to τ when β ,u > D . The following resultholds. Theorem 1.3 (Asymptotic behavior of I ∗ with respect to τ = qη ) . Let u ≥ and suppose that β ,u = max { β j,u : j = 0 , , } > D , and let E ∗ = ( S ∗ , S ∗ , S ∗ , Z ∗ , I ∗ ) T denotes the unique endemicequilibrium solution of (1.2) . The following holds.(i) I ∗ is non-increasing with respect to τ = qη and B ( β ,u − D )( µ + d )( β ,u + qγη ) ≤ I ∗ ≤ min (cid:26) B ( β ,u − D ) µβ ,u , µ D (cid:27) . (1.10) Furthermore if β ,u > min { β ,u , β ,u } or γ > then I ∗ is strictly monotone decreasing withrespect to τ = qη .(ii) If D ≥ γ (cid:80) j =1 β j,u γ j then lim qη →∞ I ∗ = 0 . (1.11) In addition, if D > γ (cid:80) j =1 β j,u γ j then there is a positive constant C > such that B ( β ,u − D )( µ + d )( β ,u + qγη ) ≤ I ∗ ≤ Cq/η ∀ qη > . (1.12) (iii) If D < γ (cid:80) j =1 β j,u γ j , then lim qη →∞ I ∗ = I ∗∞ > where I ∗∞ is the unique positive solution of the equation D = µγ (cid:80) j =1 β j,u γ j µ + β j,u I ∗∞ . Theorem 1.3 confirms that an appropriate education level of the population about the diseasewill always help to contain and minimize significantly its impact. In fact Theorem 1.3 (ii) shows I ∗ converges to 0 at the same rate as ( qη ) − does as qη → ∞ when hypothesis D > γ (cid:80) j =1 β j,u γ j relative infection weight γ (cid:80) j =1 β j,u γ j of the susceptible classes S and S is smaller than the critical death rate D .In such scenario, the size of the infected portion of the population at the endemic equilibriumwill be controlled significantly. On the other hand, according to Theorem 1.3 (iii), we see thatif the relative infection weight γ (cid:80) j =1 β j,u γ j of the susceptible classes S and S is greater thanthe critical death rate D , then the size of the infected portion of the population at the endemicequilibrium will remain uniformly bounded away from zero as a function of τ = qη .As a consequence of Theorem 1.3, the following result on the asymptotic behavior of thesusceptible population of E ∗ with respect to τ = qη can be proved. Theorem 1.4 (Asymptotic behavior of ( S ∗ , S ∗ , S ∗ ) with respect to τ = qη ) . Suppose that β ,u = max { β ,u , β ,u , β ,u } > D , and let E ∗ = ( S ∗ , S ∗ , S ∗ , Z ∗ , I ∗ ) T denotes the unique endemicequilibrium solution of (1.2) . The following holds.(i) If D = γ (cid:80) j =1 β j,u γ j , then lim qη →∞ ( S ∗ , S ∗ , S ∗ ) = (cid:16) , γ Bµγ , γ Bµγ (cid:17) . (1.14) (ii) If D > γ (cid:80) j =1 β j,u γ j , then lim qη →∞ ( S ∗ , S ∗ , S ∗ ) = (cid:18) Bµ + γZ ∗∞ , γ BZ ∗∞ µ ( µ + γZ ∗∞ ) , γ BZ ∗∞ µ ( µ + γZ ∗∞ ) (cid:19) , (1.15) where Z ∗∞ is the unique solution of the algebraic equation D = µ (cid:0) β ,u + Z ∗∞ µ (cid:80) j =1 β j,u γ j (cid:1) µ + γZ ∗∞ . (iii) If D < γ (cid:80) j =1 β j,u γ j and , by setting I ∗∞ = lim qη →∞ I ∗ > , it holds that lim qη →∞ ( S ∗ , S ∗ , S ∗ ) = (cid:18) , γ Bγ ( µ + β ,u I ∗∞ ) , γ Bγ ( µ + β ,u I ∗∞ ) (cid:19) . (1.16)We have seen from Theorem 1.1 that when max { β ,u , β ,u , β ,u } ≤ D then the disease willbe eventually eradicated from the population. We also note from Theorem 1.2 that (1.2) hasa unique endemic equilibrium E ∗ whenever β ,u > D . The following result is about the linearstability of the endemic equilibrium for β , > D . Theorem 1.5.
Let u = 0 and suppose that D ( D + qη γ ) ≥ qη (cid:88) j =1 β j, γ j . (1.17) There is β , , max > D such that if β , ∈ ( D , β , , max ) , then the endemic equilibrium solution E ∗ is linearly stable. Moreover if β , , max < ∞ , then E ∗ is not linearly stable at β , = β , , max . It is not clear whether β , , max can be taken to be infinity. We plan to further explore thisquestion in our future work.Next, we discuss whether the disease will become permanent in the population whenever β ,u > D . 9 heorem 1.6 (Permanence of disease) . Let u ≥ and suppose that β ,u = max { β ,u , β ,u , β ,u } > D . There exists a positive constant m u > such that any solution ( S ( t ) , S ( t ) , S ( t ) , Z ( t ) , I ( t )) of (1.2) with initial in X u satisfies µ ( β ,u − D ) D ( β ,u + qη γ ) ≤ lim sup t →∞ I ( t ) ≤ ( µ + d )( β ,u − D ) D β ,u (1.18) and lim inf t →∞ I ( t ) > m u . (1.19)We note that the eventual a priori upper bound ( µ + d )( β ,u −D ) D β ,u in (1.18) for positive solutionsdecreases with respect to the time delay u . This is because the function β ,u is decreasing in u and the function ( D , ∞ ) (cid:51) β (cid:55)→ ( µ + d )( β −D ) D β is increasing. Hence we conclude that increasing thetime delay has a positive effect in controlling the size of the infected population as time evolves.In fact, our numerical simulations in Section 5 confirm that the infection rates decrease as thetime delay u increases.The rest of the paper is organized as follows. Section 2 is devoted to the proof of Theorem 1.1.The proofs of Theorems 1.2, 1.3, and 1.5 are provided in Section 3 while the proof of Theorem1.6 is presented in section 4. The numerical investigation of the epidemic model is completed inSection 5. In this section we study the stability of the disease-free equilibrium. Equilibrium solutions aresolutions of the system of algebraic equations obtained by setting the right hand sides of (1.2)equal to zeros. At times throughout the rest of the manuscript, to simplify the notations, we willsuppress the dependence of the infection rates on u unless otherwise stated and write β j = β j,u foreach j = 0 , ,
2. Note that this will not cause any confusion in our presentation. The linearizedsystem of (1.2) at an equilibrium point ˜ E = ( ˜ S , ˜ S , ˜ S , ˜ Z, ˜ I ) T is ddt S = − ( µ + β ˜ I + γ ˜ Z ) S ( t ) − β ˜ S I ( t − u ) − γ ˜ S Z ( t ) ddt S j = γ j ˜ ZS − ( µ + β j ˜ I ) S j ( t ) − β j ˜ S j I ( t − u ) + γ j ˜ S ( t ) Z ( t ) , j = 1 , ddt Z = qI ( t ) − ηZ ( t ) ddt I = β ˜ IS ( t ) + β ˜ IS ( t ) + β ˜ IS ( t ) + (cid:80) j =0 β j ˜ S j I ( t − u ) − D Bµ I ( t ) (2.1)where γ = (cid:80) j =0 γ j . The characteristic equation of (2.1) is0 = P u ( λ ; ˜E ) := det( λ I − M u ( λ ; ˜E )) (2.2)where I denotes the identity matrix and M u ( λ ; ˜E ) := − B/ ˜ S − γ ˜ S − β ˜ S e − λu γ ˜ Z − ( µ + β ˜ I ) 0 γ ˜ S − β ˜ S e − λu γ ˜ Z − ( µ + β ˜ I ) γ ˜ S − β ˜ S e − λu − η qβ ˜ I β ˜ I β ˜ I (cid:80) j =0 β j ˜ S j e − λu − D Bµ . P u ( λ ; E ) = (cid:0) λ − β Bµ e − λu + D Bµ (cid:1) ( λ + µ ) ( λ + η ) . (2.3)We present the following lemma. Lemma 2.1.
The disease-free equilibrium E is linearly stable if and only if β < D .Proof. It is clear from (2.3) that λ = − µ and λ = − η are always roots of the characteristicequation P u ( λ ; E ) = 0. Thus the linear stability of E is determined by the signs of the realparts of the roots of the equation0 = λ + D Bµ − β Bµ e − λu := Q ( λ ) − β Bµ e − λu . (2.4)Observe that Q ( − D Bµ ) = 0 and Q (0) = D Bµ . Case 1. D > β . Since λ = − D Bµ < Q ( λ ) = 0 and Q (0) > β Bµ , then allroots of (2.4) have negative real part. Whence E is linearly stable if β < D . Case 2. β = D . It is clear that λ = 0 is a solution of (2.4). Hence E is not linearly stable. Case 3. D < β . Now, observe thatlim λ →∞ Q ( λ ) = + ∞ and lim sup λ →∞ β Bµ e − λu ≤ β Bµ .
Hence the intermediate value theorem guarantees that there exists a positive real number λ u satisfying the equation (2.4). Thus E is not linearly stable.The first statement of Theorem 1.1 follows from the above result. It remains to show the non-linear stability of E when β ≤ D . We first prove the following general result about solution ofthe Cauchy problem. Lemma 2.2.
Let ( S ( t ) , S ( t ) , S ( t ) , Z ( t ) , I ( t )) be a positive solution of (1.2) with initial in X u and let β = max { β , β , β } . Then I (cid:18) βI − Bµ ( β − D ) (cid:19) ≤ . (2.5) where I = lim sup t →∞ I ( t ) .Proof. Recall that (cid:88) j =0 S j ( t ) + I ( t ) ≤ Bµ ∀ t ≥ . Hence (cid:88) j =0 β j S j ≤ β (cid:88) j =0 S j ≤ β (cid:18) Bµ − I ( t ) (cid:19) , ∀ t ≥ . As a result, we obtain that ddt I ( t ) ≤ β (cid:18) Bµ − I ( t ) (cid:19) I ( t − u ) − ( µ + d ) I ( t ) , ∀ t > . < λ (cid:28)
1, there exists t λ (cid:29) I ( t − u ) ≤ I + λ ∀ t ≥ t λ . From the last two inequalities it follows that ddt I ≤ βB ( I + λ ) µ − ( µ + d + β ( I + λ )) I ( t ) ∀ t ≥ t λ , since Bµ ≥ I ( t ) for every t ≥
0. Thus, by the comparison principle for odes, we obtain that I ≤ βB ( I + λ ) µ ( µ + d + β ( I + λ ))Letting λ → I ≤ βBIµ ( µ + d + βI ) . Since, D Bµ = ( µ + d ) and µ + d + β I >
0, the last inequality is equivalent to I (cid:18) βI − Bµ ( β − D ) (cid:19) ≤ , as needed.Thanks to the previous two lemmas, we present the proof of Theorem 1.1. Proof of Theorem 1.1.
Thanks to Lemma 2.1, it remains to show the nonlinear stability of E when β = max { β , β , β } ≤ D . Let ( S ( t ) , S ( t ) , S ( t ) , Z ( t ) , I ( t )) be a positive solution of(1.2) with initials in X u . If β < D , inequality (2.5) implies that I = 0. If β = D , it followsagain from (2.5) that I = 0, which gives I = 0. Therefore, since I ( t ) ≥ t ≥ t →∞ I ( t ) = 0. Next, observe that lim sup t →∞ Z ( t ) ≤ qη lim sup t →∞ I ( t ) = 0. Hencelim t →∞ Z ( t ) = lim t →∞ I ( t ) = 0. This in turn, implies thatlim t →∞ ( S ( t ) , S ( t ) , S ( t )) = (cid:18) Bµ , , (cid:19) . As a result, it follows that lim t →∞ ( S ( t ) , S ( t ) , S ( t ) , I ( t ) , Z ( t )) = E . This completes the proof ofTheorem 1.1. In this section we shall show that if β ,u > D then (1.2) has a unique endemic equilibrium point.Moreover in the next section, we show that the disease will be permanent irrespective of theamount of information and/or education disseminated about it whenever β ,u > D . Our firstresult is about the existence and uniqueness of the endemic equilibrium. Recall that we have set β j = β j,u for each j = 0 , , G ( I, β , τ ) = β Bµ + (cid:0) β + τ γ (cid:1) I + Bτµ + (cid:0) β + τ γ (cid:1) I (cid:88) j =1 β j γ j Iµ + β j I for I ≥ , τ > , β > , (3.1)12here min { β , β , β } >
0. Note that ˜E = ( ˜ S , ˜ S , ˜ S , ˜ Z, ˜ I ) T is an equilibrium of (1.2) with ˜ I (cid:54) = 0if and only if G ( ˜ I, β , qη ) = µ + d = D Bµ . The following results hold. Lemma 3.1.
For every τ > and β ≥ it holds that sup I ≥ G ( I, β , τ ) is finite and achieved.Moreover, the function [0 , ∞ ) (cid:51) β (cid:55)→ max I ≥ G ( I, β , τ ) is strictly increasing and max I ≥ G ( I, β , τ ) > D Bµ for every β > D . Therefore, the followingquantity β τ := min (cid:26) β ≥ I ≥ G ( I, β , τ ) ≥ D Bµ (cid:27) is well definedProof. It is clear thatlim I →∞ G ( I, β , τ ) = 0 and G ( I, β , τ ) ≥ ∀ τ, I, β ≥ . Thus, since I (cid:55)→ G ( I, β , τ ) is continuous, it follows that sup I ≥ G ( I, β , τ ) is finite and achieved.Observe that ∂ β G = B (cid:0) µ + τ γI − τ I (cid:80) j =1 β j γ j Iµ + β j I (cid:1) ( µ + ( β + γτ ) I ) ≥ B (cid:0) µ + τ γI − τ I (cid:80) j =1 γ j (cid:1) ( µ + ( β + γτ ) I ) = B ( µ + τ γ I )( µ + ( β + γτ ) I ) > β ≥ τ > I ≥
0. Thus, we conclude that the function[0 , ∞ ) (cid:51) β (cid:55)→ max I ≥ G ( I, β , τ )is strictly increasing. Note that G (0 , β , τ ) = β Bµ ≥ D Bµ for every β ≥ D . The result thusfollows. Theorem 3.2.
For every τ > let β τ ≤ D be given by Lemma 3.1. The following hold. (i) If < β < β τ then the algebraic equation G ( I, β , τ ) = D Bµ has no nonnegative root. (ii) If β τ < β < D then the algebraic equation G ( I, β , τ ) = D Bµ has two positive roots I − ( β , τ ) < I + ( β , τ ) . Furthermore, ∂ I G ( I − ( β , τ ) , β , τ ) > and ∂ I G ( I + ( β , τ ) , β , τ ) < ,and the functions ( β , τ ) (cid:55)→ I ± ( β , τ ) are of class C . (iii) If β τ < D then for every β ∈ { β τ , D } there is a unique positive solution I ( β , τ ) of thealgebraic equation G ( I, β , τ ) = D Bµ . (iv) If β > D there is a unique positive root I ( β , τ ) of the algebraic equation G ( I, β , τ ) = D Bµ .The functions β (cid:55)→ I ( β , τ ) and τ (cid:55)→ I ( β , τ ) are smooth and ∂ I G ( I ( β , τ ) , β , τ ) < .Proof. ( i ) It is clear from Lemma 3.1 that max I ≥ G ( I, β , τ ) < D Bµ for every 0 < β < β τ . Hencethe result follows. 13 ii ) Let β ∈ ( β τ , D ). Then G (0 , β , τ ) = β Bµ < D Bµ < max I ≥ G ( I, β , τ ). This shows thatmax I ≥ G ( I, β , τ ) is achieved at an interior point I max . Next, observe that ∂ I G = B (cid:0) µ + b I (cid:1) − β ( β + τ γ ) + τ (cid:88) j =1 γ j β j ( µ − β j b I ) (cid:0) µ + β j I (cid:1) . (3.2)Hence ∂ I G ( I, β , τ ) = 0 ⇐⇒ τ (cid:88) j =1 γ j β j − β ( β + τ γ ) = τ (cid:88) j =1 γ j β j (2 µβ j I + ( β j I ) + β j b I ) (cid:0) µ + β j I (cid:1) . An easy computation shows that [0 , ∞ ) (cid:51) I (cid:55)→ γ j β j (2 µβ j I +( β j I ) + β j b I ) (cid:0) µ + β j I (cid:1) is strictly increasing forevery j = 1 ,
2. As a result, the equation ∂ I G = 0 has exactly a unique positive root, whichis I max . Moreover, the functions (0 , I max ) (cid:51) I (cid:55)→ G ( I, β , τ ) and ( I max , ∞ ) (cid:51) I (cid:55)→ G ( I, β , τ )are strictly increasing and decreasing respectively. The intermediate value theorem then im-plies that there exist unique elements I − ( β , τ ) ∈ (0 , I max ) and I + ( β , τ ) ∈ ( I max , ∞ ) such that G ( I ± ( β , τ ) , β , τ ) = D Bµ . Furthermore, ∂ I G ( I − ( β , τ ) , β , τ ) > ∂ I G ( I + ( β , τ ) , β , τ ) < β , τ ) (cid:55)→ I ± ( β , τ ) are ofclass C .( iii ) If D > β τ then as in case ( ii ) max I ≥ G ( I, D , τ ) is achieved at an interior point I max > . In this case, we have I − ( D , τ ) = 0 and I + ( D , τ ) ∈ ( I max , ∞ ) are the two roots of G = D Bµ .Similarly we have I ( β τ , τ ) = I max > iv ) Suppose that β > D . Note from above that D Bµ < G (0 , β , τ ) = min ≤ I ≤ I max G ( I, β , τ ).Thus, as in the above, we can employ the intermediate value theorem to infer that there is exactlya unique solution I ( β , τ ) > G ( I, β , τ ) = D Bµ , with I ( β , τ ) ∈ ( I max , ∞ ) . The next result complements Theorem 3.2 and provides necessary and sufficient condition onthe parameters for which β τ < D . Proposition 3.3.
For every τ > , let β τ be given by Lemma 3.1 and define D ,τ = γ (cid:80) j =1 γ j β j (cid:113) τγ (cid:80) j =1 β j γ j . (3.3) The function τ (cid:55)→ D ,τ is strictly increasing and D = β τ if and only if D ≥ D ,τ . Therefore if D ≥ γ (cid:88) j =1 γ j β j then D = β τ for every τ > , while if D < γ (cid:88) j =1 γ j β j there is τ ∗ > such that D = β τ for every < τ < τ ∗ and D > β τ for every τ ≥ τ ∗ . roof. It is easy to see that the function τ (cid:55)→ D ,τ is strictly increasing and that D ,τ is the onlypositive solution of the algebraic equation in β : β ( β + τ γ ) = τ (cid:88) j =1 γ j β j . Hence if D ≥ D ,τ then D ( D + τ γ ) ≥ D ,τ ( D ,τ + τ γ ) = τ (cid:80) j =1 γ j β j and it is seen from (3.2)that ∂ I G ( I, D , τ ) < ∂ I G (0 , D , τ ) ≤ ∀ I > , which implies that max I ≥ G ( I, D , τ ) = G (0 , D , τ ) = D Bµ . Whence, we deduce from Lemma 3.1that D = β τ since the map β (cid:55)→ max I ≥ G ( I, D , τ ) is strictly decreasing.On the other hand, if D < D ,τ , then using again (3.2), we obtain that ∂ I G (0 , D , τ ) >
0. Thisimplies that max I ≥ G ( I, D , τ ) > G (0 , D , τ ) = D Bµ . Hence by continuity, there is some ˜ β < D close enough such that max I ≥ G ( I, ˜ β , τ ) > D Bµ . This shows that β τ ≤ ˜ β < D . The proof ofProposition 3.3 is then complete. Proof Theorem 1.2.
For every β ,u = β > D , the result follows from Theorem 3.2 ( iv ) with τ = qη . Proof of Theorem 1.3. ( i ) Recalling the function G ( I, β , τ ) introduced in (3.1), we have ∂∂τ G = BI (cid:16)(cid:80) j =1 µγ j ( β j − β ) µ + β j I − γ β (cid:17) ( µ + β I + τ Iγ ) ∀ τ > , I ≥ . Whence, since I ∗ is uniquely determined by the equation G ( I, β , qη ) = D Bµ and β ≥ max { β , β } ,the implicit function theorem implies that ∂∂τ I ∗ = − ∂ τ G∂ I G (cid:40) = 0 if γ = 0 and β = β = β < β > min { β , β } or γ > . Note that we have used Theorem 3.2 ( iv ) to conclude that ∂ I G ( I ∗ , β , qη ) <
0. This shows themonotonicity of I ∗ with respect to qη as stated in the result. A simple computation shows that G (cid:32) µ ( β − D ) D ( β + qγη ) , β , τ (cid:33) > β Bµ + µ D ( β − D ) = D Bµ ∀ β , τ > . By (3.2), we have ∂ I G = B (cid:0) µ + b I (cid:1) − β b + τ (cid:88) j =1 γ j β j ( µ − β j b I ) (cid:0) µ + β j I (cid:1) < B (cid:0) µ + b I (cid:1) − β − τ (cid:88) j =0 γ j β + τ (cid:88) j =1 γ j β j ≤ ∀ I, τ > β ≥ max { β , β } . β = max { β , β , β } then the function (0 , ∞ ) (cid:51) I (cid:55)→ G ( I, β , qη ) is strictly decreas-ing, we conclude that the left hand side inequality of (1.10) holds. Observe that B = µ (cid:88) j =0 S ∗ j + D Bµ I ∗ + γ Z ∗ S ∗ . Hence I ∗ ≤ µ D . Thus the right hand side of (1.10) follows from (2.5).( ii ) Next we proceed by contradiction to show that (1.11) holds under hypothesis D ≥ γ (cid:80) j =1 β j γ j . We first note from the monotonicity of I ∗ ( β , qη ) with respect to qη that thereis I ∗∞ ( β ) ∈ [0 , Bµ ) such that lim τ →∞ I ∗ ( β , τ ) = I ∗∞ ( β ) . Now, we claim that I ∗∞ := I ∗∞ ( β ) = 0. If not, i.e, I ∗∞ >
0, letting qη → ∞ in the equation G ( I ∗ , β , qη ) = B D µ gives B D µ = B (cid:80) j =1 β j γ j I ∗∞ µ + β j I ∗∞ I ∗∞ γ = Bγ (cid:88) j =1 β j γ j µ + β j I ∗∞ < Bµγ (cid:88) j =1 β j γ j which is in contradiction with hypothesis D ≥ γ (cid:80) j =1 β j γ j . Thus (1.11) holds under hypothesis D ≥ γ (cid:80) j =1 β j γ j . Next, suppose that hypothesis D > γ (cid:80) j =1 β j γ j holds and chose a positiveconstant C (cid:29) Bµ + γC β + C (cid:88) j =1 β j γ j µ < D Bµ .
With this choice of the constant C we obtainlim τ →∞ G ( C τ − , β , τ ) = Bµ + γC β + C (cid:88) j =1 β j γ j µ < D Bµ .
Hence, there is τ ∗ (cid:29) I ∗ ( β , qη ) ≤ C q/η ∀ qη ≥ τ ∗ . Therefore, since I ∗ ( β , qη ) is uniformly bounded above by Bµ , we conclude that the right inequalityof (1.12) holds with C = max { C , Bτ ∗ µ } . Clearly the left inequality of (1.12) follows from (1.10).( iii ) Suppose that D < γ (cid:80) j =1 β j γ j . Let I ∗∞ ( β ) denotes the unique positive solution of thealgebraic equation D Bµ = Bγ (cid:88) j =1 β j γ j µ + β j I ∗∞ . For every 0 < ε (cid:28)
1, using the fact thatlim τ →∞ G ((1 ± ε ) I ∗∞ , β , τ ) = Bγ (cid:88) j =1 β j γ j µ + (1 ± ε ) β j I ∗∞ . D Bµ = G ( I ∗ , β , qη ) and the fact that I (cid:55)→ G ( I, β , τ ) is strictlydecreasing that there is τ ε (cid:29) − ε ) I ∗∞ < I ∗ (cid:18) β , qη (cid:19) < (1 + ε ) I ∗∞ ∀ qη ≥ τ ε . Thus we conclude that lim qη →∞ I ∗ (cid:18) β , qη (cid:19) = I ∗∞ ( β ). The proof of the Theorem is complete.The proof of Theorem 1.4 follows similar arguments as those used in the proof of Theorem1.3. So, in order to avoid repeating the same arguments, we do not include its proof. We endthis section with the proof of Theorem 1.5. Our method to prove this result is based on thebifurcation theory. Suppose u = 0 and β j = β j, , j = 0 , , E = ( S , S , S , Z, I ) T and F ( β , E ) be the vector field to the right hand side of (1.2), and we consider the parameter β asour bifurcation parameter. Note that the function ( β , E ) (cid:55)→ F ( β , E ) is of class C ∞ . Recalling(2.1), we have that D E F ( β , E ) = M (0 , E ) , ∀ β > K := D β E F ( β , E ) = − Bµ Bµ (3.4)where D E F ( β , E ) stands for the Jacobian matrix with respect to E only and D β E F ( β , E )the partial derivative of the jacobian matrix with respect to β . Moreover, it follows from (2.3)that r ( β ) := β Bµ − ( µ + d ) = Bµ ( β − D ) ∀ β > D − µB min { µ, η } is the maximal eigenvalue of D E F ( β , E ). The following lemma shows that r ( D ) is a K− simpleeigenvalue of D E F ( D , E ). We denote by { e , e , e , e , e } the canonical basis of R . Givena 5 by 5 square matrix M , we denote by N ( M ) and R ( M ) the kernel and range of the linearoperator induced by M on R , respectively. Lemma 3.4.
Let β > D − µB min { µ, η } and define η ( β ) = η + r ( β ) , µ ( β ) = µ + r ( β ) and J ( β ) := (cid:18) − Bµµ ( β ) (cid:18) β + γqη ( β ) (cid:19) , Bqγ µµ ( β ) η ( β )) , Bqγ µµ ( β ) η ( β ) , qη ( β ) , (cid:19) T . The following hold.(i) r ( D ) = 0 is a K− simple eigenvalue of D E F ( D , E ) . In fact, it holds that N ( D E F ( D , E )) = span { J ( D ) } and K J ( D ) / ∈ R ( D E F ( D , E )) = R × { } . (ii) r ( β ) is a I− simple eigenvalue of D E F ( β , E ) . More precisely, it holds that N ( D E F ( β , E ) − r ( β ) I ) = span { J ( β ) } and J ( β ) / ∈ R ( D E F ( β , E )) = R × { } . roof. ( i ) It is easily seen that with J ( D ) = (cid:16) − Bµ ( D + γqη ) , Bqγ µ η , Bqγ µ η , qη , (cid:17) T it holds that N ( D E F ( D , E )) = span { J ( D ) } and R ( D E F ( D , E )) = { E ∈ R ; < E , e > = } . Thus K J ( D ) = ( − Bµ , , , , Bµ ) T / ∈ R ( D E F ( D , E )).( ii ) It can also be verified by inspection.Next, using Lemma 3.4, Lemma 3.2, and the theory of bifurcation from simple eigenvalues [5,6],we can now present the proof of Theorem 1.5. Proof of Theorem 1.5.
Let us suppose that all the parameters q, η, µ, d, γ , γ , γ > β , β ∈ (0 ,
1) are fixed and emphasize only on the dependence of the endemic equilibrium E ∗ ( β ) withrespect to β > D − µB min { µ, η } . Thanks to Lemma 3.4 ( i ) and [5, Theorem 1.7], there existsome ε > ϕ : ( − ε, ε ) → R with ϕ (0) = 0and Φ : ( − ε, ε ) → span { J ( D ) } T with Φ(0) = 0such that F ( β ( s ) , E ( s )) = 0 for every | s | < ε , where β ( s ) = D + ϕ ( s ) and E ( s ) = E + s J ( D ) + s Φ( s ) (3.5)with span { J ( D ) } T is the orthogonal complement of span { J ( D ) } . Moreover Lemma 3.4 ( ii )and [6, Corollary 1.13 & Theorem 1.16] imply that there exist smooth functions r ∗ : ( − ε, ε ) → R and J ∗ : ( − ε, ε ) → R satisfying( i ) D E F ( β ( s ) , E ( s ))( J ∗ ( s )) = r ∗ ( s ) J ∗ ( s ) ∀ s ∈ ( − ε, ε )( ii ) J ∗ (0) = J ( D ) and r ∗ (0) = r ( D ) = 0( iii ) r ∗ ( s ) and − sϕ (cid:48) ( s ) r (cid:48) ( D ) have the same zeros and, same sign whenever r (cid:48) ( D ) (cid:54) = 0 , where r (cid:48) ( · ) is the derivative of r ( · ). From the expression r ( β ) = Bµ ( β − D ) it easily follows that r (cid:48) ( β ) = Bµ for every β . Thus it follows from ( iii ) that r ∗ ( s ) has the same sign as − sϕ (cid:48) ( s ) for | s | small. Next, we determine the sign of ϕ (cid:48) ( s ) for | s | (cid:28)
1. Observe from (3.5) that1 s < E ( s ) , e > = 1+ < Φ( s ) , e > −→ s → . This shows that I (cid:48) (0) = 1 and hence E ( s ) is an equilibrium solution of (1.2) for which I ( s ) (cid:54) = 0for 0 < | s | (cid:28)
1. As a result, recalling the function G introduced in (3.1), we must have that G ( I ( s ) , β ( s ) , qη ) = D Bµ for 0 < | s | (cid:28)
1. This is equivalent to saying that (cid:18) µ + (cid:18) β ( s ) + qη γ (cid:19) I ( s ) (cid:19) D µ − β ( s ) = qη (cid:88) j =1 β j γ i I ( s ) µ + β j I ( s ) ∀| s | (cid:28) . Taking the derivative of both sides withe respect to s yields D µ (cid:18) β (cid:48) ( s ) I ( s ) + (cid:18) β ( s ) + qη γ (cid:19) I (cid:48) ( s ) (cid:19) − β (cid:48) ( s ) = µqI (cid:48) ( s ) η (cid:88) j =1 β j γ j ( µ + β j I ( s )) (3.6)18valuating this equation at s = 0 and using the fact that β (0) = D , I (0) = 0 and I (cid:48) ( s ) = 1yield ϕ (cid:48) (0) = β (cid:48) (0) = 1 µ D (cid:18) D + qη γ (cid:19) − qη (cid:88) j =1 β j γ j . From this point, we distinguish two cases.
Case 1 . ϕ (cid:48) (0) >
0, so that there is a transcritical bifurcation at β = D . Then β ( s ) > D for 0 < s (cid:28) E ( s ) is an endemic equilibrium of (1.2) corresponding to β = β ( s ) for0 < s (cid:28)
1. Using ( iii ) we conclude that r ∗ ( s ) < ϕ (cid:48) ( s ) > < s (cid:28)
1. It is clearfrom ( i ) that r ∗ ( s ) is an eigenvalue of D E F ( β ( s ) , E ( s )). Hence, recalling that the eigenvaluesof D E F ( β , E ( β )) consist of the set {− µ, − η, r ( β ) } with r ( β ) > < ε ∗ < ε such that r ∗ ( s ) is also the largest eigenvalue of D E F ( β ( s ) , E ( s )) for s ∈ ( − ε ∗ , ε ∗ ). Therefore, E ( β ( s )) is linearly stable for every s ∈ (0 , ε ∗ ).Now, recalling Theorem 3.2, we know that the only endemic equilibrium of (1.2) correspondingto β ( s ) is E ∗ ( β ( s )). Thus we have that E ∗ ( β ( s )) = E ( s ) and s ∈ (0 , ε ∗ ). Since the function(0 , ε ∗ ) (cid:51) s (cid:55)→ β ( s ) is strictly increasing and hence invertible, then for every β ∈ ( D , β ( ε ∗ )),we have that E ∗ ( β ) is linearly stable. Case 2. ϕ (cid:48) (0) = 0, that is, we have a pitchfork bifurcation at β = D . Differentiating again(3.6) with respect to s gives D µ (cid:18) β (cid:48)(cid:48) ( s ) I ( s ) + 2 β (cid:48) ( s ) I (cid:48) ( s ) + (cid:18) β ( s ) + qη γ (cid:19) I (cid:48)(cid:48) ( s ) (cid:19) − β (cid:48)(cid:48) ( s )= µqI (cid:48)(cid:48) ( s ) η (cid:88) j =1 β j γ j ( µ + β j I ) − µq [ I (cid:48) ( s )] η (cid:88) j =1 β j γ j ( µ + β j I ( s )) Evaluating the last equation at s = 0 and using the fact that β (cid:48) (0) = 0, i.e, D (cid:0) D + qη γ (cid:1) = qη (cid:80) j =1 β j γ j , we obtain β (cid:48)(cid:48) ( s ) = 2 qηµ (cid:88) j =1 β j γ j > . Thus we obtain that ϕ (cid:48) ( s ) > < s (cid:28)
1, so we are back to case 1 and we can employ thesame arguments to deduce that E ( β ) is linearly stable for β ∈ ( D , β ( ε ∗ )).In either case we may now take β , max := sup { ˜ β > D : max { Re ( λ ) : λ ∈ σ ( D E F ( β , E ( β ))) } < ∀ D < β < ˜ β } , where σ ( D E F ( β , E ( β ))) is the spectrum of D E F ( β , E ( β )) and Re ( λ ) is the real part of thecomplex number λ . This completes the proof of Theorem 1.5. β ,u > D . In this section, we discuss the permanence of the disease when β ,u > D . Recall again that wehave set β j = β j,u for each j = 1 , , emma 4.1 (Uniform weak permanence) . Let u ≥ be fixed and suppose that β > D . Thenfor every solution ( S ( t ) , S ( t ) , S ( t ) , Z ( t ) , I ( t )) of (1.2) with initial in X u satisfying I ( s ) > forevery s ∈ [ − u, , it holds that I := lim sup t →∞ I ( t ) ≥ m := µ ( β − D ) D ( β + qγη ) > . (4.1) Proof.
We proceed by contradiction to establish the result. So, we suppose that there is apositive solution ( S ( t ) , S ( t ) , S ( t ) , Z ( t ) , I ( t )) of (1.2) with initial in X u satisfying I ( s ) > s ∈ [ − u,
0] and λ ∈ (0 ,
1) such that I := lim sup t →∞ I ( t ) < λ m. (4.2)We fix λ ∈ (0 , − λ ). Using the definition of limsup, there is t (cid:29) t ≥ t I ( t ) ≤ ( λ + λ ) m. (4.3)Recalling the equation satisfied by ddt Z ( t ), we obtain that ddt Z ( t ) ≤ ( λ + λ ) qm − ηZ ( t ) ∀ t ≥ t . Thus, if we fix λ ∈ (0 , − λ + λ ), there is some t > t + u such that Z ( t ) ≤ ( λ + λ + λ ) qmη ∀ t ≥ t . (4.4)It then follows from (4.3) and (4.4) that ddt S ( t ) ≥ B − (cid:18) µ + (cid:18) ( λ + λ ) β + ( λ + λ + λ ) qγη (cid:19) m (cid:19) S , ∀ t ≥ t . As a result, there is t > t such that S ( t ) ≥ Bµ + ( λ + λ + λ ) (cid:16) β + qγη (cid:17) m ∀ t ≥ t . (4.5)This in turn implies that ddt I ( t ) ≥ β Bµ + ( λ + λ + λ ) (cid:16) β + qγη (cid:17) m I ( t − u ) − B D µ I t ≥ t , n ≥ . (4.6)Observe from the formula of (cid:16) β + qγη (cid:17) m = µ ( β −D ) D and the fact that λ + λ + λ < β Bµ (cid:16) µ + ( λ + λ + λ ) (cid:16) β + qγη (cid:17) m (cid:17) B D > β Bµ (cid:16) µ + (cid:16) β + qγη (cid:17) m (cid:17) D B = β Bµ (cid:16) µ + µ ( β −D ) D (cid:17) B D = 1 . Case 3 of the proof of Lemma 2.1 show that the algebraic equation λ − β Bµ + ( λ + λ + λ ) (cid:16) β + qγη (cid:17) m e − λu + D Bµ = 0has a unique positive root λ >
0. Hence, the solution I ( t ) to the first-order linear delay differ-ential equation ddt I = (cid:32) β Bµ +( λ + λ + λ ) (cid:16) β + qγη (cid:17) m (cid:33) I ( t − u ) − ( µ + d ) I t > t I ( s ) = I ( s ) > s ∈ [ t − u, t ] (4.7)blows up exponentially as t → ∞ . From (4.6) and (4.7), we conclude that I ( t ) ≥ I ( t ) ∀ t ≥ t . As a result, we obtain that lim sup t →∞ I ( t ) ≥ lim t →∞ I ( t ) = ∞ . This clearly contradicts (4.2). Therefore, the statement of the Lemma must hold.We first remark that both inequalities (2.5) and (4.1) complete the proof of (1.18). Thanks toLemma 4.1 and the theory developed in [23], we can now conclude our next result on the uniformpersistence of the disease whenever β > D . We first introduce few terminologies. Note that(1.2) generates a continuous semiflow on the set˜ X u := { ( S , S , S , Z , I ( · )) ∈ X u : I ( s ) > ∀ s ∈ [ − u,
0] and S > } , which we denote by Φ t , i.e, Φ t ( S , S , S , Z , I ( · )) = ( S ( t ) , S ( t ) , S ( t ) , Z ( t ) , I t ( · )) for every t ≥ I t ( s ) := I ( t + s ) for every s ∈ [ − u,
0] and t ≥
0. We note from the Arzela-Ascoli’s Theoremthat Φ t is compact for every t > u . For every ( S , S , S , Z , I ( · )) ∈ ˜ X u , define the persistencefunction ρ (( S , S , S , Z , I ( · ))) = I (0) ∀ ( S , S , S , Z , I ( · )) ∈ ˜ X u . (4.8)In particular, ρ (Φ t ( S , S , S , Z , I ( · ))) = I t (0) = I ( t ) ∀ t ≥ S , S , S , Z , I ( · )) ∈ ˜ X u . Recall that I ( t ) > t > I (0) >
0. Hence ρ (Φ t ( S , S , S , Z , I ( · ))) > ρ (Φ ( S , S , S , Z , I ( · ))) > . Proof of Theorem 1.6.
Recall that Φ t is compact for every t > u . We also recall that (cid:88) j =0 | S j ( t ) | + (cid:107) I t (cid:107) ∞ ≤ Bµ and | Z ( t ) | ≤ qBηµ , ∀ t ≥ u, which implies that Φ is pointwise dissipative and uniformy bounded for t ≥ u . Therefore, it followsfrom [23, Theorem 2.30 Page 41] that (1.2) has a compact attractor A of ˜ X u . Next, observe that21he map ρ ◦ Φ is continuous, and Lemma 4.1 guarantee that Φ is uniformly weak ρ − persistent.Therefore, it follows from [23, Theorem 5.2 Page 126] that φ is uniformly ρ − persistent, i.e, thereexists m u > t →∞ ρ (Φ t ( S , S , S , Z , I ( · ))) > m u ∀ ( S , S , S , Z , I ( · )) ∈ ˜ X u . This completes the proof of the Theorem since ρ (Φ t ( S , S , S , I , Z )) = I ( t ) for every t > S , S , S , Z , I ( · )) ∈ ˜ X u . In this section we consider investigating the numerical studies related to the theoretical resultsoutlined in previous sections. We will consider model (1.1) with u = 3 k for k = 0 , , , , While any set of data related to compartment classes in our model could be used to illustratenumerical simulations, we decide to use the data from Uganda within the period 1992-2005. Usingsuch a real data has the advantage of testing our mathematical model, and we can estimate thedifferent infection rates with respect to the time delay. As a result, we are able to confirm thatthe infection rates decrease as functions of the time delay. The data is reported in [12] and isgiven in the following Table 2. Table 2: Data tableYear Susceptible Infected Information1997 600/1,2001999 6,700,000 606,0002000 6,597,470 573,693 700/1,2002001 7,130,000 383,000 717/1,2002003 7,462,000 544,0002005 7,636,000 548,261 778/1,200We recall that the data in Table 2 was collected from UN, UNAIDS, and the Uganda AIDSCommission for population, which included death rates, percentage of adults (ages 15 to 59)per year, growth of the adult class, the adult prevalence rates, and the percentage of adultpopulation infected. Moreover, subject matter experts, literature, and surveys were used todetermine organizational estimates of the information and education campaigns involved, ratesfor interaction of Z to susceptible classes, and efficacy of abstinence and faithfulness and the useof condom behavior. We emphasize that our model and the data collected considers only adultindividuals ages 15 to 59.It is important to note that our model (1.1) differs from the model considered by the authorsin [12] in the sense that in [12], the logistic growth was used for the Z class with a growthcoefficient that increases with the number of infectives and is multiplied by a ratio of the infectedto the living population. Our model assumes that the growth of information is proportional tothe number of infectives. Moreover, the R class in [12] only collects the number of people that diefrom the disease, while in our case the R population represents the number of susceptibles that22ill not contract the disease due to change in behavior once they get the information about thedisease. Finally, the model in [12] is a system of ordinary differential equations while our modelis a system of delay differential equations.The authors in [12] calculated values of the parameters B , γ , γ , d , and the general or naturaldeath rate µ . Also, it was assumed that β ,u and β ,u are directly proportional to β ,u , that is, β ,u = 0 . β ,u and β ,u = 0 . β ,u . They then use parameter estimation to fit the model to thedata in order to estimate together β ,u and the logistic growth rate for the Z class. Due to the similarity of the susceptible class equations and the infected group equation in ourmodel (1.1) to those of the model in [12] (the term with γ in (1.1) is not present in the modelin [12]), we choose to fix all parameters as in [12] as outlined in the previous section but we decideto fit the model (1.1) to the data in Table 2 in order to estimate the parameters β ,u , η , γ , and q . We recall the values of the fixed parameters in [12] in Table 3.Table 3: Values of fixed parameters in model (1.1)B β ,u β ,u d γ γ µ . β ,u . β ,u S (0) = 0, S (0) = 0, and R (0) = 0.According to the UN population reports, there were 18.38 million individuals in Uganda in1992, of which 32 .
1% represented the adult population ages 15 to 59. This allows us to determinethat there were 5,899,980 adult individuals for ( S (0) + I (0)) class for that year. At the sametime, it is also reported that the adult prevalence rate was about 15%, giving us 884,997 infectedadult individuals initially. That is, I (0) = 884 , S (0) = 5 , , Z (0), theauthors in [12] estimated that initially there were 240 organizations out of 1,200 (that is, 20%)that were involved in the abstinence and faithfulness, and the use of condom behaviors. Hence, Z (0) = 0 . S (0) , S (0) , S (0) , I (0) , Z (0) , R (0)) T = (5 . , , , . , . , T . We use MATLAB to run the numerical simulations. We develop and implement an algorithmto fit our model to the data. Specifically, for the delay differential system (1.1), our programcreates three nested functions within a main function. One nested function specifies the historyof the solution ( S ( t ) , S ( t ) , S ( t ) , I ( t ) , Z ( t ) , R ( t )) T for t ≤
0. In our case we use the initial23onditions for this history. The other nested function implements the delay system to evaluatethe right side of the differential equations (1.1). Finally, a third nested function uses the delaydifferential equations (DDE) solver dde23 to solve the system (1.1) then calculates and returnsthe relative sum of squared errors (SSE) between the data in Table 2 and the model solution (1.1).On the other hand, we will also be interested in the solution of system (1.1) with u = 0, i.e.,there is no delay. In this case, our program includes two nested functions within a main function.One nested function uses the ODE solver ode45 to implement and solve the system (1.1) with u = 0. The other nested function receives as input arguments the parameters to be estimatedand time t , then calculates and returns the relative sum of squared errors (SSE) between the datain Table 2 and the model solution (1.1) with u = 0. The relative SSE is calculated as follows:Model error ( β , , η, γ , q ) = (cid:88) j (cid:20) | Model( S + S + S ) j − Data( S ) j | | Data( S ) j | + | Model( I ) j − Data( I ) j | | Data( I ) j | (cid:21) + (cid:88) k (cid:20) | Model( Z ) k − Data( Z ) k | | Data( Z ) k | (cid:21) , where j = 1999 , , , , k = 1997 , , , u > u = 0, we use the MATLAB nonlinear programming solver fminsearch , which finds minimumof unconstrained multi-variable function using derivative-free method. In our case, fminsearch receives two input arguments, a handle to the nested error function and the initial guesses of theparameters to be estimated, then it returns the values of the fitted parameters.We also emphasize that in all of our simulations, we use mass action incidence terms like β ,u S ( t ) I ( t ). The standard incidence terms like β ,u S ( t ) I ( t ) /N , where N is the total populationdidn’t seem to provide a good fit between data and the model solution.Now, obtaining pre-estimates or initial guesses of the parameters to be determined is not aneasy task. In order to obtain “good” pre-estimates, we reached out to Wolfram Mathematica’s Manipulate command. This command allows, after setting up the model (1.1) with quantities tobe estimated treated as parameters, to plot both the data and the model solutions on the samegraph. Then “manipulate” them until an agreement between the data and the model solutions isachieved. We then select the parameter values for which this agreement is made and take themas initial guesses.
We will now investigate the numerical behavior of the model (1.1). To demonstrate the importanceof the information Z , we will first describe model (1.1) without such information, i.e., with Z = 0.In this case the system reduces into the simplest HIV epidemic model, the following SI model (cid:40) ddt S ( t ) = B − β ,u S ( t ) I ( t − u ) − µS ( t ) , ddt I ( t ) = β ,u S ( t ) I ( t − u ) − ( µ + d ) I ( t ) . (5.1)It is important to note that the model (5.1) does not fit HIV data well. This could possibly beexplained from the fact that for a simple SI model, the time spent in the infectious stage decreasesexponentially, which is an unrealistic assumption for HIV, for which the infectious stage is longand the duration is affected by several variations. See for instance, [15, page 131]. For this reason24e will not fit the model (5.1) to the data. However, for the purpose of illustration to show howuseful information campaigns are, we will provide numerical solutions to model (5.1). To thisend, we use parameter values for B , d , and µ in Table 3. For the choice of β ,u , we have to makesure that neither the I class nor the S class fits their respective data. For if this happens, it willresult in favoring one class over the other. Thus, β ,u = 0 . u = 6. The graphs of the twoclasses of the model together with their data are given in the first two figure windows in Figure 3.The last figure window provides only the graphs of the two classes of the model (5.1) plottedtogether. Clearly, we can observe that in the absence of information, more and more people joins time in years D a t a and m ode l s u sc ep t i b l e Susceptible DataModel Susceptible time in years D a t a and M ode l i n f e c t ed Infected DataModel infected time in years c o m pa c t m en t s o f t he m ode l Infected and Susceptible of the model solution
Model S Model I
Figure 3: The graphs of the infected and susceptible model solutions and their data for model (5.1)calculated with parameter values from Table 3 and β , in Table 4.the at-risk population, and as a consequence, the number of infected individuals increases.Now, we will fit the model (1.1) to the data in Table 2 for one time delay value u = 6 andobtain the resulting fitted model solutions over the interval [0 ,
15] representing the years from1992 to 2006. We use the Wolfram Mathematica’s
Manipulate command to obtain the initialguesses of the parameters. The reader interested in this valuable process can learn from our codeand the accompanying graph in Figure 4.We obtain the initial guesses β , = 0 . η = 0 . γ = 0 . q = 0 . u = 6Delay u β , η γ q SSE error6 0.021509 0.055251 0.252445 0.099323 1.309197Figure 5 shows the graphs of the fit between the data and the model solutions with u = 6,where, by the total susceptible population we mean the sum S + S + S .To see the effect of information campaigns on the control of the disease, we also plot in Figure 6the graphs of the solutions for all subgroups of the model obtained with time delay u = 6.The numerical results in Figure 6 suggest that thanks to the information about the disease,the number S of general susceptible people that has not yet changed their behavior before25 * This code sets up the system of delay differential equations,treating quantities to be estimated as parameters and use "Manipulate"command to obtain "good" pre - estimates of the parameters. S, M,and Z represent the susceptible, infected, and information, respectively. *) Sdata = {{
8, 6.7 } , {
9, 6.59747 } , {
10, 7.13 } , {
12, 7.462 } , {
14, 7.636 }} ;Mdata = {{
8, 0.606 } , {
9, 0.573693 } , {
10, 0.383 } , {
12, 0.544 } , {
14, 0.548261 }} ;Zdata = (cid:1)(cid:1)
6, 1 (cid:2) (cid:3) , (cid:1)
9, 7 (cid:2) (cid:3) , (cid:1)
10, 717 (cid:2) (cid:3) , (cid:1)
14, 778 (cid:2) (cid:3)(cid:3) ;Manipulate (cid:4) sol = First (cid:4)
NDSolve (cid:4)(cid:1)
S0' [ t ] == - gamma0 * S0 [ t ]* Z [ t ] - * S0 [ t ]* Z [ t ] - * S0 [ t ]* Z [ t ] - beta0 * S0 [ t ]* M [ t - ] - * S0 [ t ] , S0 [ t / ; t <= ] == [ t ] == * S0 [ t ]* Z [ t ] - * beta0 * S1 [ t ]* M [ t - ] - * S1 [ t ] ,S1 [ t / ; t <= ] == [ t ] == * S0 [ t ]* Z [ t ] - * beta0 * S2 [ t ]* M [ t - ] - * S2 [ t ] ,S2 [ t / ; t <= ] == [ t ] == beta0 * S0 [ t ]* M [ t - ] + * beta0 * S1 [ t ]* M [ t - ] + * beta0 * S2 [ t ]* M [ t - ] - (cid:5) + (cid:6) * M [ t ] , M [ t / ; t <= ] == [ t ] == q * M [ t ] - eta * Z [ t ] , Z [ t / ; t <= ] == [ t ] == gamma0 * S0 [ t ]* Z [ t ] - * R [ t ] , R [ t / ; t <= ] == (cid:3) , { S0, S1, S2, M, Z, R } , { t, 0, 15 } (cid:7)(cid:7) ;Show [ Plot [ Evaluate [{ S0 [ t ]+ S1 [ t ]+ S2 [ t ] , M [ t ] , Z [ t ]} / . sol ] , { t, 0, 15 } ,PlotLegends → { "S [ t ] ", "M [ t ] ", "Z [ t ] " }] ,ListPlot [{ Sdata, Mdata, Zdata } , PlotLegends → { "Sdata", "Mdata", "Zdata" }]] , { beta0, 0, 1, 0.001 } , { eta, 0, 1, 0.001 } , { gamma0, 0, 1, 0.001 } , { q, 0, 1, 0.001 } (cid:7) ����� ����� ��� ����� ������ ����� � ����� S [ t ] M [ t ] Z [ t ] SdataMdataZdata fithivddeDelay6.nb Figure 4: Wolfram Mathematica code using the
Manipulate
Command to pre-estimate parameterswith delay u = 6 in model (1.1). time in years D a t a and f i tt ed s u sc ep t i b l e Susceptible DataTotal Fitted Susceptible time in years D a t a and f i tt ed i n f e c t ed Infected DataFitted infected time in years D a t a and f i tt ed i n f o r m a t i on Education DataFitted Education
Figure 5: Plots of the graphs of the fit between the model solutions and the data with time delay u = 6 and total susceptible population S + S + S . time in years c o m pa c t m en t s o f t he m ode l Plot of all subgroups of the fitted model solution
Fitted S Fitted S Fitted S Fitted IFitted ZFitted R
Figure 6: The fitted model with each of its population subgroup obtained with delay u = 6.26nformation decreases from 5 million to about a million once they are aware of the information.At the same time, once they acquire knowledge about the disease, the numerical simulationsshow that more people will tend to use condoms than practicing abstinence and faithfulness.As a consequence, the number of infected people decreases considerably. The fitted model herealso shows that the number of removed people will surpass the number of susceptible practicingabstinence and faithfulness but will be smaller compared to the condom users. In light of this, it isclear that the information Z plays a central role in decreasing the number of ignorant susceptibleand the infected within this population in contrast to the numerical results from the last figurewindow in Figure 3 where there was no information education within the population.To understand the effect of the delay on the number of infected individuals, we use twoapproaches. We will fix the estimated parameter values obtained by fitting the model to the datawith u = 6 to obtain the model solutions for different other chosen time delay values. So, in thisfirst step we are assuming that the infection rates are constant with respect to the time delay.Next, we will fit the model for different time delays u = 3 k , k = 0 , , , ,
4, and compare theinfections rates β ,u obtained. This will help to check whether the rates decrease with the timedelay.For the first approach, we fix the estimated parameter values obtained from fitting the modelto the data over the interval [0 ,
15] with the time delay u = 6 as reported in Table 4 and runthe model solutions for the time delays u = 0, u = 3, u = 6, u = 9, and u = 12 with theseestimated parameter values without doing another parameter estimation over the time interval[0 , time in years I n f e c t ed Infected classes with delays u = 0, 3, 6, 9, 12
Infected with delay u=0Infected with delay u=3Infected with delay u=6Infected with delay u=9Infected with delay u=12
Figure 7: The infected model solutions for the time delays u = 0 , , , ,
12 over [0 , u = 6 over [0 , ,
15] for each delay u = 3 k , k = 0 , , , ,
4. The initial guesses usedfor all delays are the same as reported in Figure 4. At the same time, we also record the effectiverate of education dissemination τ = q/η and the corresponding SSE error in Table 5.27able 5: Estimated parameters with their corresponding errors for different delay valuesDelay u β ,u η γ q τ SSE error0 0.025801 0.041590 0.249756 0.087683 2.108258 1.3773403 0.023290 0.049984 0.252212 0.094138 1.881626 1.3507536 0.021509 0.055251 0.252445 0.099323 1.797668 1.3091979 0.020364 0.061623 0.252566 0.105675 1.714851 1.26727112 0.019605 0.063318 0.253071 0.107944 1.704781 1.238921As expected, we see from Table 5 that the infection rates decrease with respect to the timedelay. Consequently, the number on infectives would decrease as the time delay increases. As aconsequence, we expect that the effective education dissemination rate τ := qη should decrease asthe time delay increases. This is in fact confirmed by the estimates values for τ obtained in Table5. Therefore, the numerical results support our theoretical results.We conclude this section with an illustration of Theorem 1.2 on the existence and uniquenessof the endemic equilibrium and Theorem 1.6 about the permanence of the disease. Using thevalues of the model parameters described in Table 3 and in Table 4, we have that D = µ ( µ + d ) B = 0 . , µβ , − D ) D (cid:16) β , + qη γ (cid:17) = 0 . , ( µ + d )( β , − D ) D β , = 39 . . Clearly, we have that β , > D . Therefore: • Theorem 1.2 guarantees the existence and uniqueness of the endemic equilibrium E ∗ =( S ∗ , S ∗ , S ∗ , Z ∗ , I ∗ ) T = (0 . , . , . , . , . T , where S ∗ , S ∗ , S ∗ , Z ∗ aregiven by (1.5) with ˜ I = I ∗ and I ∗ is the unique positive solution of Equation (1.8). • The hypotheses of Theorem 1.6 are fulfilled and inequality (1.18) is satisfied and becomes µ ( β , − D ) D (cid:16) β , + qη γ (cid:17) = 0 . ≤ lim sup t →∞ I ( t ) ≤ . µ + d )( β , − D ) D β , , confirming one more time that, the disease will be permanent within this population. A close look at the results established in this paper suggests several directions to consider inthe future. First, note that Theorem 1.5 provides the linear stability of the endemic equilibriumfor the model without delay under the condition stated there. However, it will be interesting toestablish this stability for the model with delay. In particular, constructing a Lyapunov functionfor the endemic equilibrium solution to establish global stability would be valuable.Recalling that the time delay u > β ,u decreasesin u and that lim u →∞ β ,u = 0. In addition, due to practical application of the model, we can alsosuppose that β ,u > max { β ,u , β ,u } for every u ≥
0. Under such hypotheses it then becomesnatural to take the time delay u as the extra bifurcation parameter in the model. In this case,28y Theorem 1.1, we note that the disease-free equilibrium E is the global attractor for solutions(1.2) for every u ≥ β , ≤ D . Now, let us suppose that β , > D . Then there is a unique u ∗ > β ,u > D for every u ∈ [0 , u ∗ ) and β ,u ≤ D for every u ≥ u ∗ . As a consequenceof Theorem 1.1 we conclude that E is the global attractor for solutions of (1.2) for every u ≥ u ∗ .Moreover by Theorem 1.2, for every u ∈ [0 , u ∗ ), (1.2) has a unique endemic equilibrium E ∗ ( u ).As a result, u = u ∗ is a bifurcation point. It is then interesting to know whether the endemicequilibrium E ∗ ( u ) is stable for every u ∈ [0 , u ∗ ).Another important aspect to consider is the notion of spreading speeds of the disease. Oneway to study this problem is through the existence of traveling wave solutions with a certainspeed c and connecting the disease-free equilibrium and the endemic equilibrium solutions. Weplan to study this question and how the time delay affects the spreading speeds in our futurework. References [1] Un aids data 2018. , 2018.[2] Hiv basics. https://hivcare.org/hiv-basicshttps://hivcare.org/hiv-basics