A mathematical model of dengue transmission with memory

Accepted Manuscript A mathematical model of dengue transmission with memory Tridip Sardar, Sourav Rana, Joydev Chattopadhyay PII: DOI: Reference:

S1007-5704(14)00392-X http://dx.doi.org/10.1016/j.cnsns.2014.08.009 CNSNS 3308

To appear in:

Communications in Nonlinear Science and Numerical Simulation

Please cite this article as: Sardar, T., Rana, S., Chattopadhyay, J., A mathematical model of dengue transmission with memory, Communications in Nonlinear Science and Numerical Simulation (2014), doi: http://dx.doi.org/ 10.1016/j.cnsns.2014.08.009

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A mathematical model of dengue transmission with memory Tridip Sardara , Sourav Ranab , Joydev Chattopadhyaya,1,∗ a

Agricultural and Ecological Research Unit, Indian Statistical Institute, 203, B. T. Road, Kolkata, 700108, West Bengal, India b Department of Statistics, Visva-Bharati University, Santiniketan, West Bengal, India

Abstract We propose and analyze a new compartmental model of dengue transmission with memory between human-to-mosquito and mosquito-to-human. The memory is incorporated in the model by using a fractional differential operator. A threshold quantity R0 , similar to the basic reproduction number, is worked out. We determine the stability condition of the disease-free equilibrium (DFE) E0 with respect to the order of the fractional derivative α and R0 . We determine α dependent threshold values for R0 , below which DFE (E0 ) is always stable, above which DFE is always unstable, and at which the system exhibits a Hopf-type bifurcation. It is shown that even though R0 is less than unity, the DFE may not be always stable, and the system exhibits a Hopf-type bifurcation. Thus, making R0 < 1 for controlling the disease is no longer a sufficient condition. This result is synergistic with the concept of backward bifurcation in dengue ODE models. It is also shown that R0 > 1 may not be a sufficient condition for the persistence of the disease. For a special case, when α = 12 , we analytically verify our findings and determine the critical value of R0 in terms of some

important model parameters. Finally, we discuss about some dengue control strategies in light of the threshold quantity R0 . Keywords: Dengue transmission, Mathematical model, Fractional order differential equations

1. Introduction Dengue is a viral disease transmitted by the bite of an Aedes mosquito infected with one of the four dengue virus serotypes (DEN-1, DEN-2, DEN-3 and DEN-4) (World Health Organization, 2013a,b). Dengue can affect almost all age groups (infant to adult), and symptoms appear 3-14 days after the infected mosquito bite (World Health Organization, 2013b). A person recovers from one ∗

Corresponding author Email addresses: [email protected] (Tridip Sardar), [email protected] (Sourav Rana), [email protected] (Joydev Chattopadhyay) 1 Tel.: +91-33-25753231, Fax.: +91-33-25773049 Preprint submitted to Elsevier

August 21, 2014

of the dengue serotype having life-long immunity to that serotype but prone to infection from other three serotypes. About 12 weeks time the person becomes more susceptible to develop dengue hemorrhagic fever or dengue shock syndrome (Gubler, 1998). In recent decades, the burden of dengue increases rapidly and according to WHO there may be 50-100 million dengue cases occur each year worldwide (World Health Organization, 2013a). To understand the dengue transmission dynamics fully, it is of utmost importance to realize the blood feeding behavior of Aedes mosquito. Recent entomological studies on Aedes aegypti revealed that mosquito did not feed randomly on host blood, but they use their prior experience about a host location and a host defensiveness to select a host to feed on (McCall and Kelly, 2002; Hii et al., 1991; Chaves et al., 2010; Chilaka et al., 2012; Vinauger et al., 2011; Kelly, 2001; Takken and Verhulst, 2013). Furthermore, on human population also, the memory plays a key role in dengue transmission. In epidemic and endemic area’s awareness about dengue will lessen the contact rate between host and mosquitoes (Acharya et al., 2005; Rosenbaum et al., 1995). Thus, in dengue transmission, a future state depends upon the full history of the transmission process (Schutz and Trimper, 2004). In order to understand and control dengue infection quite a good number of mathematical (compartmental) models have been proposed and analyzed (for example, see Dietz (1975); Newton and Reiter (1992); Esteva and Yang (2005); Chowell et al. (2007); Garba et al. (2008); Thome et al. (2010)). These studies mainly focused on a deterministic integer-order compartmental epidemic model (which consists of a system of ordinary differential equation) as a key tool for their investigation. However, integer-order systems are in general memory less (Stanislavsky, 2000; Ahmed and Elgazzar, 2007). Therefore, a dengue system modelled through integer order derivatives cannot reflect the memory effect in dengue transmission. A possible generalization of the deterministic ODE dengue model would be a system that carries information about its different previous states. The behavior of a trajectory of a fractional differential operator is non-local (Oldham and Spanier, 1974; Podlubny, 1999; Hanert et al., 2011; Agarwal et al., 2013). This property of the fractional differential operator can be useful to include memory in a dynamical process. However, memory can also be included in a dynamical system using distributed delay, in which a future state depends upon the full history of the process (Gopalsamy and He, 1994; Chen et al., 2006; Sipahi et al., 2007). In a dynamical system with distributed delay, memory is not measurable in terms of any particular parameter of the model. Therefore, for a system in which memory is included using distributed delay, we may face difficulty 2

to interpret the sensitivity of the memory (how strong or weak) over the solution of the system. In this regard, Du et al. (2013) provides a justification of fractional order derivative, which can be interpreted as an index of memory i.e. if α ∈ (0, 1] is the order of the fractional derivative, then α → 0 implies the system has an ideal memory and α → 1 implies the system has no memory. Atangana and Secer (2013) also provides a justification of using fractional differential operator in modelling context by showing that fractional derivatives can provide a better account of both physical and engineering processes. The present investigation aims at proposing and studying a more accurate model of dengue transmission than those previously presented in the literature, when memory influences the dynamics. The memory is incorporated in the model by using fractional differential operator and order of the fractional derivative as an index of memory (Du et al., 2013). We follow Dokoumetzidis et al. (2010b) approach to fractionalize a dengue ODE system to consider memory only in the dengue transmission process (human-to-vector and vector-to-human), and we study some aspects of the dynamics of this model. We also discuss possible policies to fight the dengue. The rest of the paper is arranged as follows. In section 2, we formulate the dengue model with memory in the transmission process. We study some basic mathematical properties and stability of equilibrium of this formulated dengue model in section 3. The numerical investigations of our analytical findings are carried out in section 4. The paper ends with a conclusion. 2. A dengue model In dengue dynamics, only adult female mosquitoes are involved in vector-host and host-vector transmission process. In the basic dengue ODE model, we consider the total full-grown female mosquito population at time t, M (t), is subdivided into two mutually exclusive sub-populations of susceptible female mosquitoes, MS (t), and infected female mosquitoes, MI (t). We have not considered any recovered class in female mosquito population since once infected with dengue, the mosquito will remain infected with the virus for its entire lifespan (Chamberlain and Sudia, 1961). In case of human population, we assume a total human population to be constant over time, denoted by H, is subdivided into three mutually exclusive sub-populations of susceptible humans, HS (t), infected humans, HI (t), and recovered humans, HR (t). The increment in susceptible adult female vector population is due to a constant recruitment of grown-up wing female mosquitoes at a rate, ΠV . The susceptible vector population is decreased due to those mosquitoes who acquired infection by biting infected human individuals and those who 3

die naturally at a rate, µm . The infected vector population is increased by the influx of infected vectors from the susceptible class and decreased by the natural deaths at a rate, µm . The susceptible human class is increased by the recruitment of newborn at a rate, µh and decreased by those who getting infected due to bite of infected vectors and those who die naturally at a rate, µh . Infected human class is increased by the newly recruited infected individuals from the susceptible class. Infected class is reduced due to natural recovery from dengue infection and natural deaths at rates, θh and µh , respectively. Recovered human class is increased due to the influx of recovered humans from the infected class at a rate, θh and decreased due to natural deaths of recovered individuals at a rate, µh . New infection in both the infection compartments (vector and host) depends upon average mosquitoes biting rate to host, dengue transmission probabilities (host-vector and vector-host) and the number of susceptible and infective in each population. Therefore, two forces of infection that h arises due to host to vector and vector to host transmission are ( bβHm )HI and ( bβ H )MI , respectively.

Since we assume closed human population (i.e. total human population is always constant), three human compartments can be reduced to two by using the conservation relation, defined as HR (t) = H − HS (t) − HI (t). Therefore, the basic dengue compartmental ODE model is given by dMS dt

