Fractional-order modeling of visceral leishmaniasis disease transmission dynamics : strategies in eastern Sudan

Fractional calculus has become an increasingly important tool for biological process modeling. In this study, we propose a novel fractional-order model that is used to analyze visceral leishmaniase (VL) transmission dynamics, incorporating the e ff ects of interventions. By fitting the model to 22 years of cumulative VL cases, we estimate intervention e ff ectiveness and key parameters. Our results show a high disease prevalence, with reported cases accounting for less than 20% of all infections. Both human and vector control measures show limited e ff ectiveness. Future projections point to a further increase in cumulative cases.


Introduction
The SEIR model is a mathematical framework frequently employed in epidemiology to depict the dynamics of infectious diseases in a population.It comprises four main compartments denoted by the acronym SEIR, representing Susceptible, Exposed, Infected, and Recovered individuals [1][2][3].The model begins with individuals classified as susceptible to the disease, progresses to the exposed phase, where they have been exposed but are not yet infectious, moves to the infected phase with the potential for transmission, and ultimately concludes with the recovered category, assumed to confer immunity.By utilizing the SEIR model, researchers and public health authorities can gain insights into the progression of infectious diseases, estimate crucial epidemiological parameters, and assess the effectiveness of interventions such as vaccination or social distancing measures in curbing the disease's spread within a population [4][5][6].
Visceral Leishmaniasis, also referred to as Kala-Azar, is a chronic, neglected, and deadly infectious disease caused by the flagellate protozoan parasite Leishmania.It affects 2 − 4 hundred thousand individuals worldwide, with an approximately 0.1% mortality rate [7][8][9].VL is a zoonotic disease infecting humans and animals such as dogs, jackals, rats, and foxes [10].The vector that transmits the Leishmania parasite to humans and other animal reservoirs is typically the phlebotomine female sandfly, which is about one-third the size of a mosquito.The parasite is rarely transmitted through needle sharing, blood transfusions, or mother-to-child transmission [11].
Sandflies transmit Leishmania parasites between hosts (individuals and animals) through biting, enabling promastigotes (the pathogenic component) to penetrate macrophages and other mononuclear phagocytic cells under the skin [9].The pathogens then evolve into non-flagellated intracellular ovoid amastigotes with a large nucleus and a smaller kinetoplast [12].Leishmania is harmful to host cells, like macrophages, and affects the liver, spleen, bone marrow, and lymph nodes.Recognized VL symptoms include prolonged irregular fever with or without rigors, cough, loss of appetite, diarrhea, vomiting, jaundice, edema, epistaxis, and anemia.Most VL-positive individuals experience one or more of these symptoms within 3 − 6 months of infection [13].Moreover, VL may lead to post-kalaazar dermal leishmaniasis (PKDL) within six months to a year after recovery.PKDL appears as a maculopapular or nodular rash on the face, upper arms, trunks, and other body parts.People with PKDL are capable of spreading the VL disease.
An effective vaccine against VL has yet to be developed, and the main control measures are case management based on prompt diagnosis and treatment, which is vital for reducing morbidity and averting death, detection and treatment of infected animals [14], insecticide spraying for vector population control, and use of insecticide-treated bed nets.VL control efforts also include the establishment of surveillance systems for monitoring and reporting, as well as community engagement and public education [15].
VL is common in East Africa, particularly in impoverished and conflict-torn areas.Sudan is home to Phlebotomus orientalis sandflies [7], the primary carrier of the disease pathogen, and VL endemic has been established in the country since 1908 [16,17].The disease is endemic in the Darfur region of western Sudan, the Nuba Mountains of Kordofan region [18], and the vast area from the western bank of the White Nile to the Ethiopian border, which encompasses central and eastern states.The last region includes El Gadaref state the epicenter of the disease in the country [19].VL-related statistics in Sudan are alarming; it has the highest prevalence of PKDL in the world [20,21] with an incidence rate of 38/1000 per year and a case fatality rate as high as 20.5% in population-based studies.Clinical characteristics associated with VL and PKDL, such as enlarged lymph nodes, hepatomegaly, and lymphadenopathy, are more prevalent in Sudan than elsewhere [22].The disease also incurs a huge economic burden on the already poor population treatment for a single VL episode costs 53% of household income [23].
Sudanese health authorities, in partnership and collaboration with various stakeholders such as the World Health Organization, the Leishmaniasis East Africa Platform, and the Drugs for Neglected Diseases Initiative, have implemented numerous initiatives to control the disease over the last two decades.These efforts include the establishment of early diagnosis and treatment sites, the introduction of a computerized registration system, the establishment of pharmacovigilance and passive surveillance mechanisms for sandflies in villages and other remote areas, vector control strategies executed through the malaria control program, operational research, and outreach activities.Much of these efforts were focused on El Gadaref state.
However, the disease is still persistent, with an average annual number of reported cases in the last two decades exceeding 3, 500 cases (see Figure 1).The blue horizontal line represents the average number of reported cases per year.The data were obtained from the VL control program, Ministry of Health, El Gadaref state, through personal communications with officials.
Mathematical modeling has been an important aspect of investigating infectious disease transmission, and numerous models for VL transmission dynamics have been proposed.The major theme in these models is that they are built on deterministic SEIR-type models and take into account three populations that are important in VL dynamics: human beings, animal reservoirs, and vectors [24][25][26][27][28]. Elmojtaba, in collaboration with others, has exhaustively examined and established several theoretical and simulation-based results about VL transmission in Sudan from various perspectives.For example, in Ref. [29], it was shown that human treatment and vector control can effectively eliminate the disease under specific conditions, and in Ref. [30,31], the dynamics of co-infection between malaria and VL were investigated.Additionally, Refs.[27,32] applied the optimal control theory within dynamical models for VL transmission to identify the optimum combination of strategies for combating the disease.
Following in the manner of previous studies, the current research attempts to develop an SEIR-type model for VL transmission dynamics that explicitly accounts for the influence of different interventions to control the disease in the three populations.We model the effect of the interventions as percentages and assume that the reported human VL cases only represent a fraction of actual infections in the community.
Using a Bayesian approach, we fitted the proposed model to annually reported human VL cases from El Gadaref state over 22 years (2000-2021) and estimated the value of each intervention, as well as the values of the transmission rates that determine the disease's dynamic in this region.
The remainder of this article is structured as follows.In Section 2.1, we present the formulation of the proposed model, and in Section 2.2, we prove its validity from a mathematical perspective, including the expression for the basic reproduction number R 0 .In section 3 we present the analysis of the model.Section 4.1 describes the parameter estimation process, and Section 4.2 presents the sensitivity analysis of the model, where we identify the most significant parameters.A brief summary and discussion of the results are provided in Section 5.

