# A Fractional-Order Model for Zika Virus Infection with Multiple Delays

Authors: R. Rakkiyappan; V. Preethi Latha; Fathalla A. Rihan
Journal: Complexity
 (2019)

---

## Abstract

Time delays and fractional order play a vital role in biological systems with memory. In this paper, we propose an epidemic model for Zika virus infection using delay differential equations with fractional order. Multiple time delays are incorporated in the model to consider the latency of the infection in a vector and the latency of the infection in the infected host. We investigate the necessary and sufficient conditions for stability of the steady states and Hopf bifurcation with respect to three time delaysτ1, τ2, and τ3. The model undergoes a Hopf bifurcation at the threshold parameters τ1∗, τ2∗, and τ3∗. Some numerical simulations are given to show the effectiveness of obtained results. The numerical simulations confirm that combination of fractional order and time delays in the epidemic model effectively enriches the dynamics and strengthens the stability condition of the model.


---

## Body

## 1. Introduction

Zika infection is a mosquito-borne disease, transmitted to humans through the bite of an infectedAedes mosquito. It was first discovered in Uganda in 1947 in rhesus monkey. The first human cases were reported in Nigeria in 1954. Zika was thought to cause mild symptoms in humans, including mild fever, skin rashes, conjunctivitis, muscle and joint pain, and headache, which lasts for three to twelve days normally. However, the World Health Organization (WHO) has concluded that Zika virus infection during pregnancy is also a cause of congenital brain abnormalities, including microcephaly [1]. Moreover, Zika virus is a trigger of Guillain–Barre syndrome [2]. There is no doubt that mathematical modeling of Zika infection plays an important role in gaining understanding of transmission of disease and to predict the behaviour of any outbreak [3, 4].Recently, mathematical modeling of dynamics of infectious diseases, using differential equations with memory (time-delay terms or fractional orders), has attracted much attention of many researchers (see, e.g., [5] and references therein). Time delay in models of population dynamics and in particular in macroscopic models of the immune response are natural and common [6]. Naturally, time delay or memory is an unavoidable factor in dynamics of most real-life phenomena. Time delay has influence on dynamical behaviours of biological systems in various aspects. Therefore, considering time delays in the investigation of biological systems is significant in both theoretical and practical point of views. In fact, when immune system works against the non-self-cells, it may take some time (time lag) to interact with the pathogen. Therefore, time delays cannot be ignored in models for immune response. Accordingly, the analysis of dynamical properties of system with time delays is important (see [5, 7–12]). Dengue fever is analyzed in [13], using a system of four nonlinear differential equations with two time delays. In [12], the authors considered the vector-borne epidemic model with time delay. The authors intensively discussed the impact of time delay in the host-to-vector transmission term that can destabilize the system. Periodic solutions can also be raised through Hopf bifurcation.In the existing literature, most of the biological problems are studied through the integer-order mathematical modeling by using ordinary, partial, and delay differential equations [9, 10, 14]. In the last few decades, fractional-order models have been incorporated in several areas of science, engineering, applied mathematics, economics, and bioengineering [15–20]. One advantage of the fractional-order differential equation is that they provide a powerful instrument for incorporation of memory and hereditary properties of the systems as opposed to the integer-order models, where such effects are neglected or difficult to incorporate. In addition, when fitting data, the fractional models have one more degree of freedom than the integer-order model (see [21]). Based on these advantages, some authors have developed interesting applications to investigate the dynamics of such fractional-order models with systems of memory [22–26]. In [5, 22], the authors studied fractional-order cancer immune systems. In [25], a fractional-order model for HIV with nonlinear incidence has been considered and stability for various equilibrium points has also been discussed. The authors in [27] investigated the dynamics of Ebola virus with time delay and fractional order and reported that combination of time delay and fractional order can effectively enrich the dynamics and strengthen the stability condition of the infection model. Analysis and dynamics of Zika transmission have been examined by many researchers (see, e.g., [3, 28, 29]. In [3], a mathematical model for transmission of Zika virus has been proposed with control measures of Zika virus. Stability properties of the Zika infection model have been investigated in [30]. The authors in [31] have compared the Zika infection model with dengue to show effect of the virus on population. The dynamical analysis of the SIS model is studied by considering bifurcation parameters in [32]. The authors [33] have discussed absence and presence of diffusion in the Zika virus disease model. The stability analysis and Hopf bifurcation point for various generalized epidemic models have been discussed in the literature [33–35]. However, the dynamics of fractional order with multiple time-delay models for Zika virus infection has not been yet studied in mathematical epidemiology.Herein, we demonstrate that a nonlinear fractional-order differential equations model, with multiple time delays, can simulate the dynamics of Zika virus infection much more than the classical epidemic models. The application of fractional derivatives is in several cases justified because they provide a better model than integer-order derivative models do [36, 37]. One important feature of fractional derivatives is that they are nonlocal opposed to the local behaviour of integer derivatives. In this way, the next state of a fractional system depends not only upon its current state but also upon all of its historical states [38–40].Motivated by the above discussion, in this paper, we investigate the dynamics of Zika virus infection with fractional order and time delays. In Section2, we formulate the model and study the nonnegativity of the solutions. In Section 3, we investigate the asymptotic stability analysis and Hopf bifurcation properties by taking time-delay parameters as bifurcation parameters. Sufficient conditions are derived to ensure the asymptotic stability and Hopf bifurcation behaviours of the addressed model. Finally, some numerical simulations are provided with various fractional orders and time delays to demonstrate the effectiveness of our theoretical findings in Section 4. We then conclude in Section 5.Before we start analysis, we provide some useful preliminaries.


### 1.1. Preliminaries

Herein, we provide some basic definitions and properties of integration and differentiation with fractional-order (free order)α (see [41]).Definition 1.
Letα∈0,∞, the operator Iaα on L1a,b is defined by(1)Iaαft=1Γα∫att−sα−1fsds,f∈L1a,b,t∈a,b,which is called the fractional integral (or Riemann–Liouville integral) of order α, where Ia0=Id is the identity operator.Definition 2.
Letα∈0,∞ and n=α, where x=mink∈ℤ:k≥x, and the operator DRLaα is defined for f∈L1a,b by(2)DRLaαft=1Γn−αddtn∫att−sn−α−1fsds,which is called the Riemann–Liouville fractional derivative of order α.Definition 3.
Letα∈0,∞ and f is such that Ian−αfn exists, where n=α,f∈Ana,b (the set of all function f:a,b⟶ℝ provided that fn−1 be absolutely continuous), then we define the operator DCaα by(3)DCaαft=1Γn−α∫att−sn−α−1fnsds,which exists for almost everywhere x∈a,b. The operator DCaαft is called the Caputo fractional derivative of order α. In particular, when 0<α≤1, we have(4)DCaαft=1Γ1−α∫atf′st−sαds.Remark 1.
Letβ,γ∈ℝ+ and α∈0,1. Then,(i)
IfIaβ:L1⟶L1 and if ft∈L1, then IaβIaγft=Iaβ+γft(ii)
limβ⟶nIaβfx=Ianft uniformly on a,b, n=1,2,3,…, where Ia1ft=∫atfsds(iii)
limβ⟶0Iaβft=ft weakly(iv)
Ifft is absolutely continuous on a,b, then limα⟶1Daαft=dft/dt(v)
Thus,Daαft=d/dtIa1−αft (Riemann–Liouville sense) and Dαft=Ia1−αd/dtft (Caputo sense)Remark 2.
We notice that the fractional derivatives involve an integration and are nonlocal operators, which can be used for modeling systems with memory.
We should mention here that Caputo’s definition of fractional derivative is a modification of the Riemann–Liouville definition and has the advantage of dealing with initial value problems in a proper way.


## 1.1. Preliminaries

Herein, we provide some basic definitions and properties of integration and differentiation with fractional-order (free order)α (see [41]).Definition 1.
Letα∈0,∞, the operator Iaα on L1a,b is defined by(1)Iaαft=1Γα∫att−sα−1fsds,f∈L1a,b,t∈a,b,which is called the fractional integral (or Riemann–Liouville integral) of order α, where Ia0=Id is the identity operator.Definition 2.
Letα∈0,∞ and n=α, where x=mink∈ℤ:k≥x, and the operator DRLaα is defined for f∈L1a,b by(2)DRLaαft=1Γn−αddtn∫att−sn−α−1fsds,which is called the Riemann–Liouville fractional derivative of order α.Definition 3.
Letα∈0,∞ and f is such that Ian−αfn exists, where n=α,f∈Ana,b (the set of all function f:a,b⟶ℝ provided that fn−1 be absolutely continuous), then we define the operator DCaα by(3)DCaαft=1Γn−α∫att−sn−α−1fnsds,which exists for almost everywhere x∈a,b. The operator DCaαft is called the Caputo fractional derivative of order α. In particular, when 0<α≤1, we have(4)DCaαft=1Γ1−α∫atf′st−sαds.Remark 1.
Letβ,γ∈ℝ+ and α∈0,1. Then,(i)
IfIaβ:L1⟶L1 and if ft∈L1, then IaβIaγft=Iaβ+γft(ii)
limβ⟶nIaβfx=Ianft uniformly on a,b, n=1,2,3,…, where Ia1ft=∫atfsds(iii)
limβ⟶0Iaβft=ft weakly(iv)
Ifft is absolutely continuous on a,b, then limα⟶1Daαft=dft/dt(v)
Thus,Daαft=d/dtIa1−αft (Riemann–Liouville sense) and Dαft=Ia1−αd/dtft (Caputo sense)Remark 2.
We notice that the fractional derivatives involve an integration and are nonlocal operators, which can be used for modeling systems with memory.
We should mention here that Caputo’s definition of fractional derivative is a modification of the Riemann–Liouville definition and has the advantage of dealing with initial value problems in a proper way.


## 2. Model Formulation

The literature reveals that most mathematical modeling of biological systems with memory is based either on delay differential equations (DDEs) with integer-order or fractional-order differential equations without a delay. However, fractional-order calculus is more suitable, than integer-order ones, in modeling biological systems with intrinsic memory and long-range interactions such as epidemic evolution systems [42]. Modeling of such systems by fractional-order differential equations has more advantages than classical integer-order mathematical modeling, in which the effects of memory or long-range interactions are neglected. Indeed, memory effects play an essential role in the spreading of diseases. Including memory effects in the susceptible-infected-recovered (SIR) epidemic models seems very appropriate for such an investigation (see Remark 2). Herein, we investigate the impact of combining both time delays and fractional order in an epidemic model for Zika virus infection.The underlying model is governed by a system of fractional-order differential equations with multiple time delays for Zika virus infection. The model includes the dynamics of susceptible individuals,HSt, with Zika symptoms and infected portion, HIt, and recovered portion, HRt, individuals recovered from Zika, the susceptible mosquitoes, MSt, in infected mosquitoes, MIt. Thus, the total human population NHt=Hst+HIt+HRt. The overall vector (mosquito) population, at time t, is Nmt=Mst+MIt. Assume that βh is the transmission rate from humans to mosquitoes. βm is the transmission rate of Zika from the vector (mosquitoes) to humans. Natural death rate of host is denoted by dh. The recruitment rate into susceptible population is denoted by λh. Natural death rate of vector is denoted by dm. η is the recovery rate from treatment. λm is the recruitment rate into susceptible mosquito population. Also, γ is the average infectious period for humans. We use time delays in the model to consider the latency of the infection in a vector and the latency of the infection in an infected host. In our model, we consider time-delay τ1 to represent the transferring of the infection from infected mosquitoes into suspected humans. The incubation period (time delay) τ2 is incorporated to represent the time required for an individual/susceptible to become infectious, after becoming infected. τ3 is the incubation period of susceptible mosquitoes to become infectious (see Figure 1). The memory of the earlier times, which are represented by time lags, could have less effect on the present situation, as compared to more recent times. However, it is expected that long-range memory, represented by fractional order, effects decay in time more slowly than an exponential decay but can typically behave like a power-law damping function. The model then takes the following form:(5)DαHSt=λh−βhHSt−τ1MIt−τ1−βhHSt−τ2HIt−τ2−dhHS,DαHIt=βhHSt−τ1MIt−τ1+βhHSt−τ2HIt−τ2−dhHI−γHI,DαHRt=γHI−dhHR+ηHI,DαMSt=λm−βmMSt−τ3HIt−τ3−dmMs,DαMIt=βmMSt−τ3HIt−τ3−dmMI.Figure 1
Transmission and dynamics of Zika virus infection between the host (human) and the vector (mosquitoes).The initial conditions for system (5) should be provided so that HR0=HR0, HSt=ϕ1t, HIt=ϕ2t, MSt=ϕ3t, and MIt=ϕ4t, when t∈max−τi,0 for i=1,2,3., time lag, τi≥0.Remark 3.
The fractional derivativeα∈0,1 is defined by Caputo sense (4), so that introducing a convolution integral with a power-law memory kernel is useful to describe memory effects in dynamical systems. The decaying rate of the memory kernel (a time correlation function) depends on α. A lower value of α corresponds to more slowly decaying time-correlation functions (long memory). Therefore, as α⟶1, the influence of memory decreases.


### 2.1. Nonnegative Solution

Since model (5) monitors the dynamics of human populations, therefore, all the parameters are assumed to be nonnegative. Furthermore, it can be shown that all state variables of the model are nonnegative and bounded for all time t≥0 (see [42]).Lemma 1.
The closed setΩ=Hs,HI,HR,MS,MI∈ℝ+5:HS+HI+HR≤λh/dh,MS+MI≤λm/dm is positively invariant with respect to model (5).Proof.
In order to prove the nonnegativity of system (5), it is assumed that there exists a t∗>t0 such that HSt∗=0 and HSt<0 for t∈t∗,t1 where t1 is sufficiently close to t∗. If HSt=0,(6)DαHSt∗=λh.
Thus, one obtainsDαHSt>0 for all t∈t∗,t1 and DαHS>ϵHS, where ϵ>0. Hence, one derives(7)HSt>HSt∗Eαϵt−t∗α,t∈t∗,t1.
SinceHSt∗=0, one gets HSt>0,t∈t∗,t1, which contradicts the assumption. Hence, HSt>0 for any t>t0. In the same manner, we have HIt,HRt,MSt,andMIt are nonnegative.
To show that the system is bounded, we add the first three equation of System (5), and we get(8)DαHS+HI+HR=λh−dhHs−dhHI−dhHR+ηHI.
We know that all parameters value is positive, and one can obtain(9)DαHS+HI+HR≤λh−dhHS+HI+HR,DαNH≤λh−dhNH,where NH=HSt+HIt+HRt, and solving this equation, we have(10)NHt≤−λhdh+NH0Eα−dhtα+λhdh.
The solution is given byNHt=NH0Eα,1−dhtα+λhtαEα,α+1−dhtα, where Eα,β is the Mittag-Leffler function. Considering the fact that Mittag-Leffler function has an asymptotic behaviour,(11)Eα,βz∼−∑K=1ωz−KΓβ−αK+Oz−1−ω,z⟶∞,απ2<argz≤π.
One can observe thatNHt⟶λh/dh as t⟶∞. The proof of the mosquitoes (vector) population is similar to human (host) population, and we obtain NMt⟶λm/dm. Therefore, all solutions of the model with initial conditions in Ω remain bounded in the positively invariant region Ω for all t∈0,∞. The region Ω is positively invariant with respect to model (5).The equilibrium points (steady states) are obtained by settingDαHS=DαHI=DαHR=DαMS=DαMI=0, in model (5). The model has two equilibrium points: (i) disease-free equilibrium point E0=HS0,HI0,HR0,MI0,MR0=λh/dh,0,0,λm/dm,0 and (ii) endemic steady state E∗, which is(12)E∗λhβhβmλmHI∗/dmβmHI∗+dm+βhHI∗+dh,HI∗,HI∗η+γdh,λmβmHI∗+dm,βmλmHI∗dmβmHI∗+dm.Here,HI∗ is the positive root of the following equation:(13)βhλhβhβmλmHI∗/dmβmHI∗+dm+βhHI∗+dhβmλmHI∗dmβmHI∗+dm+βhλhβhβmλmHI∗/dmβmHI∗+dm+βhHI∗+dhHI∗−dhHI∗−γHI∗=0.


## 2.1. Nonnegative Solution

Since model (5) monitors the dynamics of human populations, therefore, all the parameters are assumed to be nonnegative. Furthermore, it can be shown that all state variables of the model are nonnegative and bounded for all time t≥0 (see [42]).Lemma 1.
The closed setΩ=Hs,HI,HR,MS,MI∈ℝ+5:HS+HI+HR≤λh/dh,MS+MI≤λm/dm is positively invariant with respect to model (5).Proof.
In order to prove the nonnegativity of system (5), it is assumed that there exists a t∗>t0 such that HSt∗=0 and HSt<0 for t∈t∗,t1 where t1 is sufficiently close to t∗. If HSt=0,(6)DαHSt∗=λh.
Thus, one obtainsDαHSt>0 for all t∈t∗,t1 and DαHS>ϵHS, where ϵ>0. Hence, one derives(7)HSt>HSt∗Eαϵt−t∗α,t∈t∗,t1.
SinceHSt∗=0, one gets HSt>0,t∈t∗,t1, which contradicts the assumption. Hence, HSt>0 for any t>t0. In the same manner, we have HIt,HRt,MSt,andMIt are nonnegative.
To show that the system is bounded, we add the first three equation of System (5), and we get(8)DαHS+HI+HR=λh−dhHs−dhHI−dhHR+ηHI.
We know that all parameters value is positive, and one can obtain(9)DαHS+HI+HR≤λh−dhHS+HI+HR,DαNH≤λh−dhNH,where NH=HSt+HIt+HRt, and solving this equation, we have(10)NHt≤−λhdh+NH0Eα−dhtα+λhdh.
The solution is given byNHt=NH0Eα,1−dhtα+λhtαEα,α+1−dhtα, where Eα,β is the Mittag-Leffler function. Considering the fact that Mittag-Leffler function has an asymptotic behaviour,(11)Eα,βz∼−∑K=1ωz−KΓβ−αK+Oz−1−ω,z⟶∞,απ2<argz≤π.
One can observe thatNHt⟶λh/dh as t⟶∞. The proof of the mosquitoes (vector) population is similar to human (host) population, and we obtain NMt⟶λm/dm. Therefore, all solutions of the model with initial conditions in Ω remain bounded in the positively invariant region Ω for all t∈0,∞. The region Ω is positively invariant with respect to model (5).The equilibrium points (steady states) are obtained by settingDαHS=DαHI=DαHR=DαMS=DαMI=0, in model (5). The model has two equilibrium points: (i) disease-free equilibrium point E0=HS0,HI0,HR0,MI0,MR0=λh/dh,0,0,λm/dm,0 and (ii) endemic steady state E∗, which is(12)E∗λhβhβmλmHI∗/dmβmHI∗+dm+βhHI∗+dh,HI∗,HI∗η+γdh,λmβmHI∗+dm,βmλmHI∗dmβmHI∗+dm.Here,HI∗ is the positive root of the following equation:(13)βhλhβhβmλmHI∗/dmβmHI∗+dm+βhHI∗+dhβmλmHI∗dmβmHI∗+dm+βhλhβhβmλmHI∗/dmβmHI∗+dm+βhHI∗+dhHI∗−dhHI∗−γHI∗=0.


## 3. Stability and Bifurcation Analysis

To study the stability of model (5), suppose E∗HS∗,HI∗,HR∗,MS∗,MI∗ is the steady state of the linearized system:(14)DαHSt=−βhHS∗MIt−τ1−βhMI∗HSt−τ1−βhHS∗HIt−τ2−βhHI∗HSt−τ2−dhHS,DαHIt=βhHS∗MIt−τ1+βhMI∗HSt−τ1+βhHS∗HIt−τ2+βhHI∗HSt−τ2−dhHI−γHI,DαHRt=γHI+ηHI−dhHR,DαMSt=−βmMS∗HIt−τ3−βmHI∗MSt−τ3−dmMS,DαMIt=βmMS∗HIt−τ3+βmHI∗MSt−τ3−dmMI.Taking Laplace transform [43] on both sides of the linearized system (14), we obtain(15)sαX1s=sα−1φ10+βhHS∗e−sτ1−X5s−∫−τ10e−stφ5tdt+βhMI∗e−sτ1−X1s−∫−τ10e−stφ1tdt+βhHS∗e−sτ2−X2s−∫−τ20e−stφ2tdt+βhHI∗e−sτ2−X1s−∫−τ20e−stφ1tdt−dhX1s,sαX2s=sα−1φ20+βhHS∗e−sτ1X5s+∫−τ10e−stφ5tdt+βhMI∗e−sτ1X1s+∫−τ10e−stφ1tdt+βhHS∗e−sτ2X2s+∫−τ20e−stφ2tdt+βhHI∗e−sτ2X1s+∫−τ20e−stφ1tdt−dhX2s−γX2s,sαX3s=sα−1φ30+γX2s+ηX2s−dhX3s,sαX4s=sα−1φ40+βmMS∗e−sτ3−X2s−∫−τ30e−stφ2tdt+βmHI∗e−sτ3−X4s−∫−τ30e−stφ4tdt−dmX4s,sαX5s=sα−1φ50+βmMS∗e−sτ3X2s+∫−τ30e−stφ2tdt+βmHI∗e−sτ3X4s+∫−τ30e−stφ4tdt−dmX5s,where X1s, X2s, X3s, X4s, and X5s are Laplace transforms of HS,HI, HR, MS, and MI, respectively, with X1s=ℒHSt, X2s=ℒHIt, X3s=ℒHRt, X4s=ℒMSt, and X5s=ℒMIt. Then, (15) can be written in the following matrix form as(16)ΔsX1sX2sX3sX4sX5s=k1sk2sk3sk4sk5s,in which(17)Δs=sα+a1e−sτ1+a2e−sτ2+a3a4e−sτ200a4e−sτ1−a1e−sτ1−a2e−sτ2sα−a4e−sτ2+a500−a4e−sτ10a6sα+a3000a7e−sτ30sα+a8e−sτ3+a900−a7e−sτ30−a8e−sτ3sα+a9,k1s=sα−1φ10−βhe−sτ1HS∗∫−τ10e−stφ5tdt−βhe−sτ1MI∗∫−τ10e−stφ1tdt−βhe−sτ2HS∗∫−τ20e−stφ2tdt−βhe−sτ2HI∗∫−τ20e−stφ1tdt,k2s=sα−1φ20+βhe−sτ1HS∗∫−τ10e−stφ5tdt+βhe−sτ1MI∗∫−τ10e−stφ1tdt+βhe−sτ2HS∗∫−τ20e−stφ2tdt+βhe−sτ2HI∗∫−τ20e−stφ1tdt,k3s=sα−1φ30,k4s=sα−1φ40−βmMS∗e−sτ3∫−τ30e−stφ2tdt−βme−sτ3HI∗∫−τ30e−sτ3φ4tdt,k5s=sα−1φ50+βmMS∗e−sτ3∫−τ30e−stφ2tdt+βme−sτ3HI∗∫−τ30e−sτ3φ4tdt,where a1=βhMI∗, a2=βhHI∗, a3=dh, a4=βhHS∗, a5=dh+γ, a6=−η−γ, a7=βmMS∗, a8=βmHI∗, and a9=dm and Δs is considered as the characteristic matrix of system (5) The characteristic polynomial is then(18)Ps=ℙ1s+ℙ2se−sτ1+ℙ3se−sτ2+ℙ4se−sτ3+ℙ5se−2sτ1+ℙ6se−2sτ2+ℙ7se−sτ1+τ2+ℙ8se−sτ2+τ3+ℙ9se−sτ1+τ3+ℙ10se−2sτ1−sτ3+ℙ11se−2sτ2−sτ3+ℙ12se−sτ1+τ2+τ3.The coefficientsℙis, i=1,…,12, are estimated by Mathematica and given in Appendix.Case 1.
τ
1
>
0
,
τ
2
=
0
,
and
τ
3
=
0.
Whenτ1>0,τ2=0,andτ3=0, the characteristic equation (18) becomes(19)P1s+P2se−sτ1+P3se−2sτ1=0,where(20)P1s=ℙ1s+ℙ3s+ℙ4s+ℙ6s+ℙ8s+ℙ11s=s5α+D1s4α+D2s3α+D3s2α+D4sα+D5,P2s=ℙ2s+ℙ7s+ℙ9s+ℙ12s=G1s4α+G2s3α+G3s2α+G4sα+G5,P3s=ℙ5s+ℙ10s=H1s3α+H2s2α+H3sα+H4.
Now, we prove that the characteristic equation (19) has no pure imaginary roots for any τ1>0. Assume that characteristic equation (19) has pure imaginary root, and let it be s=iξ=ξcosπ/2+isinπ/2,ξ>0. If we multiply esτ1 on both sides of equation (19), we get(21)P1sesτ1+P2s+P3se−sτ1=0.
Now, we substitute the expression ofs into (21) to have(22)A1+iℬ1esτ1+A2+iℬ2+A3+iℬ3e−sτ1=0.
The coefficientsA1,A2,andA3 and ℬ1,ℬ2,andℬ3 are real and imaginary parts of P1s,P2s, and P3s, respectively, so that(23)A1=ξ5αcos5απ2+D1ξ4αcos4απ2+D2ξ3αcos3απ2+D3ξ2αcos2απ2+D4ξαcosαπ2+D5,ℬ1=ξ5αsin5απ2+D1ξ4αsin4απ2+D2ξ3αsin3απ2+D3ξ2αsin2απ2+D4ξαsinαπ2,A2=G1ξ4αcos4απ2+G2ξ3αcos3απ2+G3ξ2αcos2απ2+G4ξαcosαπ2+G5,ℬ2=G1ξ4αsin4απ2+G2ξ3αsin3απ2+G3ξ2αsin2απ2+G4ξαsinαπ2,A3=H1ξ3αcos3απ2+H2ξ2αcos2απ2+H3ξαcosαπ2+H4,ℬ3=H1ξ3αsin3απ2+H2ξ2αsin2απ2+H3ξαsinαπ2.
Separating real and imaginary parts yields(24)A1cosξτ1−ℬ1sinξτ1=−A3cosξτ1+ℬ3sinξτ1+A2,A1sinξτ1−ℬ1cosξτ1=−ℬ3cosξτ1−A3sinξτ1+ℬ2.
It follows from (14) that(25)A12+ℬ12−A22−ℬ22−A32−ℬ32=2ℬ3A2sinξτ1+ℬ2cosξτ1+A3A2cosξτ1−ℬ2sinξτ1.
Using the fact thatcos2θ+sin2θ=1, we have sinξτ1=1−cos2ξτ1, and then (25) can be written in the following form:(26)A12+ℬ12−A22−ℬ22−A32−ℬ32−2ℬ2ℬ3+A2A3cosξτ12=21−cos2ξτ1ℬ3A2−A3ℬ22.
It can be concluded from (26) that(27)Q1cos2ξτ+Q2cosξτ+Q3=0,where(28)Q1=4A22A32+4A32ℬ22+4A22ℬ32+4ℬ22ℬ32,Q2=4ℬ2ℬ3+A2A3−A12−ℬ12+A22+ℬ22+A32+ℬ32,Q3=A12−A2−ℬ32−A3−ℬ1+ℬ2A3+ℬ2+ℬ1A12−A2+ℬ32−A3−ℬ1−ℬ2A3+ℬ1−ℬ2.
The quadratic equation (27) has roots, so we can obtain the expression of cosξτ1 and denote cosξτ1=f1ξ, where f1ξ is a function of ξ.
Substituting the expression ofcosξτ1=1−sin2ξτ1 into (27), we can get expression of sinξτ1. Assume that sinξτ2=f2ξ. Moreover, we have f12ξ+f22ξ=1. Thus, it follows from cosξτ1=f1ξ that(29)τ1=1ξarccosf1ξ+2kπ,k=0,1,2,….
We suppose thatf12ξ+f22ξ=1 has at least one positive root, and thus, the bifurcation point is defined as(30)τ1∗=minτ1k,k=0,1,2,….
Now, differentiating equation (21) with respect to τ1, we obtain(31)P1′sesτ1dsdτ1+P1esτ1τ1dsdτ1+s+P2′sdsdτ1+P3′se−sτ1dsdτ1+P3e−sτ1−τ1dsdτ1−s=0,where P1′s,P2′s, and P3′s are derivatives of P1s,P2s, and P3s, respectively. It follows that(32)dsdτ1=−sP1sesτ1+sP3se−sτ1P1′sesτ1+τ1P1sesτ1+P2′s+P3′se−sτ1−τ1P3se−sτ1=MsNs.
From (32), by some computation, we deduce that(33)Redsdτ1τ1=τ1∗,ξ=ξ0=M1N1+M2N2N12+N22,where M1,N1 and M2,N2 are the real and imaginary parts of Ms,Ns. Also, ξ0 stands for the critical value and τ1∗ denotes the bifurcation point. Here,(34)M1=A1ξ0sinξ0τ1∗+ℬ1ξ0cosξ0τ1∗+A3ξ0sinξ0τ1∗−ℬ3ξ0sinξ0τ1∗,M2=−A1ξ0cosξ0τ1∗+ℬ1ξ0sinξ0τ1∗+ξ0A3cosξ0τ1∗−ℬ3ξ0cosξ0τ1∗,N1=A1∗cosξ0τ1∗−ℬ1∗sinξ0τ1∗+τ1∗A1cosξ0τ1∗−τ1∗ℬ1sinξ0τ1∗+A2∗+A3∗cosξ0τ1∗+B3∗sinξ0τ1∗−τ1∗A3cosξ0τ1∗−τ1∗ℬ3sinξ0τ1∗,N2=A1∗sinξ0τ1∗+ℬ1∗cosξ0τ1∗+τ1∗A1sinξ0τ1∗+τ1∗ℬ1cosξ0τ1∗+ℬ2∗−A3∗sinξ0τ1∗+B3∗cosξ0τ1∗+τ1∗A3sinξ0τ1∗−τ1∗ℬ3cosξ0τ1∗,A1∗=5αξ05α−1cos5α−1π2+4αD1ξ04α−1cos4α−1π2+3αD2ξ03α−1cos3α−1π2+2αD3ξ02α−1cos2α−1π2+αD4ξ0α−1cosα−1π2,ℬ1∗=5αξ05α−1sin5α−1π2+4αD1ξ04α−1sin4α−1π2+3αD2ξ03α−1sin3α−1π2+2αD3ξ02α−1sin2α−1π2+αD4ξ0α−1sinα−1π2,A2∗=4αG1ξ04α−1cos4α−1π2+3αG2ξ03α−1cos3α−1π2+2αG3ξ02α−1cos2α−1π2+αG4ξ0α−1cosα−1π2,ℬ2∗=4αG1ξ04α−1sin4α−1π2+3αG2ξ03α−1sin3α−1π2+2αG3ξ02α−1sin2α−1π2+αG4ξ0α−1sinα−1π2,A3∗=3αH1ξ03α−1cos3α−1π2+2αH2ξ02α−1cos2α−1π2+αH3ξ0α−1cosα−1π2,ℬ3∗=3αH1ξ03α−1sin3α−1π2+2αH2ξ02α−1sin2α−1π2+αH3ξ0α−1sinα−1π2.Case 2.
τ
1
=
0
,
τ
2
>
0
,
and
τ
3
=
0.
Whenτ1=0,τ2>0,τ3=0, the characteristic equation (18) becomes(35)P4s+P5se−sτ2+P6se−2sτ2=0,where(36)P4s=ℙ1s+ℙ2s+ℙ4s+ℙ5s+ℙ9s+ℙ10s=s5α+J1s4α+J2s3α+J3s2α+J4sα+J5,P5s=ℙ3s+ℙ7s+ℙ8s+ℙ12s=L1s4α+L2s3α+L3s2α+L4sα+L5,P6s=ℙ6s+ℙ11s=R1s3α+R2s2α+R3sα+R4.
Now, we prove that the characteristic equation (35) has no pure imaginary roots for any τ2>0. Assume that characteristic equation (35) has pure imaginary root, and let it be s=iξ=ξcosπ/2+isinπ/2,ξ>0. Now, multiplying esτ2 on both sides of equation (35), we get(37)P4sesτ2+P5s+P6se−sτ2=0.
Substitute the expression ofs into (37) to have(38)A4+iℬ4esτ2+A5+iℬ5+A6+iℬ6e−sτ2=0,where A4,A5,A6 and ℬ4,ℬ5,ℬ6 are real and imaginary parts of P4s,P5s, and P6s,respectively. Here,(39)A4=ξ5αcos5απ2+J1ξ4αcos4απ2+J2ξ3αcos3απ2+J3ξ2αcos2απ2+J4ξαcosαπ2+J5,ℬ4=ξ5αsin5απ2+J1ξ4αsin4απ2+J2ξ3αsin3απ2+J3ξ2αsin2απ2+J4ξαsinαπ2,A5=L1ξ4αcos4απ2+L2ξ3αcos3απ2+L3ξ2αcos2απ2+L4ξαcosαπ2+L5,ℬ5=L1ξ4αsin4απ2+L2ξ3αsin3απ2+L3ξ2αsin2απ2+L4ξαsinαπ2,A6=R1ξ3αcos3απ2+R2ξ2αcos2απ2+R3ξαcosαπ2+R4,ℬ6=R1ξ3αsin3απ2+R2ξ2αsin2απ2+R3ξαsinαπ2.
Separating real and imaginary parts yields(40)A4cosξτ2−ℬ4sinξτ2=−A6cosξτ2+ℬ6sinξτ2+A5,A4sinξτ2−ℬ4cosξτ2=−ℬ6cosξτ2−A6sinξτ2+ℬ5.
It follows from (40) that(41)A42+ℬ42−A52−ℬ52−A62−ℬ62=2ℬ6A5sinξτ2+ℬ5cosξτ2+A6A5cosξτ2−ℬ5sinξτ2.
We know thatcos2θ+sin2θ=1; by using it, we have sinξτ2=1−cos2ξτ2, and then (41) can be written in the following form:(42)A42+ℬ42−A52−ℬ52−A62−ℬ62−2ℬ5ℬ6+A5A6cosξτ22=21−cos2ξτ2ℬ6A5−A6ℬ52.
It can be concluded from (42) that(43)Q4cos2ξτ+Q5cosξτ+Q6=0,where(44)Q4=4A52A62+4A62ℬ52+4A52ℬ62+4ℬ52ℬ62,Q5=4ℬ5ℬ6+A5A6−A42−ℬ42+A52+ℬ52+A62+ℬ62,Q6=A42−A5−ℬ62−A6−ℬ4+ℬ5A6+ℬ5+ℬ4A42−A5+ℬ62−A6−ℬ4−ℬ5A6+ℬ4−ℬ5.
As we know, the quadratic equation (43) has roots, we can obtain the expression of cosξτ2 and denote cosξτ2=f1ξ, where f1ξ is a function of ξ.
Substituting the expression ofcosξτ2=1−sin2ξτ2 into (43), we can get expression of sinξτ2. Let us denote sinξτ2=f2ξ, where f2ξ is a function with respect to ξ. Moreover, f12ξ+f22ξ=1. Thus, it follows from cosξτ2=f1ξ that(45)τ2=1ξarccosf1ξ+2kπ,k=0,1,2,….
Clearly,f12ξ+f22ξ=1 has at least one positive root. The bifurcation point is defined as(46)τ2∗=minτ2k,k=0,1,2,….
We obtain the transversality condition of the occurrence for Hopf bifurcation atτ2=τ2∗.
Differentiating equation (37) with respect to τ2 yields(47)P4′sesτ2dsdτ2+P4sesτ2τ2dsdτ2+s+P5′sdsdτ2+P6′se−sτ2dsdτ2+P6se−sτ2−τ2dsdτ2−s=0,where P4′s,P5′s, and P6′s are derivatives of P4s,P5s, and P6s, respectively. It follows that(48)dsdτ2=−sP4sesτ2+sP6se−sτ2P4′sesτ2+τ2P4sesτ2+P5′s+P6′se−sτ2−τ2P6se−sτ1=ℳsNs.
From (48), by some computation, we deduce that(49)Redsdτ2τ2=τ2∗,ξ=ξ0=ℳ1N1+ℳ2N2N12+N22,where ℳ1,N1 and ℳ2,N2 are the real and imaginary parts of ℳs,Ns. Also ξ0 stands for the critical value and τ2∗ denotes bifurcation point. Here,(50)ℳ1=A4ξ0sinξ0τ2∗+ℬ4ξ0cosξ0τ2∗+A6ξ0sinξ0τ2∗−A6ξ0sinξ0τ2∗,ℳ2=−A4ξ0cosξ0τ2∗+ℬ4ξ0sinξ0τ2∗+ξ0A6cosξ0τ2∗−ℬ6ξ0cosξ0τ2∗,N1=A4∗cosξ0τ2∗−ℬ2∗sinξ0τ2∗+τ2∗A4cosξ0τ2∗−τ2∗ℬ4sinξ0τ2∗+A5∗+A6∗cosξ0τ1∗+ℬ6∗sinξ0τ2∗−τ2∗A6cosξ0τ2∗−τ2∗ℬ6sinξ0τ2∗,N2=A4∗sinξ0τ2∗+ℬ4∗cosξ0τ2∗+τ2∗A4sinξ0τ2∗+τ2∗ℬ4cosξ0τ2∗+ℬ5∗−A6∗sinξ0τ2∗+ℬ6∗cosξ0τ2∗+τ2∗A6sinξ0τ2∗−τ2∗ℬ6cosξ0τ2∗,A4∗=5αξ05α−1cos5α−1π2+4αJ1ξ04α−1cos4α−1π2+3αJ2ξ03α−1cos3α−1π2+2αJ3ξ02α−1cos2α−1π2+αJ4ξ0α−1cosα−1π2,ℬ4∗=5αξ05α−1sin5α−1π2+4αJ1ξ04α−1sin4α−1π2+3αJ2ξ03α−1sin3α−1π2+2αJ3ξ02α−1sin2α−1π2+αJ4ξ0α−1sinα−1π2,A5∗=4αL1ξ04α−1cos4α−1π2+3αL2ξ03α−1cos3α−1π2+2αL3ξ02α−1cos2α−1π2+αL4ξ0α−1cosα−1π2,ℬ5∗=4αL1ξ04α−1sin4α−1π2+3αL2ξ03α−1sin3α−1π2+2αL3ξ02α−1sin2α−1π2+αL4ξ0α−1sinα−1π2,A6∗=3αR1ξ03α−1cos3α−1π2+2αR2ξ02α−1cos2α−1π2+αR3ξ0α−1cosα−1π2,ℬ6∗=3αR1ξ03α−1sin3α−1π2+2αR2ξ02α−1sin2α−1π2+αR3ξ0α−1sinα−1π2.Case 3.
τ
1
=
0
,
τ
2
=
0
,
and
τ
3
>
0.
Whenτ1=0,τ2=0,andτ3>0, the characteristic equation (18) becomes(51)P7s+P8se−sτ3=0,where(52)P7s=ℙ1s+ℙ2s+ℙ3s+ℙ5s+ℙ6s+ℙ7s=s5α+U1s4α+U2s3α+U3s2α+U4sα+U5,P8s=ℙ4s+ℙ8s+ℙ9s+ℙ10s+ℙ11s+ℙ12s=V1s4α+V2s3α+V3s2α+V4sα+V5.
Again, we prove that the characteristic equation (51) has no pure imaginary roots for any τ2>0. Here, we assume that characteristic equation (51) has pure imaginary root, let it be s=iξ=ξcosπ/2+isinπ/2,ξ>0. Now, we substitute the expression of s into (51), and we have(53)A7+iℬ7+A8+iℬ8e−sτ3=0,where A7,A8 and ℬ7,ℬ8 are real and imaginary parts of P7s and P8s, respectively. Here,(54)A7=ξ5αcos5απ2+U1ξ4αcos4απ2+U2ξ3αcos3απ2+U3ξ2αcos2απ2+U4ξαcosαπ2+U5,ℬ7=ξ5αsin5απ2+U1ξ4αsin4απ2+U2ξ3αsin3απ2+U3ξ2αsin2απ2+U4ξαsinαπ2,A5=V1ξ4αcos4απ2+V2ξ3αcos3απ2+V3ξ2αcos2απ2+V4ξαcosαπ2+V5,ℬ5=V1ξ4αsin4απ2+V2ξ3αsin3απ2+V3ξ2αsin2απ2+V4ξαsinαπ2.
Separation of real and imaginary parts yields(55)A8cosξτ3+ℬ8sinξτ3=−A7,−A8sinξτ3+ℬ8cosξτ3=−ℬ7.
From (14), we have(56)cosξτ3=−A7A8−ℬ7ℬ8A82+ℬ82=f¯1ξ,sinξτ3=ℬ7A8−A7ℬ8A82+ℬ82=f¯2ξ.
It is clear thatcos2θ+sin2θ=1; from (56),(57)f¯1ξ2+f¯2ξ2=1.
Hence, it follows fromcosξτ3=f¯1ξ that(58)τ3=1ξarccosf¯1ξ+2kπ,k=0,1,2,….
We suppose that (57) have at least one positive root. The bifurcation point is defined as(59)τ3∗=minτ3k,k=0,1,2,….
We obtain the transversality condition of the occurrence for Hopf bifurcation atτ3=τ3∗.
Now, differentiating equation (51) with respect to τ3, we obtain(60)P7′sdsdτ3+P8′se−sτ3dsdτ3+P8se−sτ3−τ3dsdτ3−s=0,where P7′s and P8′s are derivatives of P7s and P8s, respectively. It follows that(61)dsdτ3=−sP8sesτ3P7′s+P8′se−sτ3−τ3P8se−sτ3=MsNs.
From (61), by some computation, we deduce that(62)Redsdτ3τ3=τ3∗,ξ=ξ0=M1N1+M2N2N12+N22,where M1,N1 and M2,N2 are the real and imaginary parts of Ms,Ns. Also ξ0 stands for the critical value and τ3∗ denotes bifurcation point. Here,(63)M1=A8ξ0sinξ0τ3∗−ℬ8ξ0cosξ0τ3∗,M2=A8ξ0cosξ0τ3∗+ℬ8ξ0sinξ0τ3∗,N1=A7∗+A8∗cosξ0τ3∗+ℬ8∗sinξ0τ3∗−τ3∗A8cosξ0τ3∗+ℬ8sinξ0τ3∗,N2=ℬ7∗+ℬ8∗cosξ0τ3∗−A8∗sinξ0τ3∗−τ3∗ℬ8cosξ0τ3∗−A8sinξ0τ3∗,A7∗=5αξ05α−1cos5α−1π2+4αU1ξ04α−1cos4α−1π2+3αU2ξ03α−1cos3α−1π2+2αU3ξ02α−1cos2α−1π2+αU4ξ0α−1cosα−1π2,ℬ7∗=5αξ05α−1sin5α−1π2+4αU1ξ04α−1sin4α−1π2+3αU2ξ03α−1sin3α−1π2+2αU3ξ02α−1sin2α−1π2+αU4ξ0α−1sinα−1π2,A8∗=4αV1ξ04α−1cos4α−1π2+3αV2ξ03α−1cos3α−1π2+2αV3ξ02α−1cos2α−1π2+αV4ξ0α−1cosα−1π2,ℬ5∗=4αV1ξ04α−1sin4α−1π2+3αV2ξ03α−1sin3α−1π2+2αV3ξ02α−1sin2α−1π2+αV4ξ0α−1sinα−1π2.Case 4.
τ
1
=
τ
2
=
τ
3
=
0.
Whenτ1=τ2=τ3=0, the characteristic equation (18) becomes(64)ω5α+Z1ω4α+Z2ω3α+Z3ω2α+Z4ωα+Z5=0.
Here,(65)Z1=2a3+a5+2a9+a1−a4+a2+a8,Z2=a32+4a3a9+2a3a5+2a5a9+a92+a1a3−2a3a4+a1a5+2a1a9−2a4a9+a2a3+a2a5+2a2a9+2a3a8+a5a8+a8a9−a1a4+a2a4+a1a4−a2a4+a2a8−a4a7+a1a8−a4a8,Z3=2a3a92+a5a92+2a32a9+a32a5+4a3a5a9−a32a4+a1a3a5+2a1a3a9−4a3a4a9+2a1a5a9+a1a92−a4a92+a2a3a5+2a2a3a9+2a2a5a9+a2a92+a32a8+2a3a5a8+2a3a8a9+a5a8a9−a1a3a4−2a1a4a9+a2a3a4+2a2a4a9+a1a3a4−a2a3a4+2a1a4a9−2a2a4a9+a2a3a8+a2a5a8+a2a8a9−2a3a4a7+a1a3a8−2a3a4a8+a1a5a8−a4a7a9+a1a8a9−a4a8a9,Z4=a32a92+2a3a5a92+2a32a5a9−2a32a4a9+2a1a3a5a9+a1a3a92−2a3a4a92+a1a5a92+a2a3a92+a2a5a92+2a2a3a5a9+a32a5a8+a32a8a9+2a3a5a8a9−a1a4a92−2a1a3a4a9+2a2a3a4a9+a1a4a92−a2a4a92+a2a3a5a8+a2a3a8a9+a2a5a8a9−a32a4a7+a1a3a5a8−a32a4a8−2a3a4a7a9+a1a3a8a9−2a3a4a8a9+a1a5a8a9−a1a3a4a8−a1a4a8a9+a2a3a4a8+a2a4a8a9+a1a3a4a8−a2a3a4a8+a1a4a8a9−a2a4a8a9,Z5=a32a5a92−a32a4a92+a1a3a5a92+a2a3a5a92+a32a5a8a9−a1a3a4a92+a2a3a4a92−a2a3a4a92+a2a3a5a8a9−a32a4a7a9+a1a3a5a8a9.
From the Routh–Hurwitz criteria, if we chooseZi>0,i=1,2,3,4,5, Z1Z2Z3>Z32+Z12Z4 and Z1Z4−Z5Z1Z2Z3−Z32−Z12Z4>Z5Z1Z2−Z32+Z1Z52, then the five eigenvalues of the characteristic equation (64) have negative real parts. Hence, the steady state fractional-order system (5) is asymptotically stable when τ1=τ2=τ3=0 (without time delays).
We arrive at the following theorem.Theorem 1.
Ifα∈0,1 and an endemic equilibrium point E∗ exists for system (5), then the following results hold:(i)
Whenτ1>0,τ2=0,andτ3=0, the endemic steady state E∗ is asymptotically stable for τ1∈0,τ1∗ and the system undergoes a Hopf bifurcation at the origin at τ1=τ1∗ and the transversality condition holds, Reds/dτ1τ1=τ1∗,ξ=ξ0≠0(ii)
Whenτ1=0,τ2>0,andτ3=0, the endemic steady state E∗ is asymptotically stable for τ2∈0,τ2∗ and system undergoes a Hopf bifurcation at the origin when τ2=τ2∗ and the transversality condition holds, Reds/dτ2τ2=τ2∗,ξ=ξ0≠0(iii)
Whenτ1=0,τ2=0,andτ3>0, E∗ of is asymptotically stable for τ3∈0,τ3∗ and system (5) undergoes a Hopf bifurcation at the origin at τ3=τ3∗ and transversality condition holds, Reds/dτ3τ3=τ3∗,ξ=ξ0≠0(iv)
WhenZ1Z2Z3>Z32+Z12Z4 and Z1Z4−Z5Z1Z2Z3−Z32−Z12Z4>Z5Z1Z2−Z32+Z1Z52 holds, the endemic steady state E∗ is asymptotically stable for τ1=τ2=τ3=0Remark 4.
Theorem 1 reports the asymptotic stability of the endemic equilibrium pointE∗. The analysis can be extended to investigate the stability of infection-free equilibrium points E0 for the fractional-order model.Table 1
ξ
0 values and τ1∗ values for different fractional-order α.
Fractional orderα
Critical frequencyξ0
Bifurcation pointτ1∗
1
0.30125
1.2104
0.9
0.17811
4.5874
0.8
0.09392
11.8356
0.7
0.04573
23.2562
0.6
0.00152
744.420


## 4. Numerical Simulations and Observations

In this section, we provide some numerical simulations for system (5) to demonstrate the effectiveness of our main results. The simulations have been done by using stable implicit Euler approximation scheme, discussed in [44]. Of course, many other methods have been used for fractional-order delay differential equations such as the Adams–Bashforth–Moulton scheme [45]. The parameter values of system (5) are taken as follows:(66)λh=0.5,λm=4.58,βh=0.05,βm=0.09,η=0.01,γ=0.2,dh=0.714,dm=0.437.Case 1.
τ
1
>
0
,
τ
2
=
0
,
and
τ
3
=
0. In this case, time-delay τ1 is chosen as the bifurcation parameter. We then discuss the dynamic effect of system (5) with the above parameter values. We calculate the critical frequency ξ0 and bifurcation point τ1∗ of various fractional-order α. Figure 2 shows the numerical simulations of model (5) when τ1=10,τ2=0.0, and τ3=0.0, with different fractional orders α=0.9,0.8,and0.7 and estimated bifurcation point τ1∗=4.587,11.835,and23.256 (see Table 1). Here, τ1=10∉0,τ1∗ for the fractional-order α=0.9 and whereas τ1=10∈0,τ1∗ which satisfies the condition (i) in Theorem 1. The equilibrium E∗ of the model (5) is asymptotically stable for α=0.8,0.7. When τ1=10∉0,τ1∗, which does not satisfies the condition (i) of Theorem 1, the system undergoes a Hopf bifurcation for the functional-order α=0.9.Figure 2
State trajectories for model (5) with α=0.9,0.8,and0.7 and τ1=10,τ2=0.0,andτ3=0.0. For α=0.9 and τ1=10>τ1∗, the equilibrium point is unstable (red trajectory) for (5); however, for τ1<τ1∗ and α=0.8,0.7, it is asymptotically stable (blue and green trajectories).Table 2
ξ
0 values and τ2∗ values for different fractional-order α.
Fractional orderα
Critical frequencyξ0
Bifurcation pointτ2∗
1
0.12290
11.1762
0.9
0.09602
14.2931
0.8
0.05174
35.797
0.7
0.0363
47.963
0.6
0.0158
131.001Case 2.
τ
1
=
0
,
τ
2
>
0
,
and
τ
3
=
0. We choose time-delay τ2 as a bifurcation parameter of system (5) with parameter values:(67)λh=0.5,λm=4.58,βh=0.05,βm=0.08,η=0.05,γ=0.2,dh=0.3,dm=0.78.
We then calculate the critical frequencyξ0 and bifurcation point τ2∗ of various fractional-order α (see Table 2). Figure 3 shows the dynamics of system (5) for τ1=0.0,τ2=14, and τ3=0.0, with values of α=1,0.9, and 0.8. The corresponding bifurcation point is τ2∗=11.176,14.293,and35.797. τ2=14∈0,τ2∗ satisfies the condition (ii) of Theorem 1. Therefore, the equilibrium E∗ of the model (5) is asymptotically stable for α=0.9and0.8, which is shown in Figure 3. However, for τ2=14∉0,τ2∗, a Hopf bifurcation occurs for the functional-order α=1.Table 3
ξ
0 values and τ3∗ values for different fractional-order α.
Fractional orderα
Critical frequencyξ0
Bifurcation pointτ3∗
1
0.3777
2.255
0.9
0.2943
3.281
0.8
0.2172
4.912
0.7
0.01490
7.755
0.6
0.0920
13.404Figure 3
State trajectories for model (5) for α=1,0.9,and0.8 and τ1=0.0,τ2=14,andτ3=0.0, when α=1 and τ2=14>τ2∗; the equilibrium point is unstable (red trajectory) for (5). However, for τ2<τ2∗ with α=0.9and0.8, it is asymptotically stable (blue and green trajectories).Case 3.
τ
1
=
0
,
τ
2
=
0
,
and
τ
3
>
0. We consider time-delay τ3 as a bifurcation parameter of system (5) with parameter values:(68)λh=0.5,λm=10,βh=0.05,βm=0.4,η=0.05,γ=0.2,dh=0.714,dm=0.437.
We calculate the critical frequencyξ0 and bifurcation point τ3∗ of various fractional-order α. When τ1=0.0,τ2=0.0, and τ3=8.0, the dynamics of system (5) is shown in Figure 4 with different fractional-order α=1,0.9, and 0.8, its corresponding bifurcation points τ3∗=2.255,3.281,and4.912 (see Table 1). Here, τ3=8∉0,τ3∗ and a Hopf bifurcation occurs for the factional-order α=1,0.9,and0.8 which not satisfies the condition (iii) in Theorem 1. Therefore, the equilibrium point E∗ is of model (5) is unstable, which is shown in Figure 4.Figure 4
State trajectories for the model (5) for various values of α=1,0.9, and0.8 and τ1=0.0,τ2=0.0,andτ3=8. The equilibrium point E∗ is unstable when τ3∉0,τ3∗.Case 4.
τ
1
=
0
,
τ
2
=
0
,
and
τ
3
=
0, without time delays. Assume the parameter values:(69)λh=0.5,λm=10,βh=0.05,βm=0.4,η=0.05,γ=0.2,dh=0.714,dm=0.437.
Hence, system (5) is asymptotically stable, which is shown in Figure 5.Figure 5
State trajectories of model (5) for various values of α=0.9,0.8,and0.7, when τ1=τ2=τ3=0. The steady state of the system is asymptotically stable.


## 5. Conclusion

Fractional derivatives have the unique property of capturing the history of the variable; that is, they have short and long memory. This cannot be easily done by means of the integer-order derivatives. In this paper, we proposed a fractional-order model for Zika virus infection with multiple time delaysτ1, τ2, and τ3. We studied the asymptotic stability and Hopf bifurcation properties for the model. Time delays and fractional order play a vital role in the stability and complexity of the model. By evaluating the characteristics, some sufficient conditions have derived to ensure the asymptotic stability in terms of the fractional order and time delays. Moreover, we estimated the thresholds bifurcation parameters: τ1∗, τ2∗, and τ3∗. The transversality conditions have been obtained to confirm the existence of Hopf bifurcations for different values at the threshold parameters and particular values of fractional orders. Our findings illustrate that using the time delays as bifurcation points, one can conclude that when time delay increases, the equilibrium loses its stability and Hopf bifurcation occurs. These models can be used to understand key aspects of the viral life cycle and to predict antiviral efficacy. Finally, numerical simulations show that a combination of fractional order and time delays in the model effectively enriches the dynamics and strengthens the stability condition of the model.Including control variables in the model is desirable to determine the best strategy of treatment and control and eliminate the infection, which will be considered in future work.



---

*Source: 4178073-2019-11-14.xml*