= ΠV −

dMI dt

=

bβm H MS H I

bβm H MS H I

− µm MS ,

− µm MI ,

(2.1) dHS dt

= µh (H − HS ) −

dHI dt

=

bβh H H S MI

bβh H H S MI ,

− (µh + θh )HI .

All model parameters are assumed to be positive. Description of model (2.1) parameters and their biologically feasible ranges are given in Table 1. A first attempt to fractionalizing the dengue ODE model was proposed by Pooseh et al. (2011). Diethelm (2013) pointed out the flaws of the fractionalization of the dengue model proposed by Pooseh et al. (2011). The author first replaced all the fractional derivatives from Riemann-Liouville definition to Caputo’s definition (see, Podlubny (1999) for definition of different fractional order derivatives). The main advantage using Caputo’s fractional derivative is that the classical initial conditions can be used without encountering any problem during solvability. Furthermore, there is 4

a dimension mismatch in left-hand side and right-hand side of the proposed model of Pooseh et al. (2011). Diethelm (2013) solved this issue by changing the time-dependent parameters to their corresponding fractional order. Therefore, if we fractionalize the system (2.1) following Diethelm (2013) approach, we obtain

Dα MS = ΠαV −

D α MI

=

bα β m H MS H I

bα β m H MS H I

− µαm MS ,

− µαm MI ,

(2.2) D β HS

= µβh (H − HS ) −

D β HI

=

bβ β h H H S MI

bβ β h H H S MI ,

− (µβh + θhβ )HI .

Here 0 < α, β ≤ 1 and Dα y(t) is the Caputo fractional derivative of order α, defined as Z t 1 α p−α (p) p (t−s)p−1 y(s)ds. D y(t) := J y (t), p = [α] i.e. nearest integer value of α and J y(t) = Γ (p) 0

Note- Caputo’s fractional derivative in the vector and the human compartments are different, since dynamics of human and mosquito are not same. The main disadvantage in the Diethelm (2013) way of fractionalization is that when we consider a compartmental model with two or more compartments, an outgoing mass flux is an incoming flux to the next compartment. Thus, an outgoing mass which is a rate of some fractional order, cannot appear as an incoming mass to another compartment, as a rate of different fractional order, without violating mass balance (Dokoumetzidis et al., 2010a,b). For example in system (2.2), the infected human populations are in fractional time with dimension (time)−β . They are in interaction with susceptible mosquitoes, which are in fractional time with dimension (time)−α produces newly infected mosquitoes, which are in fractional time with dimension (time)−α . This produces an inconsistency in mass balance for the system (2.2). It is therefore, impossible to fractionalize multi-compartmental system by simply changing the integer order derivatives to fractional order in the left-hand side of the system (Dokoumetzidis et al., 2010a,b). Another issue regarding the fractionalization of a compartmental dengue model proposed by Diethelm (2013) is that, to avoid the dimension mismatch, the author transformed each rate parameter to their corresponding fractional order. There are good reasons to believe that the infection/transmission processes follow a power-law kind of dynamics, since in experimental studies on blood feeding vectors, it is pointed 5

out that memory and associated learning mechanism in mosquitoes play a crucial role in vector borne disease transmission (Chilaka et al., 2012; Vinauger et al., 2011; McCall and Kelly, 2002; Takken and Verhulst, 2013). Selecting a host location and host choice, mosquito uses its previous experience (Hii et al., 1991; Kelly, 2001; Kelly and Thompson, 2000). However, this is not the case for birth and death. Recruitment of newborns and mortality are random process and should not have a fractional order time dynamics. The trick used in the model (2.2) to avoid the dimension mismatch does not solve this issue. This problem was carefully discussed by Dokoumetzidis et al. (2010b) in the context of pharmacokinetics. To fractionalize a compartmental system Dokoumetzidis et al. (2010b) first integrate the ODE system to obtain a system of integral equations and then fractionalize these integrals by using a power-law kernel and hence obtain Riemann-Liouville integrals. Following Dokoumetzidis et al. (2010b) approach to fractionalize the dengue ODE system (2.1) we have the following system of equations

dMS dt

= ΠαV 11

dMI dt

=

(t − τ )α11 −1 bα12 βm 1−α12 (MS HI ) − µαm13 R D1−α13 (MS ), − RD H Γ (α11 )

bα12 βm 1−α12 (MS HI ) − µαm22 R D1−α22 (MI ), RD H

(2.3) dHS dt

= µαh 31 H

dHI dt

=

(t − τ )α31 −1 bα32 βh 1−α32 (HS MI ) − µαh 33 R D1−α33 (HS ), − RD H Γ (α31 )

bα32 βh 1−α32 (HS MI ) − µαh 42 R D1−α42 (HI ) − θhα43 R D1−α43 (HI ). RD H

In equation (2.3), 0 < αij ≤ 1 and R Dα y(t) represents Riemann-Liouville fractional derivative, defined as follows

dk 1 R D y(t) = Γ (k − α) dtk α

Z

0

t

(t − s)k−α−1 y(s)ds,

(k − 1 ≤ α < k).

Since birth and death in human and vector do not possess memory, we can take α11 = α13 = α22 = α31 = α33 = α42 = α43 = 1 and also we simplify our assumption that host-vector and vector-host dengue transmission follow a similar dynamics, so we take 0 < α12 = α32 = α < 1. Applying this assumptions to the system (2.3), we have the following dengue model with memory in transmission process

6

dMS dt

= ΠV −

dMI dt

=

bα βm 1−α (MS HI ) − µm MS , RD H

bα βm 1−α (MS HI ) − µm MI , RD H

(2.4) dHS dt

= µh (H − HS ) −

dHI dt

=

bα βh

H

RD

1−α

(HS MI ),

bα βh 1−α (HS MI ) − (µh + θh )HI . RD H

In the system (2.4) fractional derivatives in the right-hand side of each equation are in RiemannLiouville (RL) sense. One major draw back using RL derivatives is that the RL derivative of a constant is not zero. Therefore, we may face difficulty showing bounded ness of the variables. Furthermore, initial conditions for solving the system with RL derivative is involved with derivative of the variable which has no practical application. One way out of this problem is by fractional derivatives in Caputo’s sense which accepts usual initial conditions during solution and also Caputo’s derivative of a constant is zero. We have the following relation between RL and Caputo’s derivative of a variable

1−α A(t) = D1−α A(t) + RD

A(0)tα−1 . Γ (α)

(2.5)

Using the relation (2.5), the system (2.4) become

dMS dt

= ΠV −

dMI dt

=

C1 tα−1 bα βm 1−α , D (MS HI ) − µm MS − Γ (α) H

C1 tα−1 bα βm 1−α , D (MS HI ) − µm MI + Γ (α) H

(2.6)

where, C1 =

dHS dt

= µh (H − HS ) −

dHI dt

=

bα βh

H

D1−α (HS MI ) −

tα−1

C2 , Γ (α)

C2 tα−1 bα βh 1−α , D (HS MI ) − (µh + θh )HI + Γ (α) H

bα βh bα βm HS (0)MI (0). MS (0)HI (0) and C2 = H H

7

The model (2.6) has a type-II finite-time singularity at t = 0, which makes the system a-priori ill defined. This singularity occurs due to using Riemann-Liouville fractional derivative in (2.4), which have a singularity at t = 0 (Chikrii and Matychyn, 2011). However, for Caputo fraction derivative this singularity does not occur (Chikrii and Matychyn, 2011). Now, as 0 < α < 1, from the equation (2.5), as t → ∞, the definitions become equivalent. Thus, for the study of the behavior of solution near steady state in a dynamical process, the Riemann-Liouville definition and the Caputo’s definition must give the same results (page 81, Podlubny (1999)). This result is proven by many other authors also during the study of the stability of different equilibriums of a dynamical system (see, Qian et al. (2010); Zhang et al. (2011)). The main focus of this study is to analyze the long term behavior of the solution near steady states of the system (2.6). Therefore, we can neglect the time-dependent decay term in the righthand side of the system (2.6) and arrive at the following autonomous differential equation with Caputo fractional derivatives:

dMS dt

= ΠV −

dMI dt

=

bα βm 1−α D (MS HI ) − µm MS , H

bα βm 1−α D (MS HI ) − µm MI , H

(2.7) dHS dt

= µh (H − HS ) −

dHI dt

=

bα βh

H

D1−α (HS MI ),

bα βh 1−α D (HS MI ) − (µh + θh )HI , H

System (2.7) does not have any singularity at t = 0. We will analyze stability properties of the system (2.7) in the following section. 3. Basic Mathematical properties of the dengue model (2.7) We claim the following result. Lemma 1. The compact set