Definition 2 ([20]
).The Atangana-Baleanu integral of f (in Caputo sense) of order ℑ is defined by The Laplace transform of the Atangana-Baleanu fractional derivative of f of order ℑ in Caputo sense is given by where L is the Laplace transform operator.Definition 4 ( [29]).
The Laplace transform of the Mittag-Leffler function E ℑ (.) is given by

Model formulation
This section describes the transmission dynamics of VL using a modified version of an SEIR type model.The model includes three populations: people, vectors, and reservoirs, the total human population N h is divided into six sub-populations, namely susceptible, exposed, detected infected, undetected infected, PKDL-infected, and recovered individuals denoted, respectively, by S h , E h , I d , I u , P h , and R h .Likewise the total vector population N v is divided into two sub-populations: susceptible S v and infected I v vectors, and the total reservoir population is divided into susceptible S r , and infected I r reservoirs.At given time t, we have The transmission procedure of VL that we propose is as follows: Initially, all the three populations remain susceptible to acquiring the disease's parasite, and that they grow at constant birth rates λ h , λ v and λ r and incur natural death, regardless of their infection status, at rates µ h , µ v , and µ r for human, vector, and reservoir populations, respectively.Susceptible humans are exposed to the disease's parasite with force of infection  (t) , and further reduced by rate α r (t), where α r (t) denotes for the efficiency of the measures (if any) taken to control the reservoir population such as killing infected animals and canine-vaccinations.It is worth noting that the transmission probability from an infected vector (β 1 ) is assumed to be the same for susceptible humans and reservoirs.The chance of VL transmission from an infected host to vector is determined by the host's stage [33].Therefore, we include parameters β 2 , β 3 , β 4 , and β 5 to represent, respectively, the transmission probability per detected, undetected, and PKDL individuals, as well as an infected reservoir.Thus, the overall infection force for vector population is given by Susceptible vectors are also further reduced by rate α v (t) that corresponds to the efficiency of measures taken ABC 0 to control the vector population, such the spray of insecticides.The VL transmission dynamic described above is depicted graphically in Figure 2 and mathematically by the ten differential equations in system 1.
The following three equations are also implicit in the model's formulation: . All of the compartments, their initial values, and parameters in system 1 are non-negative.In section 3 we present tables that contain a full description of the model's parameters and their values.
The flowchart of the proposed VL transmission model, the definition of all variables, parameters,transition rates and f 1 , f 2 and f 3 are the force of infections for human, reservoir, and vector population, respectively are shown in Fig. 2 2

