domain: nlin math

text:
Rubin Fandio

Hamadjam Abboubakar
email: h.abboubakar@gmail.com

Henri Paul

Ekobena Fouda
email: hekobena@gmail.com

Anoop Kumar

Mathematical modelling and projection of Buruli ulcer transmission dynamics using classical and fractional derivatives: A case study of Cameroon

Keywords: Buruli ulcer, Fractional derivatives, Caputo derivative, Asymptotic stability, Parameter estimation, Strength number, Adams-Bashforth-Moulton scheme (ABM) 26A33, 34C60, 92C60, 92D30

The goal of this work is to derive and study a new model which traduces the transmission dynamics of the Buruli ulcer (BU) in which we replace mass action incidences with standard incidences, consider latent period, and by using both integer and fractional derivatives. We first compute the basic reproduction number denoted by R 0 and prove the asymptotic stability of the Buruli-free equilibrium whenever R 0 < 1. Then, for R 0 > 1, we prove the existence of an unique positive equilibrium point and its global stability using the general theory of Lyapunov. We perform parameter estimation to calibrate our model with real data from Cameroon, while sensitivity analysis is conducted to determine important parameters in the BU dynamics. From this model calibration, we obtain R 0 = 2.0843 which is greater than one and implies that BU ulcer is endemic in the country. To determine whether if or not the BU admits one or multiple waves, we compute the strength number denoted by A 0 . We find that A 0 ≈ 17 > 0 which means that Bu admits multiple waves. We then replace integer derivative with fractional derivative and prove the existence of equilibrium points as well as its asymptotic stability. This follows by the proof of existence and the uniqueness of solutions of the fractional model. We then construct a numerical scheme based on the Adams-Bashforth-Moulton (ABM) Method. Theoretical results are validated from numerical simulations. These last also permits to evaluate the impact of varying fractional order α in the BU dynamics.

Introduction

For several decades, mathematical modelling has been positioned as this discipline that allows us to better understand the dynamics of disease transmission and the appropriate control mechanisms that lead to successful intervention. This modelling involves mathematical models constructed for infectious diseases to study the non-linear process involved in these disease dynamics and determine the better strategy to control them. Like Examples of these mathematical models, see [START_REF] Abboubakar | Projections and fractional dynamics of the typhoid fever: A case study of mbandjock in the centre region of cameroon[END_REF][START_REF]New concept in calculus: Piecewise differential and integral operators[END_REF][START_REF] Ding | A fractional-order differential equation model of hiv infection of cd4+ t-cells[END_REF][START_REF] Singh | A fractional epidemiological model for computer viruses pertaining to a new fractional derivative[END_REF]. Most mathematical models are often constructed using classical (integer) derivatives (partial and ordinary differential equations), and in recent years, the concept of fractional derivatives which generalize classical derivatives has been considered very effective in different scientific fields. For example, modelling the typhoid fever disease with fractional derivative [START_REF] Abboubakar | Fractional dynamics of typhoid fever transmission models with mass vaccination perspectives[END_REF][START_REF] Abboubakar | Projections and fractional dynamics of the typhoid fever: A case study of mbandjock in the centre region of cameroon[END_REF], mathematical modelling of the computer virus spreading [START_REF] Singh | A fractional epidemiological model for computer viruses pertaining to a new fractional derivative[END_REF], the COVID-19 modelling [START_REF] Atangana | Modeling third waves of covid-19 spread with piecewise differential and integral operators: Turkey, spain and czechia[END_REF]. Since then, several works are appear in the literature for infectious diseases such as Malaria [START_REF] Chitnis | Bifurcation analysis of a mathematical model for malaria transmission[END_REF][START_REF] Teboh-Ewungkem | Models and proposals for malaria: a review[END_REF], Dengue [START_REF] Rodrigues | Vaccination models and optimal control strategies to dengue[END_REF], Chikungunya [START_REF] Abboubakar | Mathematical modeling and projections of a vector-borne disease with optimal control strategies: A case study of the chikungunya in chad[END_REF][START_REF] Dumont | Vector control for the chikungunya disease[END_REF][START_REF] Dumont | Mathematical studies on the sterile insect technique for the chikungunya disease and aedes albopictus[END_REF], HIV/AIDS [START_REF] Bryan | Mediational analysis in hiv/aids research: Estimating multivariate path analytic models in a structural equation modeling framework[END_REF][START_REF] Karrakchou | Optimal control and infectiology: application to an hiv/aids model[END_REF], Cholera [START_REF] Capasso | A mathematical model for the 1973 cholera epidemic in the european mediterranean region[END_REF][START_REF] Edward | A mathematical model for the dynamics of cholera with control measures[END_REF][START_REF] Sun | Transmission dynamics of cholera: Mathematical modeling and control strategies[END_REF][START_REF] Wang | Mathematical models for cholera dynamics-a review[END_REF], Tuberculosis [START_REF] Egonmwan | Analysis of a mathematical model for tuberculosis with diagnosis[END_REF][START_REF] Feng | On the role of variable latent periods in mathematical models for tuberculosis[END_REF][START_REF] Houben | The global burden of latent tuberculosis infection: a re-estimation using mathematical modelling[END_REF][START_REF] Waaler | The use of mathematical models in the study of the epidemiology of tuberculosis[END_REF], Measles [START_REF] Abboubakar | Fractional dynamics of a measles epidemic model[END_REF][START_REF] Qureshi | Modeling of measles epidemic with optimized fractional order under caputo differential operator[END_REF][START_REF] Thompson | Evolution and use of dynamic transmission models for measles and rubella risk and policy analysis[END_REF], Coronavirus (COVID-19) [START_REF] Djaoue | Mathematical modeling, analysis and numerical simulation of the covid-19 transmission with mitigation of control strategies used in cameroon[END_REF][START_REF] Khan | Mathematical modeling and analysis of covid-19: A study of new variant omicron[END_REF][START_REF] Nabi | Forecasting of covid-19 pandemic: From integer derivatives to fractional derivatives[END_REF], and Buruli ulcer [START_REF] Chu | Mathematical modeling and stability analysis of buruli ulcer in possum mammals[END_REF][START_REF] Cynthia | Modelling transmission of buruli ulcer in the central region of ghana[END_REF][START_REF] Edholm | A risk-structured mathematical model of buruli ulcer disease in ghana[END_REF][START_REF] Momoh | Modeling, optimal control of intervention strategies and cost effectiveness analysis for buruli ulcer model[END_REF][START_REF] Nyabadza | On the transmission dynamics of buruli ulcer in ghana: Insights through a mathematical model[END_REF][START_REF] Zhao | A mathematical model for the coinfection of buruli ulcer and cholera[END_REF].

Speaking of the latter, Buruli Ulcer (BU) is a disease that attacks the skin and sometimes the bones, and may also be known as necrotizing skin disease. The pathogen agent of BU is called Mycobacterium ulcerans. It is the most common mycobacterial infection after Hansen's disease (Leprosy) and tuberculosis [START_REF] Portaels | Buruli ulcer[END_REF][START_REF] Weir | Buruli ulcer: the third most common mycobacterial infection[END_REF]. Infection by Mycobacterium ulcerans causes chronic necrotizing ulcers [START_REF] Van Der Werf | Mycobacterium ulcerans disease[END_REF], implying in functional limitation, physical deformities, as well as social stigma. Long-term disabilities can be observed if the infection is not treated [START_REF] Agbenorku | Buruli-ulcer induced disability in ghana: a study at apromase in the ashanti region[END_REF][START_REF] Owusu-Sekyere | Perceptions and attitudes: The challenge of managing buruli ulcer morbidity in ghana[END_REF]. The disease is generally endemic in regions around slow-flowing or stagnant water bodies, and periodically flooded or continuously humid areas covered with a shallow water table which worldwide, has proven to be a cause of risk for infection [START_REF] Sopoh | Family relationship, water contact and occurrence of buruli ulcer in benin[END_REF][START_REF] Williamson | Distribution of mycobacterium ulcerans in buruli ulcer endemic and non-endemic aquatic sites in ghana[END_REF]. In addition, BU is age and sexes independent. The clinical evolution of BU in people co-infected with HIV is more aggressive. Indeed, complications due to the management of HIV imply a poor treatment outcome of BU, which has allowed WHO to publish a technical guide to help clinicians manage co-infection and especially children under 15 years old [START_REF]Buruli ulcer (Mycobacterium ulcerans infection)[END_REF]. However, children under the age of five are less affected to Mycobacterium ulcerans compared to older [START_REF] Bratschi | Geographic distribution, age pattern and sites of lesions in a cohort of buruli ulcer patients from the mapé basin of cameroon[END_REF][START_REF] Röltgen | Late onset of the serological response against the 18 kda small heat shock protein of mycobacterium ulcerans in children[END_REF]. Mycobacterium ulcerans is heat-sensitive, it mainly causes skin lesions. 90% of these lesions are observed on the limbs although the rest of the body can be affected [START_REF] Bratschi | Geographic distribution, age pattern and sites of lesions in a cohort of buruli ulcer patients from the mapé basin of cameroon[END_REF]. Nowadays, the treatment of BU is not limited to surgery alone [START_REF] Asiedu | Socioeconomic implications of buruli ulcer in ghana: a three-year review[END_REF], it is based on the combination of complementary treatments and antibiotics [START_REF]World Health Organisation, Treatment of Mycobacterium ulcerans disease (Buruli ulcer): Guidance for health workers[END_REF]. For health workers, a treatment guidance is available in [START_REF] Asiedu | Buruli ulcer: Mycobacterium ulcerans infection, tech. rep[END_REF]. The current recommended treatment consists of a combination of Clarithromycin and Rifampicin [START_REF]World Health Organisation, Treatment of Mycobacterium ulcerans disease (Buruli ulcer): Guidance for health workers[END_REF][START_REF]Buruli ulcer (Mycobacterium ulcerans infection)[END_REF].

