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

1.   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-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-7], ratio-dependent incidence function [8], Beddington-DeAngelis incidence function [9], and some authors considers even generalized incidence function as the researches [10-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^{\delta }$ 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$. $\Lambda$ is the constant entering flux into the susceptible populations. $\mu$ is the natural mortality rate. $\beta$ is the transmission rate of individuals with severe symptoms. $q\beta$ is the transmission rate corresponding to individuals in A-class. $\alpha, \; \gamma$ are positive constants related to the generalized Holling type III incidence function. This function is $\frac{\beta I^2}{\alpha+\gamma I+I^2}$ which increasing function in $I$ and saturates to $\beta$ 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. $\delta, \; \eta$ are recovering rates for individuals in A-class, I-class, respectively. $0 < \phi < 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.   Preliminary

Definition 2.1. The fractional Caputo's derivative with $\delta \in (0, 1)$ for $p:\mathbb{R }^+\rightarrow \mathbb{R}$ can be expressed:

$\Gamma$ is Euler's Gamma functional. We put

$\delta \in (0, 1), \; g\in C^{1}(\mathbb{R}^{n}, \mathbb{R}^{n}), \; A\in\mathbb{R}^{n\times n}, \; Dg(0) = 0.$ The local stability concept for fractional differential systems is provided through the following theorem:

Theorem 2.2. We presume that the origin is an equilibrium of (2.1), the linear stability of the origin is guaranteed if each eigenvalue denoted $\lambda$ of $A$ verifies $\left\vert \arg (\lambda)\right\vert > \frac{\delta \pi }{2},$\ and it is unstable if $\left\vert \arg (\lambda)\right\vert < \frac{\delta \pi }{2}$ for some values of $\lambda.$

Now, defining the Hopf bifurcation for the fractional order system with a parameter $\mu \in \mathbb{R}$ as

where $\delta \in (0, 1), A(\mu)\in\mathbb{R} ^{n\times n}, g(X, \mu)\in C^{1}(\mathbb{R}\times\mathbb{R}^{n}, \mathbb{R} ^{n}), Dg(0, \mu) = 0$.

It is well known that the Hopf bifurcation conditions for the system (2.2) in the case of the first order derivative are

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.

3.   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:

3.1.   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 $\left(\frac{\Lambda}{\mu}, 0, 0 \right)$. 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^* < \frac{\Lambda}{\mu}$. By substituting (3.4), (3.5) into the second equation of (3.2) we obtain:

where

We denote by $H(S)$ the right hand side of the Eq (3.6). Now, we seek for the zeros of the function $H$ for $S\in \left(0, \frac{\Lambda}{\mu}\right)$. The derivative of $H$ with respect to $S$ is

We put $H'(S) = 0$ then the following quantity has an important role in determining the number of the positive zeros of $H$:

where

Obviously, if $B_0 < b_0$ then $H$ is increasing in $S$. The existence and the uniqueness of PE is guaranteed if ($\Lambda > \Lambda_0$ and $B_0 > b_1$) or ($\lambda < \Lambda_0$ and $B_0 < b_1$), where

and

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:

$\rm (i)$ The equation $H(S) = 0$ has no roots when [$H(S_-) < 0,$ and ($\Lambda < \Lambda_0$ and $B_0 > b_1$) or ($\Lambda > \Lambda_0$ and $B_0 < b_1$)].

$\rm (ii)$ The equation $H(S) = 0$ has one root when [$H(S_-) < 0,$ and ($\Lambda > \Lambda_0$ and $B_0 > b_1$) or ($\Lambda < \Lambda_0$ and $B_0 < b_1$)], or [$H(S_+) > 0,$ and $S_+ < \frac{\Lambda}{\mu}$].

$\rm (iii)$ The equation $H(S) = 0$ has two positive roots when [$H(S_-) > 0,$ and $S_- < \frac{\Lambda}{\mu}$, and [ ($\Lambda < \Lambda_0$ and $B_0 > b_1$) or ($\Lambda > \Lambda_0$ and $B_0 < b_1$)]],

$\rm (iv)$ The equation $H(S) = 0$ has three positive roots when [$H(S_-) < 0,$ $H(S_+) > 0,$ $S_+ < \frac{\Lambda}{\mu},$ and ($\Lambda > \Lambda_0$ and $B_0 > b_1$) or ($\Lambda < \Lambda_0$ and $B_0 < b_1$)].

Now we denoted the following assumptions:

$\mathbf{(A_1)}$: [$B_0 < b_0$, and ($\Lambda > \Lambda_0$ and $B_0 > b_1$) or ($\Lambda < \Lambda_0$ and $B_0 < b_1$)] or [$B_0 > b_0$ and $H(S_-) < 0,$ and ($\Lambda > \Lambda_0$ and $B_0 > b_1$) or ($\Lambda < \Lambda_0$ and $B_0 < b_1$), ] or [$B_0 > b_0$, and $H(S_+) > 0,$ and $S_+ < \frac{\Lambda}{\mu}$].

$\mathbf{(A_2)}$: $B_0 > b_0$ and [$H(S_-) > 0,$ and $S_- < \frac{\Lambda}{\mu}$, and [($\Lambda < \Lambda_0$ and $B_0 > b_1$) or ($\Lambda > \Lambda_0$ and $B_0 < b_1$)]].

$\mathbf{(A_3)}$: [$B_0 > b_0$ and $H(S_-) > 0,$ $H(S_+) < 0,$ $S_+ < \frac{\Lambda}{\mu},$ and ($\Lambda > \Lambda_0$ and $B_0 > b_1$) or ($\Lambda < \Lambda_0$ and $B_0 < b_1$)].

The obtained results are summarized in the following theorem:

Theorem 3.1. $\rm (i)$ If $\mathbf{(A_1)}$ holds then the system (3.1) has a unique PE denoted $E^* = (S^*, A^*, I^*)$.

$\rm (ii)$ If $\mathbf{(A_2)}$ holds then the system (3.1) has a two PEs denoted $E^*_i = (S^* _i, A^* _i, I^* _i), \; \; i = 1, 2$.

$\rm (iii)$ If $\mathbf{(A_3)}$ holds then the system (3.1) has a three PEs denoted $E^*_i = (S^* _i, A^* _i, I^* _i), \; \; i = 1, 2, 3$,

where $S^*_i$ are the positive roots of the cubic Eq (3.6) verifying $S^*_i < \frac{\Lambda}{\mu}$, $A^*_i = \frac{\phi}{\delta+\mu}(\Lambda-\mu S^*_i)$ and $I^*_i = \frac{1-\phi}{\eta+\mu+m}(\Lambda-\mu S^*_i).$

Remark 3.2. Based on the above results and by taking $B_0$ as bifurcation parameter, the saddle-node bifurcation can appears at $H(S_-) = 0,$ (which it can be written as $B_0 = B_0 ^{SN_-}$) and $B_0 ^{SN_-}$ verifies either $b_0 < B_0 ^{SN_-} < b_1$ if $\Lambda > \Lambda_0$ or $\max\{b_0, b_1\} < B_0 ^{SN_-}$ if $\Lambda < \Lambda_0$,

Also, at $H(S_+) = 0,$ (which it can be written as $B_0 = B_0 ^{SN_+}$) and $B_0 ^{SN_+}$ verifies either $b_0 < B_0 ^{SN_+} < b_1$ if $\Lambda > \Lambda_0$ or $\max\{b_0, b_1\} < B_0 ^{SN_+}$ if $\Lambda < \Lambda_0$.

3.2.   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 $\lambda_1 = -\mu$, $\lambda_2 = \frac{\phi q \beta\Lambda}{\mu} -(\delta + \mu),$ and $\lambda_3 = -(\eta+\mu + m)$, hence $|\arg (\lambda_1)| = \pi > \frac{\delta \pi}{2}$, $|\arg (\lambda_3)| = \pi > \frac{\delta \pi}{2}$, and

hence the DFE is locally stable for $\Lambda < \Lambda_0$, and unstable for $\Lambda > \Lambda_0$. Transcritical bifurcation appears at $\Lambda = \Lambda_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

$\rm (i)$ $R > 0,$ $C_{2} > 0,$ $C_{0} > 0,$ $C_{2}C_{1} > C_{0}.$

$\rm (ii)$ $R < 0,$ $C_{2}\geq 0,$ $C_{1}\geq 0,$ $C_{0}\geq 0,$ and $\delta < \frac{2}{3}$.

Now let us focus on providing some sufficient conditions for obtaining Hopf bifurcation by considering $\delta$ as bifurcation parameter. This proof is inspired by results obtained in Proposition 2 in [25]. Now we set

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 \{\lambda \}| = \frac{\delta \pi}{2}$).

