Mathematical analysis of a fractional-order epidemic model with nonlinear incidence functio

1 Faculty of Exact and Computer Sciences, Mathematics Department, Hassiba Benbouali university, Chlef, Algeria 2 Laboratoire d’Analyse Non Linéaire et Mathématiques Appliquées, University of Tlemcen, Tlemcen, Algeria 3 Institute for Groundwater Studies, faculty of natural and agricultural science, University of the Free State, Bloemfontein, 9300, South Africa 4 Department of Mathematics, COMSATS University Islamabad, Abbottabad Campus, Abbottabad, 22060, Khyber Pakhtunkhwa, Pakistan 5 Research Institute for Natural Sciences, Hanyang University, Seoul 04763, Korea


Introduction
Mathematical models visualize the method of spreading infectious diseases and the potential outcome of a pandemic and it is very helpful in making a proper decision for the public health interventions. Models use basic presumptions or pooled statistics (which it can be used to approximate the real value of the parameters as example the papers [1][2][3]) along with mathematics to obtain criteria for different infectious diseases and use these parameters to determine the influence of various interventions, such as mass vaccination programs, asymptomatic effect, so on. Modeling can help determine which intervention(s) are more proper to avoid the fatality of the disease and the degree of the outbreak of this infection, also, can determine which measure is efficient, and highlighting the outcome of considering a specific measure next to predicting future growth patterns, etc [4].
In mathematical epidemiology, the crucial component that describes the speed and the manner of the spread of a contagious disease in a sample of the population is incidence function. There are numerous researchers that study the influence of this important component, as for example Holling I-III incidence function [5][6][7], ratio-dependent incidence function [8], Beddington-DeAngelis incidence function [9], and some authors considers even generalized incidence function as the researches [10][11][12][13][14][15][16][17][18][19][20]. These incidence functions highlights the method of the transmission of the disease and can model different reaction of population (as fear, caution, using protection materials) and can reflect the measures taken by government as isolation and crowding behavior (that can lead eventually to eliminate the disease in the case when the number of infected persons is high), which shows the big importance of choosing the right incidence function that describes the studied disease. There are some infectious diseases where the infection will affect some individuals and do not affect the other individuals (as COVID-19 disease), where there are some infected persons with a visualized symptoms, or even severe, these individuals have more chance to develop complications that can lead to death. The other part of infected individuals has a non visualized symptoms, this category has more resistance to infection, and developing complications that lead to death is not probable. The importance of this approach (asymptomatic and symptomatic infected person in model building) can be seen when the individuals interact with the two categories of infected persons, wherein the case of the large size of the infected class, the susceptible persons will avoid the contact with the visualized symptoms, and cannot recognize persons with mild symptoms. In terms of modeling this behavior, it is wise to choose a saturated incidence function for modeling the avoiding infected persons with severe symptoms (where we will choose generalized Holling type III incidence function [21]), and a non-saturated incidence functional for modeling the transmission of the infection for infected persons with mild symptoms (where we will choose the bilinear incidence function). Furthermore, we will investigate the effect of memory obtained by considering the fractional-order derivative in such a model. This point of view attract numerous researches as example [22]. As a result, we consider the following fractional-order system: where D δ is Caputo's derivative operator concerning t, which will be defined in the next section. S , A, I R are respectively the densities of the susceptible individuals, asymptomatic infected individuals, symptomatic infected individuals, recovered individuals at time t. Λ is the constant entering flux into the susceptible populations. µ is the natural mortality rate. β is the transmission rate of individuals with severe symptoms. qβ is the transmission rate corresponding to individuals in A-class. α, γ are positive constants related to the generalized Holling type III incidence function. This function is βI 2 α+γI+I 2 which increasing function in I and saturates to β when I goes to infinity. The reason behind choosing such as incidence function due to the awareness of the population when the number of infection cases is high Contrarily, the asymptomatic individuals will continue infecting persons without being revealed, which is the main reason behind choosing linear incidence function for this class on infected persons, where this category is characterized by non saw symptoms. δ, η are recovering rates for individuals in A-class, I-class, respectively. 0 < φ < 1 represents the probability of the newly infected person to a person with asymptomatic symptoms. m is the mortality due to this infection due to complications and severe symptoms.
The main goal of this research is to study mathematically the model (3.1). For this aim we organize our research in the following manner. In the next section, we will provide some preliminary results that include some definitions that can be very helpful in dealing with the fractional operator. In section 3 we will analyze the model (3.1) mathematically, where we will calculate the equilibria of this system and proving that the system has saddle-node bifurcation. For the ecdemic equilibrium, we will show that the system (3.1) undergoes Hopf bifurcation. The second part of the mathematical analysis is to provide a numerical scheme that can be very helpful in the graphical representation of the solution of the system (3.1), where the trapezoidal product-integration rule has been used for building this numerical scheme. As a result, some graphical representations are plotted and properly interpreted, which insures the mathematical finding. The conclusion section ends the research.
(2.1) δ ∈ (0, 1), g ∈ C 1 (R n , R n ), A ∈ R n×n , Dg(0) = 0. The local stability concept for fractional differential systems is provided through the following theorem: We presume that the origin is an equilibrium of (2.1), the linear stability of the origin is guaranteed if each eigenvalue denoted λ of A verifies arg(λ) > δπ 2 , and it is unstable if arg(λ) < δπ 2 for some values of λ.
Now, defining the Hopf bifurcation for the fractional order system with a parameter µ ∈ R as It is well known that the Hopf bifurcation conditions for the system (2.2) in the case of the first order derivative are Re But the conditions for Hopf bifurcation in the fractional order derivative is given in the following form To mention that the presence of Hopf bifurcation in the case of the first-order derivative means the possibility of having a periodic solution under some conditions on the parameters (in the case of stable periodic orbits). In the case of the presence of fractional order derivative, it is been shown that a system with fractional derivative cannot have periodic solutions. So, the proof of the existence of Hopf bifurcation in the case of the fractional derivative does not mean the existence of periodic solutions [23], where we can mention the presence of oscillations in time only.

Mathematical analysis of fractional model (3.1)
In this section, we are interested in providing a qualitative analysis of the solution of the system (3.1). At first, we can easily remark that the three first equations of (3.1) are independent of R, hence the R-equation (fourth equation of the same system) can be omitted. So the behavior of the system (3.1) can be deduced through studying the following reduced system: Based on the fact that the right hand side of (3.1) is continuously Lipcshitz we deduce the existence and uniqueness of solution. For achieving this result, we split our section into the following subsections:

Equilibria
The equilibria of the fractional system (3.1) are the solution of the following system: Obviously, the system (3.2) has always disease free equilibrium (DFE) which is Λ µ , 0, 0 . Now focusing on distinguishing the existence conditions for the positive equilibrium (PE), which is the positive solution of the system (3.2). From the first equation of (3.2) we have: substituting this result into the second equation of (3.2) we obtain Similarly, we have hence the positivity condition appears which is S * < Λ µ . By substituting (3.4), (3.5) into the second equation of (3.2) we obtain: We denote by H(S ) the right hand side of the Eq (3.6). Now, we seek for the zeros of the function We put H (S ) = 0 then the following quantity has an important role in determining the number of the positive zeros of H: Obviously, if B 0 < b 0 then H is increasing in S . The existence and the uniqueness of PE is guaranteed if (Λ > Λ 0 and B 0 > b 1 ) or (λ < Λ 0 and B 0 < b 1 ), where 2 . Now we presume that B 0 > b 0 hence the equation H (S ) = 0 has the following roots: For the goal of discussing the roots of H(S ) = 0, we have the following results: (ii) If (A 2 ) holds then the system (3.1) has a two PEs denoted E