Although BU disease has been reported worldwide, in several countries with high infection, 33 countries are in Africa, Asia, South American countries, and the Western Pacific, except for Australia, China, and Japan. According to the World Health Organization (WHO), only fourteen countries [START_REF] Butcher | Numerical methods for ordinary differential equations in the 20th century[END_REF] regularly report data on the disease. Until 2010, the number of suspected cases of UB reported each year worldwide was in the order of 5,000. The total number of BU cases then declined until 2016, when 1961 cases were reported, the lowest level on record. Subsequently, the number of cases increased again each year to 2,713 cases in 2018. And then in 2020, 1,258 cases were reported, compared to 2,271 in 2019 [START_REF]World Health Organization (WHO, Ulcère de Buruli (infection à Mycobacterium ulcerans)[END_REF]. This decrease was observed in 2020 following the onset of the infectious disease Coronavirus (COVID- [START_REF] Cynthia | Modelling transmission of buruli ulcer in the central region of ghana[END_REF], could be attributable to the effects of this infection on active case detection activities. In each endemic country, the disease is usually present in households, which unfortunately affects people living in these environments who are generally very poor and have serious difficulties in accessing quality medical care within communities [START_REF] Asiedu | Socioeconomic implications of buruli ulcer in ghana: a three-year review[END_REF]. Although there are some studies which prove that the pathogen can leave in others mammal reservoirs (animal), the transmission dynamics of BU remains unclear [START_REF] Johnson | Buruli ulcer (m. ulcerans infection): new insights, new hope for disease control[END_REF]. In Cameroon, for example, BU was first described in 1975. The 47 cases studied at the time all came from a very localized outbreak in the Nyong Valley between the cities of Ayos, and Akonolinga, in south-central Cameroon [START_REF] Alphonse | L'ulcère de Buruli au Cameroun[END_REF] but also from the different regions of the country such as the Far North, and the South-West. In addition, an epidemiological survey conducted by, among others, the Emmaus-Switzerland Leprosy Aid (ALES) in August 2001 in the Nyong basin and focusing on Buruli ulcer identified 438 cases of Buruli ulcer (active and inactive forms combined). Of these 438 cases, 97 were recorded in Ayos district and 331 cases in Akonolinga district. This survey has made it possible to classify Buruli ulcer as a health problem in Cameroon [START_REF] Alphonse | L'ulcère de Buruli au Cameroun[END_REF]. So far, health workers have been reporting suspected cases of Buruli ulcer for irregular time intervals without confirmation of diagnosis. These observations come from different regions of the country, but especially from the provinces of the Far North, South-West and Centre. In the latter province, the reported cases come mainly from the Ayos-Akonolinga area.

The main contribution in this work consists to the formulation of a new BU compartmental model including latent period, standard force of infection, with both classical and fractional derivatives in the Caputo sense. We begin by the formulation of the proposed model with integer derivative. We then prove the positivity and boundedness of solutions, compute equilibrium points and perform stability analysis of these equilibrium points in term of the basic reproduction number R 0 . After that, we use real data to calibrate the model by performing parameter estimation and obtain the value of the basic reproduction number. We also compute the strength number A 0 which the sign indicates the presence of more than one wave in the disease spread, and compute its numerical value. The sign of this last threshold permits to know if the epidemic has one or multiple wave (epidemic peaks) [START_REF] Atangana | Mathematical model of survival of fractional calculus, critics and their impact: How singular is our world?[END_REF]. To determine the important parameters in the disease dynamics, we conduct both local sensitivity analysis. We then formulate our BU fractional model in the Caputo sense and perform theoretical analysis like stability of the stationary points as well as existence and uniqueness of solutions. To validate our theoretical results, we perform several numerical simulations based on a numerical scheme constructed using the dams-Bashforth-Moulton method.

The rest of the work is organized as follows: Section 2 deals with the model description and basic results (positivity and boundedness of solutions). Section 3 is devoted to equilibrium points and stability analysis of these particular points. Parameter estimation, computation of the strength number as well as existence of multiple waves and sensitivity analysis are performed in Section 4. Preliminaries on fractional calculus are presented in Section 5. The BU model with Caputo derivative is formulated in Section 6 as well as its mathematical analysis. A numerical scheme, simulation results as well as discussions are done in Section 7. The work finish with a conclusion.

BU model and basics results

We present here the transmission dynamics of the model of Buruli ulcer disease using firstly classical derivative. In what follows, N h (t) denotes the total human population at any time positive t (t > 0). N h (t) is divided into fours subpopulation, also called compartments, namely S h (t) for susceptible humans, E b (t) for infected (but not yet infectious ) humans, I b (t) for infectious humans, and R b (t) for recovered humans. The total vector population

                                             dS h dt = π h + ϕR b -β h S h I v N h -µ h S h , dE b dt = β h S h I v N h - k 1 (µ h + γ 1 ) E b , dI b dt = γ 1 E b - k 2 (µ h + κ + η) I b , dR b dt = ηI b - k 3 (µ h + ϕ) R b , dS v dt = π v -β v S v I b N h -µ v S v , dI v dt = β v S v I b N h -µ v I v , (1)
with the following initial conditions

S h (0) = S h0 > 0, E b (0) = E b0 ≥ 0, I b (0) = I b0 ≥ 0, R b (0) = R b0 ≥ 0, S v (0) = S v0 ≥ 0, I v (0) = I v0 ≥ 0. ( 2
)
Remark 1. Model system (1) is a extension of the one from proposed by Zhao et al. in [START_REF] Zhao | A mathematical model for the coinfection of buruli ulcer and cholera[END_REF], in which we added a compartment for latent individuals, and force of infection (FOI) with standard incidences. By cons, we do not consider reservoir dynamics as done by Nyabadza and Bonyah in [START_REF] Nyabadza | On the transmission dynamics of buruli ulcer in ghana: Insights through a mathematical model[END_REF]. Others works concerning mathematical modeling and control of Buruli ulcer transmission dynamics can be found in [START_REF] Chu | Mathematical modeling and stability analysis of buruli ulcer in possum mammals[END_REF][START_REF] Cynthia | Modelling transmission of buruli ulcer in the central region of ghana[END_REF][START_REF] Edholm | A risk-structured mathematical model of buruli ulcer disease in ghana[END_REF][START_REF] Khan | Mathematical modeling and optimal control strategies of buruli ulcer in possum mammals[END_REF][START_REF] Momoh | Modeling, optimal control of intervention strategies and cost effectiveness analysis for buruli ulcer model[END_REF].

We first prove that state variables S h (t), E b (t), I b (t), R b (t), S v (t), I v (t), are positive for all t ≥ 0.

Theorem 1. ∀t > 0, a solution x(t) = (S h (t), E b (t), I b (t), R(t), S v (t), I v (t)) ′ of model (1) with initial conditions x(0) = (S h (0), E b (0), I b (0), R(0), S v (0), I v (0)) ′ ∈ R 6 + is positive. Proof. From (1), we have                  dS h dt (t) S h =0,R b ≥0 = π h + ϕR b (t > 0, dE b dt (t) E b =0,S h ,I b ,R b ,Iv≥0 = β h S h (t)I v (t) S h + I b + R b ≥ 0, dI b dt (t) I b =0,E b ≥0 = γ 1 E b (t) ≥ 0, dR b dt (t) R b =0,I b ≥0 = ηI b (t) ≥ 0, dS v dt (t) Sv=0,N h ≥0 = π v > 0, dI v dt (t) Iv=0,N h ,Sv≥0 = β v S v (t)I b (t) N h ≥ 0.
(3) So, the non-negativity of each state variable of system (1) come from the barrier theorem [START_REF] Gauthier | CIMPA lecture notes on Mathematical Epidemiology[END_REF].

Adding the first four equations of model [START_REF] Abboubakar | Fractional dynamics of a measles epidemic model[END_REF], we obtain

dS h (t) dt + dE b (t) dt + dI b (t) dt + dR b (t) dt = π h -µ h N h (S h + E b + I b + R b ) -kI b , ≤ π h -µ h N h (t)
Solving this inequality gives

N (t) ≤ π h µ h + N h (0) - π h µ h exp(-µ h t), for all t ≥ 0, with N h (0) = S h (0) + E b (0) + I b (0) + R b (0) > 0. When t -→ ∞, we get lim sup t→∞ N h (t) ≤ π h µ h .
Adding the last two equations of (1) gives

dN v (t) dt = dS v (t) dt + dI v (t) dt = π v -µ v .
Solving the above equality gives

N v (t) = π v µ v + N v (0) - π v µ v exp(-µ v t), for all t ≥ 0, When t -→ ∞, we get lim sup t→∞ N v (t) = π v µ v
. This prove that all solutions of the model [START_REF] Abboubakar | Fractional dynamics of a measles epidemic model[END_REF] lies in the following absorbing set

Ψ = (S h , E b , I b , R b , S v , I v ) ∈ R 4 + × R 2 + : N h (t) ≤ π h µ h and N v (t) = π v µ v , ∀t ≥ 0 , (4)
in which Buruli ulcer model ( 1) define a dynamical system.

Equilibrium points and asymptotic stability

Before determining the equilibrium points of the Buruli ulcer model (1), we first define the following threshold, called the "basic reproduction number," which drives the qualitative dynamics of the model.

R 0 = β h µ h π v π h µ v γ 1 (µ h + γ 1 ) 1 (κ + µ h + η) R vh β v 1 µ v R hv = R vh R hv . (5)
For model ( 1), the following result holds.

Theorem 2. Let us define the following threshold:

R 2 c := 2 1 - κγ 1 k 3 κγ 1 k 3 + µ h [(k 2 + γ 1 )ϕ + µ h (k 1 + κ) + ηk 1 ] + µ h β v k 3 γ 1 µ v {[µ h k 2 + (µ h + κ)γ 1 ] ϕ + µ h k 1 k 2 } . (6) 1. If R 0 > 1 or (R 0 = 1 and R c < 1)
, then Buruli ulcer model (1) admits two feasible equilibrium points: the BU-free equilibrium M 0 and the endemic equilibrium M 1 ; 2. If R c < R 0 < 1. then Buruli ulcer model (1) has two positive equilibrium points in addition with the BU-free equilibrium; 3. No endemic equilibrium point otherwise.

Proof. Let E b = (S * h , E * b , I * b , R * b , S * v , I * v ) ′
any arbitrary steady states of model [START_REF] Abboubakar | Fractional dynamics of a measles epidemic model[END_REF]. Setting the right-hand side of Eq. ( 1) to zero with

X * = I * v N * h , Y * = I * b N * h and N * h = π h -κI * b µ h , we obtain                π h + ϕR * b -(β h X * + µ h ) S * h = 0, β h X * S * h -k 1 E * b = 0, γ 1 E * b -k 2 I * b = 0, ηI * b -k 3 R * b = 0, π v -β v Y * S * v -µ v S * v = 0, β v S * v Y * -µ v I * v = 0. (7)
Resolution of the above system gives

                                       R * b = ηγ 1 β h X * π h [k 1 k 2 k 3 (β h X * + µ h ) -ηγ 1 β h X * ϕ] , S * h = π h + ϕR * b (β h X * + µ h ) , E * b = β h X * S * h k 1 , I * b = γ 1 E * b k 2 , S * v = π v β v Y * + µ v , I * v = β v S * v Y * µ v . (8)
Using [START_REF] Atangana | Mathematical model of survival of fractional calculus, critics and their impact: How singular is our world?[END_REF] in the expressions of X * and Y * gives

X * = πv µv βvY * βvY * + µv µ h k 1 k 2 (β h X * + µ h ) [k 1 k 2 k 3 (β h X * + µ h ) -ηγ 1 β h X * ϕ] π h k 1 k 2 (β h X * + µ h ) [k 1 k 2 k 3 (β h X * + µ h ) -ηγ 1 β h X * ϕ] -κγ 1 β h X * (π h [k 1 k 2 k 3 (β h X * + µ h ) -ηγ 1 β h X * ϕ] + ϕηγ 1 β h X * π h )
and

Y * = γ 1 (X * ) 2 ηβ h π h φ + γ 1 k 3 (X * ) 2 b h π h κ + -k 1 k 2 k 3 X * µ h -k 1 k 2 k 3 (X * ) 2 β h π h µ 2 v γ 1 X * ηβ h µ h φ -k 1 k 2 k 3 µ 2 h -k 1 k 2 k 3 X * β h µ h βvπv + (-γ 1 (X * ) 2 ηβ h π h φ -γ 1 k 3 (X * ) 2 β h π h κ + (k 1 k 2 k 3 X * µ h + k 1 k 2 k 3 (X * ) 2 b h ) π h ) βvµv .
By substituting the expression of Y * in the expression of X * , it thus follows that X * is a nonnegative solution to the following equation

X * a 2 (X * ) 2 + a 1 X * + a 0 = 0, ( 9
)
where

a 2 = R 4 0 π 3 h µ 5 v (k 2 k 3 + γ 1 (k 3 + η)) (µ v (k 2 k 3 + γ 1 (k 3 + η)) + γ 1 β v k 3 ) > 0, a 1 = R 2 0 γ 1 k 3 µ 2 h π 2 h β v µ 4 v [(µ h k 2 + µ h γ 1 + γ 1 κ) ϕ + µ h k 1 k 2 ] π v (R 2 c -R 2 0
), and

a 0 = (1 -R 0 ) (R 0 + 1) γ 2 1 µ 4 h π h (ϕ + µ h ) 2 β 2 v µ 2 v π 2 v .
Descartes' rule of signs is used to obtain each item of the theorem 2.

Stability of the BU-free equilibrium

The trivial solution of equation ( 9) is given by X * = 0, which corresponds to the BU-free

equilibrium M = π h µ h , 0, 0, 0, π v µ v , 0 ′ .
The Jacobian matrix of system (1) evaluated at M is given by:

W =            -µ h 0 0 ϕ 0 -β h 0 -k 1 0 0 0 β h 0 γ 1 -k 2 0 0 0 0 0 η -k 3 0 0 0 0 -β v π v µ h µ v π h 0 -µ v 0 0 0 β v π v µ h µ v π h 0 0 -µ v            (10)
The eigenvalues of W are w 1 = -µ h , w 2 = -k 3 , w 3 = -µ v , and those of the following sub-matrix

W 1 =    -k 1 0 β h γ 1 -k 2 0 0 β v π v µ h µ v π h -µ v    . (11)
The characteristic polynomial of W 1 is given by

P(x) = x 3 + a 2 (µ v + k 2 + k 1 ) x 2 + a 1 ((k 2 + k 1 ) µ v + k 1 k 2 ) x + a 0 1 -R 2 0 k 1 k 2 µ v . ( 12
)
Coefficients a 2 and a 1 are always positive, and a 0 > 0 ⇐⇒ R 0 < 1. It follows that all root of P have negative real part. Thus, we conclude that the condition R 0 < 1 ensures the negativity of all eigenvalues of the matrix W. We resume the above analysis as follows:

Lemma 1. The BU-free equilibrium M 0 = π h µ h , 0, 0, 0, π v µ v , 0 ′ is LAS whenever R 0 < 1. Theorem 3. If R 0 < 1, then M 0 = π h µ h , 0, 0, 0, π v µ v , 0 ′ is GAS in Ψ.
Proof. By considering the equations of the model (1) reflecting the dynamics of the infected populations, the obtained system can be rewritten as follows:

   Ėb (t) İb (t) İv (t)    = W 1   E b (t) I b (t) I v (t)   -K (S h , E b , I b , R b , S v , I v ) (13)
where W 1 is given at Eq. ( 11), and

K (S h , E b , I b , R b , S v , I v ) =       β h 1 - S h N h 0 β v I b S 0 v N 0 h - S v N h       . ( 14
)
It then follows that if

S 0 v N 0 h - S v N h ≥ 0, then K (S h , E b , I b , R b , S v , I v ) ≥ 0 R 3 . We thus have    Ėb (t) İb (t) İv (t)    ≤ W 1   E b (t) I b (t) I v (t)   . (15)
From Lemma 1, we obtained that W 1 has eigenvalues with real part. Then if R 0 < 1, BU model ( 1) is stable. We thus have (E b , I b , I v ) -→ (0, 0, 0) as t -→ ∞. Note that

W 1 = F    0 0 β h 0 0 0 0 β v π v µ h µ v π h 0    - V   k 1 0 0 -γ 1 k 2 0 0 0 µ v  
where matrices F and V are the matrices defined in [START_REF] Shuai | Global stability of infectious disease models using Lyapunov functions[END_REF][START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF] to compute the next generation matrices (NGM). Since F ≥ O R 3 * 3 and

V -1 =    1 k 1 0 0 γ 1 k 1 k 2 1 k 2 0 0 0 1 µv    ≥ O R 3 * 3 , it follows from [55, Theorem 2.1]
, that there exists a Lyapunov function for system (1) expressed as

L (S h , E b , I b , R b , S v , I v ) = w ′ V -1 (E b , I b , I v ) ′
where w ′ is the left eigenvector of V -1 Z corresponding to the eigenvalue R 0 . Thus,

dL dt = (R 0 -1) w ′ (E b , I b , I v ) -w ′ V -1 K (S h , E b , I b , R b , S v , I v ) ≤ 0. Since K (S h , E b , I b , R b , S v , I v ) ≥ 0 R 6 , then dL dt < 0 if R 0 < 1, with dL dt = 0 if and only if (E b , I b , I v ) = 0 R 3 . Thus, {M 0 } is the largest invariant set contained in (S h , E b , I b , R b , S v , I v ) ∈ R 6 + :
dL dt = 0 . From LaSalle Invariance Principle [START_REF] La Salle | The stability of dynamical systems[END_REF], every solution of (1) with initial conditions in Ψ converge to

M 0 when t -→ +∞. That is (E b , I b , I v ) -→ 0 R 3 , S -→ S 0 h , S v -→ S 0 v when t -→ +∞, which is equivalent to (S h , E b , I b , R b , S v , I v ) -→ M 0 when t -→ +∞. Thus, if R 0 < 1, then M 0 is GAS in Ψ.

Stability of the endemic equilibrium point

From item 1 of Theorem 2, model (1) admits a unique endemic equilibrium whenever R 0 > 1. The global stability of this equilibrium is given as follows.

Theorem 4. If R 0 > 1, then the unique positive equilibrium point of model (1), M 1 = (S * h , E * b , I * b , R * b , S * v , I * v ) ′
where the component is given by (8) with X * ̸ = 0 and Y * ̸ = 0, is globally asymptotically stable provided that

1 - S v I * v S * v I v 1 - Y Y * ≥ 0, 1 - S h S * h E * b E b 1 - X X * ≥ 0, S * h S h I * b I b R b R * b ≤ 1. ( 16
)
Proof. See Appendix A.

Model calibration, multiples waves and sensitivity analysis

Parameter estimations

Here, we consider data from Cameroon during the period between 2001 to 2014 [START_REF] Tabah | Buruli ulcer in cameroon: the development and impact of the national control programme[END_REF] to calibrate our BU model [START_REF] Abboubakar | Fractional dynamics of a measles epidemic model[END_REF]. One parameter, namely µ h is estimated while others are fitted using real data. The beginning date 01/01/2001 and end date 31/12/2014 correspond to t = 0 and t = 14, respectively. The parameter estimation is performed using the routine lsqcurvefit of MATLAB software [START_REF] The Mathworks | MATLAB version 9[END_REF]. Calibrate BU model ( 1) is equivalent to solving the optimization problem min

Σ ∥ I predict b -I data b ∥, (17)
where Σ = (π h , π v , ϕ, µ v , β h , β v , κ, η, γ 1 ). The results are displayed in Table 2 and Figure 1 (panel (a)). From panel (b) of Figure 1, We can be led to conclude that, it is possible that the cumulative number of cases will be decreased after (t ≈ 16). In what follows we will perform another analysis which could confirm whatever if or not the BU epidemic model (1) has only a single wave or multiple waves of epidemic peaks.

In what follows, we use the package cftool of MATLAB software [START_REF] The Mathworks | MATLAB version 9[END_REF] to predict the number of Buruli ulcer in Cameroon from 2015 to 2030. Before do this, we begin by choose the function which better fit the reported data. We obtain that the Gaussian function of order 4 given by

f (x) = 4 i=1 a i exp - x -b i c i 2 ( 18
)
where x is the year, and a 1 = 758.

Strength number and epidemic waves

The concept of "basic reproduction number " is one of most important concept in the study of epidemiological models. Indeed, it drives the qualitative behaviour of epidemiological systems. From [START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF], one must decompose the model into two component V and F , which represent the matrices of the transfer terms and infectious terms, respectively, evaluated at the disease-free equilibrium. The biological threshold R 0 is equal of the spectral radius of the matrix F V -1 , i.e. the maximal solution of

(F.V -1 -λI 3 * 3 ) = 0, ( 19
)
where I 3 * 3 is a 3 * 3 identity matrix.

To obtain the strength number [START_REF] Atangana | Mathematical model of survival of fractional calculus, critics and their impact: How singular is our world?[END_REF] of BU model (1), we need to compute the following quantity

∂ 2 ∂I b 2 β v S v I b N h = - β v S v (N h ) 2 - β v S v N h -2β v S v I b (N h ) 3 . ( 20
)
At the BU-free equilibrium point M 0 , we have

∂ 2 ∂I b 2 β v S v I b N h M 0 = -2 β v µ 2 h π v µ v π 2 h . ( 21
)
We have

F A =     0 0 β h 0 0 0 0 -2 β v π v µ 2 h π 2 h µ v 0     and V =   µ h + γ 1 0 0 -γ 1 µ h + k + η 0 0 0 µ v   .

By solving det((F

A .V -1 ) -λI 3 * 3 ) = 0, (22)
we thus obtain the strength number which is the spectral radius of F A V -1 A [START_REF] Atangana | Mathematical model of survival of fractional calculus, critics and their impact: How singular is our world?[END_REF]:

A 0 = ρ(F A V -1 A ) = 2γ 1 µ 2 h β v π 2 v π 2 h µ 2 v > 0. ( 23
)
Using parameter values of Table 2, we obtain A 0 ≈ 17. Since A 0 > 0, it follows that the Bu model can have multiple waves [START_REF] Atangana | Mathematical model of survival of fractional calculus, critics and their impact: How singular is our world?[END_REF]. This is depicted in figures 4 where it is clear that the BU will have other waves (at least 2).

Sensitivity analysis

To determine important parameter in the BU dynamics, we perform here a local sensitivity analysis (LSA). To this aim, we compute for each model parameter, its sensitivity index, which is obtained by applying the following formula [START_REF] Chitnis | Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model[END_REF]

ϖ z = z R 0 ∂R 0 ∂z , (24)
where z represents each model parameter (z

∈ {π h , π v , ϕ, µ h , µ v , κ, β h , β v , η, γ 1 }).
Let us evaluate the sensitivity of R 0 to each parameters. We thus obtain:

ϖ β h = 1 2 > 0, ϖ βv = 1 2 > 0, ϖ πv = 1 2 > 0, ϖ µv = -1 < 0, ϖ π h = - 1 2 < 0, ϖ k = - κ 2(µ h + κ + η) < 0, ϖ η = - κ 2(µ h + κ + η) < 0, ϖ γ 1 = µ h 2(µ + γ 1 )
> 0,

ϖ µ h = -µ 2 h + γ 1 (κ + η) 2 , ϖ ϕ = 0 (25)
Hence, R 0 increasing with the increase of the following parameters β h , β v , π v , and γ 1 . It decreasing with the increase of µ v , π h , κ and η. The sign of ϖ µ h depends of parameter values.

Using parameter values in Table 2, the sensitivity indices of each model parameters are depicted in Table 3 and Figure 5. Note that increase µ v of 10% permits to decrease the value of R 0 to 10%, while increase the value of β h of 10% (respectively β v and π v ) permits to increase the value of R 0 to 5%. The sensitivity index of ϕ is equal to zero since R 0 does not depend of this model parameter.

a) R 0 = f (β h , β v ) (b) R 0 = f (µ v , β h ) (c) R 0 = f (µ v , β v ) (d) R 0 = f (π v , β h ) (e) R 0 = f (π v , β v ) (f) R 0 = f (π h , β h ) (g) R 0 = f (π h , β v ) (h) R 0 = f (π h , π v )

Preliminaries on fractional calculus

The beauty of fractional calculus is its ability to accurately capture the exact behaviour of many complex models in science, engineering, and finance [START_REF] Atangana | Mathematical model of survival of fractional calculus, critics and their impact: How singular is our world?[END_REF][START_REF] Machado | Recent history of fractional calculus[END_REF][START_REF] Richard | New analytical modelling of fractional generalized kuramoto-sivashinky equation via atangana-baleanu operator and j-transform method[END_REF]. Before formulate the fractional model of Buruli ulcer, we introduce here some useful definitions and results.

Definition 1 ([49]). Let ν ∈ R * + , f ∈ C ([a, b]
) and a < t < b. The Riemann-Liouville fractional integral of f of order ν is defined by:

I ν a f (t) := 1 Γ(ν) t a f (τ )(t -τ ) ν-1 dτ, ( 26
)
where Γ is the gamma function defined as follows: Γ(τ ) := ∞ 0 σ τ -1 exp(-σ)dσ for τ > 0.

Definition 2 ([49]

). Given a function f : [a, b] → R of class C m , the Caputo fractional derivative of f of order ν is defined by :

C t 0 D ν t f (t) = 1 Γ(m -ν) t t 0 (t -θ) m-ν-1 f (m) (θ)dθ (27) if ν / ∈ N, m = [ν] + 1.

Definition 3 ([46]

). The generalized Caputo-type fractional derivative of order ν is defined as follows

C D ν,ρ a f (s) (t) = ρ ν-m+1 Γ(m -ν) ξ a (s ρ -τ ρ ) m-ν-1 τ 1-ρ d dτ m f (τ )dτ, s > a, (28)
where ρ > 0, a ≥ 0, and m -

1 < ν ≤ m Proposition 1 (Linearity). Let f, g : [a, b] → R be such that C t 0 D ν t f (t) and C t 0 D ν t g(t) exist almost everywhere and let ξ 1 , ξ 2 ∈ R. Then C t 0 D ν t (c 1 f (t) + c 2 g(t)) exists almost everywhere with C t 0 D ν t (ξ 1 f (t) + ξ 2 g(t)) = ξ 1 × C t 0 D ν t f (t) + ξ 2 × C t 0 D ν t g(t). Lemma 2. Suppose that f (t) ∈ C([a, b]) and C t 0 D ν t f (t) ∈ C([a, b]) for 0 < ν ≤ 1. Then f (t) = f (x) + 1 Γ(ν) C t 0 D ν t f (τ )(t -t 0 ) ν for 0 ≤ τ ≤ t; ∀t ∈ [a, b].

The fractional model and its analysis

The proposed BU fractional model in the Caputo sense is given as follows:

                               χ α-1 × C D α t S h = π h + ϕR b -β h S h I v N h -µ h S h , χ α-1 × C D α t E b = β h S h I v N h -(µ h + γ 1 ) E b , χ α-1 × C D α t I b = γ 1 E b -(µ h + κ + η) I b , χ α-1 × C D α t R b = ηI b -(µ h + ϕ) R b , χ α-1 × C D α t S v = π v -β v S v I b N h -µ v S v , χ α-1 × C D α t I v = β v S v I b N h -µ v I v . (29)
where C D α t denotes the Caputo fractional derivative of order α, and χ is use to balance the units [START_REF] Ullah | An efficient numerical technique for a new fractional tuberculosis model with nonsingular derivative operator[END_REF].

Without lost of generalities, the fractional model ( 29) define a dynamical system in the following absorbing set

Ψ α = (S h , E b , I b , R b , S v , I v ) ∈ R 4 + × R 2 + : N h (t) ≤ π α h µ α h and N v (t) = π α v µ α v , ∀t ≥ 0 , (30)
in which solutions of (29) are bounded.

In what follows, we set

k α 1 = µ α h + γ α 1 , k α 2 = µ α h + κ α + η α , k α 3 = µ α h + ϕ α .

Equilibrium points and stability analysis

For the fractional model ( 29), we define the fractional basic reproductive number

R 0 := γ α 1 µ α h β α h β α v π α v (µ α h + γ α 1 ) (κ α + µ α h + η α ) π α h (µ α v ) 2 . (31)
As for the model (1), BU model ( 29) admits always a BU-free equilibrium (DFE), M = (S 0 , 0, 0, 0, S 0 v , 0)

′ where S 0 h = π α h µ α h and S 0 v = π α v µ α v .
The following result, obtained as a similar manner as which of Theorem 2, holds:

Theorem 5. Let us define the following threshold:

R 2 c : = 2 1 - κ α γ α 1 (µ α h + ϕ α ) κ α γ α 1 (µ α h + ϕ α ) + µ α h [(k α 2 + γ α 1 )ϕ α + µ α h (k α 1 + κ α ) + η α k α 1 ] + µ α h β α v k α 3 γ α 1 µ α v {[µ α h k α 2 + (µ α h + κ α )γ α 1 ] ϕ α + µ α h k α 1 k α 2 } , (32)
1. If R 0 > 1 or (R 0 = 1 and R c < 1), then Buruli ulcer model (29) admits two feasible equilibrium points: the BU-free equilibrium M 0 and the endemic equilibrium M 1 ; 2. If R c < R 0 < 1. then Buruli ulcer model [START_REF] Feng | On the role of variable latent periods in mathematical models for tuberculosis[END_REF] has two positive equilibrium points in addition with the BU-free equilibrium; 3. No endemic equilibrium point otherwise.

We thus claim what follows (we omit the proof which is similar to the one from Theorem 1):

Lemma 3. If R 0 < 1, then the BU-free equilibrium M of the fractional model (29) is LAS in Ψ α .
We also claim the following result: Theorem 6. The BU-free equilibrium M 0 of the fractional model (29) is GAS in Ψ α whenever R 0 < 1.

Proof. The proof of Theorem 6 follows the proof of Theorem 3. It is sufficient to consider the following Lyapunov function

L (S h , E b , I b , R b , S v , I v ) = w ′ V -1 (E b , I b , I v ) ′
where w ′ is the left eigenvector of V -1 Z corresponding to the eigenvalue R 0 , and applying the LaSalle Invariance Principle [START_REF] La Salle | The stability of dynamical systems[END_REF], to conclude that (E b , I b , I

v ) -→ 0 R 3 , S -→ S 0 h , S v -→ S 0 v when t -→ +∞, i.e. (S h , E b , I b , R b , S v , I v ) -→ M 0 = (S 0 h , 0, 0, 0, S 0 v , 0) ′ when t -→ +∞. Thus, if R 0 < 1, then the BU-free equilibrium M 0 is GAS in Ψ α .
From [START_REF] Li | Stability analysis of a fractionalorder linear system described by the caputo-fabrizio derivative[END_REF], it follows that the fractional model [START_REF] Feng | On the role of variable latent periods in mathematical models for tuberculosis[END_REF], as the ODE model ( 1), has a unique endemic equilibrium. Thank to Boukhouima et al. [START_REF] Boukhouima | Lyapunov functions for fractional-order systems in biology: Methods and applications[END_REF]Corollary 2], the Lyapunov function Z given by Eq. ( 44) can also use to prove the GAS of the unique endemic equilibrium of the the fractional model [START_REF] Abboubakar | Fractional dynamics of a measles epidemic model[END_REF]. Thus the following result holds.

Theorem 7. If R 0 > 1, then M 1 of the fractional model (29) is GAS in Ψ α whenever 1 - S v I * v S * v I v 1 - Y Y * ≥ 0, 1 - S h S * h E * b E b 1 - X X * ≥ 0, S * h S h I * b I b R b R * b ≤ 1. (33)

Existence and uniqueness of solutions of the BU fractional model

Let us rewrite the system (29) in the following form:

C D α t x(t) = d(t, x(t)), x l (0) = x l 0 , l = 1, 2..., 6, (34)
where 0

< α ≤ 1, d(t, x) = (d 1 , d 2 , d 3 , d 4 , d 5 , d 6 ) ′ , x = (S h , E b , I b , R b , S v , I v ) ′ and x(0) = (S h (0), E b (0), I b (0), R b (0), S v (0), I v (0)) ′ .
Here, d l , l = 1, 2, ...6 are right-hand side of system (29), i.e.,

d 1 = π α h + ϕ α R b -β α h S h I v N h -µ α h S h , etc.
The function d(t, x) : R × R 6 → R 6 defines a vector field.

Let us consider the well-defined function Z(t) : R + → R + differentiable on (0, t max ) such as:

Z(t) = 6 i=1 d i (t). ( 35
)
We obtain

Z(t) = π α h + ϕ α R b -β α h S h I v N h -µ α h S h + β α h S h I v N h -(µ α h + γ α 1 ) E b + γ α 1 E b -(µ α h + κ α + η α ) I b + η α I b -(µ α h + ϕ α ) R b + π α v -β α v S v I b N h -µ α v S v + β α v S v I b N h -µ α v I v =⇒ Z(t) = (π α h + π α v ) -µ α h S h + µ α h E b -(µ α h + κ α ) I b -µ α h R b -µ α v S v -µ α v I v =⇒ Z(t) ≤ (π α h + π α v ) -a 0 I b -(µ α h S h + µ α h E b + µ α h I b + µ α h R b + µ α v S v + µ α v I v ) , =⇒ Z(t) ≤ (π α h + π α v ) -a 1 (S h + E b + I b + R b + S v + I v ) , ( 36
)
where a 0 = κ α and a 1 = min {µ α h , µ α v }. Clearly, a 0 > 0. From [START_REF] La Salle | The stability of dynamical systems[END_REF] we have

d(t) ≤ (π α h + π α v ) -a 1 x(t), ∀t ∈ (0, t max ). ( 37
) Setting J = [t 0 -ϵ, t 0 + ϵ], B = {x ∈ R 6 : ∥ x -x 0 ∥≤ ζ} and K = {(t, x) ∈ R × R 6 : t ∈ R, x ∈ R 6 }
with ϵ > 0 and ζ > 0, and considering that d : K -→ R 6 satisfies all conditions of [39, Theorem 2.1], we thus claim what follows:

Theorem 8. Assume that there exists

Q(t) ∈ L 6 (J) such that ∥ d(t, x) -d(t, y) ∥≤ Q(t) ∥ x -y ∥, ( 38
)
for t ∈ J and g, z ∈ B. The IVP (34) admits a unique solution on [t 0 -ϵ, t 0 + ϵ] with ϵ > 0.

Proof. Note that d l , l = 1, 2, ..., 6, are continuous on B, measurable on J and bounded ∀t ∈ [t 0 -ϵ, t 0 + ϵ]. From Eq. ( 37), it follows that d(t, x) satisfies all the conditions of [39, Theorem 2.1] with m(t) = (π α h + π α v ) -a 1 x(t). Thus, there exists a solution of the fractional model ( 29) in (0, t max ).

Using Eq. ( 37), we have

∥ d(t, x) -d(t, y) ∥≤ Q(t) π α h + π α v a 1 ∥ x -y ∥, (39)
with a 1 = min {µ α h , µ α v }. This ends the proof.

Numerical scheme

In general, exact solutions of fractional disease models are not always available. Hence an iterative solution of model ( 29) is proposed using the Adams-Bashforth technique [START_REF] Nabi | Forecasting of covid-19 pandemic: From integer derivatives to fractional derivatives[END_REF][START_REF] Owolabi | Analysis and application of new fractional adams-bashforth scheme with caputo-fabrizio derivative[END_REF]. This method is generally based on the discretization of the independent variable and includes modifications of the integer order. The advantage of the Adams-Bashforth method is that it uses only one additional function evaluation per step yet achieves high-order accuracy [START_REF] Butcher | Numerical methods for ordinary differential equations in the 20th century[END_REF]. By considering the following fractional differential equation [START_REF] Li | On the fractional adams method[END_REF] C D α t h(t) = g(t, q(t)),

0 ≤ t ≤ T q k (0) = q k 0 , k = 0, 1, 2..., m -1, where m = [α], (40)
which is equivalent to the following Volterra integral equation:

q(t) = m-1 . k=0 q k t k k! + 1 Γ (α) t 0 (t -φ) α-1 q(φ, g(φ))dφ, (41)
and using the Diethelm and Neville method [START_REF] Diethelm | Analysis of fractional differential equations[END_REF] based on the well know Adams-Bashforth-Moulton algorithm [START_REF] Diethelm | An algorithm for the numerical solution of differential equations of fractional order[END_REF], we obtain the following scheme of the BU fractional model ( 29):

Setting α ∈ (0, 1], 0 ≤ t ≤ T , q = T N , t n = n.q, n = 0, 1, 2..., N ∈ N, and S h = S, E b = E, I b = I, R b = R, S v = V , I v = W
as state variables, the solution of the fractional model ( 29) is given by

                                                                         S n+1 = S 0 + q α Γ (α + 2) (π α h + ϕ α R p n+1 -β α h S p n+1 W p n+1 N h -µ α h S p n+1 ) + q ι Γ (α + 2) n . j=0 a j,n+1 (π α h + ϕ α R j -β α h S j W j N h -µ α h S j ), E n+1 = E 0 + q α Γ (α + 2) (β α h S p n+1 W p n+1 N h -(µ α h + γ α 1 ) E p n+1 ) + q ι Γ (α + 2) n . j=0 a j,n+1 (β α h S j W j N h -(µ α h + γ α 1 ) E j ), I n+1 = I 0 + q α Γ (α + 2) (γ α 1 E p n+1 -(µ α h + κ α + η α ) I p n+1 ) + q ι Γ (α + 2) n . j=0 a j,n+1 (γ α 1 E j -(µ α h + κ α + η α ) I j ), R n+1 = R 0 + q α Γ (α + 2) (η α I p n+1 -(µ α h + ϕ α ) R p n+1 ) + q ι Γ (α + 2) n . j=0 a j,n+1 (η α I j -(µ α h + ϕ α ) R j ), V n+1 = V 0 + q α Γ (α + 2) (π α v -β α v V p n+1 I p n+1 N h -µ α v V p n+1 ) + q ι Γ (α + 2) n . j=0 a j,n+1 (π α v -β α v V j I j N h -µ α v V j ), W n+1 = W 0 + q α Γ (α + 2) (β α v V p n+1 I p n+1 N h -µ α v W p n+1 ) + q ι Γ (α + 2) n . j=0 a j,n+1 (β α v V j I j N h -µ α v W j ), (42)
where

S p n+1 = S 0 + 1 Γ (ι) n . j=0 b j,n+1 (π α h + ϕ α R j -β α h S j W j N h -µ α h S j ), E p n+1 = E 0 + 1 Γ (ι) n . j=0 b j,n+1 (β α h S j W j N h -(µ α h + γ α 1 ) E j ), I p n+1 = I 0 + 1 Γ (ι) n . j=0 b j,n+1 (γ α 1 E j -(µ α h + κ α + η α ) I j ), R p n+1 = R 0 + 1 Γ (ι) n . j=0 b j,n+1 (η α I j -(µ α h + ϕ α ) R j ), V p n+1 = V 0 + 1 Γ (ι) n . j=0 b j,n+1 (π α v -β α v V j I j N h -µ α v V j ), W p n+1 = W 0 + 1 Γ (ι) n . j=0 b j,n+1 (β α v V j I j N h -µ α v W j ). (43)
and

a j,1+n =      n 1+p -(-ι + n)(1 + n), if j = 0, (2 + n -j) 1+ι -2(n + 1 -j) 1+ι + (-j + n) 1+ι , if 1 ≤ j ≤ n, 1, if j = 1 + n. b j,1+n = q ι ι ((1 + n -j) ι -(-j + n) ι ), 0 ≤ j ≤ n.

Numerical results and discussions

The ODE model

Let us consider the parameter values as listed in Table 2. We begin by illustrating the result of Lemma 1 (resp. Lemme 3) and Theorem 3 (resp. Theorem 4). Note that for estimated parameters, we have 1 < R 0 = 2.0843 > R c = 0.8795. The numerical values of coefficients of Eq. ( 9) are a 2 = 261.9718984431078 > 0, a 1 = -36.51673116398468 < 0, and a 0 = -2.414751617064077 < 0. It then follows that Eq. ( 9) admits X * = 0 and the unique positive solution X * = 0.1883345. Thus, model (1) admits as equilibrium points M 0 = (2090, 0, 0, 0, 69971, 0) ′ and M 1 = (888, 262, 128, 521, 292620, 339) > 0 R 6 + . Now, setting β h = 0.1 such that R 0 = 0.744080618507066 < R c = 0.8795695327518487, the coefficients of Eq. ( 9) are a 2 = 33.38675690832024 > 0, a 1 = 0.2867204474332617 > 0, and a 0 = 0.3222820573018668 > 0. It then follows that the unique nonnegative solution of Eq. ( 9) is X * = 0 which correspond to the disease-free equilibrium M 0 = (2090, 0, 0, 0, 69971, 0) ′ . This validate item 3 of Theorem 2.

Setting β h = 0.18 such that R c = 0.8795695327518487 < R 0 = 0.9982889062331324 < 1, the coefficients of Eq. ( 9) are a 2 = 60.09616243497641 > 0, a 1 = -0.5230213257082926 < 0, and a 0 = 0.002468871467002343 > 0. The discriminant of Eq. ( 9) is given by ∆ = a 2 1 -4a 2 a 0 = -319.9275 × 10 -3 < 0. It then follows that the unique nonnegative solution of Eq. ( 9) is X * = 0 which correspond to the disease-free equilibrium M 0 = (2090, 0, 0, 0, 69971, 0) ′ . This validate the fact that backward bifurcation does not occurs in model [START_REF] Abboubakar | Fractional dynamics of a measles epidemic model[END_REF].

Setting β v = 0.1806175784054825 such that R 0 = 1, the coefficients of Eq. ( 9) are a 2 = 60.30235183593316 > 0, a 1 = -0.5328650411460003 < 0, and a 0 = 0. It then follows that the nonnegative solutions of Eq. ( 9) are X * = 0 and X * = 0.0088366 which correspond to the disease-free equilibrium M 0 = (2090, 0, 0, 0, 69971, 0) ′ and the unique endemic equilibrium M 1 = (3164, [START_REF]New concept in calculus: Piecewise differential and integral operators[END_REF][START_REF] Abboubakar | Projections and fractional dynamics of the typhoid fever: A case study of mbandjock in the centre region of cameroon[END_REF][START_REF] Diethelm | An algorithm for the numerical solution of differential equations of fractional order[END_REF]292931,[START_REF] Egonmwan | Analysis of a mathematical model for tuberculosis with diagnosis[END_REF], respectively. This validate the fact that R 0 ≥ 1 implies that the DFE M 0 becomes unstable, and the emergence of a unique endemic equilibrium M 1 which is GAS.

The above analysis are displayed in figures 7-8, which display the time series of model ( 2) with different initial conditions. On figure 7 trajectories of infected states tend to zero whenever R 0 < 1, while for R 0 > 1, trajectories of infected states tend to the endemic equilibrium [START_REF] Atangana | Mathematical model of survival of fractional calculus, critics and their impact: How singular is our world?[END_REF].

The fractional model

Figure 9 shows the influence of varying the fractional-order parameter α between 0.6 and 1 on the model dynamics. The blue curve in each of these figures represents the numerical results of model ( 29) when the fractional order is equal to 1. From the results of Figure 9, it follows that the variation of the fractional parameter has a great impact on the quantitative dynamics of the model. Indeed, panels (b) and (c), the classes of human infected peak after 3 years and decrease according to the decrease of the fractional parameter α. We also note that the total numbers of susceptible humans decreases rapidly according to the increase of the fractional parameter (panel (a)), while the total numbers of susceptible vectors increases when the the fractional parameter α increases (panel (f)).

Figure 10 validates Theorem 6. It is clear that varying the fractional order parameter α does not influence the model dynamics whenever R 0 < 1. Indeed, whatever the value of α, the infected compartments tend to zero asymptotically whenever R 0 < 1. This validate the fact that the Buruli-free equilibrium of the fractional model is GAS whenever R 0 < 1.

Conclusion

In this work, we formulated and studied a compartmental model for the Buruli-Ulcer transmission dynamics using classical and fractional (in the Caputo sense) derivatives. We first formulate and study the model with integer derivative by proving the existence of equilibrium points, determining the basic reproduction number R 0 , and proving the asymptotic stability of the Buruli-free equilibrium whenever R 0 < 1. Then, for R 0 > 1, we also proved the global stability of the unique endemic equilibrium point using the general theory of Lyapunov.

We calibrated our model by performing parameter estimation using real data from Cameroon, while sensitivity analysis is conducted to determine important parameters in the disease dynamics. It followed from the estimated value of the basic reproduction number R 0 is equal to 2.0843, which implies that the system is in an endemic state since R 0 is greater than unity. To determine the number of waves (epidemic peaks), we computed the strength number A 0 . We found that A 0 = 17 > 0 which means that BU epidemic has multiple waves. Numerically, we found that BU has at least 2 epidemic peaks. To determine parameters which play an important rule in the model dynamics, we performed sensitivity analysis of R 0 by computing sensitivity index of each model parameters. The results of this sensitivity analysis showed that the mortality rate of vector µ v , the transmission probabilities β h and β v , as well as the recruitment rate of vectors are the important parameters which govern the BU dynamics. 2 such that R 0 = 0.9982889062331324 < 1. In this case the Buruli-free equilibrium M 0 = (2090, 0, 0, 0, 69971, 0) ′ is globally asymptotically stable. 2 such that R 0 = 2.0843 > 1.

In this case the Buruli-free equilibrium is unstable and the unique endemic equilibrium point is globally asymptotically stable. By replacing integer derivatives with Caputo derivative, we obtained our Buruli ulcer fractional model. As for the classical model, we performed stability analysis of equilibrium points, as well as existence and uniqueness of solution. Based on Adams-Bashforth-Moulton (ABM) method, we constructed a numerical scheme. We thus conducted several numerical simulations to validate theoretical results and also, see the role playing by the variation of the fractional parameter α in the disease dynamics. Indeed, simulation results confirmed the theoretical results (Theorems 6 and 7) which claimed that the disease-free equilibrium (respectively the endemic equilibrium point) of both models are GAS whenever R 0 ≤ 1 (respectively R 0 > 1). Also, we found that, in quantitative point of view, fractional model allows more flexibility than the classical one. The epidemic peaks depend of the variation of the fractional-order α.

The direct perspectives of this work consist to: (1) use another fractional derivatives like Caputo-Fabrizio, and Atangana-Baleanu derivatives and compare numerically the obtained results, and (2) implement optimal control which permits to determine which control strategy is better to control of BU disease in Cameroon.

where e i for i = 1, ...7 are positive constants to be determined latter, its time-derivative gives

Ż(S h , E b , I b , R b , S v , I v ) = e 1 1 - S * S Ṡh + e 2 1 - E * b E b Ėb + e 3 1 - I * b I b İb + e 4 1 - I * b I b İb + e 5 1 - R * b R b Ṙb + e 6 1 - S * v S v Ṡv + e 7 1 - I * v I v İv = e 1 1 - S * S [π h + ϕR b -β h S h X -µ h S h ] + e 2 1 - E * b E b [β h S h X -k 1 E b ] + e 3 1 - I * b I b [γ 1 E b -k 2 I b ] + e 4 1 - I * b I b [γ 1 E b -k 2 I b ] + e 5 1 - R * b R b [ηI b -k 3 R b ] + e 6 1 - S * v S v [π v -β v S v Y -µ v S v ] + e 7 1 - I * v I v [β v S v Y -µ v I v ]
At equilibrium, we have the following relations

               π h = -ϕR * b + β h S * h X * + µ h S * h , β h S * h X * = k 1 E * b , γ 1 E * b = k 2 I * b , ηI * b = k 3 R * b , π v = β v S * v Y * + µ v S * v , β v S * v Y * = µ v I * v . (45)
We thus obtain

Ż(S h , E b , I b , R b , S v , I v ) = e 1 µ h S * h 2 - S h S * h - S * h S h + e 6 β v S * v Y * S * v I * v 2 - S v S * v - S * v S v + R *
S h S * h - S * h S h + β v S * v Y * S * v I * v 2 - S v S * v - S * v S v + β v S * v Y * 2 + Y Y * - S * v S v - I v I * v - S v Y S * v Y * I * v I v + 2µ h E * b 2 + X X * - S * h S h - S h X S * h X * E * b E b - E b E * b + 2k 2 I * b 3 + X X * - S * h S h - I * b I b E b E * b - I b I * b - S h X S * h X * E * b E b + 2ϕR * b I b I * b 1 - R * b R b + S * h S h 1 - R b R * b
Since the arithmetic means is greater than the geometric means, we conclude that the first two terms of the right-hand side of above expression are negative. It then follows that Ż(S h , E b , I b , S v , I v ) ≤ 0 if the following conditions hold:

2 + Y Y * - S * v S v - I v I * v - S v Y S * v Y * I * v I v ≤ 0, 2 + X X * - S * h S h - E b E * b - S h X S * h X * E * b E b ≤ 0,
I b I * b 1 - R * b R b + S * h S h 1 - R b R * b ≤ 0. ( 46
)
Using the function f (ξ) := 1 -ξ + ln(ξ) which is negative for all ξ > 0 and equal to zero whenever ξ = 1, we obtain Thanks to the LaSalle invariance principle [START_REF] La Salle | The stability of dynamical systems[END_REF], we conclude that if R 0 > 1, then M 1 is GAS in Ψ. This ends the proof.

2 + Y Y * - S * v S v - I v I * v - S v Y S * v Y * I * v I v = 2 + Y Y * - S * v S v - I v I * v -1 - S v I * v S * v I v 1 - Y Y * -1 + S v I * v S * v I v + Y Y * = -1 - S v I * v S * v I v 1 - Y Y * + 1 - S * v S v + 1 - I v I * v + 1 - S v I * v S * v I v ≤ -1 - S v I * v S * v I v 1 - Y Y * ≤ 0 2 + X X * - S * h S h - E b E * b - S h X S * h X * E * b E b = 2 + X X * - S * h S h - E b E * b -1 - X X * 1 - S h S * h E * b E b -1 + S h S * h E * b E b + X X * = -1 - X X * 1 - S h S * h E * b E b + 1 - S * h S h + 1 - E b E * b + 1 - S h S * h E * b E b ≤ -1 - S h S * h E * b E b 1 - X X * ≤ 0, 3 + X X * - S * h S h - I b I * b - I * b I b E b E * b - S h X S * h X * E * b E b = 3 + X X * - S * h S h - I b I * b - I * b I b E b E * b -1 - X X * 1 - S h S * h E * b E b -1 + S h S * h E * b E b + X X * = -1 - X X * 1 - S h S * h E * b E b + 1 - S * h S h + 1 - I b I * b + 1 - I * b I b E b E * b + 1 - S h S * h E * b E b ≤ -1 - X X * 1 - S h S * h E * b E b ≤ 0.

Note that

I b I * b 1 - R * b R b + S * h S h 1 - R b R * b < 0 ⇔ S *

  Figure 1: (a) Real data versus model simulation; (b) Forecasting of Buruli Ulcer in the next fifteen (15) years (from 2015 to 2030).

  3; b 1 = 2004; c 1 = 0.6758; a 2 = -391.8; b 2 = 2014; c 2 = 3.04; a 3 = 64.25; b 3 = 2006; c 3 = 0.3512; a 4 = 552.9; b 4 = 2017 and c 4 = 11.5.The results are depicted in Figures2 and 3

Figure 2 :

 2 Figure 2: Gaussian fitting of order 4 (red line) using the reported cases of Buruli Ulcer (BU) in Cameroon between 2001 and 2014 (blue dots).

Figure 3 :

 3 Figure 3: Gaussian fitting of order 4 (red line) using the reported cases of Buruli Ulcer (BU) in Cameroon between 2001 and 2014, and prediction of of BU in the next sixteen (16) years (blue line).

Figure 4 :

 4 Figure 4: Expecting number of cases with multiple waves.

Figure 5 :

 5 Figure 5: Sensitivity indices for R 0 against model parameters.

Figures 6

 6 Figures 6 depict how the basic reproduction number R 0 varies in term of most important model parameters.

(

Figure 6 : 3 -

 63 Figure 6: 3-D plots of the basic reproduction number R 0 in term of most important model parameters.

Figure 7 :

 7 Figure 7: Time-series of model (2) with β v = 0.18 and the other parameter values as in Table2such that R 0 = 0.9982889062331324 < 1. In this case the Buruli-free equilibrium M 0 = (2090, 0, 0, 0, 69971, 0)′ is globally asymptotically stable.

Figure 8 :

 8 Figure 8: Time-series of model (2) with the parameter values as in Table2such that R 0 = 2.0843 > 1.In this case the Buruli-free equilibrium is unstable and the unique endemic equilibrium point is globally asymptotically stable.

Figure 9 :Figure 10 :

 910 Figure 9: Time-series of fractional model[START_REF] Feng | On the role of variable latent periods in mathematical models for tuberculosis[END_REF] showing the impact of varying the fractional-order parameter α on model dynamics. The fractional-order being α = 1, α = 0.9, α = 0.8, α = 0.7, and α = 0.6.

b e 5 k 3 -e 1 ϕ + e 5 k 3 + 1 X+ β v S * v Y * e 6 + e 7 + e 6 Y Y * -e 6 S 7 -e 6 )

 33166676 β h S * h X * e 1 + e 2 + e S v Y S * v Y * -e 2 µ h E b .Setting e 3 = e 4 = e 6 = 1, and e 1 = e 2 = (e 3 + e 4 ), e 5 = (e 3 + e 4 ) ϕ k 3 and e 6 = e 7 , We obtainŻ(S h , E b , I b , R b , S v , I v ) = 2µ h S * h 2 -

  (S h , E b , I b , S v , I v ) ≤ 0. So that lim t→∞ (S h (t), E b (t), I b (t), R b (t), S v (t), I v (t)) → M 1 := (S * h , E * b , I * b , R * b , S * v , I * v ) .

Table 1 :

 1 Description of state variables and model parameters. h is the human life mean expectancy while µ v is the mortality rate of vectors; κ is the disease-induced death rate; the transmission rate of infection between susceptible humans and infectious vectors is denoted by β h while β v denoted the transmission rate of infection between infectious humans and susceptible vectors; γ 1 is the transition rate between E b and I b ; η 1 is the recovered rate of infectious humans. We resume states variables and parameter descriptions in Table1. We thus express our new Buruli ulcer model using ordinary derivatives as follows:

	State variables Description
	S h	Susceptible humans
	E b	Infected humans in latent stage
	I b	Infectious population
	R b	Recovered population
	S v	Susceptible vectors
	I v	Infected vectors
	Parameters	Description
	π h	Recruitment rate of susceptible humans
	π v	Recruitment rate of water bugs
	ϕ	Rate of recovered people who loosed immunity after
		recovered from Buruli ulcer
	µ h	Natural death rate of human population
	µ v	Mortality rate of the vector population
	β h	Transmission rate of Mycobacterium ulcerans
		from infectious vector to a susceptible human
	β v	Transmission rate of Mycobacterium ulcerans
		from infectious human to a susceptible vector
	κ	Disease-induced death rate
	η	Recovery rate of Buruli Ulcer infected people
	γ 1	Contact rate of Mycobacterium Ulcer with the population
	is denoted by N v (t) and includes the susceptible vector population noted S v (t) and the
	infectious vector population noted I v (t), at any time t. So, N h (t) = S h (t)+E b (t)+I b (t)+R b (t)
	and N v (t) = S v (t) + I v (t). The following model parameters, which are nonnegative, have
	described as follows: π h represents the recruitment rate of population; π v denoted recruitment
	rate of water bugs while ϕ denotes the rate of recovering people who have loosed immunity
	after recovering from Buruli ulcer; 1/µ

Table 2 :

 2 Estimated parameters. The corresponding value of the basic reproduction number is R 0 = 2.0843.

	Parameters	Value	Source Parameters	Value	Source
	µ h	1/55.5	[2, 4]	β h	0.784684 Fitted
	k	0.203789	Fitted	β v	0.036985 Fitted
	π h	58.660064	Fitted	η	0.760527 Fitted
	π v	172813.870585 Fitted	γ 1	0.482389 Fitted
	ϕ	0.169881	Fitted	µ v	0.589891 Fitted

Table 3 :

 3 Sensitivity indices of model parameters.

	Parameters Sensitivity index Parameters Sensitivity index
	π h	-0.5000	β h	0.5000
	π v	0.5000	β v	0.5000
	ϕ	0	κ	-0.1037
	µ h	0.4728	η	-0.3871
	µ v	-1	γ 1	0.0180

Acknowledgements

HA and HPEF acknowledge the grants and support of the Cameroon Ministry of Higher Education through the initiative for the modernization of research in Cameroon's Higher Education.

Declaration of Competing Interest

The authors declare that they have no known competing interests.

Author Contributions

Conceptualization: HA & HPEF; Methodology: RF, HA, AK; Formal analysis and investigation: RF & HA; Writing-original draft preparation: RF & HA; Writing-review and editing: AH, HPEF, AK; Funding acquisition: HA; Resources: HA, HPEF, AK; Supervision: HA & HPEF.

A. Proof of Theorem 4

Considering the function