# The Dynamics of Epidemic Model with Two Types of Infectious Diseases and Vertical Transmission

Authors: Raid Kamel Naji; Reem Mudar  Hussien
Journal: Journal of Applied Mathematics
 (2016)

---

## Abstract

An epidemic model that describes the dynamics of the spread of infectious diseases is proposed. Two different types of infectious diseases that spread through both horizontal and vertical transmission in the host population are considered. The basic reproduction numberR0 is determined. The local and the global stability of all possible equilibrium points are achieved. The local bifurcation analysis and Hopf bifurcation analysis for the four-dimensional epidemic model are studied. Numerical simulations are used to confirm our obtained analytical results.


---

## Body

## 1. Introduction

Mathematical models can be defined as a method of emulating real life situations with mathematical equations to expect their future behavior. In epidemiology, mathematical models play role as a tool in analyzing the spread and control of infectious diseases. Although one of the most famous principles of ecology is the competitive exclusion principle that stipulates “two species competing for the same resources cannot coexist indefinitely with the same ecological niche” [1, 2], Volterra was the first scientist who used the mathematical modeling and showed that the indefinite coexistence of two or more species limited by the same resource is impossible [3]. Moreover, Ackleh and Allen [4] were the first who used the competitive exclusion principle of the infectious disease with different levels in single host population.It is well known that one of the most useful parameters concerning infectious diseases is called basic reproduction number. It can be specific to each strain of an epidemic model. In fact the basic reproduction number of the model is defined as the maximum reproduction numbers of other strains [5–7]. Diekmann et al. [8] had studied epidemic models with one strain, while Martcheva in [9] studied the SIS-type of disease with multistrain. However, Ackleh and Allen [10] studied SIR-type of disease with n strain and vertical transmission.Keeping the above in view, in our proposed model two strains with two different types of infectious diseases are considered. Accordingly two different reproduction numbers are obtained and then competitive exclusion principle is presented. It is assumed that two different types of diseases transmission, say horizontal and vertical transmission, are used too. The horizontal transmission occurs by direct contact between infected and susceptible individuals, while vertical transmission occurs when the parasite is transmitted from parent to offspring [11–13]. The incidence of an epidemiological model is defined as the rate at which susceptible becomes infectious. Different types of incidence rates are introduced into literatures [14–17]. Finally two types of incidence rates, say bilinear mass action and nonlinear type, are used with the horizontal and vertical transmission, respectively. The local and global stability for all possible equilibria are carried out with the help of Lyapunov function and LaSalle’s invariant principle [18]. An application of Sotomayor theorem [19, 20] for local bifurcations is used to study the occurrence of local bifurcations near the equilibria. The Hopf bifurcation [21, 22] conditions are derived. Finally, numerical simulations are used to confirm our obtained analytical results and specify the control set of parameters.


## 2. Model Formulation

Consider a real world system consisting of a host populationN(t) that is divided into four compartments: St which represents the number of susceptible individuals at time t; I1t and I2t that represent the number of infected individuals at time t for SIRS-type of disease and SIS-type of disease, respectively; finally R(t) that represents the number of recovered individuals at time t, thus Nt=St+I1(t)+I2(t)+R(t). Now in order to formulate the dynamics of the above system mathematically, the following assumptions have been adopted:(1)
There is a constant number of the host populations entering to the system with recruitment rateΛ>0.(2)
There is a vertical transmission of both of the diseases; that is, the infectious host gives birth to a new infected host of rates0≤p1≤1 and 0≤p2≤1 for the diseases I1 and I2, respectively. Consequently p1I1 and p2I2 individuals enter into infected compartments I1 and I2, respectively, and the same quantities are disappearing from recruitment in the susceptible compartment.(3)
The diseases are transmitted by contact, according to the mass action law, between the individuals in theS-compartment and those in Ii(i=1,2) compartments with nonlinear incidence rate for I1 that is given by β1SI1/(1+I1), in which β1>0 represents the infection force rate while 1/(1+I1) represents the inhibition effect of the crowding effect of the infected individuals, and linear incidence rate for I2 that is given by β2SI2, where β2>0 represents the infection rate.(4)
The individuals in theI1 compartment are facing death due to the disease with infection death rate α1≥0. They recover from disease and get immunity with a recovery rate δ>0.(5)
The individuals in theI2 compartment are facing death due to the disease with infection death rate α2≥0. They also recover from the disease but return back to be susceptible with recovery rate γ>0.(6)
The individuals in theR compartment are losing the immunity from the I1 disease and return back to be susceptible again with losing immunity rate 0≤η<1.(7)
There is a natural death rateμ>0 for the individuals in the host population. Finally, it is assumed that both the diseases cannot be transmitted to the same individual simultaneously.According to these assumptions the dynamics of the above real world system can be represented mathematically by the following set of differential equations:(1)dSdt=Λ-β1I11+I1+β2I2S+γ-p2I2-μS-p1I1+ηR,dI1dt=β1SI11+I1-μ+α1+δ-p1I1,dI2dt=β2SI2-μ+α2+γ-p2I2,dRdt=δI1-η+μRwith the initial condition S0>0, I10>0, I20>0, and R0>0. Moreover to insure that the recruitment Λ in the susceptible compartment is always positive the following hypotheses are assumed to be holding always: (2)δ≥p1,γ≥p2.Theorem 1.
The closed setΩ=S,I1,I2,R∈R+4:N≤Λ/μ is positively invariant and attracting with respect to model (1).Proof.
LetSt,I1t,I2t,R(t) be any solution of system (1) with any given initial condition. Then by adding all the equations in system (1) we obtain that(3)dNdt=Λ-μS-μ+α1I1-μ+α2I2-μR≤Λ-μN.Thus, from standard comparison theorem [20], we obtain(4)Nt≤N0e-μt+Λμ1-e-μt.Consequently it is easy to verify that(5)Nt≤Λμ,when0≤Λμ.Thus, Ω is positively invariant. Further, when N0>Λ/μ, then either the solution enters Ω in finite time, or N(t) approaches Λ/μ as t→∞. Hence, Ω is attracting (i.e., all solutions in R+4 eventually approach, enter, or stay in Ω).Therefore, the system of equations given in model (1) is mathematically well-posed and epidemiologically reasonable, since all the variables remain nonnegative ∀t≥0. Further since the equations of model (1) are continuous and have continuously partial derivatives then they are Lipschitzian. In addition to that from Theorem 1, model (1) is uniformly bounded. Therefore the solution of it exists and is unique. Hence, from now onward it is sufficient to consider the dynamics of model (1) in Ω.


