Loading [MathJax]/jax/output/SVG/jax.js
Research article Special Issues

Fractional order SIR model with generalized incidence rate

  • The purpose of this work is to present the dynamics of a fractional SIR model with generalized incidence rate using two differential derivatives, that are the Caputo and the AtanganaBaleanu. Firstly, we formulate the proposed model in Caputo sense and carried out the basic mathematical analysis such as positivity, basic reproduction number, local and global dynamics of the disease free and endemic case. The disease free equilibrium is locally and globally asymptotically stable when the basic reproduction number is less than 1. We also show that the SIR fractional model is locally and globally asymptotically stable at endemic state when the basic reproduction number greater than unity. Further, to analyze the dynamics of Caputo model we provide some numerical results by considering various incidence rates and with different fractional order parameters. Then, we reformulate the same model using the Atangana-Baleanu operator having nonsingular and non-local kernel. We prove the existence and uniqueness of the Atangana-Baleanu SIR model using PicardLindelof approach. Furthermore, we present an iterative scheme of the model and provide graphical results for various values of the fractional order α. From the graphical interpretation we conclude that the Atangana-Baleanu derivative is more prominent and provides biologically more feasible results than Caputo operator.

    Citation: Muhammad Altaf Khan, Muhammad Ismail, Saif Ullah, Muhammad Farhan. Fractional order SIR model with generalized incidence rate[J]. AIMS Mathematics, 2020, 5(3): 1856-1880. doi: 10.3934/math.2020124

    Related Papers:

    [1] Ahmad A. Abubaker, Khaled Matarneh, Mohammad Faisal Khan, Suha B. Al-Shaikh, Mustafa Kamal . Study of quantum calculus for a new subclass of q-starlike bi-univalent functions connected with vertical strip domain. AIMS Mathematics, 2024, 9(5): 11789-11804. doi: 10.3934/math.2024577
    [2] Ekram E. Ali, Georgia Irina Oros, Rabha M. El-Ashwah, Abeer M. Albalahi . Applications of fuzzy differential subordination theory on analytic p -valent functions connected with q-calculus operator. AIMS Mathematics, 2024, 9(8): 21239-21254. doi: 10.3934/math.20241031
    [3] Ebrahim Amini, Mojtaba Fardi, Shrideh Al-Omari, Rania Saadeh . Certain differential subordination results for univalent functions associated with q-Salagean operators. AIMS Mathematics, 2023, 8(7): 15892-15906. doi: 10.3934/math.2023811
    [4] Shujaat Ali Shah, Ekram Elsayed Ali, Adriana Cătaș, Abeer M. Albalahi . On fuzzy differential subordination associated with q-difference operator. AIMS Mathematics, 2023, 8(3): 6642-6650. doi: 10.3934/math.2023336
    [5] Muhammad Sabil Ur Rehman, Qazi Zahoor Ahmad, H. M. Srivastava, Nazar Khan, Maslina Darus, Bilal Khan . Applications of higher-order q-derivatives to the subclass of q-starlike functions associated with the Janowski functions. AIMS Mathematics, 2021, 6(2): 1110-1125. doi: 10.3934/math.2021067
    [6] Syed Ghoos Ali Shah, Shahbaz Khan, Saqib Hussain, Maslina Darus . q-Noor integral operator associated with starlike functions and q-conic domains. AIMS Mathematics, 2022, 7(6): 10842-10859. doi: 10.3934/math.2022606
    [7] Aisha M. Alqahtani, Rashid Murtaza, Saba Akmal, Adnan, Ilyas Khan . Generalized q-convex functions characterized by q-calculus. AIMS Mathematics, 2023, 8(4): 9385-9399. doi: 10.3934/math.2023472
    [8] Georgia Irina Oros, Gheorghe Oros, Daniela Andrada Bardac-Vlada . Certain geometric properties of the fractional integral of the Bessel function of the first kind. AIMS Mathematics, 2024, 9(3): 7095-7110. doi: 10.3934/math.2024346
    [9] Ekram E. Ali, Miguel Vivas-Cortez, Rabha M. El-Ashwah . New results about fuzzy γ-convex functions connected with the q-analogue multiplier-Noor integral operator. AIMS Mathematics, 2024, 9(3): 5451-5465. doi: 10.3934/math.2024263
    [10] Ekram E. Ali, Nicoleta Breaz, Rabha M. El-Ashwah . Subordinations and superordinations studies using q-difference operator. AIMS Mathematics, 2024, 9(7): 18143-18162. doi: 10.3934/math.2024886
  • The purpose of this work is to present the dynamics of a fractional SIR model with generalized incidence rate using two differential derivatives, that are the Caputo and the AtanganaBaleanu. Firstly, we formulate the proposed model in Caputo sense and carried out the basic mathematical analysis such as positivity, basic reproduction number, local and global dynamics of the disease free and endemic case. The disease free equilibrium is locally and globally asymptotically stable when the basic reproduction number is less than 1. We also show that the SIR fractional model is locally and globally asymptotically stable at endemic state when the basic reproduction number greater than unity. Further, to analyze the dynamics of Caputo model we provide some numerical results by considering various incidence rates and with different fractional order parameters. Then, we reformulate the same model using the Atangana-Baleanu operator having nonsingular and non-local kernel. We prove the existence and uniqueness of the Atangana-Baleanu SIR model using PicardLindelof approach. Furthermore, we present an iterative scheme of the model and provide graphical results for various values of the fractional order α. From the graphical interpretation we conclude that the Atangana-Baleanu derivative is more prominent and provides biologically more feasible results than Caputo operator.


    Shortly after the first local transmissions of the coronavirus disease 2019 (COVID-19) in the Philippines were confirmed, the government imposed social distancing policies, termed community quarantines, which were largely implemented by the police and military [1,2]. By March 30, 2020, the country only had six laboratories that accommodated up to 1000 tests daily [3]. Contact tracing began slowly due to an insufficient number of contact tracers [4]. Testing capacity was increased to about 35000 tests daily by the end of September 2020 [5]. From April to September 2020, the weekly positivity rate ranged between 4.5 and 28%, higher than the 5% threshold set by the World Health Organization during this time [6,7].

    The Philippines' capital region, Metro Manila, comprised around 12% of the country's population in 2020 [8]. Metro Manila was placed under Enhanced Community Quarantine (ECQ) on March 16, 2020 [9]. Under ECQ, movement was restricted to essential goods and services. Public transportation and mass gatherings were suspended [2]. People were encouraged to work from home and businesses were advised to do transactions online [1]. After two months, ECQ was replaced by the Modified Enhanced Community Quarantine (MECQ), a transition phase before easing further to General Community Quarantine (GCQ). On June 1, 2020, Metro Manila was placed under GCQ, wherein public transportation and other establishments, except those for leisure, were allowed to operate [2]. A surge in the number of cases occurred from July to August 2020 and consequently, Metro Manila was again placed under MECQ. During this time, the utilization of ICU, isolation, and ward beds in Metro Manila reached 77%, 74%, and 84%, respectively, placing most facilities on critical or high-risk, and prompting 80 medical societies, representing 80000 doctors and a million nurses, to demand a 'timeout' [10,11]. On August 19, 2020, Metro Manila returned to GCQ [12]. By the end of September 2020, 53% of the 309303 total confirmed cases in the Philippines belonged to Metro Manila [5].

    Because of the lack of vaccines and limited antiviral therapies during the early phase of the COVID-19 pandemic, NPIs such as wearing masks, school and workplace closures, and travel restrictions were crucial disease control measures. In the Philippines, compliance with policies was not only prompted by public health campaigns, but also driven by uncertainty and anxiety about the disease, and fear of getting reprimanded by the authorities [13,14,15]. Some of those who got infected suffered stigma and were blamed for not following the protocols [16,17]. A study among low-income households in the Philippines done in the early phase of the pandemic reported that 66% of respondents who might experience symptoms considered staying at home instead of seeking medical attention [13].

    Mathematical modeling has been extensively used to understand the dynamics of COVID-19 transmission [18,19,20,21,22,23,24]. Non-pharmaceutical interventions, behavior change, and underreporting of cases have been incorporated into mathematical models of COVID-19 [25,26,27,28,29,30,31,32]. In this study, we aim to investigate the effects of reporting and behavior on the timing and magnitude of the peak of COVID-19 infections in Metro Manila from March to September 2020. We extend the SEIQR model by Kim et al., which includes a compartment for behavior-changed susceptible individuals [28], by adding an unreported compartment to account for individuals who were undetected due to inadequate testing and tracing, or unwillingness to be detected. We provide a detailed mathematical analysis of the model showing the existence, nonnegativity, boundedness of solutions, computation of the threshold parameter, and two forms of disease-free equilibria (DFEs). We show the global stability of one of the DFEs by utilizing a suitable Lyapunov function and analyzing the asymptotic behavior of the states. Furthermore, a bi-objective optimization problem is formulated that allows the easing from ECQ to GCQ at an earlier time and minimizes the number of additional beds necessary to ensure sufficient capacity in healthcare facilities. The multiple optimal solutions (Pareto optimal solutions) of the bi-objective problem offer trade-off solutions in which the user, a policymaker or healthcare authority, can use in decision-making.

    The number of cumulative confirmed cases from March 8 to September 30, 2020, was obtained from the Philippine Department of Health (DOH) data drop [6]. These data were used to estimate the rates of transmission, behavior change, and reporting. The data on COVID-19 bed capacity and occupancy rate in Metro Manila were gathered from the number of occupied and available isolation, ward, and ICU beds from the weekly DOH bulletin from June 20 to September 26, 2020 [33]. The total population of Metro Manila was set to 13484462, based on the 2020 census data of the Philippine Statistics Authority [8].

    The model we present is an extension of the model in [28] wherein an unreported compartment is added to represent the undetected or unreported COVID-19 cases in the early phase of the pandemic in the Philippines. The constant total population N is divided into eight mutually exclusive compartments: susceptible (S), behavior-changed susceptible (SF), exposed (E), reported infectious (I), unreported infectious (Iu), isolated (Q), recovered (R), and deaths (D). A schematic diagram of the model is shown in Figure 1.

    Figure 1.  COVID-19 transmission model that incorporates behavior change and unreported cases. Susceptible (S) may change behavior (SF) and vice versa at rates βF or μ. These classes can be exposed (E) to the virus and become infectious (I,Iu) in 1/κ days on average. Transmission rate β is reduced by a factor δ for a behavior-changed SF. Reporting ratio is denoted by ρ. Confirmed cases are isolated (Q) in 1/α days and recover (R) 1/γ days on average or die (D). The average fatality rate is denoted by f. Those in Iu recover 1/η days on average.

    Assuming a local, prevalence-based spread of fear of the disease and following the study of Perra et al. [34], the transition rate of a susceptible to a behavior-changed susceptible is given by βFQ/(NQD). This means that a susceptible individual is more likely to change behavior as the number of confirmed cases among one's contacts increases. The movement back to S is assumed to be influenced by the number of recoveries and susceptible individuals without behavior change [34]. As the recoveries and susceptible individuals increase among the contacts of a behavior-changed susceptible, the more likely the individual to exit the SF class and resume regular social behavior. The parameter μ represents the rate of the easing of behavior and its value is assumed to be 1/14 [28].

    Susceptible individuals (S and SF) move to the exposed class upon contact with infectious individuals (I and Iu) at a rate β. The transmission rate for the behavior-changed susceptible class is assumed to be reduced by a factor δ. The reporting ratio ρ partitions the exposed class to reported I and unreported Iu classes. Assuming that an individual becomes infectious 2 days before symptom onset [35], the mean incubation period of the original virus strain is 6 days [36], and the mean duration between symptom onset and the first medical consultation in the Philippines during this time was 6.75 days [37], we set the mean latent period (1/κ) to 4 days and the mean infectious period of the reported cases (1/α) to 8.75 days. From isolation, individuals recover 1/γ days on average or die. The average fatality rate, which is the ratio of daily deaths to daily active cases and denoted by f, is set to 1.9% [38]. Those in the unreported class are assumed to have less severe symptoms and move to the recovered class 1/η days on average. The following system of differential equations describe the model used in the study:

    dSdt=βSI+Iu˜N+μSFS+R˜NβFSQ˜N,dSFdt=δβSFI+Iu˜NμSFS+R˜N+βFSQ˜N,dEdt=δβSFI+Iu˜N+βSI+Iu˜NκE,dIdt=ρκEαI,dIudt=(1ρ)κEηIu,dQdt=αIγQ,dRdt=(1f)γQ+ηIu,dDdt=fγQ,˜N=NQD,N=S+SF+E+I+Iu+Q+R+D, (2.1)

    where S(0)0, SF(0)0, E(0)0, I(0)0, Iu(0)0, Q(0)0, R(0)0, and D(0)0. We introduce the term ˜N because those in Q and D are assumed to be isolated and have no contact with the rest of the population. All the parameters are assumed to be positive and ρ,f(0,1]. In the numerical simulations, we set the initial population of the infectious I(0), exposed E(0), and unreported Iu(0) classes equal to the number of cases 1/α, 1/α+1/κ, and 10×I0 days from March 8, 2020, respectively. The initial number of isolated individuals Q(0) was the number of cases at the start of the estimation period. The initial susceptible population S(0) was computed by getting the difference between the total population and E(0), Iu(0), I(0), and Q(0). The rest of the state variables were initially set to zero. The model parameters and initial values of the state variables are shown in Table 1.

    Table 1.  List of model parameters, their values, and references. The subscript denotes the estimation period; 1 if from March 8 to May 31, 2020 or 2 if from June 1 to September 30, 2020.
    Symbol Description (unit) Value Ref.
    β1,β2 Transmission rate of COVID-19 in 0.199, 0.361 Estimated
    Periods 1 and 2 (1/day)
    βF,1,βF,2 Transmission rate of awareness or fear of 471.057, 68.783 Estimated
    the disease in Periods 1 and 2 (1/day)
    μ Behavior change ease rate (1/day) 1/14 Assumed
    δ Transmission reduction factor for 0.202 Estimated
    behavior-changed individuals
    1/κ Mean latent period (day) 4 [35,36]
    1/α Mean infectious period of reported cases (day) 8.75 [35,37]
    1/γ Mean recovery period (day) 16 [37]
    f Mean fatality rate 1.9% [38]
    ρ1,ρ2 Reporting ratio in Periods 1 and 2 0.289, 0.866 Estimated
    1/η Mean infectious period of unreported cases (day) 10 [39]
    S(0) Initial susceptible population 13483232 Calculated
    SF(0) Initial behavior-changed susceptible population 0 Assumed
    E(0) Initial exposed population 139 Assumed
    I(0) Initial reported infectious population 99 Calculated
    Iu(0) Initial unreported infectious population 990 Assumed
    Q(0) Initial isolated population 2 [6]
    R(0) Initial recovered population 0 Assumed
    D(0) Initial deaths population 0 Assumed

     | Show Table
    DownLoad: CSV

    The values of the transmission rate, reporting ratio, and reduction factor were estimated from the cumulative cases data in Metro Manila from March 8 to September 30, 2020. We divide the estimation period into two: Period 1 is from March 8 to May 31, and period 2 is from June 1 to September 30, 2020. Metro Manila was mostly under ECQ during period 1, while it was mostly under GCQ during period 2. It was during period 2 that the first major epidemic wave in the Philippines occurred. Since the intensity of NPIs and the behavior of the population during ECQ and GCQ vary, the values for the transmission rates (β and βF) and reporting ratio (ρ) in periods 1 and 2 are assumed to be different. The reduction in transmission (δ) for the behavior-changed susceptible class is assumed to have the same value in the two periods. We denote the transmission rates and reporting ratio for periods 1 or 2 by the subscripts 1 or 2, respectively.

    Denote by p=[β1,β2,βF,1,βF,2,δ,ρ1,ρ2] the vector of parameters to be estimated and set the domain as Γ=[0,1]×[0,1]×[0,103]×[0,103]×[0,1]×[0,1]×[0,1]. Define the objective functional J:ΓR7R as

    J(p)=i(αI(ti;p)˜C(ti))2, (2.2)

    where ˜C(ti) is the cumulative reported cases on day ti based on the data and I(ti;p) is the solution I of (1) at t=ti using the parameter vector p. To estimate the unknown parameters, we calculate the minimizer p of the objective function (2.2) using a trust-region method based on the interior-reflective Newton method [40,41]. Each iteration of this approach solves a large linear system using preconditioned conjugate gradient method [42]. We utilize the Matlab built-in function lsqcurvefit to implement this numerical optimization technique. The program requires the bounds, which we set based on the domain Γ, and initial guess, which we set to pinit=[0.3,0.3,285,285,0.2,0.25,0.25].

    Sensitivity analysis is a numerical technique that is widely used in identifying and ranking critical parameters for a given model output [43]. A parameter is said to be influential to an output if small perturbations of its value lead to significant changes in the model output. In this work, we use the Partial Rank Correlation Coefficient (PRCC) method paired with the Latin Hypercube Sampling (LHS) technique. LHS-PRCC is a global sensitivity analysis technique and is widely used in infectious disease modeling [44,45,46,47,48]. We investigate the sensitivity of the ten parameters χ:=[β,ρ,α,δ,βF,η,μ,κ,γ,f][0.01,1]×[0.01,1]×[0.01,0.5]×[0.01,1]×[0.01,103]×[1/40,1/7]×[0.05,1]×[0.1,0.5]×[1/40,1/7]×[0.01,0.1]. To consider every infection, we use the cumulative number of infected individuals,

    F(t;χ)=t0κE(σ;χ)dσ, (2.3)

    as the model output. For each parameter χj, where j{1,2,,10}, we sample 10000 values {Xij}10000i=1, all following a uniform distribution, using LHS. Hence, XR10000×10 with each row representing a sampled combination of χ. For brevity, we denote Xi as the ith row of X and Xj as the jth column of X. Denote by YR10000 the output vector whose ith component is calculated by evaluating F(t;Xi). The correlation coefficient (CC) is a metric to quantify the linear association between input and output. If the data is ranked-transformed first before calculating the CC, the result is a Spearman or rank correlation coefficient (RCC). Ranking is done by sorting both Xj and Y in descending order, and the integer rank values XRj and YR, respectively, are stored. Partial correlation coefficient (PCC) is used to characterize the linear relationship between Xj and Y after the linear effects on Y of the other inputs are discounted [43]. The partial rank correlation coefficient (PRCC) between Xj and Y is the PCC between XRj and YR. The PRCC is computed numerically using the Matlab built-in function partialcorr, with the 'type' set to 'Spearman'. Note that the output Y is time-dependent as seen in (2.3). In our simulations, the PRCC values are calculated at five time points: April 19, May 31, August 2, October 4, and November 1, 2020.

    Parameter bootstrapping is a statistical technique to quantify uncertainty and construct confidence intervals of estimated parameters. In this study, we utilize the algorithm introduced in [49], where large samples of synthetic data sets using the estimated model parameters are generated assuming a certain probability distribution structure. In our simulations, parameters are re-estimated from 10000 (number of realizations) synthetic data sets each generated by assuming a Poisson error structure. The mean, standard deviation, and 95% confidence intervals of the re-estimated parameters are determined. This technique is summarized in Algorithm 1.

    Algorithm 1 Uncertainty analysis.
    1: Input: The estimated parameter vector p, the number of realizations M
    2: Calculate αI(ti;p)) at each time point ti by solving (1).
    3: for j=1:M do
    4:   Generate the noisy data ˆCj(ti) by adding Poisson noise to αI(ti;p)).
    5:   Using the Matlab built-in function lsqcurvefit, calculate the minimizer ˆpj of
              ˆJ(p)=i(αI(ti;p)ˆCj(ti))2.
    6: end for
    7: Perform statistical analysis on the re-estimated parameter vectors {ˆpj}Mj=1.
    8: Output: Histograms to display the empirical distributions of the estimates and the corresponding mean, standard deviation, and 95% confidence interval.

     | Show Table
    DownLoad: CSV

    Using the bed occupancy data from the DOH [33], we calculated that an average of 16% of the active cases Q(t) occupied COVID-19 beds from June to September 2020. Fitting the weekly data on available beds, a linear function representing 75% of the bed capacity was obtained. The DOH categorizes a facility as high risk if the bed occupancy is 70% to 85%, and critical if bed occupancy is greater than 85%. From July 18 to August 8, 2020, most facilities in Metro Manila were on critical or high-risk, with a combined COVID-19 bed occupancy (isolation, ward, and ICU beds) exceeding 75% of the capacity. Here, we propose an optimization approach to determine the number of additional beds per day so that the number of cases requiring beds does not exceed 75% of the total bed capacity and if it is possible to transition from ECQ to GCQ earlier than June 1, 2020.

    We denote by QH(t) the number of active cases requiring beds and H(t) the linear, time-dependent, data-fitted bed capacity. We consider two objectives:

    (1) Minimize the number of additional hospital beds needed per day (ω) so that the 75% capacity is not reached;

    (2) Determine an earlier timing (τ) of easing from ECQ to GCQ.

    The problem can be formulated as a bi-objective constrained optimization problem expressed as follows:

    min[ωτ] (2.4)

    such that

    QH(t)H(t;ω,τ):=0.75[ω(tτ)+H0] for all t, (2.5)

    where QH(t)=0.16Q(t), H0 is the baseline number of beds at time τ based on the data, and Q(t) is solved from (2.1) by setting τ as the day when GCQ started. Since this is a bi-objective optimization problem, the solution is not unique but a pareto optimal set. We solve (2.4) and (2.5) using Genetic Algorithm, which has found a growing number of applications in various fields of science and engineering [50,51,52]. In particular, we implement the Matlab built-in function gamultiobj, which is based on a variant of Non-dominated Sorting Genetic Algorithm Ⅱ (NSGA-Ⅱ) [53,54].

    In the following theorems, we establish the existence, nonnegativity, and boundedness of the solutions to model (2.1) for t0. Moreover, we derive the disease-free equilibria and a threshold parameter for the given model.

    Theorem 1 (Existence and uniqueness of solutions). There exists a unique solution to (2.1) for t0.

    Proof. Model (2.1) can be written as

    ˙x=F(t,x(t)), x(0)=x0, (3.1)

    where the column vector x(t)=(S(t),SF(t),E(t),I(t),Iu(t),Q(t),R(t),D(t))TR8 defines a mapping from [0,+] to R8. Moreover, let F(t,x(t))=(F1(t,x(t)),,F8(t,x(t)))T be a vector-valued function from R8 to R8 equal to the right-hand side of the equations in (2.1). It can be verified that F is continuous and each of the first-order partial derivatives of F with respect to x is continuous in its domain. By the existence and uniqueness theorem in [55], system (2.1) has a unique solution for t0.

    Theorem 2 (Nonnegativity and boundedness of solutions). Solutions to system (2.1) with nonnegative initial conditions will remain nonnegative for all t>0. Moreover, the solutions are bounded.

    Proof. Using the notation in (3.1), assume that x00. First, we show by contradiction that S(t)0 for all t>0. Following the approach in [56,57], suppose that there exists t1>0 such that S(t1)=0, S(t1)<0, SF(t)>0, E(t)>0, I(t)>0, Iu(t)>0, Q(t)>0, R(t)>0, and D(t)>0 for t(0,t1). Then, from the first equation of system (2.1),

    S(t1)=μSF(t1)R(t1)˜N(t1)>0,

    which is a contradiction to the assumption S(t1)<0. Here we choose t1 such that SF(t1)>0 and R(t1)>0. Hence, S(t)0 for all t0. Similarly, we can show that SF, I, and Iu remain nonnegative for t0.

    Next, since I(t)0 for all t0 then Q(t)=αIγQγQ and hence, Q(t)Q(0)eγt0. It also follows that D(t)=fγQ0 and thus, D(t)0 for all t0. In a similar approach, it can be shown that E(t) and R(t) are nonnegative for all t0.

    Since N is constant and S,SF,E,I,Iu,Q,R,D0, then from the last equation of (2.1), S,SF,E,I,Iu,Q,R,DN.

    By equating the right-hand side of (2.1) to zero, the model has DFEs of the following forms:

    E1=[θN,0,0,0,0,0,ϕN,(1θϕ)N] and E2=[0,N,0,0,0,0,0,0] (3.2)

    for any 0θ1 and 0ϕ1 such that θ+ϕ1. Using the next-generation method and following the notations in [58], the matrix F of new infections and the matrix V representing the flow of individuals between compartments are given by

    F=[0ββ000000]andV=[κ00κρα0κ(1ρ)0η]. (3.3)

    The basic reproduction number R0 is calculated by finding the spectral radius of the matrix FV1 and is given by

    R0=ρβα+(1ρ)βη. (3.4)

    Moreover, the reproductive number R(t), which is the average number of secondary infections from an individual during one's infectious period, is expressed as

    R(t)=βρα(δSF(t)+S(t)N)+β(1ρ)η(δSF(t)+S(t)N).

    In the following theorems, we analyze the stability of the DFEs of the model.

    Theorem 3 (Local asymptotic stability). E1 is locally asymptotically stable if R0<1 and E2 is unstable.

    Proof. We calculate the eigenvalues λ of the Jacobian matrix derived from system (2.1) at the DFEs. For E1, the corresponding characteristic polynomial is given by

    C1(λ)=λ3(λ+γ)(λ+μ)P1(λ),

    where

    P1(λ)=λ3+a1λ2+a2λ+a3,a1=α+η+κ,a2=ϵθ(αη+ακβκ+ηκ)+ϵϕ(αη+ακ+ηκ),a3=ϵϕαηκ+ϵθ(αηκαβκ+ραβκρβηκ),

    and ϵ=1/(ϕ+θ). It is clear from C1(λ) that five of its roots are nonpositive. Thus, it is left for us to show that all of the roots of the polynomial P1(λ) are negative or have negative real parts. By the Routh-Hurwitz criteria, it is sufficient to show that a1>0, a3>0, and a1a2a3>0 [59].

    Suppose R0<1. Since all of the model parameters are positive, a1>0. Next, we can rewrite a3 as

    a3=ϵϕαηκ+ϵθαηκ(1R0). (3.5)

    Since R0<1, then a3>0. Note that

    a1a2a3=A0+ϵθ(A1+A2+A3+A4), (3.6)

    where

    A0=ϵϕ(α+η+κ)(αη+ακ+ηκ)ϵϕαηκ,A1=α2η+αηκ+αη2,A2=ακ(αρβ),A3=ηκ(η(1ρ)β),A4=κ2(αβ+η).

    It suffices to show that A0,A1,A2,A3,A4>0. By expanding A0, its last term will be cancelled out and so A0>0. Clearly, A1>0. Define

    R0,I:=ρβαandR0,Iu:=(1ρ)βη

    so that R0=R0,I+R0,Iu. Since 0<R0<1, then R0,I<1 and R0,Iu<1. Hence,

    αρβ>0andη(1ρ)β>0.

    Using the above-mentioned inequalities, it follows that A2 and A3 are both positive. Finally,

    A4=κ2(αβ+η)=κ2(A2ακ+A3ηκ)>0. (3.7)

    Therefore, E1 is stable.

    For E2, the characteristic polynomial is given by

    C2(λ)=λ3(λ+γ)(λμ)P2(λ)

    for some third-degree polynomial P2(λ). Observe that one of the eigenvalues is μ>0. Thus, E2 is unstable.

    Remark. If θ=0, then a3 in (3.5) and a1a2a3 in (3.6) are both positive regardless of the value of R0, which means that [0,0,0,0,0,0,ϕN,(1ϕ)N)] is stable. Meanwhile if ϕ=0, then [θN,0,0,0,0,0,0,(1θ)N)] is unstable if R0>1 since a3<0. By some algebraic manipulation, E1 is unstable if ϕ<θ(R01).

    In the next theorem, the global stability of the DFE E1 was studied. Lyapunov stability was established using Lyapunov's direct method and asymptotic stability was shown following the definition in [60].

    Theorem 4 (Global asymptotic stability of E1). The DFE E1 of the model (2.1) is globally asymptotically stable in the domain

    Ω:={(S,SF,E,I,Iu,Q,R,D)R8:0S,SF,E,I,Iu,Q,R,DN},

    whenever R0<1.

    Proof. Consider the Lyapunov function

    V(x)=R0E+βαI+βηIu, (3.8)

    where x is the vector of state variables as in (3.1) and R0=βρα+β(1ρ)η is the basic reproduction number. Note that V is continuous and differentiable in Ω, V(E1)=0, and V(x)>0 for xΩ, xE1.

    Next, we obtain ˙V(x)=R0E(t)+βαI(t)+βηIu(t). Note that ˙V(E1)=0. Substituting the expressions for E(t),I(t), and Iu(t) from the model and simplifying, we get

    ˙V(x)=R0δβSFI+Iu˜N+R0βSI+Iu˜Nβ(I+Iu)=β(I+Iu)(R0δSF˜N+R0S˜N1)β(I+Iu)(R01),

    since δSF+S˜N1. It follows that ˙V(x)0 for all xΩ if R0<1. Therefore, E1 is Lyapunov stable [60]. To show that the DFE E1 is globally asymptotically stable, we need to show that xE1 as t [60].

    First, we show that E,I,Iu0 by showing that V0 as t. Since V is decreasing and nonnegative whenever R0<1, then by the monotone convergence theorem, V converges to, say, ξ as t. We claim that ξ=0. Suppose otherwise, that is, ξ>0. Then ξV(x(t))V(x(0)). Define the set U:={xΩ:0<ξV(x(t))V(x(0))}. Clearly, U is compact and E1U. Since ˙V is continuous, then by the extreme value theorem, there exists ζ>0 such that

    supxU˙V=ζ<0.

    Because ˙V(x(t))ζ for all t, we get

    V(x(T))=V(x(0))+T0˙V(x(t))dtV(x(0))ζT.

    If we take T>V(x(0))/ζ, then V(x(T))<0, which is a contradiction. Therefore, V0 as t, which implies that E(t),I(t),Iu(t)0.

    Next, since R(t),D(t)0 then R(t) and D(t) are both increasing. Moreover, since R(t) and D(t) are bounded, then by the monotone convergence theorem, there exist R,D[0,N] such that R(t)R and D(t)D.

    Let ε1>0. Since I(t)0, then there exists T1>0 such that Q(t)+γQ(t)=αI<ε1 for all t>T1. Letting ε10, we have Q(t)γQ(t). Equivalently, Q(t)Q(T1)eγt. Therefore, as t, Q(t)0.

    By removing some negative terms and noting that S˜N, it follows from the second equation in (2.1) that

    SF(t)μSFR˜N+βFQ. (3.9)

    Let ε2>0. Because Q(t)0, then there exists T2>0 such that

    Q(t)<ε22βF (3.10)

    for all t>T2. Also, since R(t)R and D(t)D, then there exists T3>0 so that

    |μR(t)ND(t)ψ|<ε22N

    for all t>T3 and ψ=μRND. It follows that

    ψε22N<μR(t)ND(t). (3.11)

    Taking T=max{T2,T3}, (3.9)–(3.11) imply that

    SF(t)+(ψε22N)SF(t)<ε22 (3.12)

    for all t>T. Because SF(t)N, (3.12) becomes

    SF(t)+ψSF(t)<ε22+ε22NSF(t)ε22+ε22NN=ε2.

    Letting ε20, we get SF(t)+ψSF(t)0 as t. Thus, SF(t)SF(T)eψt and SF(t)0.

    Finally, S=N(SF+E+I+Iu+Q+R+D) implies S(t)NRD=:S.

    Because S+R+D=N, then there exist θ,ϕ[0,1] with θ+ϕ1 so that E1=[θN,0,0,0,0,0,ϕN,(1θϕ)N]. Thus, if R0<1, all solutions of (2.1) converge to the DFE of the form E1.

    The best model fit for the cumulative and daily cases is shown as the black curves in Figure 2. The red circles represent the data points. The estimated transmission rates β1 and βF,1 in period 1 were 0.199 and 471.057, respectively. In period 2, the transmission rate of the disease β2 increased to 0.361, while the transmission rate of awareness or fear of the disease βF,2 dropped to 67.783. The reporting ratio ρ1 in period 1 was estimated at 28%, which increased to 86% in period 2. The reduction in transmission induced by behavior change δ was estimated at 0.202. The parameter estimates are given in Table 1.

    Figure 2.  The best model fit to the cumulative cases from March 8 to September 30, 2020. The black curves are the plots of the cumulative cases (top) and daily cases (bottom) using the model and parameter estimates. The vertical dashed lines mark the end of periods 1 and 2. Metro Manila was under ECQ from March 8 to May 15, 2020 (dark gray), MECQ from May 16 to May 31 and August 4 to 18, 2020 (gray), and GCQ from June 1 to August 3 and August 19 to September 30, 2020 (light gray). The red circles represent the data points and the blue curve is the reproductive number R(t).

    The reproductive number R(t) is shown as the blue curve in Figure 2. Initially, R(t) was at 1.9 then decreased to 0.9 by the end of period 1. During period 2, R(t) remained above 1 from early June until mid-August, with a peak of up to 1.5 in mid-July 2020.

    Using the cumulative number of infected individuals as the model output, the results of the sensitivity analysis showed that β (range: 0.82 to 0.92) and δ (range: 0.52 to 0.68) have the highest PRCC values, followed by ρ (range: 0.49 to 0.44) and the parameters for the infectious periods, α (range: 0.48 to 0.46) and η (range: 0.45 to 0.28). PRCC values of κ declined over time, with its highest value at 0.436. The rest of the parameters have small magnitudes of PRCC. Figure 3 illustrates the results. Moreover, results of the parameter bootstrapping showed that the re-estimated values of β,βF,δ and ρ follow a normal distribution and the mean values of the estimates all fall within their respective 95% confidence intervals. Figure 4 shows the distributions of the parameter re-estimates. The mean, standard deviation, and confidence interval are also indicated in the figure.

    Figure 3.  The PRCC of the model parameters with respect to the cumulative number of infected individuals. The bars represent the PRCC values on April 19, May 31, August 2, October 4, and November 1, 2020.
    Figure 4.  The distribution of the re-estimates of β1,β2,βF,1,βF,2,δ,ρ1, and ρ2 using the parameter bootstrapping method. The mean, standard deviation (SD), and 95% confidence interval (CI) are also shown.

    The solid curves in the upper panel of Figure 5 are the plots of the susceptible class S, while the dashed lines are the behavior-changed susceptible class SF. Here we investigate what happens if people changed their behavior one (orange), two (yellow), three (purple), or four (green) weeks earlier. At the beginning of period 2, we calculated that the proportion of SF with respect to the total susceptible population was 88% (SF: 11849000; S: 1592900). To incorporate early behavior change, we scale the value of βF to yield the same proportion of SF one to four weeks before the start of GCQ. The black curves in Figure 5 represent the plots of the model using the parameters in Table 1.

    Figure 5.  The dynamics of the susceptible population (top panel; S solid and SF dashed curves), daily (lower left panel), and cumulative cases (lower right panel) if the population changed their behavior one (orange), two (yellow), three (purple), or four (green) weeks earlier than the start of period 2.

    During period 1, we observe a switch in the populations of SF and S. By the end of period 1, SF comprises the majority (88%) of the susceptible classes. The impact of early behavior change is seen in the daily and cumulative cases in period 2, shown in the bottom panels of Figure 5. As behavior changed one, two, three, or four weeks earlier, the cumulative cases decreased to 140468, 115573, 93041, or 73328 from 163191 (model, black), respectively. These translate to reductions of 14%, 29%, 43%, or 55% in cumulative cases by September 30, 2020. The peak of the daily cases also reduced to 2148, 1870, 1618, and 1392 from 2408 (model, black), and the timing was delayed from one to four weeks.

    When the community quarantine was relaxed to GCQ, we observe that the size of SF declined, while S increased. Around mid-July, when the number of daily cases were increasing (R(t)>1), the behavior of the susceptible classes switched again. From that point on until the peak of the first big wave in Metro Manila (August 14 according to the model), the proportion of SF among the susceptible classes increased from 59% to 89%.

    The upper panels in Figure 6 show the effect of varying the values of μ and βF on the cumulative cases and timing of the peak of infections. Higher values of βF and lower factors of the behavior change ease rate μ in period 2 result in notable reductions in the number of cases and delay in the occurrence of the peak. For instance, if in period 2 we set βF,2 as 3 times βF,1 and μ reduced by 90%, then the cumulative cases by the end of September 2020 would have been approximately 30000 and the peak would have occurred around July 25 (140 days from March 8). The blue area on the heatmap for peak timing indicates that the big wave in period 2 did not occur until September 2020.

    Figure 6.  The effect of varying the behavior parameters (μ and βF) in period 2 and reporting ratio in period 1 (ρ1) to the timing of the peak and cumulative cases by September 30, 2020.

    Finally, the bottom panel in Figure 6 shows the effect of increasing the reporting ratio ρ1 in period 1 on the cumulative cases by September 30, 2020. Only slight differences in the reported cumulative cases (blue bars; range: 11817 to 13421) were observed if ρ1 was increased by factors of 1.5, 2, 2.5, or 3. On the other hand, cumulative cases, including the unreported (red bars), can be reduced by 29%, 45%, 54%, or 61%, if ρ1 was increased by factors of 1.5, 2, 2.5, or 3, respectively.

    Panel (A) in Figure 7 shows the fitted bed capacity H(t) (dashed curve) from the data (red circles) and the number of cases requiring beds QH(t) from the model. Note that the red circles depict 75% of the total COVID-19 bed capacity during this time and QH(t) is 16% of Q(t), representing the average number of active cases that occupy beds. By calculating the slope of H(t), denoted by ωdata, an estimated 33 beds were added per day from June 21 to September 30, 2022, in Metro Manila. Notably, QH(t)>H(t) during the second MECQ, when healthcare workers demanded a 'timeout'. The peak of QH(t) was calculated at 5307 cases.

    Figure 7.  The Pareto optimal solutions of the bi-objective optimization problem. (A) The black curve QH(t) is the number of cases requiring beds, calculated as 16% of the reported active cases Q(t) and the black dashed line H(t) is the bed capacity obtained by fitting the data (red circles) using linear regression. (B) Pareto optimal set of (2.4). (C) Plots of QH(tω,τ) (cases requiring beds) and the optimal hospital bed capacity H(t;ω,τ) corresponding to the three Pareto optimal solutions colored blue, purple, and yellow. (D) Peaks of QH(t;ω,τ) corresponding to the Pareto optimal solutions, compared to the model.

    The optimal solutions of the bi-objective optimization problem in (2.4) form a Pareto optimal set illustrated in Figure 7 Panel (B). The circles are all optimal solutions depicting different policies. The blue-colored optimal solution corresponds to the earliest easing to GCQ on May 2, 2020, and requires 131 additional beds per day to ensure that the bed capacity is adequate (up to 75% occupancy) during the surge in cases following the lifting of restrictions (blue curves in Panel (C)). On the other hand, the yellow-colored optimal solution has the least number of additional beds per day (47 beds) and the latest start of GCQ (on May 28, 2020). This solution has a delayed and lower peak of infections compared to the blue Pareto optimal solution (see Panel (C)). Compared to the curves in Panel (A), the number of cases requiring beds QH(t;ω,τ) shown in Panel (C) is below the optimal bed capacity H(t;ω,τ). Hence, constraint (2.5) of the optimization problem is satisfied. Panel (D) shows the peak sizes of the epidemic waves corresponding to the various Pareto optimal solutions. The smallest peak size (4807 cases, purple) is the optimal solution with GCQ starting on May 20, 2020, and with at least 56 additional beds.

    Using the estimated parameter values, we observe that the model captures the trend of the daily and cumulative data from March until November 2020. The model shows a small peak in the number of daily cases (204 cases) and a slow increase in cumulative cases during period 1. A much higher peak (2408 cases) around mid-August 2020 and a sharp rise in cumulative cases are seen from the model during period 2. A delay of about one month between the drop in Rt and the decline in daily cases was also observed. Results of the parameter bootstrapping suggest good reliability of the estimated parameters. Moreover, sensitivity analysis showed that the transmission rate β was the most sensitive parameter with respect to the number of cumulative infections. A higher reporting ratio ρ or shorter mean infectious period of reported cases 1/α reduces the cumulative infections. These results suggest that intensifying testing and tracing efforts can effectively reduce new infections. The average latent period (1/κ), which has the effect of delaying infection, becomes less sensitive to the cumulative number of infections as the epidemic progresses, while the mean infectious period of unreported cases (1/η) becomes more sensitive as the epidemic progressed.

    Reporting was low in the early pandemic phase, possibly resulting from low testing capacity, slow contact tracing, uncertainty, and lack of knowledge about the disease and protocols, or fear of social stigma. This changed during period 2, where the estimated reporting ratio went up three times. These results are consistent with a study on time-varying under-reporting estimates in various countries, including the Philippines, during the same period [61]. The impact of the community quarantines imposed by the government is reflected in the reduction parameter δ and transmission rates β and βF. The estimated 20% reduction in transmission for the behavior-changed class compares with a mathematical model of COVID-19 transmission in the Philippines which showed that the minimum health standards reduced the probability of transmission per contact by 1327% [62]. As the community quarantine was relaxed from period 1 to 2, the transmission rate of the disease (β) increased from period 1 to 2, while the rate of behavior change or transmission rate of awareness or fear of COVID-19 (βF) decreased from period 1 to 2.

    Results in Figures 5 and 6 emphasize the importance of early public health campaigns and positive behavior changes (e.g., mask-wearing, improved hygiene practices, and social distancing) on reducing and delaying the peak of infections. Although fear and stigma can influence behavior changes [63,64], these can also affect reporting and negatively impact disease control. We see in Figure 6 that as reporting increased, total cumulative cases including the unreported decreased significantly.

    Without vaccines and antiviral therapy, control of epidemic diseases relies on effective NPIs and the management of healthcare systems. The bi-objective optimization approach can be used as a decision support tool because of the multiple optimal solutions provided by the method, wherein policymakers can choose depending on how much priority is given to minimizing the number of additional beds or earlier easing of restrictions. Although the method is applied to the Philippines, the optimization approach can also be used by other cities or countries by adapting location-specific epidemiological parameters. For countries with limited resources, the solutions corresponding to later easing of restrictions and a smaller number of additional beds may be better options. On the other hand, for those that can provide sufficient additional beds, the approach can be used to identify the optimal timing of adjusting NPIs. Results in Figure 7 suggest that if Metro Manila eased to GCQ on June 1, 2020, at least 47 beds per day should have been prepared so that the bed occupancy in the capital did not reach critical or high-risk, and MECQ was not needed to be reimposed. The blue solutions in Figure 7 prioritizes the earlier timing (τ) of easing protocols over the number of additional beds (ω). With this policy, GCQ could have been started 30 days earlier. However, this requires 131 additional beds per day, which is 4 times ωdata. On the other hand, the yellow solutions correspond to implementing GCQ on May 29, 2020. This would require 47 additional beds per day, which is still more than ωdata. In Panel (D), the policy in purple has the lowest peak among all the Pareto solutions. For this policy, even though GCQ starts on May 20 (12 days earlier than what happened), the peak of cases (purple curve in Panel (C)) was 500 less than the peak from the model (black curve in (A)). This policy would have required 56 beds per day, which is almost double than ωdata.

    In this work, we used an SEIQR model that considers behavior change and underreporting to study the spread of COVID-19 during the early phase of the pandemic in Metro Manila, Philippines. Behavior change can be influenced by awareness or fear of the disease, and willingness to observe NPIs such as social distancing and mask-wearing. It was incorporated into the model by introducing a two-compartment susceptible population: one for the behavior-changed population and the other for those with regular behavior. The probability of getting infected is reduced for the behavior-changed susceptible class. Due to limited testing and tracing, or negative attitudes of people towards seeking healthcare, a compartment for the unreported cases was also added.

    The results of this study highlight the importance of early behavior change and a high reporting rate in reducing the number of cases and delaying the peak of infections. These can be done by intensifying case surveillance and public health campaigns promoting compliance with NPIs, seeking healthcare, and discouraging social stigma. Moreover, this study provides an optimization approach that quantifies the additional bed requirement when policies are eased. The approach can be helpful in planning strategies that address strengthening or easing of policies, especially during the early phase when NPIs were the only control measures. Although the study focused on the transmission of COVID-19 in the Philippines, the proposed model is general enough that it can be applied to any city or country. The optimization problem can also be applied to other disease outbreaks by adjusting key epidemiological parameters.

    The model is best suited to represent the early phase of an epidemic. So, for a future work, the model can be extended to describe succeeding epidemic waves and incorporate vaccination, variants, and NPIs. Since the simulations did not consider the economic impact of NPIs, a model that maximizes the economic output of a city or country while minimizing the number of infections is another work that can be pursued. Moreover, one can also modify the model using fractional derivatives to account for memory effects [65,66]. Recently, numerous works on the theoretical analysis and applications of fractional differential equations have been done [67,68,69,70,71,72]. To the best of our knowledge, this approach has not yet been applied to a behavior change model similar to the one utilized in this study. Although a compelling research direction, this generalization requires an intensive investigation and demands a separate study.

    This paper is supported by the Korea National Research Foundation (NRF) grant funded by the Korean government (MEST) (NRF-2021M3E5E308120711). This paper is also supported by the Korea National Research Foundation (NRF) grant funded by the Korean government (MEST) (NRF-2021R1A2C100448711). The authors acknowledge Dr. Peter Julian Cayton of the UP Resilience Institute for his assistance with the data collection.

    All authors declare no conflicts of interest in this paper.



    [1] S. Ullah, M. A. Khan, J. F. Gez-Aguilar, Mathematical formulation of hepatitis B virus with optimal control analysis, Optim. Contr. Appl. Meth., 40 (2019), 529-544. doi: 10.1002/oca.2493
    [2] E. Bonyah, M. A. Khan, K. O. Okosun, et al. A theoretical model for Zika virus transmission, PLoS ONE, 12 (2017), 1-26.
    [3] S. Ullah, M. A. Khan, M. Farooq, et al. Modeling and analysis of Tuberculosis (TB) in Khyber Pakhtunkhwa, Pakistan, Math. Comput. Simulat., 165 (2019), 181-199. doi: 10.1016/j.matcom.2019.03.012
    [4] D. P. Ahokpossi, A. Atangana and D. P. Vermeulen, Modelling groundwater fractal flow with fractional differentiation via Mittag effler law, Eur. Phy. J. Plus, 132 (2017), 165.
    [5] M. A. Khan, Y. Khan and S. Islam, Complex dynamics of an SEIR epidemic model with saturated incidence rate and treatment, Phy. A., 493 (2018), 210-227. doi: 10.1016/j.physa.2017.10.038
    [6] M. Caputo and M. Fabricio, A New Definition of Fractional Derivative without Singular Kernel, Progr. Fract. Differ. Appl., 1 (2015), 73-85.
    [7] A. Atangana and D. Baleanu, New Fractional Derivatives with Nonlocal and Non-Singular Kernel: Theory and Application to Heat Transfer Model, Therm. Sci., 20 (2016), 763-769. doi: 10.2298/TSCI160111018A
    [8] M. A. Khan, S. Ullah, M. Farhan, The dynamics of Zika virus with Caputo fractional derivative, AIMS Math., 4 (2019), 134-146. doi: 10.3934/Math.2019.1.134
    [9] S. Ullah, M. A. Khan, M. Farooq, A new fractional model for tuberculosis with relapse via Atangana-Baleanu derivative, Chaos, Solitons and Fractals, 116 (2018), 227-238. doi: 10.1016/j.chaos.2018.09.039
    [10] S. Qureshi, A. Yusuf, Modeling chickenpox disease with fractional derivatives: From caputo to atangana-baleanu, Chaos, Solitons and Fractals, 122 (2019), 111-118. doi: 10.1016/j.chaos.2019.03.020
    [11] S. Qureshi, E. Bonyah, A. A. Shaikh, Classical and contemporary fractional operators for modeling diarrhea transmission dynamics under real statistical data, Phy. A, 535 (2019), 1-22.
    [12] Z. Wang, Y. K. Xie, J. Lu, et al. Stability and bifurcation of a delayed generalized fractional-order prey redator model with interspecific competition, Appl. Math. Comput., 347 (2019), 360-369.
    [13] X. Wang, Z. Wang, X. Huang, et al. Dynamic Analysis of a Delayed Fractional-Order SIR Model with Saturated Incidence and Treatment Functions, Int. J. Bifurcat. Chaos, 28 (2018), 1850180.
    [14] X. Wang, Z. Wang, H. Shen, Dynamical analysis of a discrete-time SIS epidemic model on complex networks, Appl. Math. Lett., 94 (2019), 292-299. doi: 10.1016/j.aml.2019.03.011
    [15] X. Wang, Z. Wang, J. Xia, Stability and bifurcation control of a delayed fractional-order ecoepidemiological model with incommensurate orders, J. Franklin I., 356 (2019), 8278-8295. doi: 10.1016/j.jfranklin.2019.07.028
    [16] M. A. Khan, K. Shah, Y. Khan, et al. Mathematical modeling approach to the transmission dynamics of pine wilt disease with saturated incidence rate, Int. J. Biomath., 11 (2018), 1850035.
    [17] F. Rao, P. S. Mandal and Y. Kang, Complicated endemics of an SIRS model with a generalized incidence under preventive vaccination and treatment controls, App. Math. Mod., 67 (2019), 38-61. doi: 10.1016/j.apm.2018.10.016
    [18] M. A. Khan, M. Farhan, S. Islam, et al. Modeling the transmission dynamics of avian influenza with saturation and psychological effect, Discrete and Continuous Dynamical Systems-S, 12 (2019), 455-474. doi: 10.3934/dcdss.2019030
    [19] W. O. Kermack, A. G. MCKendrick, A contribution to the mathematical theory of epidemics, Proc. Roy. Soc. Lond. A, 115 (1927), 700-721. doi: 10.1098/rspa.1927.0118
    [20] V. Capasso and G. Serio, A generalization of the Kermack-Mckendrick deterministic epidemic model, Math. Biosci., 42 (1978), 41-61.
    [21] A. Korobeinikov and P. K Maini, A Lyapunov function and global properties for SIR and SEIR epidemiological models with nonlinear incidence, Math. Biosci. Eng., 1 (2004), 57-60. doi: 10.3934/mbe.2004.1.57
    [22] I. Podlubny, Fractional differential equations: an introduction to fractional derivatives, fractional differential equations, to methods of their solution and some of their applications, Elsevier, 1999.
    [23] S. G. Samko, A. A. Kilbas, O. I. Marichev, Fractional integrals and derivatives: theory and applications, 1993.
    [24] H. Delavari, D. Baleanu, J. Sadati, Stability analysis of Caputo fractional-order nonlinear systems revisited, Nonlinear Dyn., 67 (2012), 2433-2439. doi: 10.1007/s11071-011-0157-5
    [25] C. Vargas-De-Leon, Volterra-type Lyapunov functions for fractional order epidemic systems, Commun. Nonlinear Sci. Numer. Simul., 24 (2015), 75-85. doi: 10.1016/j.cnsns.2014.12.013
    [26] J. Li, Y. Yang, Y. Xiao, et al. A class of Lyapunov functions and the global stability of some epidemic models with nonlinear incidence, J. Appl. Anal. Comput., 6 (2016), 38-46.
    [27] J. J. Wang, J. Z. Zhang and Z. Jin, Analysis of an SIR model with bilinear incidence rate, Nonlinear Analysis: Real World Applications, 11 (2010), 2390-2402. doi: 10.1016/j.nonrwa.2009.07.012
    [28] M. A. Khan, Q. Badshah, S. Islam, et al. Global dynamics of SEIRS epidemic model with non-linear generalized incidences and preventive vaccination, Adv. Diff. Equ., 2015 (2015), 88.
    [29] M. A. Khan, Y. Khan and S. Islam, Complex dynamics of an SEIR epidemic model with saturated incidence rate and treatment, Phy. A., 493 (2018), 210-227. doi: 10.1016/j.physa.2017.10.038
    [30] X. Liu and L. Yang, Stability analysis of an SEIQV epidemic model with saturated incidence rate, Nonlinear Analysis: Real World Applications, 13 (2012), 2671-2679. doi: 10.1016/j.nonrwa.2012.03.010
    [31] J. R. Beddington, Mutual interference between parasites or predators and its effect on searching efficiency, The Journal of Animal Ecology, (1975), 331-341.
    [32] D. L. DeAngelis, R. A. Goldstein and R. V. Oeill, A model for tropic interaction, Ecology, 56 (1975), 881-892. doi: 10.2307/1936298
    [33] K. Hattaf, M. Mahrouf, J. Adnani, et al. Qualitative analysis of a stochastic epidemic model with specific functional response cand temporary immunity, Phys. A, 490 (2018), 591-600. doi: 10.1016/j.physa.2017.08.043
    [34] Z. M. Odibat and N. T. Shawagfeh, Generalized taylor's formula, Appl. Math. Comput., 186 (2007), 286-293.
    [35] W. Lin, Global existence theory and chaos control of fractional differential equations, J. Math. Anal. Appli., 332 (2007), 709-726. doi: 10.1016/j.jmaa.2006.10.040
    [36] M. Toufik, A. Atangana, New numerical approximation of fractional derivative with non-local and non-singular kernel: application to chaotic models, Eur. Phys. J. Plus, 132 (2017), 444.
  • Reader Comments
  • © 2020 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(6054) PDF downloads(777) Cited by(36)

Figures and Tables

Figures(9)  /  Tables(1)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog