ON THE STABILITY ANALYSIS OF A FRACTIONAL ORDER EPIDEMIC MODEL INCLUDING THE GENERAL FORMS OF NONLINEAR INCIDENCE AND TREATMENT FUNCTION

. In this paper, we propose to study a SEIR model of fractional order with an incidence and a treatment function. The incidence and treatment functions included in the model are general nonlinear functions that satisfy some meaningful biological hypotheses. Under these hypotheses, it is shown that the disease free equilibrium point of the proposed model is locally and globally asymptotically stable when the reproduction number R 0 is smaller than 1. When R 0 > 1, it is established that the endemic equilibrium of the studied system is uniformly asymptotically stable. Finally, some numerical simulations are provided to illustrate the theory.


Introduction
Studying the spread process of infectious diseases has been a very important and popular topic since outbreaks have serious impacts on the economy, daily lives and the future.For finding intervention strategies or treatments and reducing the deaths understanding this mechanism is very important.Mathematical models help us understand the dynamics of epidemiological diseases and talk about the future of the epidemics.Until today, a vast number of mathematical models have been developed for diseases such as rabies, measles, malaria, chickenpox, tuberculosis, cancer, HIV/AIDS and COVID-19 [1, 7-9, 16, 18, 26, 36-38, 40, 46].When modeling disease transmission, compartmental models such as SIR (Suscepted-Infected-Recovered), SIS, SEIR and SEIS models are mostly used in the literature (for more detail, see [8]).The general idea behind these compartmental models is to divide the total population into compartments and describe the transfer from one compartment to 286 E. KARAO GLU another under some meaningful assumptions [21].For example, for an SIR model, the population is divided into three compartments: susceptible individuals S who are not infected yet, infected individuals I who are infectious and can spread the disease to the susceptible class and recovered (removed) individuals R who recovered the infection and already gained the immunity.
In epidemic models, incidence functions are used in the transmission progress between the susceptible population and the infected or exposed population.In literature, there are different incidence rates and using different incidence rates can affect the dynamic behaviour of the system [27,31].The most common incidence rates are the bilinear incidence rate βSI [45,50], where β is the average number of contacts per individual per day, the standard incidence rate βSI/N [14,25] where N = S + I + R, the saturated incidence rate βSI/(1 + α 1 S + α 2 I) [16,19] where α 1 , α 2 are positive constants.Even though it is very hard to fit real data values in infectious disease transmission, nonlinearity is inevitable in the incidence rates.These nonlinear incidence rates seem more realistic because they may include saturation effects, heterogeneous mixing populations, environmental factors, media effects or behavioural changes of individuals, etc. [23].In 2005, Korobeinikov and Maini [24] studied the stability properties of infectious disease models with a general, arbitrary nonlinear incidence rate f (S, I, N ) and obtained the global stability by constructing a Lyapunov function under a more specific incidence rate of the form g(I)h(S).Following this work, Korobeinikov studied the global dynamics of infectious disease models with nonlinear incidence rates in [22] and [23], respectively.In 2014, Li et al. [27] considered an SIR epidemic model with a nonlinear arbitrary incidence function f (S, I) and they improved their model by incorporating a time delay representing the latent period.Recently, in [48], the authors analyzed the stability of a fractional order SEIR model with general incidence rate F (S)G(I) and in [20], the authors studied the local and global stability of the disease free and endemic equilibrium points of a fractional SIR model with a general incidence function f (S, I).
Treatment is an important strategy to reduce the number of infected people during outbreaks.Vaccination can be thought one form of treatment to protect against infection before the outbreak [6].One of the most efficient methods of treatment is, of course, hospitalisation, but sometimes the resources of countries are not adequate as in the COVID-19 pandemic.In mathematical models, to reflect and analyse the effect of treatment, scientists have incorporated another treatment class into the SIR model [6,42,47], or they used some treatment functions in the models [11,15,28].In 1991, Anderson and May [4] proposed that the treatment function is proportional to the number of infectious people.Following this work, some different modified treatment functions are used in order to reflect the treatment capacity of communities (For further details, see [11,28]).
In recent decades, fractional differential equations (FDEs) have gained great attention in the stability analysis of ecological or epidemiological models.Besides FDEs being generalizations of classical differential equations, using FDEs in epidemiology helps to model in a more realistic way since they reflect history (memory) [16,39].There are different kinds of fractional operators in the literature and these operators may answer distinct real world problems [32].Moreover, because of the memory effect, FDEs will be more suitable for epidemic models.For example, in 2020, Naik et al. [32], proposed and studied the stability of a FDE system that models COVID-19 pandemic with Atangana-Baleanu or Caputo derivative.They divided the total population into eight groups such as suspected, exposed, symptomatic (infected), asymptomatic, quarantined, treated classes etc. and used the real data from Pakistan.In 2020, Yavuz and Sene [49] studied the stability analysis of a fractional predator-prey model with a harvesting rate.In another paper, Naik et al. [33] has established the global dynamics of a fractional order model for the transmission of HIV epidemic with optimal control in 2020.In 2023, Joshi et al. [16] studied COVID-19 pandemic with Atangana-Baleanu derivative.They used an SIR model with nonlinear Beddington-DeAngelis infection rate and Holling type II treatment rate in their paper.Recently, a vast number of papers containing FDEs for different research areas have been published in the literature [5,16,17,32,34,41].
Motivated by the aforementioned works [20,23,27], in this paper, we have proposed an SEIR model including FDEs with a general nonlinear incidence function f (S, I) and a treatment function T (I): where D α t u represents Caputo fractional derivative of the function u with the following initial conditions: (2) In this model (1), S, E, I and R denote the suscepted, exposed, infected, and recovered individuals, respectively.To the best of the author's knowledge, a fractional SEIR model with a general incidence function and treatment function has not been studied yet.We have chosen this model with exposed individuals compartment, as most infectious diseases have an incubation period.Before explaining parameters, we need to emphasize that the originality of this paper comes from the choice of the incidence and treatment functions, f (S, I) and T (I) functions in model (1), respectively.These functions have not been determined specifically so that depending on the studied disease, one may choose his/her function according to the spread of the disease and treatment type.Moreover, considering these functions in a general way increases the complexity of the proofs.
In model (1), parameter λ is the recruitment rate which represents the total change in the population and assumed as a positive number, β is the rate at which exposed individuals become infectious (incubation rate), µ is the natural death rate, θ is the death rate depending on the infection, and r is the recovery rate of exposed individuals.The function T (I) represents the general treatment function.In the model, it is assumed that unsuccessfully treated infectious individuals re-enter the exposed compartment proportional to parameter p and the parameter q denotes the fraction of infectious individuals whose treatments are successful (p = 1 − q).The flowchart for the model ( 1) is given in Figure 1.The paper is organised as follows: In Section 2, the definition of Caputo fractional derivative is presented and some lemmas are given for the proofs needed for stability analyses.In Section 3, the properties of the incidence function and the treatment function are analysed and the positivity of the solution of system (1) is proved.After that, in Section 4, the equilibrium points of system (1) are determined and the global stability analysis of disease free equilibrium point and the uniform asymptotic stability of endemic equilibrium point are established.Following these theorems, some numerical simulations are carried out to show some examples in Section 5. Finally, we finish this paper with a conclusion part in Section 6.