## 3. Equilibrium Points and Basic Reproduction Number

Model (1) has four equilibrium points that are obtained by setting the right hand sides of this model equal to zero. The first equilibrium point is the disease-free equilibrium (DFE) point that is denoted by E0=(S0,0,0,0) with S0=Λ/μ. Moreover the basic reproduction number of model (1), which is denoted by R0, is the maximum eigenvalue of the next generation matrix (i.e., the maximum of the reproduction numbers, those computed of each disease). That is,(6)R0=max⁡R1,R2.Here R1=β1S0+p1/(μ+α1+δ) and R2=β2S0+p2/μ+α2+γ.The other three equilibrium points can be described as follows.The first disease-free equilibrium point, which is located in the boundarySI2-plane, is denoted by E1=S¯,0,I¯2,0 where(7)S¯=S0T2,I¯2=Λμ+α21-1T2,and here T2=β2S0/(μ+α2+γ-p2). Clearly E1 exists uniquely in the interior of SI2-plane provided that(8)T2>1.The second disease-free equilibrium point that is located in the boundary SI1R-space is given by E2=S~,I~1,0,R~ where(9)S~=S01+I~1T1;I~1=Λ1-1/T1μ+α1+δ-ηδ/η+μ+Λ/T1,R~=I~1η+μ,and here T1=β1S0/μ+α1+δ-p1. Obviously E2 exists uniquely in the interior of positive octant of SI1R-space provided that(10)T1>1.Finally, the endemic equilibrium point, which is denoted by E3=(S∗,I1∗,I2∗,R∗), where(11)S∗=S0T2;I1∗=T1T2-1;R∗=δT1/T2-1η+μ,I2∗=Λμ+α21-1T2-T1/T2-1μ+α2μ+α1+δ-ηδη+μ,exists uniquely in the interior of Ω provided that the following conditions hold:(12)T1>T2>1,Λ>ΛT2+I1∗μ+α1+δ-ηδη+μ.Keeping the above in view, it is easy to verify with the help of condition (2) that(13)Ti>1⁡Ti<1⟺Ri>1⁡Ri<1,i=1,2.Then directly we obtain Ti>1Ti<1⇔R0>1(R0<1). Consequently, Ti represent the threshold parameters for the existence of the last three equilibrium points of model (1). Moreover, it is well known that the basic reproduction number (R0) is representing the average number of secondary infections that occur from one infected individual in contact with susceptible individuals. Therefore if R0<1, then each infected individual in the entire period of infectivity will produce less than one infected individual on average, which shows the disease will be wiped out of the population. However, if R0>1, then each infected individual in the entire infection period having contact with susceptible individuals will produce more than one infected individual; this leads to the disease invading the susceptible population.


## 4. Local Stability Analysis