Local stability and Hopf bifurcation
In this section, we are interested in determining the stability of the equilibria. For achieving this result we calculate the Jacobian matrix for an arbitrary equilibrium E = (S , A, I), which is expressed in the following manner: At the DFE the Jacobian matrix (3.8) becomes: Based on the results mentioned in the second section of the paper, we can highlight that the Hopf bifurcation appears if the following condition holds It is easy to see that the eigenvalues of (3.9) are λ 1 = −µ, λ 2 = φqβΛ µ − (δ + µ), and λ 3 = −(η + µ + m), hence | arg(λ 1 )| = π > δπ 2 , | arg(λ 3 )| = π > δπ 2 , and hence the DFE is locally stable for Λ < Λ 0 , and unstable for Λ > Λ 0 . Transcritical bifurcation appears at Λ = Λ 0 . Now, we Study the stability of the PE(s), where we denote E * = (S * , A * , I * ) for one of the PEs discussed in theorem 3.1. The Jacobian matrix at this equilibrium can be expressed in the following form: The characteristic equation is where . Now defining the quantities: Making use of the Routh-Hurwitz criteria for the fractional order derivative (see [24]) we get the following stability conditions: Theorem 3.3. The stability of the PEs is guaranteed if one of the following conditions holds: If one or both of the following conditions may apply, then the equilibrium(s) E 2 (E i , i = 1, 2), is (are) stable 3 . Now let us focus on providing some sufficient conditions for obtaining Hopf bifurcation by considering δ as bifurcation parameter. This proof is inspired by results obtained in Proposition 2 in [25]. Now we set d = C 2 2 − 4C 1 , hence we consider two cases: Case 1: If d > 0 then the characteristic Eq (3.11) has three real roots hence we cannot have the Hopf bifurcation condition (this means that we cannot have | arg{λ}| = δπ 2 ). Case 2: If d < 0 then the characteristic Eq (3.11) has one real root λ 1 = −a, the other two roots are λ 2,3 = a 1 ± ia 2 . Hence, the characteristic Eq (3.11) becomes: using the change of variable Y − C 2 3 we get the Eq Y 3 + pY + q = 0, (3.12) where p = C 1 − C 2 2 3 , and q = C 1 27 (2C 2 2 − 9C 1 ) + C 0 . Using Cardano's method the real roots of (3.12) is where r = −q+ √ q 2 + 4 27 p 3 2 1 3 , and e = −q− √ q 2 + 4 27 p 3 2 1 3 . The complex roots are: Y 2 = jr +je, and Y 3 = j 2 r +j 2 e, with j = exp i 2π 3 , hence the roots of (3.11) are by taking into count the obtained roots, then the characteristic Eq (3.11) can be rewritten as the form: Hence we get: 3(r+e)+2C 2 , n = 1, 2. Hence, for 2 π tan −1 3 √ 3(e−r) 3(r+e)+2C 2 > δ then | arg{λ n }| > δπ 2 , n = 1, 2 which means that PE is stable, and unstable for 2 π tan −1 3 √ 3(e−r) 3(r+e)+2C 2 < δ. Hopf bifurcation occurs at 2 π tan −1 3 The obtained results are summarized in the following theorem: (iii) If d < 0 and C 0 > 0 then the system (3.1) undergo Hopf bifurcation at δ = 2 π tan −1 3 3(r+e)+2C 2 .

Numerical scheme using trapezoidal product-integration rule
Now investigating the following fractional initial-value problem: D δ X(t) = P(t, X(t)). (3.14) By applying the fundamental theorem of fractional calculus obtained in (3.14), we get Using t = t n = n in (3.15), we arrives at: Using the first order Lagrange interpolation we can approximate function P(t, U(t)) as: where X i = X(t i ). By substituting (3.17) into (3.16) we get (for more details see [26,27]) where, , n = 1, 2, . . . .

(3.19)
We apply the achieved numerical scheme (3.18) to (3.1) we obtain: with P 1 (S , A, I) = Λ − S µ + βI 2 α+γI+I 2 + qβA , P 2 (S , A, I) = φS βI 2        In this case we get the stability of the equilibrium E * 1 = (0.2122, 0.18237, 0.84927), where we obtained the persistence of the disease. In this case we get the stability of the equilibrium E 0 = (2.1948, 0, 0), which highlights the extinction of the disease. Remark 3.5. To mention that the parameters values used in this research are used only for discussing the mathematical finding through this research, and it can be used for predicting the disease in the case of availability of information's on the studied disease.

Discussion and conclusions
We dealt in this research with a new mathematical model that considers two different classes of infected persons, namely, symptomatic infected class, and asymptomatic infected class with a nonlinear incidence function. The main goal is to study the bifurcating solution of the system (3.1) in the presence of the fractional derivative operator. In fact, it is obtained that the system (3.1) has important different types of bifurcation, as saddle-node bifurcation that appears at B 0 = B S B 0 . This type of bifurcation is known as well by backward bifurcation (which consists the coexistence of two equilibria as in Figure 2) for epidemiological models (see [28,29]). Further, it is obtained that the system (3.1) can undergo Hopf bifurcation at the endemic equilibria, where the possibility of having this kind of bifurcation has been shown using Cardon's method for solving third order equation. As a result, we obtained the aspects highlighted in Theorem 3.4. Next, for plotting the solutions of the system (3.1) we build our numerical scheme using the trapezoidal product-integration rule. The method of proving Hopf bifurcation can be applied for other paper as [30][31][32][33], where they studied the stability of the positive equilibria only, which it can be very helpful in showing the existence of Hopf bifurcation for a fractional system with three equations. In future works, including random walk (or stochastic effect) into the model (3.1) can generate important behavior, and can be considered as a good subject of interest.