Preliminaries on the Caputo Fractional Calculus
We begin by introducing the definition of Caputo fractional derivative.
To prove the nonnegativity of the solutions of model (1) we will need the following lemma.

Basic Properties of the Model
In model (1), we assume that the functions S, E, I and R and their Caputo fractional derivatives are continuous when t > 0.
The general nonlinear incidence function f (S, I) and the treatment function T (I) are considered positive, continuously differentiable functions and they satisfy the following hypotheses: H1) f (S, I) > 0, f (0, I) = 0, f (S, 0) = 0 for all S, I > 0. These conditions are consistent with biological assumptions and in accordance with literature (see [11,20,23]).For example, for (H1 ),we can think that there will be no transmission when there are no susceptible or infected people.For (H2 ), we understand that the incidence function is a monotonically growing function for all S, I > 0. In the absence of an infected person, susceptible individuals will become stagnant and transmission will begin to increase in the case of an infected person.In addition, in the absence of an infected person, there is no need for treatment.When there is an increase in the rate of transmission, that is, the number of infected people increases, we need to apply more treatment strategies (H5 ).Now, let R 4 + = {(x 1 , x 2 , x 3 , x 4 ) ∈ R 4 : x i ≥ 0 for all i = 1, 2, 3, 4}.We will prove the existence, uniqueness and positivity of the solutions with the following theorem.
Theorem 1.There exists a unique solution of the model ( 1) with initial conditions under the hypothesis (H1) is satisfied.Moreover, the solution will remain in R 4 + for all t ≥ 0.
Proof.From Theorem 3.1 and Remark 3.2 of [29], one can see the existence and uniqueness of the solution of system (1) with the initial conditions (2).Now, we will prove that R 4  + is a positively invariant region.For this, let We can make a similar discussion as in Theorem 2 of [3] and with the help of Corollary 1, one can observe that the solution will remain in R 4  + for all t ≥ 0. □