.2. Mathematical analysis
From the three equations in system 2, it is obvious to note that in the absence of the disease (δ 1 = 0 ,δ 2 = 0) Visceral Leishmaniasis induced death rate for individuals" indicates that the death rate from Visceral Leishmaniasis (VL) for individuals is zero.This means that no individuals died from VL during the specified period or within the specified population.and in the long run (t → ∞), the total number of human population is determined by natural growth and death rates, N h = λ h µ h .Some of the reservoir and vector population control measures are implemented in other health programs, such as pesticide spraying in malaria program, and may remain in place even after the VL disease has eradicated.Therefore, both the total number of reservoir and vector populations are governed by natural causes and the effects of such interventions.In practice, however, to the best of our knowledge, there is no other program for controlling animal populations in the region of interest.Though mathematical analysis allows for it to be different, it is reasonable to set α r = 0. Thus, in the long run we have , and the disease free equilibrium coordinates of system (1), denoted by E 0 , is Considering the fact that system (1) calculate number of living things, the sensible region to study it is which is bounded and positively-invariant domain for model (1) and given non-negative initial conditions in R 10 + in Section 3.2 .Using the next generation operator [34], we obtained the basic reproduction number -the average number of susceptible individuals infected by a single infected case -R 0 of model (1) as stated in Equation (8).We also discussed and proved the disease free equilibrium conditioned on R 0 as well as the stability of system (1) in Section 3.2.

Analysis of the model
In this section, we seek to qualitatively study the dynamical properties of the VL disease model (1) in the main text, which quantifies disease spread or extinction in a population.

Positivity and boundedness
Theorem 1. (positivity) if the initial conditions are nonnegative then the corresponding solution set S h (t), E h (t), I d (t), I u (t), P h (t), R h (t), S r (t), I r (t), S v (t), I v (t) , of system (1) is non-negative for all t > 0.Moreover, and the region is positively-invariant for the model (1) any non-negative initial conditions in R 10 + .
Proof:.By substituting S h = 0, E h = 0, I d = 0, I u = 0, P h = 0, R h = 0, S r = 0, I r = 0, S v = 0, and I v = 0 in the first, second, third, ..., and the tenth equation of system (1), respectively, we have All of these equations are non-negative, and since the initial conditions are also non-negative, then using Proposition B.7 in [35] we can conclude that S h (t) ≥ 0, E h (t) ≥ 0, I d (t) ≥ 0, I u (t) ≥ 0, P h (t) ≥ 0, R h (t) ≥ 0, S r (t) ≥ 0, I r (t) ≥ 0, S v (t) ≥ 0, and I v (t) ≥ 0, for all t > 0. Now, we prove that the solution is bounded on interval [0, ∞) by using the standard comparison theory [36] and considering the positivity of the solution on [0, ∞) established for the equations of system (1) Therefore, N h (t), N r (t) and N v (t) are bounded on the domain [0, ∞), and any solution of model ( 1) is positively invariant in the region given by □ We denote the disease free equilibrium coordinates of system ( 1) by E 0 , and it is given by To obtain the above result, we assume α r = 0 in the long run.