In this section, the local stability analyses of all possible equilibrium points of model (1) are discussed by determining the Jacobian matrix with their eigenvalues. Now the general Jacobian matrix of model (1) can be written:(14)
J
=
-
β
1
I
1
1
+
I
1
-
β
2
I
2
-
μ
-
β
1
S
1
+
I
1
2
-
p
1
-
β
2
S
+
γ
-
p
2
η
β
1
I
1
1
+
I
1
-
β
1
S
I
1
1
+
I
1
2
+
f
S
,
I
1
0
0
β
2
I
2
0
β
2
S
-
μ
+
α
2
+
γ
-
p
2
0
0
δ
0
-
η
+
μ
,where fS,I1=β1S/(1+I1)-μ+α1+δ-p1. Therefore the local stability results near the above equilibrium points can be presented in the following theorems.Theorem 2.
The disease-free equilibriumE0=(S0,0,0,0) is locally asymptotically stable when R0<1 and unstable for R0>1.Proof.
The characteristic equation of the Jacobian matrix of model (1) at the disease-free equilibrium can be written as(15)λ+μλ-β1S0+μ+α1+δ-p1λ-β2S0+μ+α2+γ-p2λ+η+μ=0.So, if R0<1, then according to (6), (15) has four negative real roots (eigenvalues). Hence, the DFE is locally asymptotically stable. Further, for R0>1 (15) has at least one positive eigenvalue and then the DFE is a saddle point.Theorem 3.
The first disease-free equilibrium pointE1=(S¯,0,I¯2,0) of model (1) is locally asymptotically stable provided that(16)T2>1>T1.Proof.
The characteristic equation of the Jacobian matrix of model (1) at E1 can be written as(17)λ--β1S¯+μ+α1+δ-p1λ-+η+μλ-2-Tr1λ-+D1=0;here Tr1=-(β2I¯2+μ)<0 and D1=β2I¯2(β2S¯-γ-p1)>0 due to condition (16). Hence both the eigenvalues λ-S and λ-I2, which describe the dynamics in the S-direction and I2-direction, respectively, have negative real parts. Moreover, from (17), the eigenvalue in the I1-direction can be written as(18)λ-I1=β1S¯-μ+α1+δ-p1=1-T2T1β1S¯.Thus under the given condition (16), we have λ-I1<0, while λ-R=-(η+μ) is always negative. Hence, E1 is locally asymptotically stable.Theorem 4.
The second disease-free equilibrium pointE2=(S~,I~1,0,R~) of model (1) is locally asymptotically stable provided that(19a)T1>1>T21+I~1,(19b)P~2>δ,where P~2 is given in the proof.Proof.
The characteristic equation of the Jacobian matrix of model (1) at E2 can be written as(20)λ~-β2S~+μ+α2+γ-p2λ~3+Aλ~2+Bλ~+C=0;here(21)A=P~1+P~2+η+μ,B=P~1+P~2η+μ+P~1P~2+P~3,C=η+μP~1P~2+P~3-P~4with(22)P~1=β1I~11+I~1+μ>0,P~2=β1S~I~11+I~12>0,P~3=β1I~11+I~1β1I~11+I~12+p1>0,P~4=ηδβ1I~11+I~1>0.Clearly, the eigenvalue λ~I2 in the I2-direction can be written as(23)λ~I2=β2S~-μ+α2+γ-p2=1-T1T21+I~1β2S~,and thus λ~I2<0 under the condition (19a). In addition from (20) we have A>0 always, while C can be written as(24)C=p~5p~2δ+μp~5p~2δ+p~2β2I2∗η+μ+p~2μη+μ+p~3η+μ+p~4p~2+η+μ-p~5.Hence, C>0 provided that the sufficient condition (19b) holds. Further it is easy to verify that(25)AB-C=P~1+P~2B+η+μ2+P~4>0.Hence, due to the Routh-Hurwitz criterion the third-degree polynomial term in (20) has roots (eigenvalues) with negative real parts. Hence E2 is locally asymptotically stable.Theorem 5.
The endemic equilibrium pointE3=(S∗,I1∗,I2∗,R∗) of model (1) is locally asymptotically stable provided that(26a)T1>T2>1,(26b)Q̿>q5>δ,where Q̿ and q5 are given in the proof.Proof.
The characteristic equation of the Jacobian matrix of model (1) at E3 can be written as(27)λ∗4+A1λ∗3+B1λ∗2+C1λ∗+D1=0.Here(28)A1=q1+q2+η+μ>0,B1=η+μq1+q2+q1q2+q3+q4>0,C1=q1q2+q3η+μ+q4q2+η+μ-q5,D1=q4q2η+μ>0with(29)q1=β1I1∗1+I1∗+β2I2∗+μ;q2=β1S∗I1∗1+I1∗2;q3=β1I1∗1+I1∗β1S∗1+I1∗2+p1,q4=β2I2∗β2S∗-γ-p2;q5=ηδβ1I1∗1+I1∗.Obviously, qi>0, i=1,2,3,5, while q4 is positive under condition (26a). Now, by using the values of qi and the sufficient condition (26b), then straightforward computation gives (30)C1=Q̿+q4q2+η+μ-q5>0,and here Q̿=q1q2+q3(η+μ). Moreover we have(31)A1B1C1=Q+q1+q22η+μ+q4q1+q1q2+q3+η+μ2q1+q2Q-q5,C12+A12D1=Q-q52+q1+q2+η+μ2q4q2η+μ,where Q=Q̿+q4q2+η+μ. Therefore we obtain that(32)A1B1C1-C12-A12D1=FQ̿+q4Q̿q1+q2+q4q1A1η+μ2+q42q1q2+η+μ+q4q2q1+q2q1q2+q3-Fq5+C1q5.Here F=q1+q22η+μ+q1q2+q3+η+μ2(q1+q2)+q4q1.
Hence, according to condition (26b) it is easy to verify that A1B1C1-C12-A12D1>0. Therefore, all the coefficients of (27) are positive and A1B1C1-C12-A12D1>0. Hence, due to the Routh-Hurwitz criterion all the eigenvalues (λS∗,λI1∗,λI2∗andλR∗) of the Jacobian matrix near the endemic equilibrium point E3 have negative real parts. Thus, the proof is complete.


## 5. Global Stability Analysis