Equilibrium Points and Their Stability
In this section, we study the stability of the equilibrium points of system (1).Since the right hand sides of the first three equations of model (1) do not include R(t) we will deal with the first three variables S, E and I. System (1) has two possible equilibrium points.There is always a disease-free equilibrium point E 0 = (S 0 , 0, 0) where S 0 = λ µ provided that (H1) and (H5) is satisfied.
Proof.For finding an equilibrium point, let D α t S(t) = 0, D α t E(t) = 0 and D α t I(t) = 0. Now, we need to solve the following system λ − f (S, I) − µS = 0, It is easy to see that E 0 = (S 0 , 0, 0) where S 0 = λ µ is always disease-free equilibrium point of system (1).From system (6), we have and where F (I) is defined as in Lemma 3. By using the equations ( 7) and ( 8) one can obtain S as in the following form On the other hand, we need to find a positive root for the equation f (S, I) = ηE − pT (I).Let us define the continuous function Clearly, H(0) = 0 under the hypotheses (H1) and (H5).By using Lemma 3, one can calculate a positive I 0 value such that λβ − ηF (I 0 ) + βpT (I 0 ) = 0.At this positive I 0 value, H(I 0 ) = −λ < 0 is obtained.
Moreover, if we look at the derivative of the function H(I), we can see that and which implies that there exist some I * ∈ (0, I 0 ) such that H(I * ) = 0. Also, at this positive ).For uniqueness of the endemic equilibrium point, assume that we have another positive equilibrium point Ĩ.On the other hand, one can observe H ′ (I * ) < 0 which says that at every root of the function H(I), the function is strictly decreasing.But, this will be a contradiction.This completes the proof.Next, we will prove the local and global asymptotic stability of the disease-free equilibrium point E 0 .□ Theorem 3. The disease-free equilibrium point E 0 = (S 0 , 0, 0) is locally asymptotically stable if R 0 < 1.
Proof.The disease-free equilibrium point E 0 is locally asymptotically stable if all of the eigenvalues λ i of the Jacobian matrix evaluated at the equilibrium point E 0 , that is J(E 0 ), satisfy the Matignon's conditions [30], that is, The Jacobian matrix J(E 0 ) can be evaluated as The eigenvalues of the Jacobian matrix J(E 0 ) are λ 1 = −µ and where σ = θ + µ + T ′ (0) and γ = β(f I (S 0 , 0) + pT ′ (0)).It is easy to see that all the eigenvalues are real and if ησ − γ > 0 all the eigenvalues λ i will be negative.Hence, |arg(λ i )| > απ 2 .Also, note that the inequality ησ − γ > 0 implies that R 0 < 1.Therefore, if R 0 < 1 the disease-free equilibrium point E 0 is locally asymptotically stable and if R 0 > 1 the trivial equilibrium E 0 becomes unstable.□ Theorem 4. Let lim The disease-free equilibrium point E 0 = (S 0 , 0, 0) is globally asymptotically stable if R 0 < 1.
Proof.To establish the global stability of the disease-free equilibrium point, consider the following Lyapunov function Calculating the time fractional order derivative in Caputo sense of both sides of Eq. 16 and using the hypothesis of this theorem and (H4) we get Therefore, the equilibrium point E 0 is globally asymptotically stable.This completes the proof.□ and and Proof.We consider the following Lyapunov function where Function L is defined, continuous and positive definite for all S(t) > 0, E(t) > 0 and I(t) > 0. With the help of Lemma 2, we have Using the equations in system (1), one has After some computations, we obtain Moreover, we get By Theorem hypotheses, where strict equality holds when S = S * , E = E * .Moreover, is satisfied under the assumptions of the theorem.On the other hand, for all S, I > 0, since the arithmetic mean-geometric mean inequality is satisfied.So, by Theorem 3.1 in [10], the positive endemic equilibrium point Σ * = (S * , E * , I * ) of system ( 1) is uniformly asymptotically stable.□