The basic reproduction number
The basic reproduction number denoted by R 0 is the expected value of secondary infections rate per time unit.There are six infected classes in system (1): E h , I d , I u , P h , I r , and I v .Therefore, system (1) is reduced to the following new system: We assume that the average detection rate is constant ϵ(t) with the value ϵ for all time t and that the non-pharmaceutical interventions 2) in vector from [37], as follows where We can verify that system (3.2) satisfies the conditions (A1)-(A5) in [34].From the above functions, it is easy to conclude that: (A1) if x i ≥ 0, then F i (x), V + i (x), and V − i (x) ≥ 0, for i = 1, 2, 3, 4, 5, 6. i.e., each compartment is directed transfer of living things and thus they are all non-negative; recall that ϵ and α ∈ [0, 1].(A2) if x i = 0, then V − i (x) = 0, for i = 1, 2, 3, 4, 5, 6. i.e., removal of living things is not from an empty compartment.(A3) F i = 0 for i {1, 2, 3, 4, 5, 6}.i.e., the incidence of infection for uninfected compartments (S h , R h , S r , and S v ) is zero, see the main text.(A4) if x ∈ E 0 , then F i (x) = 0 and V + i (x) = 0, for i = 1, 2, 3, 4, 5, 6. (A5) if F(x) = 0, then the Jacobian matrix of system (3.2) is given by The eigenvalues of D f (x 0 ) are: and each one has negative real part. where The reproduction number of the system equation ( 2), R 0 , is given by the Letting

Stability analysis of the model
Theorem 2. System (1) at disease-free equilibrium point E d f e is locally asymptotically stable if R 0 < 1 and unstable if R 0 > 1 .
Proof:.The Jacobian matrix [34] with respect to the system of equation ( 1) At E d f e , is given by where .
We can clearly see that four eigenvalues are negative, namely, Z 1 = −µ h , Z 2 = −A 9 , and Z 3 = −A 10 , Z 4 = −A 11 .For the rest of eigenvalues, we consider the following reduced matrix [38]: After simplification, we get: The eigenvalues of J 1 (E d f e ) take the following form: The last eigenvalue also has the negative real part since implying that Z 10 < 0,then R 0 < 1, so the system of equation ( 1) is locally asymptotically stable for R 0 < 1, and implying that R 0 > (µ r + α r0 )(µ v + α v0 ), then Z 10 > 0, we obtain that R 0 > 1, hence the system ( 1) is unstable for R 0 > 1. □ Table 1 contains an extensive discussion of the model's parameters, including their values and ranges.For human, reservoir, and vector populations, the table includes detection rates, intervention efficiencies, transmission probabilities, incubation durations, death rates, treatment rates, natural recovery rates, progression rates, and constant growth and death rates.The parameter values are estimated or collected from a variety of sources, and they play critical roles in creating the model's dynamics.
Table 2 focuses on initially value estimations, which are critical for establishing the model's initial circumstances.These numbers show the starting numbers of different human and reservoir classifications, such as vulnerable, exposed, undiscovered infected, PKDL, and immune individuals, as well as the size of the human, reservoir, and vector populations.The table also includes the estimated reproduction number, designated as R 0 , which offers information about the disease's potential spread.
Finally, Table 3 displays the Partial Rank Correlation Coefficient (PRCC) and accompanying p-values for each parameter, including I d (undetected infected human class) and mathcalR 0 (estimated reproduction number).The PRCC calculates the strength and direction of the linear correlation between a parameter and the target variable while accounting for the influence of other variables.The p-values represent the correlation's significance.Parameters with high PRCC values and low pvalues have larger impacts on the dynamics of I d and mathcalR 0 in the model, making them crucial for understanding the disease's behavior.Parameters with low PRCC values and high p-values, on the other hand, have little influence.