ΠV ; HS + HI ≤ H} µm is positively invariant and every forward solution of the system (2.7) is ultimately bounded.

Ω = {(MS , MI , HS , HI ) ∈ R4+ : MS + MI ≤

Proof: Adding first two equations and last two equations of the system (2.7) respectively, we have, 8

dM dt

= ΠV − µm M,

(3.1)

= µh H − µh N − θh HI ,

(3.2)

and

dN dt

6 µh H − µh N, where, M = MS + MI and N = HS + HI . Solving the equation (3.1) we have,

M (t) =

  ΠV ΠV e−µm t . + M (0) − µm µm

(3.3)

Again, from (3.2) using comparison theorem (Lakshmikantham et al., 1989) we have,

N (t) 6 H + (N (0) − H) e−µh t .

(3.4)

ΠV ΠV and N (t) 6 H. and N (0) 6 H then we have, M (t) 6 µm µm Therefore, the compact set Ω is positively invariant with respect to the dengue system (2.7). From (3.3) and (3.4), if M (0) 6

Any solutions with initial conditions within Ω will stays in Ω for all t and therefore bounded. ΠV and N (0) > H, from (3.3), we have, Now, for any initial condition outside Ω i.e. M (0) > µm ΠV and therefore M (t) is ultimately every forward solution of the equation (3.3) converges to µm bounded. Therefore, there exist t1 > 0 and a positive constant L1 such that ∀ t > t1 , M (t) 6 L1 .

Again, for any initial condition N (0) > H, from (3.4), there exists T > t1 > 0 such that ∀ t > T , N (t) 6 H. Now, assume L = max{L1 , H}, therefore, we have, ∀ t > T , M (t) 6 L and N (t) 6 L. Thus, every forward solution of the system (2.7) is ultimately bounded.  3.1. Existence and stability of equilibria of the dengue model (2.7) Fractional derivatives does not satisfy Leibniz rule (Tarasov, 2013), therefore, we cannot infer stability of the equilibrium of the system (2.7) directly from the equation (2.7). We apply the following transformation:

9

X

= D α MS +

bα βm MS H I , H

Y

= D α MI −

bα βm MS H I , H

(3.5) Z

= D α HS +

bα βh

H

H S MI ,

bα βh H S MI . H Using these transform variables in equation (3.5), the dengue system (2.7) is now equivalent to W

= D α HI −

the following fractional order compartmental system:

D1−α X = ΠV − µm MS ,

D1−α Y = −µm MI ,

D1−α Z = µh (H − HS ),

D1−α W = −(µh + θh )HI ,

bα βm MS H I , H bα βh H S MI . D α HI = W + H

bα βm MS H I , H bα βh H S MI , D α HS = Z − H

D α MI = Y +

D α MS = X −

(3.6) Since system (2.7) is equivalent to the system (3.6), it is enough to study the stability properties of the system (3.6). System (3.6) has only one equilibrium point E0 = (X ∗ , Y ∗ , Z ∗ , W ∗ , MS∗ , MI∗ , HS∗ , HI∗ ) = (0, 0, 0, 0,

ΠV , 0, H, 0), µm

which corresponds to the disease-free equilibrium of the system (2.7). To study the stability of the equilibrium point E0 ,  0 0 0   0 0 0    0 0 0    0 0 0  JE0 =   1 0 0    0 1 0     0 0 1  0 0 0

we construct the Jacobian matrix J ofthe system (3.6) at E0 , we have 0 −µm 0 0 0    0 0 −µm 0 0    0 0 0 −µh 0   0 0 0 0 −(µh + θh )    bα β m Π V  0 0 0 0 − Hµm    bα β m ΠV 0 0 0 0  Hµm   0 0 −bα βh 0 0   1 0 bα βh 0 0

10

The transformed dengue system (3.6) is a multi-order fractional compartmental system (except for α =

1 2 ).

Although α can be rational or irrational or even a complex number but there is

no existing method in literature for analyzing stability of a multi-order fractional compartmental system with the order of the fractional derivatives which being either irrational or complex number. M , where M, N > 0, N > M , and Therefore, we assume that 0 < α < 1 is rational i.e. α = N gcd (M, N ) = 1.

Stability of the disease-free equilibrium E0 is determined by the roots of the following characteristic equation of JE0 :

△(JE0 − diag([λ(N −M ) , λ(N −M ) , λ(N −M ) , λ(N −M ) , λ(M ) , λ(M ) , λ(M ) , λ(M ) ])) = 0,

(3.7)

where, △ and diag denote the determinant and the diagonal matrix, respectively. The stability of the disease-free equilibrium point E0 can be determine using the following result by Tavazoei and Haeri (2008): Lemma 2. (Theorem 2, Tavazoei and Haeri (2008)) Let 0 < pq11 =: α 6= β := pq22 < 1, where pi , qi ∈ Z+ , gcd(pi , qi ) = 1 and qi > pi , for i = 1, 2, and let L := lcm(q1 , q2 ). Then disease-free π for all roots λ of the equilibrium E0 of the system (3.6) is asymptotically stable if | arg(λ)|
> 2L Lα Lα Lα Lα Lβ Lβ Lβ Lβ following equation △ JE0 − diag[λ λ λ λ λ λ λ λ ] = 0.

We claim the following result.

Theorem 1. The disease-free equilibrium E0 of the system (3.6) is locally asymptotically stable if all roots λ of the equation:   2α 2N N 2N (1−α) b βm βh ΠV =0 λ + λ (µm + µh + θh ) + µm (µh + θh ) − λ Hµm

π , ∀j = 1, 2, ..., 2N , and is unstable if at least one root λi of the above 2N π . If there exist a pair of roots (λp1 , λp2 ) of the above polynomial equation satisfies | arg (λi )| < 2N π , for some α∗ = M polynomial equation settled in the critical region i.e. | arg (λpi )| = N ∈ (0, 1), 2N ∗ then system (3.6) undergoes a Hopf-type bifurcation at α .

satisfy | arg (λj )| >

Proof: Expanding the characteristic equation (3.7) in factorized form, we have the following equation:

N

λ + µm




N

λ + µh




  2α  2N N 2N (1−α) b βm βh ΠV = 0. λ + λ (µm + µh + θh ) + µm (µh + θh ) − λ Hµm (3.8) 11

Using lemma 2, the disease-free equilibrium (E0 ) is asymptotically stable if all roots λ of the π . equation (3.8) satisfy | arg (λ)| > 2N Arguments of the roots of the first factor in equation (3.8) has the following form π π 2kπ π ∀k = > , k = 0, 1, 2, ..., (N − 1). Therefore, | arg (λk )| ≥ + arg (λk ) = 2N N N N 0, 1, 2, ..., (N − 1). π , where, k = Thus, all roots of the first factor in equation (3.8) satisfy | arg (λk )| > 2N π , 0, 1, 2, ..., (N − 1). Similarly, all roots of the second factor in equation (3.8) satisfy | arg (λl )| > 2N where, l = 0, 1, 2, ..., (N − 1).

Therefore, stability of E0 of the system (3.6) depends upon the nature of the roots of the remaining factor in the characteristic equation (3.8), which given as follows:

 b2α βm βh ΠV = 0. (3.9) λ + λ (µm + µh + θh ) + µm (µh + θh ) − λ Hµm π , ∀j = 1, 2, ..., 2N , then Therefore, if all roots λj of the equation (3.9) satisfy | arg (λj )| > 2N the disease-free equilibrium (E0 ) of the system (3.6) is locally asymptotically stable and unstable π . if at least one root (λi ) of the equation (3.9) satisfies | arg (λi )| < 2N We now study the existence of Hopf-type bifurcation for the system (3.6). The criterion for 2N

2N (1−α)

N



Hopf-type bifurcation for the system (3.6) is that if for some root λp of the equation (3.9) satisfies π , for some α∗ = M | arg (λp )| = N ∈ (0, 1). Then system (3.6) under-goes a Hopf-type bifurcation 2N at α∗ (Min and Xing, 2012).

To find a sufficient conditions for the existence of a Hopf-type bifurcation for the system (3.6), iπ

let us assume that there exists a root λ of the equation (3.9) has the form λ = re 2N , where r is √ real and positive and i = −1. Substituting λ in equation (3.9) we have following two equation in

terms of the variable r

r

2N

+r

2N (1−α)



b2α βm βh ΠV Hµm

N

r (µm + µh + θh ) − r

2N (1−α)



cos [π(1 − α)] − µm (µh + θh ) = 0,



b2α βm βh ΠV Hµm



sin [π(1 − α)] = 0.

These two equations in r can be transform to the following polynomial equation in r

r

4N

−r

4N (1−α)



b2α βm βh ΠV Hµm

2

+ r2N [µ2m + (µh + θh )2 ] + µ2m (µh + θh )2 = 0. 12

(3.10)

Using Descartes’ rule of signs, the polynomial (3.10) has either two or zero positive root. If the iπ

polynomial (3.10) has two positive roots, which are conjugate root of the form λ = re± 2N of the

equation (3.9) for some α∗ =

M N

∈ (0, 1). Then system (3.6) undergoes a Hopf-type bifurcation at

α∗ .  Remark 1: We have mentioned the possible bifurcations of the solution of the dengue system (3.6) as Hopf-type bifurcation. One of the basic difference between dynamical properties of integer-order systems and fractional-order systems is that the limit set of a trajectory of an integer order system as the limit cycle of this system is a solution for this system. However, for fractional-order system, the limit set of a trajectory may not be a solution for this system (see, Tavazoei et al. (2009b); Abdelouahab et al. (2012)). Tavazoei et al. (2009a) showed that there are no periodic orbits in fractional order systems, and Tavazoei (2010) gave an example where the solutions of a fractional-order system are not periodic but converge to periodic signals. Remark 2: An interesting fact to note down is that whether we can get back the result of the classical integer order dengue system (α = 1) from our findings. At α = 1, equation (3.9) become

   2 b βm βh ΠV = 0. (3.11) λ + λ(µm + µh + θh ) + µm (µh + θh ) 1 − Hµ2m (µh + θh ) s b2 βm βh ΠV , which is dimensionless. It is epidemioWe define a threshold quantity R0 = Hµ2m (µh + θh ) logically known as the basic reproduction number, which measure the potential for the disease to 2

spread in a population. From (3.11), if R0 < 1, then all roots of the equation (3.11) have negative real part i.e. Re(λj ) < 0, for j = 1, 2. Therefore, the disease-free equilibrium (E0 ) is locally asymptotically stable. Again, if R0 > 1, then one root (λp ) of the equation (3.11) satisfies Re(λp ) > 0, therefore, E0 is unstable and disease persist in the population. Also, from equation (3.10) in Theorem 1, for α = 1 and R0 > 1, the equation (3.10) does not have any positive root. Therefore, when α = 1 and R0 > 1, the system (2.7) does not under-go a Hopf-type bifurcation. Actually, when α = 1 and R0 > 1, the system (2.7) will have an endemic equilibrium E1 = (MS∗ , MI∗ , HS∗ , HI∗ ), which is locally asymptotically stable. Therefore, for α = 1, we get back the classical epidemiological result from our system. If R0 < 1, the disease-free equilibrium (E0 ) is locally asymptotically stable and if R0 > 1, then E0 is unstable and disease persist in the population (Diekmann et al., 1990; Diekmann and Heesterbeek, 2000; van den Driessche and Watmough, 2002).

13

Remark 3: For α ∈ (0, 1), we can define a threshold quantity R0 , similar to the basic reproduction number for classical integer order dengue system and it is given as:

R0 =

s

b2α βm βh ΠV . Hµ2m (µh + θh )

(3.12)

Threshold quantity R0 defined in equation (3.12) is of dimension [time]1−α , hence it is not a dimensionless quantity. Therefore, we cannot call R0 as the basic reproduction number, which is a dimensionless quantity. R0 defined in (3.12) for fractional order system (3.6) is a memory dependent threshold quantity, as R0 ∝ bα . Therefore, the smaller α, the smaller bα , hence less effective the average bite rate per mosquito per month. 3.2. A special case when α =

1 2

When α = 12 , the multi-order fractional compartmental dengue model (3.6) becomes a single-

order fractional compartmental system. For this kind of system, we can easily infer our main stability result in Theorem 1 in terms of the threshold quantity R0 defined in equation (3.12). We claim the following result. Proposition 1. When the fractional order derivative is α = 12 , the critical value of R0 is determined q 1 + µ1m . If R0 < R0c , then the disease-free equilibrium (E0 ) of the system (3.6) is as R0c = (µh +θ h) locally asymptotically stable, and is unstable when R0 exceeds this critical value. When R0 = R0c , then E0 is stable and the system (3.6) under-goes a Hopf-type bifurcation.

Proof: Substituting α =

1 2

in equation (3.9) and using the expression of the threshold quantity

R0 defined in equation (3.12), we have the following equation:   λ4 + λ2 (µm + µh + θh ) − µm (µh + θh )R02 + µm (µh + θh ) = 0.

(3.13)

Putting λ2 = Z in equation (3.13) we have:

  Z 2 + Z (µm + µh + θh ) − µm (µh + θh )R02 + µm (µh + θh ) = 0.

Roots of the above equation are: Z1,2

−{(µm + µh + θh ) − µm (µh + θh )R02 } ± = 2

√

D

,

where,  2 D = (µm + µh + θh ) − µm (µh + θh )R02 − 4µm (µh + θh ). s 1 1 . + Now, we define R0c = (µh + θh ) µm

14

(3.14)

If D > 0 and R0 < R0c , then both the roots Z1 and Z2 in (3.14) are real and negative. Therefore, roots of the equation (3.13) satisfy | arg (λj )| =

> π4 , for all j = 1, 2, 3, 4.

π 2

For D < 0 and R0 < R0c , we assume, D = −D1 (D1 > 0) and A = (µm + µh + θh ) − µm√ (µh + D1 ) θh )R02 (> 0). Therefore, argument of the roots Z1 and Z2 in (3.14) are arg (Z1 ) = π − arctan ( A √ D1 ). Thus, arguments of the four roots of the equation (3.13) are and arg (Z2 ) = −π + arctan ( A given as follows: √ D1 π 1 ), arg(λ1 ) = − arctan ( A 2 2

√ D1 3π 1 ), − arctan ( arg(λ2 ) = A 2 2

arg(λ3 ) = −

Since A > 0 and

√

(3.15)

√

D1 π 1 ), + arctan ( A 2 2

√ D1 π 1 ). arg(λ4 ) = + arctan ( A 2 2

D1 > 0, thus, 0 < arctan (

the equation (3.13) satisfy | arg(λj )| >

π 4

√

D1 A )

<

π 2.

Therefore, from (3.15 ) all roots λ of

for j = 1, 2, 3, 4.

Therefore, using Theorem 1, we have, for α =

1 2

and if R0 < R0c , then the disease-free equilibrium

(E0 ) of the system (3.6) is locally asymptotically stable. Now, when R0 = R0c , then argument of the roots Z1 and Z2 in (3.14) are arg (Zj ) = ± π2 , for

j = 1, 2. Therefore, using Theorem 1, we have, two roots (λp ) of the equation (3.13) falls in the critical region i.e. | arg (λp )| =

π 4

and other two roots of the equation (3.13) satisfy | arg (λ)| >

π 4.

Using Matignon (1996) stability criterion for fractional order systems we have, equilibrium E0 of the system (3.6) is stable when R0 = R0c . For, R0 > R0c , and D > 0, then all roots in equation (3.14) are positive. Therefore, two roots of the equation (3.13) have | arg (λ)| = 0 <

| arg (λ)| = π >

π 4.

π 4

and other two roots of the equation (3.13) have

Thus, using Theorem 1, when R0 > R0c , D > 0, and α =

the equation (3.13) exceed the critical region | arg (λ)| =

π 4.

1 2,

two roots of

Therefore, system (3.6) under-goes a

Hopf-type bifurcation (Min and Xing, 2012). Now, if R0 > R0c , and D < 0, then all roots in equation (3.14) are of the form: Z1,2 = √ √ A1 ± i D1 , i = −1, where, A1 = µm (µh + θh )R02 − (µm + µh + θh ) > 0 and D1 = −D√ > 0. D1 ). Therefore, arguments of the roots in equation (3.14) has the form arg (Z1,2 ) = ± arctan ( A1 Therefore, arguments of the four roots of the equation (3.13) are given as follows:

15

√ D1 1 ), arg(λ1 ) = arctan ( A1 2

√ D1 1 ) − π, arg(λ2 ) = arctan ( A1 2

(3.16)

√

D1 1 ), arg(λ3 ) = − arctan ( A1 2

√ D1 1 ) + π. arg(λ4 ) = − arctan ( A1 2

√

π D1 ) < , therefore, from (3.16), we have, two roots of the equation (3.13) 2 A1 crosses the stability margin | arg (λ)| = π4 , when R0 > R0c , and D < 0. Therefore, using Theorem 1 Since 0 < arctan (

the dengue system (3.6) under-goes a Hopf-type bifurcation.  Remark 4: If we choose, µm = 2.14 (month)−1 , µh = 0.0012 (month)−1 and θh = 4.15 (month)−1 , then from Proposition 1, the critical value of threshold quantity R0 below which dengue will be eradicated is become 0.8415. Therefore, in case of α =

1 2,

making R0 below unity is no longer

a sufficient condition for eradication of the disease, although it is necessary. Again, if we choose, µm = 0.7 (month)−1 , µh = 0.0012 (month)−1 and θh = 2.15 (month)−1 , then from Proposition 1, the critical value of threshold quantity R0 below which dengue will be eradicated is become 1.376. Therefore, in case of α = 12 , making R0 greater than unity is no longer a sufficient condition for

the persistence of dengue in the population, however it is necessary. Therefore, for α =

1 2,

our

dengue model (3.6) does not satisfy the classical epidemiological theory (Diekmann et al., 1990; Diekmann and Heesterbeek, 2000; van den Driessche and Watmough, 2002). This nature of the solution of our dengue system (3.6) in case of α = 12 , is quite alike that of the backward bifurca-

tion phenomena in integer order dengue model (Garba et al., 2008). However, in case of backward bifurcation, R0 > 1 is always a sufficient condition for the persistence of the disease (Garba et al., 2008). 4. Numerical study Using Theorem 1, we draw numerically the stability region of the disease-free equilibrium (E0 ) of the system (3.6) and the region where Hopf-type bifurcation occurs with respect to the threshold quantity R0 , defined in equation (3.12) and α ∈ (0, 1), see Figure 1. Fixed parameters in Figure 1 16

are µm , µh and θh as in the legend. The blue region in Figure 1 is the boundary of the two region. The red line in Figure 1 denotes R0 = 1. From Figure 1, we find that up to the value of R0 = 0.3 (approximate), the disease-free equilibrium (E0 ) of the system (3.6) is locally asymptotically stable irrespective of any value of α ∈ (0, 1). Crossing this threshold value of R0 , E0 becomes unstable and the system (3.6) under-goes a Hopftype bifurcation depending upon values of α ∈ (0, 1). Therefore, R0 < 1 may not be a sufficient condition for eradication of the disease for our dengue model (3.6), when α ∈ (0, 1). Again, from Figure 1, we find for R0 > 1.15 (approximate), E0 is always unstable and system (3.6) under-goes a Hopf-type bifurcation irrespective of any value of α ∈ (0, 1). However, when 1 < R0 < 1.15, E0 becomes stable depending upon the value of the memory parameter α. Therefore, R0 > 1, may not be a sufficient condition for persistent of the disease in population for the dengue model (3.6), when α ∈ (0, 1). It is interesting to note that for α ∈ (0, 1), the dengue model (3.6) in general does not satisfy the classical epidemiological theory, which is R0 < 1, the disease-free equilibrium (E0 ) is always stable and R0 > 1, E0 is always unstable (Diekmann et al., 1990; Diekmann and Heesterbeek, 2000; van den Driessche and Watmough, 2002). To verify Theorem 1 numerically (consequently for the model (2.6)), we draw time series of the solution of the dengue model (2.6) using numerical scheme given in Appendix 1. Figure 2(a) shows a Hopf-type bifurcation for the solution of the dengue model (2.6) when the threshold quantity R0 , defined in equation (3.12) is greater than unity. Parameters are as in the legend Figure 2(a). Using parameter set from the legend Figure 2(a), R0 = 2.5981 > 1 and the equation (3.10) in Theorem 1 has two positive roots. Therefore, α∗ falls in the critical region and the system (2.7) (consequently the system (2.6)) under-goes a Hopf-type bifurcation. Figure 2(b) shows a Hopf-type bifurcation for the solution of the dengue model (2.6) when the threshold quantity R0 , defined in equation (3.12) is less than unity. Parameters are as in the legend Figure 2(b). Using parameter set from the legend Figure 2(b), R0 = 0.9696 < 1 and the condition in Theorem 1 is satisfied. Therefore, system (2.7) (consequently the system (2.6)) under-goes a Hopf-type bifurcation. Asymptotic stability of the disease-free equilibrium (E0 ) of the dengue model (2.7) is depicted in Figure 3. Parameters are as in the legend Figure 3. Using parameter set from the legend Figure 3, the threshold quantity in equation (3.12) becomes, R0 = 0.4532, and E0 = (28662.42, 0, 30000, 0). The condition for the asymptotic stability of E0 in Theorem 1 is also satisfied for these set of parameters. 17

5. Conclusion and discussion We have formulated and analyzed a new kind of dengue model with memory in vector-host and host-vector transmission process. Memory is incorporated in the transmission process by using fractional differential operator. We have determined a threshold quantity R0 , similar to the basic reproduction number. Stability of the disease-free equilibrium E0 is determined with respect to R0 and the order of the fractional derivative α. We have determined a α dependent threshold values for R0 , below which E0 is locally asymptotically stable, above which E0 is always unstable and at which the system under-goes a Hopf-type bifurcation. Our results suggest that classical epidemiological requirement of making R0 less than unity to control the disease is no longer a sufficient condition, although it is necessary (Diekmann et al., 1990; Diekmann and Heesterbeek, 2000; van den Driessche and Watmough, 2002). This result is similar to the concept of backward bifurcation for classical dengue ODE model, where a stable disease-free equilibrium coexists with a stable endemic equilibrium when R0 < 1 (Garba et al., 2008). It is also shown that for our dengue model (2.7), R0 > 1 may not be a sufficient condition for the persistence of dengue in the population. For a particular case, when α = 12 , analytically we have determined the critical value

of R0 in terms of some important model parameters. In this case, R0 below this critical value, the disease-free equilibrium E0 of the dengue model (2.7) is always stable, above which E0 is always unstable, and at which the system exhibits a Hopf-type bifurcation. The threshold quantity R0 defined in equation (3.12) is memory dependent, as R0 ∝ bα . Only, memory dependent parameter in the expression of R0 in equation (3.12) is the average monthly biting rate of mosquito. Increasing memory (α → 0) will reduce the average monthly biting rate of mosquito. Also as the average monthly biting rate of mosquito increases, R0 becomes more sensitive to the transmission memory α. From the stability region in Figure 1, it is clear that the range of the set {R0 (α) : α ∈ (0, 1]} for which disease persists in the population will increase as the transmission memory parameter α decreases. An immediate conclusion is that as the transmission memory increases (α → 0) the probability of dengue persistence in a population will increase. This result is in agreement with the experimental findings that indicate the increase in memory and learning in vectors will increase the dengue transmission (McCall and Kelly, 2002; Chilaka et al., 2012; Vinauger et al., 2011; Chaves et al., 2010; Takken and Verhulst, 2013). It would be interesting to discuss the possible policies to fight the dengue in terms of the threshold quantity R0 . Mosquito control using insecticide, use of antibiotic when epidemic occurs and increasing human awareness to use of mosquito net or other mosquito repeller to reduce average 18

monthly biting rate of mosquito are the existing intervention strategies to control dengue (Gratz, 1991; World Health Organization, 2013a). It is interesting to discuss about how we can use these existing control strategies in an effective way to eradicate dengue from the population using the result from our model. From the expression of R0 in equation (3.12), the only parameter which is sensitive to transmission memory (α) is the average monthly biting rate of mosquito, b. Therefore, increasing human awareness in an effective way (proper distribution of mosquito net and other mosquito repeller) to reduce the parameter b (with no influence on the other parameters) in such a way so that the interval (inf {R0 (α) : α ∈ (0, 1]}, sup {R0 (α) : α ∈ (0, 1]}) becomes narrow. Reducing the average monthly biting rate of mosquito will reduce the sensitiveness of R0 to the value of α i.e. R0 becomes less sensitive to transmission memory. Now, we can apply other two interventions (increasing natural mortality rate of mosquito using insecticides and increasing natural recovery rate of human by using antibiotic when a person suffers from dengue fever) in an effective way such that sup {R0 (α) : α ∈ (0, 1]} becomes less than the critical threshold value of R0 (α) determined by the dengue model (2.7). In such case, dengue will be eradicated. Therefore, effective use of all three combinations of intervention, and the results obtained from the analysis of the model may reduce dengue burden from a population. Currently, there is no effective vaccination against dengue epidemic (Whitehorn and Farrar, 2010), but a number of candidate vaccines (including live attenuated mono- and tetra-valent formulation, inactivated whole virus vaccines, and recombinant subunit vaccines) are undergoing various

phases of clinical trials (Blaney Jr et al., 2005; Center For Disease Control, 2007; Center For Vaccine Developm 2007; Dung et al., 1999). Therefore, it would be interesting to discuss about the impact of a vaccine on the trend of R0 . We believe that a good vaccine should focus on minimizing the parameter θh (with no influence on the other parameters) and therefore, make the quantity inf {R0 (α) : α ∈ (0, 1]} higher, and as a result a good policy of prevention can be used. However, increasing only the quantity inf {R0 (α) : α ∈ (0, 1]} may not help since R0 becomes very sensitive to the value of α i.e. range of the interval (inf {R0 (α) : α ∈ (0, 1]}, sup {R0 (α) : α ∈ (0, 1]}) may be very high depending on the value of the average monthly biting rate of mosquito. As discussed earlier, increasing human awareness can play an important role in this respect, which will make the interval (inf {R0 (α) : α ∈ (0, 1]}, sup {R0 (α) : α ∈ (0, 1]}) narrower, hence R0 becomes less sensitive to transmission memory α. Therefore, vaccination which focus on minimizing human recovery rate (θh ) with combination of human awareness (using mosquito net or other mosquito repeller) may reduce the dengue burden from a population. 19

The present investigation has some limitations and may be extended from different perspectives. In this paper, we assume that vector-host and host-vector transmission have a same memory (the fractional parameter α value fixed to be same). However, it will be more realistic to consider vector-host and host-vector transmission follow different dynamics i.e. α12 may not be same as α32 in equation (2.3). Moreover, most of the earlier studies on dengue dynamics consider only the adult stages of the female mosquitoes ignoring the previous aquatic phases of the mosquitoes (see, Feng and Velasco-Hern´andez (1997); Esteva and Vargas (1998, 1999, 2003); Chowell et al. (2007); Garba et al. (2008) etc). This assumption also stands in our study. Considering the dynamics of the aquatic phases of mosquito population with varying human population size may provide interesting dynamics. We leave these problems for our next study. Acknowledgement The authors are grateful to the Editor Prof. Stefano Ruffo and the learned reviewers for their comments and suggestions on the earlier version of the paper. The comments immensely improve the standard of the Manuscript. Special thanks to the learned third reviewer of the manuscript for providing a wonderful review. Tridip Sardar is supported by the research fellowship from Council of Scientific and Industrial Research (CSIR), Government of India. The Funder had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. References Abdelouahab, M., Hamri, N.-E., Wang, J., 2012. Hopf bifurcation and chaos in fractional-order modified hybrid optical system. Nonlinear Dynam. 69, 275–284. Acharya, A., Goswami, K., Srinath, S., Goswami, A., 2005. Awareness about dengue syndrome and related preventive practices amongst residents of an urban resettlement colony of south Delhi. J. Vect. Borne Dis. 42, 122–27. Agarwal, R., Ntouyas, S., Ahmad, B., Alhothuali, M., 2013. Existence of solutions for integrodifferential equations of fractional order with nonlocal three-point fractional boundary conditions. Adv. Difference Equ. 128, 1–9. Ahmed, E., Elgazzar, A., 2007. On fractional order differential equations model for nonlocal epidemics. Physica A 379(2), 607–14. 20

Andraud, M., Hens, N., Marais, C., Beutels, P., 2012. Dynamic epidemiological models for dengue transmission: a systematic review of structural approaches. PLoS ONE 7(11):e49085. Atangana, A., Secer, A., 2013. A note on fractional order derivatives and table of fractional derivatives of some special functions. Abstr. Appl. Anal. 2013, 1–8. Bartley, L., Donnelly, C., Garnett, G., 2002. The seasonal pattern of dengue in endemic areas: mathematical models of mechanisms. Trans. R. Soc. Trop. Med. Hyg. 96(4), 387–397. Blaney Jr, J., Matro, J., Murphy, B., Whitehead, S., 2005. Recombinant, liveattenuated tetravalent dengue virus vaccine formulations induce a balanced, broad, and protective neutralizing anitibody response against each of the four serotypes in rhesus monkeys. J. Virol. 79(9), 5516–5528.

Center For Disease Control, 2007. Dengue fact sheet. http://www.cdc.gov/ncidod/dvbid/dengue/resources [Online; accessed 21-May-2014]. Center

For

Vaccine

Development,

2007.

Live

attenuated

tetravalent

den

vaccine.

http://www.denguevaccines.org/live-attenuated-vaccines, [Online; accessed 21-May2014]. Chamberlain, R., Sudia, W., 1961. Mechanism of transmission of viruses by mosquitoes. Ann. Rev. Entomol. 6, 371–390. Chaves, L., Harrington, L., Keogh, C., Nguyen, A., Kitron, U., 2010. Blood feeding patterns of mosquitoes: random or structured? Front. Zool. 7:3, 1–11. Chen, A., Huang, L., Liu, Z., Cao, J., 2006. Periodic bidirectional associative memory neural networks with distributed delays. J. Math. Anal. Appl. 317(1), 80–102. Chikrii, A., Matychyn, I., 2011. Riemann-Liouville, Caputo, and Sequential fractional derivatives in differential games. Vol. 11. Birkhuser Boston. Chilaka, N., Perkins, E., Tripet, F., 2012. Visual and olfactory associative learning in the malaria vector Anopheles gambiae sensu stricto. Malaria J. 11:27, 1–11. Chowell, G., Diaz-Due¯ nas, P., Miller, J., Alcazar-Velazco, A., Hyman, J., et al., 2007. Estimation of the reproduction number of dengue fever from spatial epidemic data. Math. Biosci. 208(2), 571–589.

21

Diekmann, O., Heesterbeek, J., 2000. Mathematical epidemiology of infectious diseases: model building, analysis and interpretation. Wiley, New York. Diekmann, O., Heesterbeek, J., Metz, J., 1990. On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations. J. Math. Biol. 28, 365–382. Diethelm, K., 2013. A fractional calculus based model for the simulation of an outbreak of dengue fever. Nonlinear Dynam. 71(4), 613–619. Dietz, K., 1975. Transmission and control of arbovirus diseases. In: Ludwig, D., Cooke, K. (Eds.), Epidemiology. SIAM, Philadelphia, pp. 104–121. Dokoumetzidis, A., Magin, R., Macheras, P., 2010a. A commentary on fractionalization of multicompartmental models. J. Pharmacokinet. Pharmacodyn. 37, 203–207. Dokoumetzidis, A., Magin, R., Macheras, P., 2010b. Fractional kinetics in multi-compartmental systems. J. Pharmacokinet. Pharmacodyn. 37, 507–524. Du, M., Wang, Z., Hu, H., 2013. Measuring memory with the order of fractional derivative. Sci. Rep. 3, 1–3. Dung, N., Day, N., Tam, D., Loan, H., Chau, H., et al., 1999. Fluid replacement in dengue shock syndrome: a randomized, double-blind comparison of four intravenous-fluid regimens. Clin. Infect. Dis. 29(4), 787–94. Esteva, L., Vargas, C., 1998. Analysis of a dengue disease transmission model. Math. Biosci. 150(2), 131–151. Esteva, L., Vargas, C., 1999. A model for dengue disease with variable human population. J. Math. Biol. 38(3), 220–240. Esteva, L., Vargas, C., 2003. Coexistence of different serotypes of dengue virus. J. Math. Biol. 46(1), 31–47. Esteva, L., Yang, H., 2005. Mathematical model to assess the control of Aedes aegypti mosquitoes by the sterile insect technique. Math. Biosci. 198, 132–147.

22

Feng, Z., Velasco-Hern´ andez, J., 1997. Competitive exclusion in a vector-host model for the dengue fever. J. Math. Biol. 35(5), 523–544. Garba, S., Gumel, A., Bakar, A., 2008. Backward bifurcations in dengue transmission dynamics. Math. Biosci. 215, 11–25. Gopalsamy, K., He, X.-Z., 1994. Delay-independent stability in bidirectional associative memory networks. IEEE Trans. Neural Networks 5(6), 998–1002. Gratz, N., 1991. Emergency control of Aedes aegypti as a disease vector in urban areas. J. Am. Mosq. Control Assoc. 7(3), 353–65. Gubler, D., 1998. Dengue and dengue hemorrhagic fever. Clin. Microbiol. Rev. 11(3), 480–496. Hanert, E., Schumacher, E., Deleersnijder, E., 2011. Front dynamics in fractional-order epidemic models. J. Theor. Biol. 279(1), 9–16. Hii, J., Chew, M., Sang, V., Munstermann, L., Tan, S., et al., 1991. Population genetic analysis of host seeking and resting behaviors in the malaria vector, Anopheles balabacensis (Diptera: Culicidae). J. Med. Entomol. 28(5), 675–84. Kelly, D., 2001. Why are some people bitten more than others? Trends Parasitol. 17(12), 578–81. Kelly, D., Thompson, C., 2000. Epidemiology and optimal foraging: modelling the ideal free distribution of insect vectors. Parasitology 120, 319–27. Lakshmikantham, V., Leela, S., Martynyuk, A., 1989. Stability analysis of nonlinear systems. Marcel Dekker Inc, New York and Basel. Matignon, D., 1996. Stability results for fractional differential equations with applications to control processing. In: Proceedings of the multiconference on computational engineering in systems and application IMICS. Vol. 2. IEEE-SMC, Lile, France, pp. 963–968. McCall, P., Kelly, D., 2002. Learning and memory in disease vectors. Trends Parasitol. 18(10), 429–33. Min, X., Xing, Z., 2012. Hopf-type bifurcation and synchronization of a fractional order Van der Pol oscillator. In: Proceedings of the 31st Chinese control conference. Hefei, China, pp. 193–198. Newton, E., Reiter, P., 1992. A model of the transmission of dengue fever with an evolution of the impact of ultra-low volume (ulv) insecticide applications on dengue epidemics. Am. J. Trop. Med. Hyg. 47(6), 709–20.

23

Oldham, K., Spanier, J., 1974. The fractional calculus: theory and applications of differentiation and integration to arbitrary order. Academic Press, New York. Pinho, S., Ferreira, C., Esteva, L., Barreto, F., Morato e Silva, V., et al., 2010. Modelling the dynamics of dengue real epidemics. Phil. Trans. R. Soc. A 368, 5679–93. Podlubny, I., 1999. Fractional differential equations. Academic Press, San Diego, California. Pooseh, S., Rodrigues, H., Torres, D., 2011. Fractional derivatives in dengue epidemics. In: Simos, T., Psihoyios, G., Tsitouras, C., Anastassi, Z. (Eds.), Numerical Analysis and Applied Mathematics, ICNAAM. American Institute of Physics, Melville, pp. 739–742. Qian, D., Li, C., Agarwal, R., Wong, P., 2010. Stability analysis of fractional differential system with Riemann-Liouville derivative. Math. Comput. Model. 52, 862–874. Rosen, L., Roseboom, L., Gubler, D., Lien, J., Chaniotis, B., 1985. Comparative susceptibility of mosquito species and strains to oral and parenteral infection with dengue and japanese encephalitis viruses. Am. J. Trop. Med. Hyg. 34(3), 603–15. Rosenbaum, J., Nathan, M., Ragoonanansingh, R., Rawlins, S., Gayle, C., et al., 1995. Community participation in dengue prevention and control: a survey of knowledge, attitudes, and practice in Trinidad and Tobago. Am. J. Trop. Med. Hyg. 53(2), 111–117. Schutz, G., Trimper, S., 2004. Elephants can always remember: exact long-range memory effects in a non-Markovian random walk. Phys. Rev. E 70, 045101. Sheppard, P., MacDonald, W., Tonn, R., Grab, B., 1969. The dynamics of an adult population of Aedes aegypti in relation to dengue haemorrhagic fever in Bangkok. J. Anim. Ecol. 38, 661–702. Sipahi, R., Atay, F., Niculescu, S.-l., 2007. Stability of traffic flow behavior with distributed delays modeling the memory effects of the drivers. SIAM J. Appl. Math. 68(3), 738–759. Southwood, T., Murdie, G., Yasuno, M., Tonn, R., Reader, P., 1972. Studies on the life budget of Aedes aegypti in Wat Samphaya, Thailand. Bull. World Health Organ. 46(2), 211–26. Stanislavsky, A., 2000. Memory effects and macroscopic manifestation of randomness. Phys. Rev. E 61(5), 4752–4759. Takken, W., Verhulst, N., 2013. Host preferences of blood-feeding mosquitoes. Annu. Rev. Entomol. 58, 433–53. Tarasov, V., 2013. No violation of the Leibniz rule. No fractional derivative. Commun. Nonlinear. Sci. Numer. Simulat. 18, 2945–2948. Tavazoei, M., 2010. A note on fractional-order derivatives of periodic functions. Automatica 46(5), 945–948. Tavazoei, M., Haeri, M., 2008. Chaotic attractors in incommensurate fractional order systems. Physica D 237(20), 2628 – 2637. Tavazoei, M., Haeri, M., Attari, M., 2009a. A proof for non existence of periodic solutions in time invariant fractional order systems. Automatica 45(8), 1886 – 1890.

24

Tavazoei, M., Haeri, M., Attari, M., Bolouki, S., Siami, M., 2009b. More details on analysis of fractional-order Van der Pol oscillator. J. of Vib. and Control 15(6), 803 – 819. Thome, R., Yang, H., Esteva, L., 2010. Optimal control of Aedes aegypti mosquito by the sterile insect technique and insecticide. Math. Biosci. 223(1), 12–23. van den Driessche, P., Watmough, J., 2002. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 180, 29–48. Vinauger, C., Buratti, L., Lazzari, C., 2011. Learning the way to blood: first evidence of dual olfactory conditioning in a bloodsucking insect, Rhodnius prolixus. I. Appetitive learning. J. Exp. Biol. 214, 3032–38. Watts, D., Burke, D., Harrison, B., Whitmire, R., Nisalak, A., 1987. Effect of temperature on the vector efficiency of Aedes aegypti for dengue 2 virus. Am. J. Trop. Med. Hyg. 36(1), 143–52. Whitehorn, J., Farrar, J., 2010. Dengue. Br. Med. Bull. 95, 161–173. World Health Organization, 2013a. Dengue and http://www.who.int/mediacentre/factsheets/fs117/en/index.html, 22-October-2013].

severe [Online;

World Health Organization, 2013b. Health topics http://www.who.int/topics/dengue/en/, [Online; accessed 22-October-2013].

dengue. accessed (dengue).

Zhang, F., Changpin, L., Chen, Y., 2011. Asymptotical stability of nonlinear fractional differential system with Caputo derivative. Int. J. Differ. Equ. 2011, 1–12.

25

Figures and Tables

2 Hopf−type bifurcation region

1.5

R0

R

0

Boundary of the two region

= 1

1

Asymptotic stability region of the disease−free equilibrium (E )

0.5

0

0

0

0.5 α

0.99

Figure 1: Stability region of the disease-free equilibrium (E0 ), and Hopf-type bifurcation region of the model (3.6) with respect to the parameters α and R0 (defined in equation (3.12)). The blue region is the boundary of the two region. The red line denotes R0 = 1. Fixed parameters used during plot are taken as µm = 2.74 (month)−1 , µh = 0.0012 (month)−1 and θh = 4.14 (month)−1 .

26

4

x 10

6 Infected Mosquito

Susceptible Mosquito

5

1.6 1.4 1.2 1

200

100

0

4 2 0

300

Infected Human

Susceptible Human

2

2

200

100

0

200

300

200

300

4

x 10

4

0

100

0

4

6

x 10

1.5 1 0.5 0

300

x 10

0

Time (Month)

100 Time (Month)

(a) x 10

15000 Infected Mosquito

Susceptible Mosquito

4

9.5 9 8.5 8

50

0

10000 5000 0

100

0

50

100

50

100

4

x 10

6000 Infected Human

Susceptible Human

4 3.5 3 2.5

0

50

100

Time (Month)

4000 2000 0

0

Time (Month)

(b) Figure 2: Existence of a Hopf-type bifurcation for the solution of the dengue model (2.6). Figure 2(a) represents a Hopf-type bifurcation when the threshold quantity R0 defined in equation (3.12) is greater than unity. Parameter values used to draw the Figure 2(a) are α = 0.62 > α∗ = 0.61, µm = 0.8 (month)−1 , µh = 0.0012 (month)−1 , θh = 2.14 (month)−1 , ΠV = 120000 (month)−1 , b = 10 (month)−1 , βm = 0.51, βh = 0.51 and H = 56000. Initial conditions used in simulation of Figure 2(a) are MS (0) = 150000, MI (0) = 50, HS (0) = 54400 and HI (0) = 750. With these parameters, R0 = 2.59813 > 1 and the condition of Hopf-type bifurcation in Theorem 1 is satisfied. Figure 2(b) represents a Hopf-type bifurcation when the threshold quantity R0 defined in equation (3.12) is less than unity. Parameter values and the initial conditions used to simulate Figure 2(b) are taken as: α = 0.53 > α∗ = 0.52, µm = 1.5 (month)−1 , µh = 0.21 (month)−1 , θh = 9 (month)−1 , ΠV = 141000 (month)−1 , b = 18.2 (month)−1 , βm = 0.52, βh = 0.52, H = 40000, MS (0) = 93000, MI (0) = 1000, HS (0) = 39500 and HI (0) = 500. With these parameters, R0 = 0.9696 < 1 and the condition of Hopf-type bifurcation in Theorem 1 is satisfied.

27

4

5 Infected Mosquito

Susceptible Mosquito

x 10 2.8662

2.8657 0

10

20

30

2 0

40

0

10

20

30

40

20

30

40

29,998

Infected Human

Susceptible Human

0.8 30000

0

10

20

30

0.6 0.4 0.2 0

40

Time (Month)

0

10

Time (Month)

Figure 3: Asymptotic stability of the disease-free equilibrium (E0 ) of the dengue model (2.7). Parameter values and initial conditions used for simulation are α = 0.86, µm = 3.14 (month)−1 , µh = 0.0012 (month)−1 , θh = 4.15 (month)−1 , ΠV = 90000 (month)−1 , b = 9 (month)−1 , βm = 0.5, βh = 0.128, H = 30000, MS (0) = 28657, MI (0) = 5, HS (0) = 30000 and HI (0) = 0. With these parameters, R0 = 0.4532 < 1, E0 = (28662.42, 0, 30000, 0) and the condition for the asymptotic stability of E0 in Theorem 1 is satisfied. Red line in two susceptible compartments (mosquito and human) representing the corresponding values of E0 .

28

Table 1: Description of the dengue model (2.1) parameters and their possible feasible ranges.

Parameter

Biological meaning

Range of values

Reference

ΠV

recruitment rate of adult female vector population

(12000 − 150000) month−1 (**)

b

average bite per mosquito per month

(9 − 30) month−1

βm

transmission probability from human to vector

(0.5 − 1)

µm

death rate of vector population

(0.6 − 7.5) month−1

µ−1 h

average host life expectancy = (host recruitment rate)−1

(50 − 75) (year)

(Andraud et al., 2012; Esteva and Vargas, 1998; Newton and Reiter, 1992; Garba et al., 2008; Feng and Velasco-Hern´andez, 1997; Diethelm, 2013) (Pinho et al., 2010; Newton and Reiter, 1992) (Andraud et al., 2012; Bartley et al., 2002; Newton and Reiter, 1992; Rosen et al., 1985; Watts et al., 1987) (Andraud et al., 2012; Bartley et al., 2002; Newton and Reiter, 1992; Sheppard et al., 1969; Southwood et al., 1972) (Andraud et al., 2012; Diethelm, 2013)

βh

transmission probability from vector to human

(0.1 − 1)

θh

dengue recovery rate in human population

(2.14 − 10) month−1

** The range for the vector recruitment rate (ΠV ) was used in most of the modelling studies.

29

(Andraud et al., 2012; Bartley et al., 2002; Newton and Reiter, 1992; Rosen et al., 1985; Watts et al., 1987) (Andraud et al., 2012; Pinho et al., 2010; Newton and Reiter, 1992)

Appendix1: Numerical scheme for solving dengue model (2.6) Following is the numerical scheme for solving dengue model (2.6)

MS (i) = MS (i − 1) + h{ΠV − µm MS (i − 1) −

MI (i) = MI (i − 1) + h{−µm MI (i − 1) +

C1 ((i − 1)h)(α−1)

Γ (α)

C1 ((i − 1)h)(α−1)

HS (i) = HS (i − 1) + h{µh (H − HS (i − 1)) −

Γ (α)

+

HI (i) = HI (i − 1) + h{−(µh + θh )HI (i − 1) +

C2 ((i − 1)h)

Γ (α)

bα β m

H

−

−(1−α)

h

H

bα β m

C2 ((i − 1)h)(α−1)

Γ (α)

−

i X

[MS (i)HI (i) +

(1−α)

cj

MS (i − j)HI (i − j)]},

j=1

−(1−α)

h

[MS (i)HI (i) +

i X

(1−α)

cj

MS (i − j)HI (i − j)]},

j=1

bα β h

H

(α−1)

−(1−α)

h

[HS (i)MI (i) +

i X

b βh

H

HS (i − j)MI (i − j)]},

j=1

α

+

(5.1) (1−α)

cj

−(1−α)

h

[HS (i)MI (i) +

i X

(1−α)

cj

HS (i − j)MI (i − j)]}.

j=1

Where, (MS (0), MI (0), HS (0), HI (0)) is the initial condition, time step size is h and i =   1, 2, ..., N , for N = Th . The model is simulate up to the time T . The parameter cαj is defined as (α)

c0

(α)

= 1 and cj

= (1 −

1+α (α) j )cj−1 ,

j ≥ 1. C1 and C2 are defined in equation (2.6).

Appendix2: A Matlab code for solving dengue model (2.6) A Matlab code for solving the dengue model (2.6) is given as follows:

function [T, y]=Dengue memory new(param, Y0, TSim) % time step size: h=0.05;

% number of iterations: n=round((TSim)/h)+1; %orders of fractional derivative alpha =param(1); beta =1−alpha; %Parameters of the dengue model Pi v=param(2);b=param(3);beta m =param(4);mu m =param(5); beta h =param(6);theta h=param(7);mu h = param(8);

%preallocation size M S =zeros(n,1);M I =zeros(n,1);H S =zeros(n,1);H I =zeros(n,1);c1=zeros(n,1); y =zeros(n,4);

%initial condition M S(1) = Y0(1); M I(1) =Y0(2); H S(1)= Y0(3); H I(1) = Y0(4); H=H S(1)+H I(1);

% binomial coefficients calculation: cp1=1; for j=1:n c1(j)=(1−(1+beta)/j)*cp1; cp1=c1(j);

30

end

for i=2:n % some commands before run the solver Options = optimset('Display', 'notify'); Options = optimset(Options,'TolX', 1e−100); Options = optimset(Options,'MaxIter', 3000); Options = optimset(Options,'MaxFunEvals', 10000);

% Solving Dengue model X1 = fminsearchbnd(@(X) (X(1)−M S(i−1)−h*(Pi v−mu m * M S(i−1)−... ((bˆ(alpha)* beta m)/H)*((M S(1)* H I(1))/gamma(alpha))*((i−1)*h)ˆ(−beta)−... ((bˆ(alpha)* beta m)/H)*hˆ(−beta)*(X(1).*X(4)+memo(M S. * H I,c1,i)))).ˆ2+... (X(2)−M I(i−1)−h*(−mu m * M I(i−1)+... ((bˆ(alpha)* beta m)/H)*((M S(1)* H I(1))/gamma(alpha))*((i−1)*h)ˆ(−beta)+... ((bˆ(alpha)* beta m)/H)*hˆ(−beta)*(X(1).*X(4)+memo(M S. * H I,c1,i)))).ˆ2+... (X(3)−H S(i−1)−h*(mu h *(H−H S(i−1))−... ((bˆ(alpha)* beta h)/H)*((H S(1)* M I(1))/gamma(alpha))*((i−1)*h)ˆ(−beta)−... ((bˆ(alpha)* beta h)/H)*hˆ(−beta)*(X(3).*X(2)+memo(H S. * M I,c1,i)))).ˆ2+... (X(4)−H I(i−1)−h*(−(mu h+theta h)* H I(i−1)+... ((bˆ(alpha)* beta h)/H)*((H S(1)* M I(1))/gamma(alpha))*((i−1)*h)ˆ(−beta) +... ((bˆ(alpha)* beta h)/H)*hˆ(−beta)*(X(3).*X(2)+... memo(H S. * M I,c1,i)))).ˆ2,... [M S(i−1);M I(i−1);H S(i−1);H I(i−1)],[0;0;0;0],[],Options);

M S(i)=X1(1);M I(i)=X1(2);H S(i)=X1(3);H I(i)=X1(4);

end for j=1:n y(j,1)=M y(j,2)=M y(j,3)=H y(j,4)=H

S(j); I(j); S(j); I(j);

end T=0:h:TSim;

The memo.m function is given as follows:

function [yo] = memo(r, c, k) % temp = 0; for j=1:k−1 temp = temp + c(j)*r(k−j); end yo = temp; end

The optimization function f minsearchbnd.m is freely available in MATLAB central Mathworks website. 31

• We formulate a dengue model with memory in transmission process. • Memory is incorporated in the model using fractional differential operator. • Under certain conditions classical epidemic result does not holds for our model. • This result is synergistic with the backward bifurcation concept in epidemiology.