This section deals with the global stability of the equilibrium points of model (1) using Lyapunov methods with LaSalle’s invariant principle. The obtained results are presented in the following theorems.Theorem 6.
Assume thatDFEE0=(S0,0,0,0) of model (1) is locally asymptotically stable; then it is global asymptotically stable in Ω.Proof.
ConsiderV^:Ω→R that is defined by(33)V^S,I1,I2,R=∑i=12Ii.Computing the derivative of this positive semidefinite function with respect to time along the solution of model (1) and then simplifying the resulting terms give(34)dV^dt=β1S1+I1-μ+α1+δ-p1I1+β2S-μ+α2+γ-p2I2.Since the solution of model (1) is bounded by S0=Λ/μ as t→∞, (35)dV^dt≤β1S0-μ+α1+δ-p1I1+β2S0-μ+α2+γ-p2I2=μ+α1+δ-p1T1-1I1+μ+α2+γ-p2T2-1I2.Since Ti<1,i=1,2, due to the local stability condition of E0 then dV^/dt<0. Also we have that dV^/dt=0 on the set S,I1,I2,R∈Ω:I1=I2=0, so dV^/dt is negative semidefinite and hence according to Lyapunov first theorem E0 is globally stable point. Now, since on this set we have (36)dSdt=Λ-μS+ηR=0,if and only if S=S0,R=0, thus the largest invariant set contained in this set is reduced to the disease-free equilibrium point E0. Hence according to LaSalle’s invariant principle [18], E0 is attractive point and hence it is globally asymptotically stable in Ω.Theorem 7.
Assume that the first disease-free equilibrium pointE1=(S¯,0,I¯2,0) is locally asymptotically stable; then it is global asymptotically stable in Ω provided that (37)I2>I¯2.Proof.
Consider thatL:Ω→R that is defined by(38)LS,I1,I2,R=12S-S¯+I1+I2-I¯2+R2+2μ+α1β1I1+2μ+α2β2I2-I¯2-I¯2ln⁡I2I¯2+2μ+α12δR2.Clearly L is continuous and positive definite function. Now by taking the derivative of L with respect to time along the solution of model (1), we get after simplifying the resulting terms that(39)dLdt=-μS-S¯2-μR2-2μRS-S¯-2μ+α1S-S¯I1-2μ+α2S-S¯I2-I¯2-2μ+α1+α2I2-I¯2I1-2μ+α1I1R-μ+α2I2-I¯22-2μ+α2I2-I¯2R-μ+α1I12,dLdt≤-μS-S¯+R2-μI1+I2-I¯22-2μI2-I¯2R+2μμ+α1+δ-p1β1T1T2-1I1-μη+μδR2.Hence according to local stability condition (16) along with the sufficient condition (37) it obtains that dL/dt is negative definite function. Thus due to Lyapunov second theorem E1 is global asymptotically stable in Ω.Theorem 8.
Assume that the second disease-free equilibrium pointE2=S~,I~1,0,R~ of model (1) is locally asymptotically stable; then it is global asymptotically stable in Ω if (40a)q122<q11q22,(40b)q142<q11q44,(40c)q232<q22q33,(40d)q342<q33q44,where (41)q11=μ;q22=μ+α1+S~1+I~11+I1;q33=μ+α2;q44=2μ+α1η+μ2δ+μ;q12=11+I1-2μ+α1;q14=2μ;q23=2μ+α1+α2;q34=2μ+α2.Proof.
Consider the functionW:Ω→R that is defined by(42)WS,I1,I2,R=12S-S~+I1-I~1+I2+R-R~2+1β1I1-I~1-I~1ln⁡I1I~1+2μ+α2β2I2+2μ+α12δR-R~2.Clearly W is continuous and positive definite function. Now by taking the derivative of W with respect to time along the solution of model (1), we get after simplifying the resulting terms that(43)dWdt=-q11S-S~2+q12S-S~I1-I~1-q14S-S~R-R~-q22I1-I~12-q23I1-I~1I2-q33I22-q34R-R~I2-q44R-R~2+2μμ+α2+γ-p2β2T21+I~1T1-1I2.Now by using the given conditions (40a)–(40d) we get that(44)dWdt≤-q112S-S~-q222I1-I~12-q222I1-I~1+q332I22-q112S-S~+q442R-R~2-q442R-R~+q332I22+2μμ+α2+γ-p2β2T21+I~1T1-1I2.Hence according to local stability condition (19a) it obtains that dW/dt is negative definite function. Thus due to Lyapunov second theorem E1 is global asymptotically stable in Ω.Theorem 9.
Assume that the endemic equilibrium pointE3=(S∗,I1∗,I2∗,R∗) of model (1) is locally asymptotically stable; then it is global asymptotically stable in Ω if (45a)α142<α11α44,(45b)α122<α11α22,(45c)α242<α22α44,where(46)α11=β1I11+I1+β2I2+μ,α22=β1S∗,α44=η+μ,α12=β11+I1-β1S∗1+I11+I1∗+p1,α24=δ,α14=η.Proof.
Consider the functionV:Ω→R that is defined by(47)VS,I1,I2,R=S-S∗22+I1-I1∗-I1∗ln⁡I1I1∗+2μ+α2β2I2-I2∗-I2∗ln⁡I2I2∗+R-R∗22.Clearly the function V is continuous and positive definite function. By taking the derivative of V with respect to time along the solution of model (1), we get after simplifying the resulting terms that (48)dVdt=-α11S-S∗2-α22I1-I1∗2-α44R-R∗2+α12S-S∗I1-I1∗+α14R-R∗S-S∗+α24I1-I1∗R-R∗.Then by using the given conditions (45a)–(45c) we obtain that(49)dVdt<-α11S-S∗+α22I1-I1∗2-α11S-S∗+α44R-R∗2-α44R-R∗+α22I1-I1∗2.Hence, dV/dt is negative semidefinite, and dV/dt=0 on the set S,I1,I2,R∈Ω:S=S∗,I1=I1∗,I2>0,R=R∗, so according to Lyapunov first theorem E3 is globally stable point. Further, since on this set we have(50)dSdt=Λ-β1I1∗1+I1∗+β2I2S∗+γ-p2I2-μS∗-p1I1∗+ηR∗=0,if and only if I2=I2∗, then the largest compact invariant set contained in this set is reduced to the endemic equilibrium point E3. Hence according to LaSalle’s invariant principle [18], E3 is attractive point and hence it is globally asymptotically stable in Ω.


## 6. Bifurcation Analysis