Numerical Simulations
In this section, we perform some numerical simulations to observe the results obtained in Section 4. We study system (1) for different values of the noninteger order derivative α.For the numerical simulations, the Adams-Bashforth-Moulton scheme [12] is used in Matlab.The following Table 1 shows the parameters that are used in system (1).The assumed values p and q have been chosen randomly.The value of r has taken as small due to the recovery rate of an exposed person should be small.The death rate is again randomly selected that one of the 5 patients died from the disease.The value of the parameter w is also considered as small and variable in our numerical examples.Incidence function parameter 0.01 [19] As a first example, we have chosen bilinear incidence function f (S, I) = wSI and the treatment function where w, r 1 and r 2 are positive parameters.
incidence function is mostly common in literature and it indicates that the rate of transmission increases as the number of infected people increases, that is, as the connection between the susceptible population and the infected population increases.On the other hand, the chosen treatment function is a monotonically increasing function for I > 0. On the next page, the readers can see another numerical example with a saturated treatment function T (I) = r 1 I 1 + r 2 I .
To obtain the value of R 0 smaller than one, first we choose the transmission coefficient w = 5 × 10 −6 and find R 0 = 0.1704 < 1.In this case, the diseasefree equilibrium point E 0 = (S 0 , 0, 0) where S 0 = 55500 is globally asymptotically stable which is depicted in Figure 2 and 3.        .The parameters γ, α 1 and α 2 can be found in Table 1.As in the paper [16], the parameters γ, α 1 and α 2 can be thought as the transmission rate of the disease, a measure of inhibition for the susceptible population, and a measure of inhibition for the infected population, respectively.In this example, using these parameters, the positive equilibrium point is evaluated as Σ * = (S * , E * , I * ) = (53221, 1278.8,3.26) with R 0 = 6.1298 > 1.The trajectories can be seen in Figure 6 and 7.     Following this second example, we also wondered the effect of the treatment success.As we expect, if the rate of successfully treated individuals q increase, we observe a decrease in the number of the exposed and infected individuals in Figure 11.
Finally, we need to mention our observation about the parameter α.If we look carefully at all the figures, we will see that the system will be stable over a longer period of time if the parameter α decreases.

Concluding Remarks
In this paper, we have introduced a fractional order SEIR model with a general incidence function f (S, I) and a general treatment function T (I).By analysing the equilibrium points of system (1), we have shown that the disease free equilibrium point is locally asymptotically stable if the basic reproduction number R 0 < 1 and globally asymptotically stable if the inequality (15) is satisfied when R 0 < 1.We have also constructed a Lyapunov function and found that the endemic   1 and changing parameters p and q and the initial conditions S(0) = 50000, E(0) = 2000, I(0) = 0 and R(0) = 0.
equilibrium point is uniformly asymptotically stable under the conditions ( 18)- (20).Moreover, in our model, based on the treatment model in [43], unsuccessfully treated infectious individuals re-enter the exposed compartment proportional to parameter p. Changing these parameters p and q, we also have a chance to compare treatment success.
To the best of the author's knowledge, a fractional SEIR model with a general incidence function and treatment function has not been studied yet.In 2017, Elkhaiar et al. [11] studied the stability analysis of an ordinary differential equation system of the SEIR model with treatment.In that paper, they used a general incidence function and a general treatment function T * (I).In our paper, if the parameters are changed as µ = d, p = 0, q = 1, r = 0, β = σ, θ = 0 and T (I) = γI + T * (I), then system (1) becomes to the SEIR model studied in [11] when α = 1.Another paper including a stability analysis of a fractional order SEIR model is published in 2020 by Yang et.al. [48].In our system, if the parameters are chosen as λ = Λ α , µ = d α , r = 0, β = σ α , θ = 0 and qT (I) = γ α I, f (S, I) = β α F (S)G(I) then system (1) turns into the system that studied in [48].In 2018, analysis of a fractional order SEIR model with treatment is established by Almeida [2].If the parameters in our system (1) is chosen as λ = bN (N is assumed as fixed population size), µ = b, r = 0, β = σ, p = 0, q = 1, θ = 0 and f (S, I) = βIS N and T (I) = (µ + q)I (µ and q are the parameters used in [2]) then again system (1) turns into that system used in [2].Finally, the last example studied in [13] is related to outbreaks of influenza A(H1N1).In [13], the authors proposed a fractional order SEIR model to explain and understand the outbreaks of ifluenza A(H1N1).They used real data values and tested and simulated these data values for their model and chose the best fitted order α of fractional differentiation.In system (1), if the parameters are changed as λ = µ α , µ = µ α , β = Ω α , p = q = θ = 0, q = 1, f (S, I) = β α SI and T (I) = ρ α I, then one can obtain the system studied in this paper.As a conclusion, the proposed model in this paper is rather general and as pointed out in [13], in the applications, one can choose the best fitted order of α depending on the real data values.If we look carefully at all the figures, we will see that the system will be stable over a longer period of time if the parameter α decreases.

Declaration of Competing Interests
The author has no competing interests to declare.
and T (I) is an increasing function and E * = F (I * ) β > 0. This guarantees the existence of a positive endemic equilibruim point Σ * = (S * , E * , I *

Figure 9 .
Figure 9. Trajectories for system (1) with parameters given in Table 1 and the initial conditions S(0) = 50000, E(0) = 2000, I(0) = 0 and R(0) = 0. Additionally, let us take the saturated treatment function, T (I) = r 1 I 1 + r 2 I in model (1) with the Beddington-DeAngelis infection function.In this case, the parameters r 1 and r 2 represent the treatment rate of disease and the limitation in treatment availability.This treatment function has a horizontal asymptote which

Table 1 .
Parameter values used in numerical simulations