Case 2: If $d < 0$ then the characteristic Eq (3.11) has one real root $\lambda_1 = -a$, the other two roots are $\lambda_{2, 3} = a_1 \pm i a_2$. Hence, the characteristic Eq (3.11) becomes: using the change of variable $Y-\frac{C_2}{3}$ we get the Eq

where $p = C_1 - \frac{C_2 ^2}{3}$, and $q = \frac{C_1}{27}(2 C_2 ^2 - 9 C_1) + C_0$. Using Cardano's method the real roots of (3.12) is

where $r = \left(\frac{-q + \sqrt{q ^2+\frac{4}{27}p^3}}{2} \right)^{\frac{1}{3}}$, and $e = \left(\frac{-q - \sqrt{q ^2+\frac{4}{27}p^3}}{2} \right)^{\frac{1}{3}}$. The complex roots are:

with $j = \exp ^{i \frac{2 \pi }{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:

this implies $\lambda_1 < 0$, hence $|\arg\{ \lambda _1 \}| > \frac{\delta \pi}{2}$.

Besides, $|\arg\{ \lambda _n \}| = \left|\tan ^{-1} \left(\frac{3\sqrt{3}(e-r)}{3(r+e)+2 C_2} \right) \right|, \; \; n = 1, 2.$ Hence, for $\frac{2}{\pi}\left|\tan ^{-1} \left(\frac{3\sqrt{3}(e-r)}{3(r+e)+2 C_2} \right) \right| > \delta$ then $|\arg\{ \lambda _n \}| > \frac{\delta \pi}{2}, \; \; n = 1, 2$ which means that PE is stable, and unstable for $\frac{2}{\pi}\left|\tan ^{-1} \left(\frac{3\sqrt{3}(e-r)}{3(r+e)+2 C_2} \right) \right| < \delta$. Hopf bifurcation occurs at $\frac{2}{\pi}\left|\tan ^{-1} \left(\frac{3\sqrt{3}(e-r)}{3(r+e)+2 C_2} \right) \right| = \delta$.

The obtained results are summarized in the following theorem:

Theorem 3.4. Assume that one of the assumptions $(A_1), (A_2), (A_3)$ holds, then we have the following results

$\rm (i)$ If $d > 0$ then the system (3.1) cannot undergo Hopf bifurcation.

$\rm (ii)$ If $d < 0$ and $C_0 < 0$ then the system (3.1) cannot undergo Hopf bifurcation.

$\rm (iii)$ If $d < 0$ and $C_0 > 0$ then the system (3.1) undergo Hopf bifurcation at $\delta = \frac{2}{\pi}\left|\tan ^{-1} \left(\frac{3\sqrt{3}(e-r)}{3(r+e)+2 C_2} \right) \right|$.

3.3.   Numerical scheme using trapezoidal product-integration rule

Now investigating the following fractional initial-value problem:

By applying the fundamental theorem of fractional calculus obtained in (3.14), we get

Using $t = t_{n} = n\hslash$ 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,

We apply the achieved numerical scheme (3.18) to (3.1) we obtain:

with

3.4.   Graphical representations

Figure 1: In this figure we consider the following set of parameters

where for (A) we have considered $\alpha = 0.05$, and we obtained $\Lambda_0 = 0.011,$ $B_0 = 0.0511$, $B_1 = 0.3912$, $B_2 = -0.3726$, hence we get $S_- = 0.5095 < \frac{\Lambda}{\mu} = 2$. In this case we get the second axiom in Theorem 3.1.

For (B) we considered $\alpha = 0.2$ hence we get $\Lambda_0 = 0.011,$ $B_0 = 0.2011$, $B_1 = 0.4594$, $B_2 = -0.3726$, hence we get $S_- = 0.6191 < \frac{\Lambda}{\mu} = 2$. In this case we get the non existence of PEs.

Now we consider the following set of parameters

where for (C) we have considered $q = 4.5, \; \; \phi = 0.1$ and we obtained $\Lambda_0 = 0.014,$ $B_0 = 0.04$, $B_1 = 0.7828$, $B_2 = -0.7622$, hence $H'(S) > 0$ which verifies the first axiom in Theorem 3.1. In this case we get the second axiom in Theorem 3.1.

For (D) we have considered $q = 1.5, \; \; \phi = 0.55$ and we obtained $\Lambda_0 = 0.0121,$ $B_0 = 0.08$, $B_1 = 1.255$, $B_2 = -1.0531$, hence we obtain the results obtained in the third axiom in Theorem 3.1. In this case we get the second axiom in Theorem 3.1.

Figure 2: In this figure we consider the following set of parameters

where the system (3.1) undergo saddle-node bifurcation at $B_0 = B_0 ^{SD_-} = 0.6776$.

Figure 3: We consider the following set of parameters

and the initial conditions

Figure 4: In this figure we consider the following set of parameters

where for (A) we have considered, and we obtained $\Lambda_0 = 0.011,$ $B_0 = 0.0511$, $B_1 = 0.3912$, $B_2 = -0.3726$, and the initial conditions

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.

Figure 5: In this figure we consider the following set of parameters

and the initial conditions

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.

4.   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_0 ^{SB}$. 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-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.

Conflict of interest

The authors declare that there is no known conflict of interest.