Mathematics in Applied Sciences and Engineering https://ojs.lib.uwo.ca/mase Volume 1, Number 3, Septembe 2020, pp.207-223 https://doi.org/10.5206/mase/10745 DYNAMICAL ANALYSIS OF A FRACTIONAL-ORDER MODEL INCORPORATING FEAR IN THE DISEASE TRANSMISSION RATE OF SEVERE INFECTIOUS DISEASES LIKE COVID-19 CHANDAN MAJI AND DEBASIS MUKHERJEE Abstract. This paper deals with a fractional-order three-dimensional compartmental model with fear effect. We have investigated whether fear can play an important role or not to spread and control the infectious diseases like COVID-19, SARS etc. in a bounded region. The basic results on uniqueness, non-negativity and boundedness of the solution of the system are investigated. Stability analysis ensures that the disease-free equilibrium point is locally asymptotically stable if carrying capacity exceeds a certain threshold value. We have also derived the conditions for which endemic equilibrium is globally asymptotically stable that means the disease persists in the system. Numerical simulation suggests that the fear factor has an important role which is observed through Hopf-bifurcation. 1. Introduction Over the past few years, human civilization had repeatedly been plagued by attacks of various Coro- naviruses. They are a broad family of viruses, some of which can cause several human diseases, varying from cold to SARS (Severe Acute Respiratory Syndrome). Two important coronavirus disease out- breaks have already occurred in the past two decades : SARS in 2003 [15] and MERS (Middle East Respiratory Syndrome) in 2012 [7]. By the end of 2019, human civilization is once again infected by the coronavirus which has not been previously identified in humans. On February 11, 2020 the world Health Organization (WHO) declared the official name coronavirus disease 2019. It is abbreviated as COVID-19. The most common symptoms of this disease are fever, tiredness and dry cough. More scarcely, the disease can be serious and fatal. Aged person and individuals with other medical con- straints (such as asthma, diabetes or heart disease) may be vulnerable to becoming seriously ill. It is observed in December 2019, SARS COV-2 has been found as the causal factor in a series of critical cases of pneumonia originating in Wuhan Hubei Province China [48]. Since then it spreads rest of the world. As per report of WHO(2020), 210 countries/territories had confirmed cases of this disease. With the extremely high infection rate and high mortality, individuals naturally began worrying about the COVID-19. One distinctive nature of infectious disease related with other conditions is fear. As the disease outbreak is ongoing, a wave of fear is developed in the society [28]. The fear and worry are obvious as peoples concern their health. No one wants to get infected with a virus that has a relatively high risk of death [23]. Fear is presently connected with its transmission rate [32, 4, 41]. As the out- break of COVID-19 spreads to more countries and death toll rises, the uncertainty of what lies ahead is concerning. Though present treatment on COVID-19 throughout the world has mainly confined into infection control, effective vaccine and recovery rate [14], the psychosocial aspect is neglected. As countries worldwide have to work on diminishing the transmission rate of COVID-19, they should pay Received by the editors 28 May 2020; accepted 24 August 2020; online first 30 August 2020. 2000 Mathematics Subject Classification. 34K20, 92B05. Key words and phrases. COVID-19; fear effect; fractional order; stability; bifurcation. 207 208 C. MAJI AND D. MUKHERJEE attention on individuals fears to gain a society free of COVID-19. From scientific point of view, reactions to fear are normal and potentially beneficial. Fear of the present disease outbreak is understandable, not to mention almost universal. With the fact in mind, it is reasonable to develop a mathematical model that incorporate fear factor in the disease transmission rate to control the COVID-19. Already a lot of works have been done mathematically and numerically to give an efficient prediction on COVID 19 outbreaks [46, 34, 49, 39, 35, 44, 43]. The study of fractional order differential equations has received much interest to the researchers during the past few decades due to their ability to provide a good description of certain non-linear phenomena [22]. It is an extension of classical calculus that generalizes the order of derivatives and integrals to a non-integer order. Fractional calculus was first proposed by Leibnitz and Hospitals in 1965 [36]. The fractional-order models are more realistic than integer-order model as the system has nonlocal charac- teristic which is absent in integer-order model and it has greater degree of freedom. Another reason for considering the fractional order system is to address memory which exists in most biological systems but such effects are in fact neglected in integer-order system. Apart from, using fractional order differential equations can help us to reduce the errors arising from the neglected parameters in modelling real life phenomena [12]. However, due to progress of fractional calculus many researchers in different fields such as biology, physics, engineering, finance , medicine considered fractional calculus to develop their problems [10, 17, 38, 11, 16, 29, 30, 5, 18]. By using fractional-order derivative a lot of work has done in epidemiology [20, 8, 21, 40, 45]. Motivated from the above literature survey, we proposed a fractional order SEI model of COVID-19 outbreaks. The paper is organized as follows. In Section 2, the model is proposed. Basic definitions are pre- sented in Section 3. Existence, uniqueness, non-negativity and boundedness of solutions are shown in Section 4. Equilibria and their local stability are discussed in Section 5. Global stability of endemic equilibrium point is derived in Section 6. Hopf bifurcation phenomena are demonstrated in Section 7. Numerical simulations are carried out in Section 8. A brief discussion follows in Section 9. 2. Model Formulation In epidemic models, the bilinear incidence rate βSI is frequently used. Recent study [28] indicates that the fear effect will reduce the disease transmission rate, so we modify βSI by multiplying a factor f(α,I) which leads to f(α,I)βSI. Here, the parameter α represents the level of fear. For biological justification of α, I and f(α,I), it is appropriate to consider the following : f(0,I) = 1, f(α, 0) = 1, lim α→∞ f(α,I) = 0, lim I→∞ f(α,I) = 0, ∂f(α,I) ∂α < 0, ∂f(α,I) ∂I < 0. For convenience of analysis, we assume the following form for the fear effect f(α,I) = 1 1 + αI . Here, when there is no infected individuals, there is no reduction in the susceptible individuals due to the fear factor i.e f(α, 0) = 1. Based on above assumption, in this present paper, we formulate a three dimensional compartmen- tal model with fear effect with the help of fractional-order Caputo-type derivative which is given as DYNAMICAL ANALYSIS OF A FRACTIONAL-ORDER MODEL INCORPORATING FEAR 209 follows: cDµS(t) = b̂S ( 1 − S k̂ ) − β̂SI 1 + α̂I , cDµE(t) = β̂SI 1 + α̂I − (ĉ + d̂)E, (2.1) cDµI(t) = ĉE − γ̂I with initial conditions S(0) = S0 ≥ 0,E(0) ≥ 0 and I(0) ≥ 0, where µ ∈ (0, 1) and cDµ is the standard Caputo differentiation. Here, S(t) is the total density of the susceptible individuals, E(t) is the number of individuals to the infected but not infectious and I(t) denotes the infected individuals who are infected. b̂ is the net per capita growth rate of the susceptible individuals and k̂ is the environmental carrying capacity. α̂ represents the level of fear among the individuals which have the controlling effect not to spread the disease. In this classical endemic model the transmission coefficient for the disease is denoted by β̂. ĉ is the rate per unit time (day) that infected individuals become infectious. d̂ is the natural death rate of the exposed individuals. The infected individuals removed at a rate of γ̂, which include natural death of the infected population and the recovery rate of the hospitalized infectious individuals. The system (2.1) has some defects as regard to the time dimension because the right hand side ex- pressions have dimension (time)−1, whereas the left hand side expressions have dimension (time)−µ. The corrected form of system (2.1) is as follows: cDµS(t) = b̂µS ( 1 − S k̂ ) − β̂µSI 1 + α̂I , cDµE(t) = β̂µSI 1 + α̂I − (ĉµ + d̂µ)E, (2.2) cDµI(t) = ĉµE − γ̂µI with initial conditions S(0) = S0 ≥ 0,E(0) ≥ 0 and I(0) ≥ 0, where µ ∈ (0, 1) and cDµ is the standard Caputo differentiation. For convenience, we redefine the parameters as follows: b = b̂µ,k = k̂,β = β̂µ,α = α̂,c = ĉµ,d = d̂µ,γ = γ̂µ. Therefore, the modified system is as follows: cDµS(t) = bS ( 1 − S k ) − βSI 1 + αI , cDµE(t) = βSI 1 + αI − (c + d)E, (2.3) cDµI(t) = cE −γI with S(0) = S0 ≥ 0,E(0) ≥ 0 and I(0) ≥ 0. 3. Basic Definitions Fractional calculus is a powerful tool for mathematical modeling and it has a wild application in differ- ent field of sciences. Throughout this paper, we use a Caputo fractional-order derivative as the initial conditions of fractional differential equations with Caputo derivatives consider on the identical form 210 C. MAJI AND D. MUKHERJEE as for integer-order ones, which can be used in modelling and analysis. In this section, some basic definitions for fractional calculus have been presented. Definition 3.1. [37] The Riemann-Liouville fractional integral operator of order µ of a continuous function f ∈ L1[0,a], t ∈ [0,a] is presented as Iµf(t) = 1 Γ(µ) ∫ a 0 (t− τ)µ−1f(τ)dτ, where Γ(µ) is the Euler’s Gamma function. Definition 3.2. [37] The definition of Caputo’s fractional derivative of order µ for a function f ∈ Cn([0, +∞],R) is defined by cDµf(t) = 1 Γ(n−µ) ∫ a 0 (t− τ)n−µ−1f(n)(τ)dτ, where Γ(.) is the Euler’s Gamma function and the operator cDµ is known as ”Caputo differential op- erator of order µ”. t ≥ 0 and n is the positive integer such that n− 1 < µ < n,n ∈ N. Particularly, when 0 < µ < 1, cDµf(t) = 1 Γ(n−µ) ∫ t 0 f′(τ) (t− τ)µ dτ. Riemann-Liouville (R-L) was first introduced the idea of fractional derivative but in R-L fractional differential equation, initial value is usually taken in the form of fractional derivative, which is not appropriate in real sense whereas in Caputo fractional derivative, the derivative is not defined locally at time t, it depends on the total effects of the so called n-order integer derivative on the interval [0,s]. Thus it is reasonable to consider the variation of a system in which the instantaneous change rate depends on the past rate, which is known as ”memory effect”[33]. 4. Main Results In this section we shall discuss about existence, uniqueness, non-negativity and boundedness of the solutions for fractional order system (2.3). 4.1. Existence and Uniqueness. Before we prove the existence and uniqueness of the solution of system (2.3), we need the following Lemma. Lemma 4.1. [26] Define the system cDµx(t) = f(t,x), t > 0 (4.1) with initial condition x0, where µ ∈ (0, 1],f : [0,∞) × Ω → Rn, Ω ∈ Rn, then there exists a unique solution of (2.3) whenever f(t,x) follows locally Lipschitz condition with respect to x on [0,∞) × Ω. Theorem 4.2. For any non-negative initial conditions the fractional order system (2.3) admits a unique solution. Proof. Existence and uniqueness of system (2.3) will be shown in the region ∆ × (0,T] where ∆ = {(S,E,I) ∈ R3 : max(|S|, |E|, |I| ≤ M)}. Now, we follow the approach used in [27].We de- note X = (S,E,I) and X̄ = (S̄, Ē, Ī). Consider a mapping DYNAMICAL ANALYSIS OF A FRACTIONAL-ORDER MODEL INCORPORATING FEAR 211 H(X) = (H1(X),H2(X),H3(X)) and H1(X) = bS ( 1 − S k ) − βSE 1 + αI H2(X) = βSE 1 + αI − (c + d)E, H3(X) = cE −γI. (4.2) For X,X̄ ∈ D, it follows from equation (4.2) that ‖H(X) −H(X̄)‖ = |H1(X) −H1(X̄)| + |H2(X) −H2(X̄)| + |H3(X) −H3(X̄)| = ∣∣∣∣bS ( 1− S k ) − βSI 1+αI −bS̄ ( 1− S̄ k ) + βS̄Ī 1+αĪ ∣∣∣∣+ ∣∣∣∣ βSI1+αI −(c+d)E− βS̄Ī1+αĪ + (c+d)Ē ∣∣∣∣+ ∣∣∣∣cE−γI−cĒ +γĪ ∣∣∣∣ = ∣∣∣∣b(S−S̄)− bk (S2−S̄2)−β ( SI 1+αI − S̄Ī 1+αĪ )∣∣∣∣+ ∣∣∣∣β ( SI 1+αI − S̄Ī 1+αĪ ) −(c+d)(E−Ē) ∣∣∣∣+ ∣∣∣∣c(E−Ē)−γ(I−Ī) ∣∣∣∣ ≤ ( b + 2bM k + 2β α )|S − S̄| + (2M + γ)|I − Ī| + (2c + d)|E − Ē| ≤ L||X − X̄||, where L = max{b + 2β α + 2bM k , 2c + d, 2M + γ}. Hence, Lipschitz condition is satisfied for H(X). Thus there exist a unique solution X(t) of system (2.3), follows from Lemma 4.1. � 4.2. Non-negativity and boundedness. Theorem 4.3. All the solutions of system (2.3) which are initiating in R3+ are uniformly bounded within a region Π = {(S,E,I) ∈ R3+ : V ≤ k(b+λ)2 4λb + �,� > 0}. Proof. Here we follow an approach which is used in [27]. Define a function V (t) = S(t) + E(t) + I(t). Then, cDµV (t) =c DµS(t) +c DµE(t) +c DµI(t) Now for any positive number λ, we calculate cDµV (t) + λV (t) = bS(1 − S k ) − βSE 1 + αI + βSE 1 + αI − (c + d)E + cE −γI + λ(S + E + I) = − b k S2 + (b + λ)S + (λ−γ)I + (λ−d)E = − b k (S − k(b + λ) 2b )2 + k(b + λ)2 4b + (λ−γ)I + (λ−d)E ≤ k(b + λ)2 4b , where λ = min{d,γ}. Applying the standard comparison theorem for fractional order in Chol et al. [9], we get V (t) ≤ V (0)Gµ(−λ(t)µ) + ( k(b + λ)2 4b ) tµGµ,µ+1(−λ(t)µ) 212 C. MAJI AND D. MUKHERJEE where Gµ is the Mittag-Leffler function. So application of Lemma 5 and Corollary 6 in [9] yields V (t) ≤ k(b + λ)2 4bλ ,t →∞. Therefore, all solution of fractional order system (2.3) which are initiating in R3+ will enter the region Π = {(S,E,I) ∈ R3+ : V ≤ k(b + λ)2 4λb + �,� > 0}. (4.3) � Theorem 4.4. All solutions of system (2.3) which start in R3+ are nonnegative in nature. Proof. From first equation of system (2.3), we obtain cDµS(t) = bS(1 − S k ) − βSI 1 + µI (4.4) Again from equation (4.3), we get S + E + I ≤ k(b + λ)2 4λb = k1(say). (4.5) So from equation (4.4) and (4.5), we have, cDµS(t) ≥ bS ( 1 − k1 k ) −βS = S ( b− bk1 k −β ) = k2S, where k2 = b−β − bk1k . Now according to the standard comparison theorem for fractional order in [9] and the positivity of Mittag-Leffler function Gµ,1(t) > 0 for any µ ∈ (0, 1) [47] it follows that S(t) ≥ S0Gµ,1(qtµ) =⇒ S(t) ≥ 0. From second equation of system (2.3), we obtain cDµE(t) = βSI 1 + αI − (c + d)E ≥−(c + d)E. Therefore, E(t) ≥ E0Gµ,1(−(c + d)tµ) =⇒ E ≥ 0. Again from the third equation of system (2.3), cDµI(t) = cE −γI ≥−γI. So, I(t) ≥ I0Gµ,1(−τtµ) =⇒ I ≥ 0. Hence all solution of system (2.3) are non-negative. � 5. Equilibria of the fractional order system and their stability To evaluate the equilibrium points of system (2.3), let cDµS(t) = 0, cDµE(t) = 0, cDµI(t) = 0. Then we obtain the following equilibrium points: E0(0, 0, 0), E1(k, 0, 0) and E ∗(S∗,E∗,I∗) where E∗ = γI∗ c , S∗ = γ(c + d) cβ (1 + αI∗) and I∗ is a positive root of the equation α2bγ(c + d)I2 + {2bγ(c + d)α− bαkcβ + kβ2c}I + b{γ(c + d) −kcβ} = 0. DYNAMICAL ANALYSIS OF A FRACTIONAL-ORDER MODEL INCORPORATING FEAR 213 Now,we want to check the stability analysis of the above equilibria based on the standard linearization technique by using the Jacobian matrix. The Jacobian matrix of system (2.3) around any point (S,E,I) is given by J(S,E,I) =   b− 2bS k − βI 1+αI 0 − βS (1+αI)2 βI 1+αI −(c + d) βS (1+αI)2 0 c −γ   . Theorem 5.1. The population free equilibrium point E0 of system (2.3) is always unstable while disease free equilibrium point E1 is stable if β < γ(c+d) kc . Proof. According to the Mittag-Leffler function [31], the equilibrium point Ei of system (2.3) is locally stable if all the eigenvalues λi of J(Ei) satisfy |arg(λi)| > µπ2 . The Jacobian matrix of system (2.3) at the equilibrium point E0 is given by J(E0) =   b 0 00 −(c + d) 0 0 c −γ   . The eigenvalues corresponding to the equilibrium point E0 are: λ1 = b > 0,λ2 = −(c + d) < 0, and λ3 = −γ < 0. We observed that |arg(λ1)| = 0 < µπ2 , |arg(λ2)| = π > µπ 2 , |arg(λ3)| = π > απ2 . Hence, the equilibrium E0 is always a saddle point. Again, The Jacobian matrix of system (2.3) at the equilibrium point E1 is given by: J(E1) =   −b 0 −βk0 −(c + d) βk 0 c −γ   . The eigenvalues corresponding to the equilibrium point E1 are: λ1 = −b > 0 and other two λ2,λ3 are obtained by solving the characteristic equation λ2 + λ(c + d + γ) + γ(c + d) −βkc = 0. (5.1) The eigenvalues corresponding to the equation (5.1) are λi = (c + d + γ) ± √ (c + d + γ)2 − 4{(c + d)γ −βkc} 2 , i = 2, 3. Now we see that |arg(λ1)| = π > µπ2 ; and if β < γ(c+d) kc then λi < 0, i = 2, 3 and |arg(λ2,3)| = π > µπ2 . Hence, the equilibrium E1 is locally asymptotically stable. � To analyze the stability of the endemic equilibrium point E∗, we compute the Jacobian matrix of system (2.3) around E∗ and is given by J(E∗) =   − bS∗ k 0 − βS ∗ (1+αI∗)2 βI∗ 1+αI∗ −(c + d) βS ∗ (1+αI∗)2 0 c −γ   . The eigenvalues of J(E∗) are the roots of the following characteristic equation P(λ) = λ3 + p1λ 2 + p2λ + p3 = 0 (5.2) 214 C. MAJI AND D. MUKHERJEE where p1 = bS∗ k + c + d + γ, p2 = r(c + d) − cβS∗ (1 + αI∗)2 + bS∗ k (c + d + γ), p3 = r(c + d) − cβS∗ (1 + αI∗)2 . Let D(P) be the discriminant of the cubic polynomial P(λ), which can be written as D(P) =   1 p1 p2 p3 0 0 1 p1 p2 p3 3 2p1 p2 0 0 0 3 2p1 p2 0 0 0 3 2p1 p2   = 18p1p2p3 + (p1p2) 2 − 4p3p21 − 4p 3 2 − 27p 2 3. Then, we have the following results by [2]. Proposition 5.2. Suppose β > γ(c+d) kc . Then the equilibrium E∗ of system (2.3) is asymptotically stable if one of the following conditions are satisfied. (i) D(P) > 0,p1 > 0,p3 > 0 and p1p2 > p3; (ii) D(P) < 0,p1 ≥ 0,p2 ≥ 0,p3 > 0 and µ < 23 ; (iii) D(P) < 0,p1 > 0,p2 > 0,p1p2 = p3 and for all µ ∈ (0, 1). 6. Global Stability In this section we present global stability of endemic equilibrium point E∗. Before stating our theo- rem we define the matrix A as follows: A =   −b k − β 2α − β 2(1+αI∗)2 β 2α (c + d) −1 2 ( βS∗ α(1+αI∗) + c ) β 2(1+αI∗) −1 2 ( βS∗ α(1+αI∗) + c ) γ   . Theorem 6.1. The endemic equilibrium point E∗ of system (2.3) is globally asymptotically stable if 4α2b(c + d) > kβ2 and det A > 0. Proof. Consider the following positive definite function about E∗ V (S,E,I) = S −S∗ −S∗ ln S S∗ + 1 2 (E −E∗)2 + 1 2 (I − I∗)2. We compute the µ order derivative of V (S,E,I) along the solution of system (2.3) with the help of Lemma 3.1 in [42] and Lemma 1 in [1]. Thus we have cDµV (S,E,I) ≤ ( 1 − S ∗ S ) cDµS + (E −E∗)cDµE + (I − I∗)cDµI = (S −S∗){b(1 − S k ) − βI 1+αI } + (E −E∗){ βSI 1+αI − (c + d)E} + (I − I∗){cE −γI} DYNAMICAL ANALYSIS OF A FRACTIONAL-ORDER MODEL INCORPORATING FEAR 215 = (S −S∗){βI(S−S ∗) 1+αI + βS∗(I−I∗) (1+αI)(1+αI∗) − (c + d)(E −E∗)} + (I − I∗){c(E −E∗) −γ(I − I∗)} ≤−b k (S −S∗)2 + β 1+αI∗ |S −S∗||I − I∗| + β α |E −E∗||S −S∗| + ( βS∗ α(1+αI∗) + c ) |E −E∗||I − I∗| −(c + d)(E −E∗)2 −γ(I − I∗)2. Consequently, cDµV (S,E,I) ≤ 0 when A is positive definite. The result follows by the application of Lemma 4.6 in Huo et al. [20]. � 7. Hopf-bifurcation In this section, we discuss about the Hopf-bifurcation analysis of system (2.3). Define, a function with respect to µ by m(µ) = µπ 2 − min 1≤i≤3 |arg(λi)|. Theorem 7.1. ([24]) (Existence of Hopf bifurcation) When bifurcation parameter µ passes through the critical value µ∗ ∈ (0, 1), fractional-order system (2.3) undergoes a Hopf bifurcation at the endemic equilibrium point E∗, if the following conditions are satisfied: (i) the corresponding characteristic equation (5.2) of system (2.3) has a pair of complex conjugates λ1,2 = θ + iω (where θ > 0) and one negative real root λ3; (ii) m(µ∗) = µ ∗π 2 − min1≤i≤3 |arg(λi)| = 0; (iii) dm(µ) dµ |µ=µ∗ 6= 0. We now present the conditions for the existence of a Hopf bifurcation at the endemic equilibrium point E∗ as the order of derivative passes a critical value. Theorem 7.2. Suppose the characteristic equation (5.2) of system (2.3) has two complex conjugate eigenvalues λ1,2 = θ + iω. Then the fractional-order system (2.3) undergoes a Hopf bifurcation at the endemic equilibrium point E∗ when µ passes through the critical value µ∗ = 2 π arctan(ω θ ). Proof. From the given assumptions and m(µ) = µπ 2 − min1≤i≤3 |arg(λi)|, we get, m(µ∗) = µ∗π 2 − min 1≤i≤3 |arg(λi)| = µ∗π 2 − arctan ( ω θ ) = arctan ( ω θ ) − arctan ( ω θ ) = 0. Furthermore, dm(µ) dµ |µ=µ∗ = π 2 6= 0. Therefore, from Theorem 7.1, we conclude that system (2.3) undergoes a Hopf bifurcation at E∗ when bifurcation parameter µ crosses a critical value µ∗. � 8. Numerical Simulation In this section we present some numerical simulation to check the dynamics of the fractional order system. Although there are different type of numerical method to solve nonlinear fractional differential equations [13, 25, 6], but Adams type predictor corrector method is more appropriate and useful to solve the dynamical behaviour of the solutions of fractional differential equations. Here we have considered a 216 C. MAJI AND D. MUKHERJEE 0 50 100 150 200 t 0.6 0.8 1 1.2 1.4 1.6 1.8 2 P o p u la ti o n S E I Figure 1. Unstable behaviour of system (1) for fractional order derivative µ = 1 and parameter values α = 0.125,β = 0.75,c = d = γ = 1 3 ,b = 0.764158,k = 7.8382146. 0 10 20 30 40 50 60 70 80 90 100 t 1 1.2 1.4 1.6 1.8 2 2.2 P op ul at io n S E I Figure 2. Phase diagram for solution of system (2.3) with µ = 1 and other parameter set β = 0.75,c = d = γ = 1 3 ,b = 0.764158,k = 7.8382146. time series plot and phase diagram for α = 0.5 . set of parametric values α = 0.125,β = 0.75,c = d = γ = 1 3 ,b = 0.764158,k = 7.8382146. With the help of these parameters from Figure 1 we observe that when µ = 1 the solution of system (2.3) is unstable in nature. Next, we plot the solution of system (2.3) in Figure 2 by choosing the same set of parameters except for α = 0.5. Here we observed that the system is stable in nature. So for integer order system fear effect α has an interesting role. If we increase the value of α, the system becomes stable that means fear effect can stabilizes the system. Influence of fear effect on different population is given in Figure 3. Again the bifurcation diagram of system (2.3) with respect to the parameter α has been drawn in Figure 4. DYNAMICAL ANALYSIS OF A FRACTIONAL-ORDER MODEL INCORPORATING FEAR 217 Figure 3. Influence of fear effect on each population Now to see the effect of fractional-order µ on each population we plotted a diagram in Figure 5. For the above set of parameters, from Figure 6 we have seen that for integer-order system µ = 1, system (2.3) is unstable when α = 0.125. Again, for fractional-order derivative µ = 0.99 and µ = 0.98 the system shows unstable behaviour but the system changes its stability when µ = 0.96 and µ = 0.92 and it becomes stable (Figure 6(c),6(d)). Thus, from the above figures we conclude that fractional order derivative may change the system dynamical behaviour from unstable to stable. Hence, fractional-order derivative has an important role on the system stability of our considered system and it may improve system stability. A bifurcation diagram with respect to the fractional-order µ is given in Figure 7. 218 C. MAJI AND D. MUKHERJEE Figure 4. Bifurcation diagram of system (2.3) with respect to α when µ = 1 and other parameters are β = 0.75,c = d = γ = 1 3 ,b = 0.764158,k = 7.8382146. 9. Discussion Recently, Harper et al.[19] investigated psychological predictors of behavior change and fear response to the COVID-19 pandemic 2020. In [3], the authors discussed that the fear of COVID-19 Scale is a seven-item uni-dimensional scale with robust psychometric properties. They also concluded that this approach is convincing and appropriate in examining fear of COVID-19 among the peoples. Wang et. al. [46] investigated a time-dependent mathematical model of COVID-19 to focus on the effects of medical resources on transmission of COVID-19 but fear effect and fractional-order is not addressed in their work. At present there is no proper estimate about how long this disease persist, the number of DYNAMICAL ANALYSIS OF A FRACTIONAL-ORDER MODEL INCORPORATING FEAR 219 0 10 20 30 40 50 time (t) 0.4 0.6 0.8 1 1.2 1.4 1.6 S u s c e p ti b le =0.9 =0.8 =0.7 =0.6 (a) 0 10 20 30 40 50 time(t) 0.4 0.6 0.8 1 1.2 1.4 E x p o s e d =0.9 =0.8 =0.7 =0.6 (b) 0 10 20 30 40 50 time(t) 0.6 0.8 1 1.2 1.4 1.6 In fe c te d =0.9 =0.8 =0.7 =0.6 (c) Figure 5. Time series plot of each population for parameter values α = 0.125,β = 0.75,c = d = γ = 1 3 ,b = 0.764158,k = 7.8382146 with different values of µ individuals worldwide who will be infected and how long human lives will be affected. Extensive re- search is going on to control the spread of the disease in different way . So above observation motivate us to find out a way to combat the disease. The significance of communicable disease behavior induces scientists to develop a mathematical model that can examine the spread procedure, rule, and direction etc. Its purpose is that, according to the aspect of infectious disease, the model is suitably formed, the relevant parameters are selected and reasonable variables are chosen. We have investigated the three component epidemiological fractional system which considers the fear effect in the disease transmission rate of coronavirus disease. The 220 C. MAJI AND D. MUKHERJEE Figure 6. Phase diagram for solution of system (1) when α = 0.125 and other pa- rameters are β = 0.75,c = d = γ = 1 3 ,b = 0.764158,k = 7.8382146. with (a) µ = 0.99 (b) µ = 0.98 (c) µ = 0.96 and (d) µ = 0.92. dynamical behavior of the given system is studied. The classical first order time derivative is mod- elled with the Caputo fractional derivative of order µ ∈ (0, 1]. Mathematical analysis of existence and uniqueness of solutions are shown. Basic stability properties of disease free and endemic equilibrium points are examined. We observed that as long as disease transmission rate remains below a certain threshold value (β < γ(c+d) kc ) the disease free equilibrium point is locally asymptotically stable. If the disease transmission rate increases then endemic equilibrium point appears. Stability property of en- demic equilibrium point is described in Theorem 5. For disease eradication the conditions of Theorem 5 should be avoided. From our analytical and numerical studies indicate that a Hopf bifurcation due to variation of the fractional order µ ∈ (0, 1]. From Figure 1, we observed that fear has an important role on the dynamics of our system. When the level of fear is very low the system exhibits unstable be- haviour while increasing the value of fear stabilizes the system. We have plotted a bifurcation diagram choosing α as a bifurcation parameter in Figure 4. From this figure we have seen that when the value of α in the range 0.05 ≤ α < 0.127, the system is unstable but in the range 0.125 ≤ α ≤ 0.2 the system DYNAMICAL ANALYSIS OF A FRACTIONAL-ORDER MODEL INCORPORATING FEAR 221 Figure 7. Bifurcation diagram of system (2.3) with respect to µ when α = 0.1251 and other parameters are β = 0.75,c = d = γ = 1 3 ,b = 0.764158,k = 7.8382146. is stable. Therefore we conclude that the level of fear may decrease disease transmission rate and the epidemic may be under control. It is also noted that the fractional order model is more stable than the integer order model. Acknowledgement: The authors are grateful to the Editor and anonymous reviewers for their valuable comments and suggestions for improving the paper. 222 C. MAJI AND D. MUKHERJEE References [1] N. Aguila-Camacho, M. A. Durate-Mermol and J. a. Galleges, Lyapunov functions for fractional order systems, Commun. Nonlinear. Sci. Numer. Simul. 19(2014), 2951-2957. [2] E. Ahmed, A. El-Sayed and H. EL-Saka, Equilibrium points, stability and numerical solutions of fractional-order predator-prey and rabies models, J. Math. Anal. Appl. 325 (2007), 542-553. [3] D. K. Ahorsu, C-Y. Lin, V. Imani, M. Saffari, M. D. Griffiths and A. H. Pakpour, The fear of COVID-19 Scale: Development and initial validation, Int. J. of Mental Health and Addtion. 2020. Doi. org/10.1007/s11469-020-00270-8. [4] A. Apisarnthanarak et al., Impact of anxiety and fear for COVID-19 toward infection control practices among Thai healthcare workers, Infection Control & Hospital Epidemiology, (2020) 1-2, doi:10.1017/ice.2020.280. [5] K. Assaleh and W. M. Ahmad, Modelling of speech signals using fractional calculus, In: 9th International Symposium on Signal Processing and its Applications (2007) (ISSPA 2007). [6] S. Bhalekar and V. Daftardar-Gejji, A predictor-corrector scheme for solving nonlinear delay differential equations of fractional order, J. Fract. Calc. Appl. 1(5) (2011), 1-9. [7] S. Cauchemez, C. Fraser and M. D. Van Kerkhove, et al., Middle East respiratory syndrome coronavirus: quantifi- cation of the extent of the epidemic, surveillance biases, and transmissibility, Lancet. Infect. Dis 14 (2014), 50-56. [8] R. Chinnathambi and F. A. Rihan, Stability of fractional-order prey-predator system with time-delay and Monod- Haldane functional response, Nonlinear Dynamics 92 (2018), 1637-1648. [9] S. K. Choi, B. Kang and N. Koo, Stability for Caputo fractional differential systems, Abstr. Appl. Anal. 2014 (2014), 1-6. [10] K. S. Cole, Electric conductance of biological systems, In: Cold. Spring Harbor Symposium on Quantitative Biology, (1993) 107-116. [11] L. Debnath, Recent applications of fractional calculus to science and engineering, Int. J. Math. Sci. 54 (2003), 3413-3442. [12] E. Demirci, A. Unal and N. Ozalp, A fractional order SEIR model with density dependent death rate, Hacettepe J. Math. Stat. 40 (2011), 287-295. [13] K. Diethem, N. J. Ford and A. D. Freed, A predictor corrector approach for the numerical solution of fractional differential equations, Nonlinear. Dynamics 29 (2002), 23-22. [14] L. Dong and S. Hu, J. Gao, Discovery drugs to treat coronavirus disease 2019 (COVID-19), Drug Discoveries and therapeutics 14 (2020), 58-60. [15] C. A. Donnely, A. C. GhaniG and M Leung et al., Epidemiological determinants of spread of causal agent of severe acute respiratory syndrome in Hong Kong, Lancet 361(2003), 1761-1766. [16] S. El-Sayed, Nonlinear functional differential equations of arbitrary orders, Nonlinear. Anal. Theory. Methods. Appl. 33(2), (1998) 181-186. [17] A. El-Sayed, A. El-Mesiry and H. EL-Saka, On the fractional order logistic equation, Appl. Math. Lwtt. 20(7) (2007), 817-823. [18] Y. Fedri, Some applications of fractional order calculus to design digital filters for biomedical signal processing, J. Mech. Med. Biol. 12(2) (2012), 13. [19] C. A. Harper, L. Satchell, D. Fido and R. Latzman, Functional fear predicts public health compliance in the COVID- 19 pandemic-2020, Int J Ment Health Addiction (2020): https://doi.org/10.1007/s11469-020-00281-5 [20] J. Huo, H. Zhao and L. Zhu, The effect of vaccines on backward bifurcation in a fractional order HIV model, Nonlinear Anal. 26 (2015), 289-305. [21] H. Kheiri and M. Jafari, Stability analysis of a fractional order model for the HIV/AIDS epidemic in a patchy environment, J. Comput. Appl. Math. 346 (2019), 323-339. [22] A. A. Kilbas, H. M. Srivastava and J. J. Trujilo, Theory and Applications of Fractional Differential Equations, Elsevier, 2006. [23] T. Kobayshi, S. M. Jung, N. M. Linton, R. Kinoshita, K. Hayashi and T. Miyama et al., Communicating the risk of death from novel coronavirus disease (CPVID-19), J. Clin. Med. 9 (2020) 580. [24] X. Li and R. Wu, Hopf bifurcation analysis of a new commensurate fractional-order hyperchaotic system, Nonlinear Dynamics. 78(2014), 279-288. [25] C. Li and C. Tao, On the fractional Adams method, Comput. Math. Appl. 58(8) (2009), 1573-1588. [26] Y. Li, Y. Chen and I. Podlubny, Stability of fractional-order non linear dynamic systems: Lyapunov direct method and generalized Mittag-Leffler stability, Comput. Math. Appl. 59 (2010), 1810-1821. [27] H. Lili, Z Long and H. Cheng, et al., Dynamical analysis of a fractional-order predator-prey model incorporating a prey refuge, J. Appl. Math. Comput. 54 (2016), 435449. [28] C. Y. Lin, Social reaction toward the 2019 novel coronavirus (COVID-19), Social Health and Behavior 3 (2020) 1-2. DYNAMICAL ANALYSIS OF A FRACTIONAL-ORDER MODEL INCORPORATING FEAR 223 [29] J. A. T. Machado, Entropy analysis of integer and fractional dynamical systems, Nonlinear Dynamics 62(1-2) (2010), 371-378. [30] J. A. T. Machado and A. M. S. F. Galhano, Fractional order inductive phenomenona based on the skin effect, Nonlinear Dynamics 68(1-2) (2012), 107-115. [31] D. Matignon, Stability results in fractional differential equation with applications to control processing. In: Pro- ceedings of the Multiconference on Computational Engineering in Systems and Application IMICS.IEEE-SMC. Life, France. 1996; Vol.2: 963-968. [32] G. Mertens et al., Fear of the coronavirus COVID-19): Predictors in an online study conducted in March 2020, Journal of Anxiety Disorders, (2020) , doi.org/10.1016/j.janxdis.2020.102258. [33] Y. S. Mishura, Stochastic calculus for Fractional Brownian Motion and Related Process. vol. 1929 for Lectures Notes in Mathematics, Springer, Berlin, Germany, 2008. [34] K. Muniz-Rodriguez, G. Chowell, C. H. Cheung, D. Jia, P.-Y. Lai, Y. Lee, M. Liu, S. K. Ofori, K. M. Roosa, L. Simonsen, I. C.H. Fung, Epidemic doubling time of the 2019 novel coronavirus outbreak by province in mainland china, medRxiv, 2020. [35] H. Nishiura, N. M. Linton, and A. R. Akhmetzhanov, Serial interval of novel coronavirus (2019-ncov) infections, medRxiv, 2020. [36] I. Petras, Fractional-order nonlinear systems: Modeling, analysis and simulation, Berlin, Springer-Verlag. (2011) doi:10.1007/978-3-642-18101-6. [37] I. Podlubny, Fractional Differential Equations, vol.198 of Mathematics in Science and Engineering, Academic Press, San Diego, Galif, USA, 1999. [38] F. A. Rihan and D. H. Abdelrahman, Delay differential model for tumor-immune dynamics with HIV infection of CD+ T-cells, Int. J. Comput. Math. 90(3) (2013) 594-614. [39] S. Sanche, Y. T. Lin, C. Xu, E. Romero-Severson, N. Hengartner, and R. Ke. The novel coronavirus, 2019-ncov, is highly contagious and more infectious than initially estimated, medRxiv, 2020. [40] C. J. Silva and D. F. M. Torres, Stability of a fractional HIV/AIDS model, Math. Comput. Simulat. 164 (2019), 180-190. [41] J. M. Shultz, B. M. Althouse and F. Baingana et al., Fear factor: The unseen perils of the Ebola outbreak, Bull. At. Sci. 72(5) (2016), 304-310. doi:10.1080/00963402.2016.1216515. [42] C. Vargas-De-Leon, Volterra- type Lyapunov function for fractional order epidemic systems, Commun. Nonlinear. Sci. Numer. Simul. 24 (2015), 75-85. [43] V. Volpert, M. Banerjee and S. Petrovskii, On a quarantine model of coronavirus infection and data analysis, Math. Model. Nat. Phenom. 15 (2020): 24, DOI: 10.1051/mmnp/2020006 [44] A. O. Victor, Mathematical Predictions for COVID-19 as a Global Pandemic, 2020; medRxiv preprint. [45] Z. Wang, K. Xie, J. Lu and Y. Li, Stability and bifurcation of a delayed generalized fractional-order prey-predator model with interspecific competition, Appl. Math. Comp. 347(C) (2019), 360-369. [46] L. Wang, J. Wang, H. Zhao, Y. Shi, K. Wang, P. Wu and L. Shi, Modelling and assessing the effects of medi- cal resources on transmission of novel coronavirus (COVID-19) in Wuhan, China, Mathematical Biosciences and Engineering 17(4) (2020), 2936-2949. [47] Z. Wei, Q. Li and J. Che, Initial value problems for fractional differential equations involving Riemann-Liouville sequential derivative, J. Math. Appl. 367 (2010), 260-72. [48] P. Wu, X. Hao and E. H. Y Lau et al., Real-time tentative assessment of the epidemiological characteristics of novel coronavirus infections in Wuhan, China, as at 22 January 2020, Euro. Surveill 25 (2020) 2000044. [49] Y. Yang, Q. Lu, M. Liu, Y. Wang, A. Zhang, N. Jalali, N. Dean, I. Longini, M. Elizabeth Halloran, B. Xu, X. Zhang, L. Wang, W. Liu, and L. Fang. Epidemiological and clinical features of the 2019 novel coronavirus outbreak in china, medRxiv, 2020. Department of Mathematics, Vivekananda College, Thakurpukur, Kolkata-700063, India Current address: same E-mail address: chandanmaji.ju@gmail.com Department of Mathematics, Vivekananda College, Thakurpukur, Kolkata-700063, India E-mail address: mukherjee1961@gmail.com