In this section the local bifurcations near the equilibrium points of model (1) are investigated as shown in the following theorems with the help of Sotomayor theorem [20]. Note that model (1) can be rewritten in a vector form dX/dt=f(X), where X=S,I1,I2,RT and f=f1,f2,f3,f4T with fi,i=1,2,3,4, are given in the right hand side of model (1). Moreover, straightforward computation gives that the general second derivative of the Jacobian matrix (14) can be written:(51)
D
2
f
X
,
β
U
,
U
=
-
2
β
1
ζ
1
ζ
2
1
+
I
1
2
-
2
β
2
ζ
1
ζ
3
+
2
β
1
S
ζ
2
2
1
+
I
1
3
2
β
1
ζ
1
ζ
2
1
+
I
1
2
-
2
β
1
S
ζ
2
2
1
+
I
1
3
2
β
2
ζ
1
ζ
3
0
,where β is any bifurcation parameter and U=ζ1,ζ2,ζ3,ζ4T is any eigenvector.Theorem 10.
Assume thatT2<1; then as T1 passes through the value T1=1, model (1) near the disease-free equilibrium E0 has(1)
no saddle-node bifurcation;(2)
a transcritical bifurcation;(3)
no pitchfork bifurcation.Proof.
SinceT1=β1S0/(μ+α1+δ-p1)=1; then β1=(μ+α1+δ-p1)/S0=β1∗. Now straightforward computation shows that the Jacobian matrix of model (1) at E0 with β1=β1∗ has zero eigenvalue (λI1=0) and can be written as follows:(52)J0=-μ-β1∗S0-p1-β2S0+γ-p2η000000β2S0-μ+α2+γ-p200δ0-η+μ.Let V=v1,v2,v3,v4T be the eigenvector corresponding to λI1=0. Thus J0V=0 gives(53)V=av2v20bv2,where v2 is any nonzero real number, a=ηδ/η+μ-μ+α1+δ/μ, and b=δ/η+μ.
Similarly,W=w1,w2,w3,w4T represents the eigenvector corresponding to eigenvalue λI1=0 of J0T. Hence J0TW=0 gives that(54)W=0w200,and here w2 is any nonzero real number. Now, since (55)dfdβ1=fβ1X,β1=-I1S1+I1,I1S1+I1,0,0T,thus fβ1E0,β1∗=(0,0,0,0)T, which gives WTfβ1E0,β1∗=0.
Thus, according to Sotomayor’s theorem for local bifurcation, model (1) has no saddle-node bifurcation near DFE at β1=β1∗.
Now since(56)Dfβ1E0,β1∗=0-S0000S00000000000,then, WTDfβ1E0,β1∗V=S0v2w2≠0. Now, by substituting E0 and β1∗ in (51) we get(57)D2fE0,β1∗·V,V=-2β1∗v1v2-2β2v1v3+2β1∗S0v222β1∗v1v2-2β1∗S0v222β2v1v30.Therefore,(58)WTD2fE0,β1∗·V,V=2β1∗v1v2-2β1∗S0v22w2=2β1∗v22a-S0w2≠0.So, according to Sotomayor’s theorem model (1) has a transcritical bifurcation at E0 with parameter β1=β1∗≈(T1=1), while the pitchfork bifurcation cannot occur.Note that similar results as those of Theorem10 are obtained at T2=1 or β2=β2∗=μ+α2+γ-p2/S0.Theorem 11.
Assume thatT1=T2; then model (1) near the first disease-free equilibrium point E1=(S¯,0,I¯2,0) has(1)
no saddle-node bifurcation;(2)
a transcritical bifurcation;(3)
no pitchfork bifurcation.Proof.
There are two cases; in the first case, it is assumed thatT1=T2=1; then straightforward computation shows that JE1=JE0; that is, I-2=0. So (by Theorem 10) model (1) has no bifurcation and then the proof is complete.
Now in the second case it is assumed thatT1=T2≠1 or equivalently β1=β1∗∗=β2(μ+α1+δ-p1)/(μ+α2+γ-p2). So straightforward computation shows that the Jacobian matrix of model (1) at E1 with β1=β1∗∗ has zero eigenvalue (λI1=0) and can be written as follows: (59)J1=-β2I-2-μ-β1∗∗S¯-p1-β2S¯+γ-p2η0000β2I-20β2S¯-μ+α2+γ-p200δ0-η+μ.Let X=x1,x2,x3,x4T be the eigenvector corresponding to λI1=0, which satisfies J1X=0, so we get(60)X=0x2a1x2bx2,where x2 is any nonzero real number, a1=-β1S-+p1-ηδ/η+μ/β1S--(γ-p2), and b=δ/η+μ.
Similarly the eigenvectorY=y1,y2,y3,y4T that is corresponding to the eigenvalue λI1=0 of J1T satisfies J1TY=0, so we get(61)Y=0y200,and here y2 is any nonzero real number. Now, since(62)dfdβ1=fβ1X,β1=-I1S1+I1,I1S1+I1,0,0Tthen by substituting the values of E1 and β1∗∗ we obtain that fβ1E1,β1∗∗=(0,0,0,0)T and hence we get that YTfβ1E1,β1∗∗=0.
Thus according to Sotomayor’s theorem for local bifurcation, model (1) has no saddle-node bifurcation near E1 at β1=β1∗∗. Now since(63)Dfβ1E1,β1=0-S¯000S¯0000000000,then YTDfβ1E1,β1∗∗X=S-x2y2≠0. Thus by substituting E1 and β1∗∗ in (51) we get(64)D2fE1,β1∗∗·X,X=-2β1∗∗x1x2-2β2x1x3+2β1∗∗S¯x222β1∗∗x1x2-2β1∗∗S¯x222β2x1x30.Therefore, YTD2fE1,β1∗∗·X,X=-2β1∗∗S-(x2)2y2≠0.
So, model (1) has a transcritical bifurcation at E1 with parameter β1=β1∗∗≈(T1=T2=1), while the pitchfork bifurcation cannot occur and hence the proof is complete.Theorem 12.
Assume that condition (19b) holds and let T1=T2(1+I~1); then model (1) near the second disease-free equilibrium point E2=(S~,I~1,0,R~) undergoes(1)
no saddle-node bifurcation;(2)
a transcritical bifurcation;(3)
no pitchfork bifurcation.Proof.
FromT1=T2(1+I~1) it is obtained that(65)β2=β2∗∗=β1μ+α2+γ-p2μ+α1+δ-p11+I~1,and then straightforward computation shows that the Jacobian matrix of model (1) at E2 with β2=β2∗∗ has zero eigenvalue (λI2=0) and can be written as follows:(66)
J
2
=
-
β
1
I
~
1
1
+
I
~
1
-
μ
-
β
1
S
~
1
+
I
~
1
2
-
p
1
-
β
2
∗
∗
S
~
+
γ
-
p
2
η
β
1
I
~
1
1
+
I
~
1
-
β
1
S
~
I
~
1
1
+
I
~
1
2
0
0
0
0
0
0
0
δ
0
-
η
+
μ
.Further the eigenvector L=l1,l2,l3,l4T that is corresponding to λI2=0 satisfies J2L=0, so we get(67)
L
=
S
~
1
+
I
~
1
l
2
l
2
-
E
β
2
∗
∗
S
~
-
γ
-
p
2
l
2
b
l
2
,where l2 is any nonzero real number and E=(β1S~/1+I~12)I~1+1+μS~/1+I~1+p1-ηb.
Similarly the eigenvectorK=k1,k2,k3,k4T that is corresponding to eigenvalue λI2=0 of J2T satisfies J2TK=0, so we get(68)K=00k30and here k3 is any nonzero real number. Now, since df/dβ2=fβ2X,β2=-I2S,I2S,0,0T, therefore fβ2E2,β2∗∗=(0,0,0,0)T, which yields KTfβ2E2,β2∗∗=0. Consequently according to Sotomayor’s theorem for local bifurcation, model (1) has no saddle-node bifurcation near E2 at β2=β2∗∗.
Now since(69)Dfβ2E2,β2=00-S~0000000S~00000then KTDfβ2E2,β2∗∗L=S~l3k3≠0. Now, by substituting E2 and β2∗∗ in (51) we get(70)
D
2
f
E
2
,
β
2
∗
∗
·
L
,
L
=
-
β
1
1
+
I
~
1
2
l
1
l
2
-
2
β
2
∗
∗
l
1
l
3
-
β
1
l
1
l
2
1
+
I
~
1
2
+
2
β
1
S
l
2
2
1
+
I
~
1
3
β
1
l
2
1
+
I
~
1
2
l
1
l
2
+
β
1
1
+
I
~
1
2
l
1
l
2
-
2
β
1
S
l
2
2
1
+
I
~
1
3
2
β
2
∗
∗
l
1
l
3
0
.Therefore, KTD2fE2,β2∗∗·L,L=2β2∗∗l1l3k3≠0. Thus model (1) undergoes a transcritical bifurcation at E2 with parameter β2=β2∗∗≈(T1=T21+I~1), while the pitchfork bifurcation cannot occur.Moreover, the following results are obtained too:(1)
AlthoughT1=T2≠1 gives β2=β`2∗∗, model (1) does not undergo any of the above types of bifurcation near the equilibrium point E1 with parameter β2=β`2∗∗.(2)
AlthoughT1=T21+I~1 gives β1=β`1∗∗, model (1) does not undergo any of the above types of bifurcation near the equilibrium point E2 with parameter β1=β`1∗∗.(3)
The determinant of the Jacobian matrix atE3, say J(E3), cannot be zero and hence it has no real zero eigenvalue. So there is no bifurcation near E3.Keeping the above in view, in the following theorem we detect of the possibility of having Hopf bifurcation.Theorem 13.
Assume that condition (26a) holds and let the following conditions be satisfied. Then model (1) undergoes Hopf bifurcation around the endemic equilibrium point when the parameter β2 crosses a critical positive value β`2∗,(71a)N2>0,N4<0,(71b)σ1>2σ2,(71c)δ<q5<q4q2,q4q1,(71d)q2>1,and here σ1=q12+q22+η+μ2 and σ2=q1+q2η+μ+2q1q2+2q3+2q4, while N2 and N4 are given in the proof.Proof.
It is well known that, in order for Hopf bifurcation in four-dimensional systems to occur, the following conditions should be satisfied [21, 22]:(1)
The characteristic equation given in (28) has two real and negative eigenvalues and two complex eigenvalues, say, λ∗β2=τ1β2±τ2β2.(2)
τ1β`2∗=0.(3)
d/dβ2τ1β2|β2=β`2∗≠0 (The transversality condition).Accordingly the first two points are satisfied if and only if (72)Δ2β`2∗=A1B1C1-C12-A12D1=0,A13-4Δ1>0,while the third condition holds provided that (73)ddβ2τ1β2β2=β`2∗=-Ψβ2Θβ2+Φβ2Γβ2Ψβ22+Φβ22≠0.That means Ψβ2Θβ2+Φβ2Γβ2≠0.
Here(74)Ψβ2=4τ13-12τ1τ22+3A1τ12-τ22+2B1τ1+C1,Φβ2=12τ12τ2-4τ23+6A1τ1τ2+2B1τ2,Θβ2=A1˙τ13-3τ1τ22+B1˙τ12-τ22+C˙1τ1+D1˙,Γβ2=A1˙3τ12τ2-τ23+2B1˙τ1τ2+C˙1τ2.Now, straightforward computation shows the condition (75)Δ2β`2∗=0givesthatN1β23+N2β22+N3β2+N4=0,where(76)N1=q2η+μη+μ+q2+η+μμ+α2η+μ+2q2+μ+α2I2∗3>0,N2=η+μq2η+μ3q6+2+η+μ+q3η+μ+2μ+α2+2q2+μ+α2q62η+μ+μ+α2+q2η+μ+4q6+q2+q223q6+q2+μ+α2η+μ2-q5+μ+α2q221+q2+q2q3+q6μ+α2-q2q5I2∗2,N3=η+μ2q23+q32q6+2q2+η+μ+3q6q2-q5+q6+q2q6μ+α2+η+μq6q223q6+2q2+2q6+q22q2q3+q2η+μ2+μ+α21+q6η+μ2+2q6+q2q6q2+q3-2q6q5-q2q5+q2μ+α2q6q21+q2+q6+q2q3+1-q22q6+q2+q3q5I2∗,N4=η+μ2q6q2q62+q6q2+q22+q3q6q6+2q2+q6+q2η+μq6q2+q3+q3q22-q2+q6+q2η+μq62q22+2q6q2q3+q32+η+μq3-q6q2-q6η+μq6+η+μ-q22η+μ-q6+q2q6q2+q3q5-q52with q1,q2,q3,q4, and q5 given in (28) and q6=q1-β2I2∗. Clearly, from condition (71a) there is unique positive root, say, β2=β`2∗. Consequently by using Δ2=0 in the characteristic equation and then doing some algebraic computation we get four roots,(77)λ1,2β`2∗=±iC1A1=τ2β`2∗,λ3,4=12-A1±A12-4Δ1A1.Now, it is easy to verify that λ3 and λ4 are real and negative provided that (71b).
Further forβ2∈(β`2∗-ϵ,β`2∗+ϵ) the general form of complex eigenvalues can be written as(78)λ1=τ1β2+iτ2β2,λ2=τ1β2-iτ2β2.Substituting λβ2=τ1β2+iτ2(β2) into characteristic equation and after that calculating the derivative with respect to β2 and then comparing the real and imaginary parts give that(79)Ψβ2τ1˙β2-Φβ2τ2˙β2=-Θβ2,Φβ2τ1˙β2+Ψβ2τ2˙β2=-Γβ2.Moreover, by solving the above linear system for τ1˙(β2) and τ2˙(β2) then we get that(80)Ψβ2Θβ2+Φβ2Γβ2=q1q2q3μ+α2+2η+μ3+μ+α2η+μη+μ+3q2+3q1-1+q2q3q2μ+α2+η+μμ+α2+2q2+2q1-1+q1η+μμ+α2η+μ+2q3+η+μ2+q4+η+μ3q2μ+α2+q22+q3+q4+q5μ+α22q2+η+μ+q1q3η+μ2+q2η+μμ+α2q2q2-1+η+μ2q2+2q1-1+q1μ+α2q1q22+q23+q2q4-q5+q1η+μ2q1q21+q2+2q23+q2q4-q5+q2q4-q5q1q2+q22+η+μ23q1q22+q23+q1q4-q5.Thus it is easy to verify that Ψβ2Θβ2+Φβ2Γβ2≠0 provided that (71c) and (71d). Thus, the proof is complete.


## 7. Numerical Simulations

In this section, the global dynamics of model (1) is investigated numerically for different sets of initial values and different sets of parameters values. The objectives of such investigation are to determine the effect of varying the parameters values and confirm our obtained results. It is observed that, for the following biologically feasible set of hypothetical parameters values(81)Λ=20,β1=0.75,β2=0.1,γ=0.75,p1=0.01,p2=0.01μ=0.3,η=0.5,α1=0.1,α2=0.6,δ=0.7,the solution of model (1) approaches asymptotically to the endemic equilibrium point E3=(16.4,10.28,9.18,8.99) as shown in Figure 1, started from different sets of initial points.Figure 1
Globally asymptotically stable positive equilibrium point of model (1) for the parameters set (81), started from different sets of initial point.
(a)
(b)
(c)
(d)Clearly Figure1 confirms our obtained analytical results regarding existence of a globally asymptotically stable positive equilibrium point when the parameters values are satisfying R1>R2>1. On the other hand, model (1) for the following set of hypothetical data approaches asymptotically to the DFE as shown in Figure 2:(82)Λ=20,β1=0.05,β2=0.05,γ=0.6,p1=0.01,p2=0.03,η=0.5,μ=0.9,α1=0.1,α2=0.1,δ=0.3.It is easy to verify that for the data (82) we have R0=0.86<1, and the solution approaches to E0=(22.22,0,0,0).Figure 2
Time series of the solution of model (1) that approaches asymptotically to DFE for the data (82).Now in order to investigate the effect of varying one parameter value at a time on the dynamical behavior of model (1), the following results are observed.(i)
Varying of the parameters values(Λ,γ,p2,η,α1,α2,δ) does not affect the dynamical behavior of model (1); that is, the system still approaches to coexistence equilibrium point.(ii)
For the data (81) with β1≤0.015, the solution of model (1) approaches asymptotically to E1=(S,0,I2,0) in the interior of positive quadrant of SI2-plane with R1<1<R2. However for 0.015<β1<0.066 the solution of the model still approaches asymptotically to E1 even when 1<R1<R2. Finally when β1≥0.066 the solution of model (1) approaches to coexistence equilibrium point with 1<R2<R1 as shown in Figure 3.(iii)
Similar results are obtained in case of varying the parameterβ2 keeping the rest of parameters values in (81) fixed. In fact for β2≤0.024, we have R2<1<R1, and the solution of model (1) approaches asymptotically to E2=(S,I1,0,R). However for 0.024<β2<0.06, we have 1<R2<R1, and it is observed that the solution of model (1) approaches asymptotically to E2=(S,I1,0,R) too, while when β2≥0.06 the solution of model (1) approaches asymptotically to E3=(S,I1,I2,R) as shown in Figure 4.(iv)
Now decreasing the parameterμ, keeping the rest of the parameters values in (81) fixed, gives similar dynamical behavior as that of varying β2. Further it is observed that when μ>0.9, for which R2<1<R1, the solution of model (1) approaches asymptotically to E2=(S,I1,0,R); however when 0.57<μ<0.9, we have 1<R2<R1, and the solution of model (1) still approaches to E2=(S,I1,0,R). Finally, for μ≤0.57 that satisfies 1<R2<R1, the solution of model (1) approaches asymptotically to E3=(S,I1,I2,R). Clearly this confirmed our obtained existence conditions (10) and (12) as well as stability conditions of these points.(v)
Finally, varying the parameterp1, keeping the rest of the parameters values in (81) fixed, showed that, for p1≥0.6, that satisfies 1<R2<R1, the solution of model (1) approaches asymptotically to E2=(S,I1,0,R); however when p1<0.6, which satisfies 1<R2<R1 too, the solution of model (1) approaches asymptotically to E3=(S,I1,I2,R) as shown in Figure 5.Figure 3
Time series of the solution of model (1) for the data given by (81). (a) For β1=0.015 the model approaches to E1=(16.4,0,16.75,0). (b) For β1=0.06 the model approaches to the same point. (c) For β1=0.3 the model approaches to E3=(16.4,3.51,14.1,3.07).
(a)
(b)
(c)Figure 4
Time series of the solution of model (1) for the data given by (81). (a) For β2=0.024 the model approaches to E2=(27.33,17.80,0,15.5). (b) For β2=0.059 the model approaches to the same point. (c) For β2=0.065 the model approaches to E3=(25.2,16.3,1.76,14.3).
(a)
(b)
(c)Figure 5
Time series of the solution of model (1) for the data given by (81). (a) For p1=0.6 the model approaches to E2=(15.9,22.9,0,20.08). (b) For p1=0.5 the model approaches to E3=(16.4,19.5,2.40,17.06).
(a)
(b)


## 8. Conclusion

In this paper, we proposed and analyzed an epidemic model involving vertical and horizontal transmission of infection with nonlinear incidence rate. It is assumed that the rates of infectionsp1, p2 are less than the recovery rates δ and γ, respectively. According to the diseases in model (1) the population is divided into four subclasses: susceptible individuals that are represented by S(t), infected individuals for SIRS-type of disease that are represented by I1(t), infected individuals for SIS-type of disease that are represented by I2(t), and recovery individuals that are denoted by R(t). The boundedness and invariant of the model are discussed. The basic reproduction number of the model and the associated threshold parameter values, namely, Ti,i=1,2, are determined. It is observed that if the basic reproduction number is less than unity then the diseases are eradicated from the model. The competitive exclusion principle occurred in model (1) such that only the second disease-free equilibrium point appeared in case of R2<1<R1. However only the first disease-free equilibrium point appeared in case of R1<1<R2. Finally the coexistence of both the diseases occurred in case of 1<R2<R1 and the sufficient condition (12) holds. The dynamical behavior of model (1) has been investigated locally as well as globally using Routh-Hurwitz criterion and Lyapunov function, respectively. The local bifurcations of model (1) and the Hopf bifurcation around the endemic equilibrium point are studied. Finally to understand the effect of varying each parameter on the global dynamics of system (1) and to confirm our obtained analytical results, model (1) has been solved numerically and the following results are obtained for the set of hypothetical parameters values given by (81).(1)
Model (1) approaches asymptotically to a globally asymptotically stable point E3=(16.4,10.28,9.18,8.99).(2)
Varying one of the parameters values(Λ,γ,p2,η,α1,α2,δ) at a time keeping other parameters fixed has no effect on the dynamical behavior of the model.(3)
As the infection rate of the first disease (β1) decreases keeping other parameters fixed as in (81) the solution of model (1) approaches asymptotically to the equilibrium point E1=(16.4,0,16.75,0). However in case of increasing this parameter the model is still globally asymptotically stable in the interior of R+4.(4)
As the infection rate of the second disease (β2) decreases keeping other parameters fixed as in (81) the solution of model (1) approaches asymptotically to the equilibrium point E2=(27.33,17.80,0,15.5). However in case of increasing this parameter the model is still globally asymptotically stable in the interior of R+4.(5)
As the mortality rate (μ) increases keeping other parameters fixed as in (81) the solution of model (1) approaches asymptotically to the equilibrium point E2=(13.8,5.16,0,2.5) and when μ decreases the model is still globally asymptotically stable in the interior of R+4. Further, it is observed that p1 has the same effect as μ on the dynamical behavior of model (1).(6)
For the parameter set given in (82) the solution of model (1) approaches asymptotically to DFE.(7)
Finally, although for our selected parameters values model (1) does not undergo periodic dynamics, the model still has possibility to have periodic dynamics for other sets of parameters, especially Hopf bifurcation existing analytically.



---

*Source: 4907964-2016-01-05.xml*