Parameters estimation
Approximate values for a subset of parameters in model ( 1) that pertain to the characteristics of VL disease, such as the average incubation period, are known and have been obtained from relevant literature (refer to Table 1).However, certain parameters that represent the rates of interactions (β 1 , β 2 , β 3 , β 4 , and β 5 ), efficacy of interventions (ϵ, α h , α r , and α v ), as well as the initial values of N h , E h , I u , P h , R h , N r , N v , I r , and I v remain unknown beforehand.We collectively denote these unknown parameters as θ, i.e.
The exact total population of El Gadaref state in 2000 is less certain, but according to [44,45] it stood around one million.We model the initial number of human population around this numbers as N h (0) = K h × 10 6 , where K h ∼ u(0.8, 1.4), and set I d (0) to be the cumulative reported VL cases in the year 2000 (see Figure 3).As in the work of [27] we modeled the initial value of susceptible human as a fraction of N h (0), i.e., S h (0 , respectively, where K rh and K vh represent the total number of reservoirs and vectors per human, respectively.The initial values for the other compartments are considered to be fractions of the compartments to which they belong (see Table 2).Information about the time points where the interventions went into effect is unavailable to us.Therefore, we assumed that the average detection rate is constant ϵ(t) with the value ϵ for all time t.Likewise, for the non-pharmaceutical interventions for protecting human population α h (t) = α h , reservoir control interventions α r (t) = α r , and vector control interventions The cumulative number of reported human VL cases at time t, represented by y, can be estimated from system (1) as follows The yearly cumulative number of reported VL cases are count data and can be considered as random number of occurrences in a year, and thus can be modeled by negative binomial (NB) distribution.Therefore, the joint posterior distribution of the parameters vector θ and hyper parameter ϕ for dispersion is then given by where L represents the number of years, y(t) is obtained from 11, and P(θ) represents the density of joint prior distribution of all the parameters in vector θ, P(ϕ) is the density of prior distribution for dispersion parameter ϕ.We assumed a flat uninformative prior distribution (either u(0, 1) or positively truncated normal distribution N(0, 1)) for each parameter in θ, and ϕ −1 ∼ Gamma(0.1, 0.1).
The following graph Figure 3a displays the annual cumulative number of visceral leishmaniasis (VL) cases that have been officially recorded for the state of El Gadaref from 2000 to 2021.The red curve depicts the outcomes simulated by a model, while the blue dots reflect the actual reported cases in the data.
The Figure 3b also offers a forecast for the total number of reported VL cases in El Gadaref state during the subsequent five years, from 2022 to 2026.This forecast, which is based on model simulations, provides insightful information on anticipated trends in the number of VL cases in the region throughout the given time period.
We fit Eq.( 12) to the cumulative number of reported VL cases in El Gadaref State from 2000 to 2021(L = 22), using Markov Chain Monte Carlo (MCMC) approach.Specifically, we implemented the No-U-Turn sampler [46], an advanced MCMC algorithm, within the Turing [47] ecosystem for Bayesian computing paradigm in Julia programming lan- guage [48].No-U-Turn algorithm converges fast to the target distribution, avoids repeated random walks, and automatically adapts the learning rate [45].We run the program for 15, 000 iterations; discarded the samples in the first 5, 000 iterations as burn-in, and calculated point estimate and 95% incredible interval for each parameter in θ from samples in the remaining iterations (see Table 1).Figure 3a compares the results from model We further used the estimated parameter values and projected the cumulative VL cases that will be reported in El Gadaref state for 5 years, from 2022-2027.The cumulative cases of VL will be steadily increasing (see Figure 3b).

Sensitivity analysis
We performe sensitivity analysis of model ( 1) using combination of Latin Hypercube Sampling (LHS) and PRCC.Such analysis allows for the identification of parameters that have an important effect on the model's dynamic [49].LHS is a stratified sampling technique in which the range of interest for each parameter is partitioned into N equidistant intervals and random samples are drawn from each interval.The PRCC is an effective technique for exploring nonlinear, monotonic relationship between input parameters and model output [49].PRCC results quantify the relationship between a parameter and an outcome.A positive result implies increasing the parameter increases the outcome, whereas a negative result indicates that increasing the parameter decreases the outcome.We assume a uniform distribution for each parameter (or the corresponding fraction of each parameter, see Section 4.1) in θ, randomly generated 3000 samples of each parameter, and calculated PRCC results for the total number of reported cases I d and the basic reproduction number R 0 .Figure 4 shows the results of the analysis for each input parameter, at 95% confidence level.From Figure 4 a, it can be noted that the parameters that have significant positive effect on I d are ϵ, β 1 , β 5 , N v , and S h ; whereas α h and N r have significant effect impact on I d .The highest positive influence on I d comes from ϵ (with correlation coefficient r = 0.562 and p-value = 1.59e −5 , see Table 3 in Section 3), which intuitively understandable -the stronger the reporting rate, the larger the number of detected cases.The highest negative influence on I d comes from α h (r = −0.192and p-value = 0.000283), which is also intuitive; the stronger the interventions for human population protection, the lesser the reported cases in the community.
However, very few parameters have significant effect on R 0 (see Figure 4 b), which include β 1 , N v , α r , α v and N r .In particular β 1 (r = 0.21 and p-value =4.92e −215 ) and N v (r = 0.37 and p-value = 4.92e −215 ) have positive correlation with R 0 (r = 0.21 and p-value =9.75e −71 ), whereas α r (r = −0.255and p-value = 5.77e −34 ), α v (r = −0.0613and p-value = 0.00153) and N r (r = −0.459and p-value = 1.01e −135 ) have negative correlation with R 0 .Interestingly, we found two observations that are inline with the conclusions of other studies: increase in β 1 , the interaction rate between sandflies and hosts (humans and reservoirs), is and important factor that increase the basic reproduction number R 0 [40], and decrease in N r , the total size of reservoir population, by culling or otherwise, leads to an increase of VL cases among human population [29,50].

Discussion and conclusion
Visceral Leishmaniasis (VL) is a major public health concern that has affected millions of people worldwide.The disease has been endemic in Sudan for several decades, and various control strategies have been established for contain it.Much of the efforts were concentrated in the eastern state of El Gadaref, the disease's epicenter in the country, where around 80% of cases occurred.The efficacy of these interventions have yet to be thoroughly examined.
In this study, we constructed a mathematical model for VL transmission dynamics that includes three population at the disease interplay: humans, animals, and vectors.The Model is built on an SEIR type model and incorporates the effects of interventions in a systematic manner.We proved the validity of the model mathematically and derived formula for its basic reproduction number R 0 (see Equation 8).Using Partial Rank Correlation Coefficient method coupled with Latin Hypercube Sampling technique we also performed sensitivity analysis on the model, and identified the important parameters that have, either positive and negative, significant influence on the number of reported human VL cases I d and R 0 (see Figure 4).
We fitted model 1 to the yearly cumulative data of reported human VL cases in EL Gadaref state for 22 years (2000−−2021), estimated the model parameters, and compared the model predictions with the actual cumulative reported VL cases (see Section 4.1 and Figure 3a).According to the model, the disease is quite prevalent among the local population -the basic reproduction number is quite larger than one.The reported cases represents fewer than 20% of the total number of infected patients, with infected reservoirs (β 5 ) and humans in the PKDL stage (β 4 ) being the primary sources of disease propagation.If the existing level of the measures remains the same, future projections show that the yearly cumulative cases will continue to rise (see Figure 3b).The model provided a very limiting value for reservoir control intervention efficacy (α r < 15%).To some extent, this is correct; only a few efforts have been undertaken to detect, treat, or cull infected animals.Furthermore, the model predicted that the efficacy of both human and vector measures would be less than 50% (see Table 1 in SI file).
In general, we reach the conclusion that existing VL control strategies are inadequate and inefficient, and that much has to be improved in order to meet Sudan's National Leishmaniasis Control Program (NLCP) goal of eliminating VL as a public health problem by 2030 [51].We regarded the yearly cumulative reported VL as random numbers for over-dispersed count data, thus we modeled them using the negative binomial distribution and calibrated our computations accordingly.We presented the average values for the efficacy of the interventions over the span of 22 years because we did not have the exact time points when these interventions were implemented, suspended or resumed.We would have better model resolutions if we had this information.On the future, we are planning to solve and analyze nonlinear fractional models by new methods [52][53][54].

Figure 2 .
Figure 2. Diagram of the proposed VL transmission model.

Figure 3 .Figure 4 .
Figure 3. Yearly Cumulative Number of Reported VL Cases in El Gadaref State (2000 − 2021) and Projection for (2022 − 2026) (1) simulation to the actual data of cumulative reported VL cases in El Gadaref.The model results are in accordance with the actual accumulated data, only a few data point were missed by the model.The estimated value for R 0 = 21.04 (95% CI: 18.28 − 36.23), which is far larger than one, and thus VL will remain in the region for the foreseeable future, if the level of exiting measures remains the same.The model estimated a limiting value for the average case report ϵ = 0.187 (95% CI: xx − xx).α r = 0.xx (95% CI: xx − xx), α h = 0.xx (95% CI: xx − xx).α v = 0.xx (95% CI: xx − xx).
where β 1 is the probability of infection transmission to host (human or animal) after being bitten by an infected vector, and α h (t) ∈ [0, 1] represents the efficiency of non-pharmaceutical interventions made for protecting human population, e.g., the use of treated bed-nets, at time t.If these interventions were 100% efficient -that is α h (t) = 1, no new incidences of VL would occurred.Exposed individuals remain in E h compartment for1ω days of year, after which they Figure 1.Yearly counts of reported human VL cases for El Gadaref state, Eastern Sudan, from 2000 to 2021.demise naturally at rate µ h , or become detected infectious I d at rate ϵ(t), or undetected infectious I u at rate 1 − ϵ(t).The parameter ϵ(t) represents the efficacy of concentrated surveillance activities for case detection-e.g., capacity building and early diagnosis-that take place at time t.Without surveillance systems, all infected individuals would go undetected.Individuals in I d compartment are subject to treatment, recover from the disease due to the treatment at rate σ 1 , develop PKDL at rate ρ 1 , or die at rates δ 1 and µ h as result of VL infection and naturally, respectively.Undetected infectious individuals, I u , are not recorded by the surveillance systems nor eligible for treatment.
However, they recover naturally at rate γ 1 owning to the individual's immunity, develop PKDL at rate ρ 2 , and die at rates δ 2 and µ h as result of VL infection and natural mortality, respectively.The I u class includes asymptomatic, sub-clinical cases, and individual who are covered by the health programs.Individuals in PKDL class, P h , recover because of treatment and naturally at rates of σ 2 and γ 2 , respectively, and die naturally at rate µ h .Recovered individuals, R h , die naturally at rate µ h and become susceptible again at an average rate of η attributable to diminishing VL-induced immunity.Susceptible reservoirs acquire VL infection from infected vectors by force

Table 1 .
Description of the model's parameters.

Table 2 .
Estimates for the initial values.

Table 3 .
Partial rank correlation coefficient r and p-value of each parameter with both I d and R 0 .