Optimal combinations of control strategies and cost-effective analysis for visceral leishmaniasis disease transmission

Visceral leishmaniasis (VL) is a deadly neglected tropical disease that poses a serious problem in various countries all over the world. Implementation of various intervention strategies fail in controlling the spread of this disease due to issues of parasite drug resistance and resistance of sandfly vectors to insecticide sprays. Due to this, policy makers need to develop novel strategies or resort to a combination of multiple intervention strategies to control the spread of the disease. To address this issue, we propose an extensive SIR-type model for anthroponotic visceral leishmaniasis transmission with seasonal fluctuations modeled in the form of periodic sandfly biting rate. Fitting the model for real data reported in South Sudan, we estimate the model parameters and compare the model predictions with known VL cases. Using optimal control theory, we study the effects of popular control strategies namely, drug-based treatment of symptomatic and PKDL-infected individuals, insecticide treated bednets and spray of insecticides on the dynamics of infected human and vector populations. We propose that the strategies remain ineffective in curbing the disease individually, as opposed to the use of optimal combinations of the mentioned strategies. Testing the model for different optimal combinations while considering periodic seasonal fluctuations, we find that the optimal combination of treatment of individuals and insecticide sprays perform well in controlling the disease for the time period of intervention introduced. Performing a cost-effective analysis we identify that the same strategy also proves to be efficacious and cost-effective. Finally, we suggest that our model would be helpful for policy makers to predict the best intervention strategies for specific time periods and their appropriate implementation for elimination of visceral leishmaniasis.


Introduction
Leishmaniasis is the world's second largest parasitic killer and is caused by protozoan parasites belonging to the Leishmania genus. There are four different clinical manifestations of the disease − visceral leishmaniasis (VL), the post-kala-azar dermal leishmaniasis (PKDL), cutaneous leishmaniasis (CL) and cutaneous leishmaniasis with involvement of lesions of the mucous membranes, which is also called mucocutaneous leishmaniasis (MCL). Visceral leishmaniasis (VL), sometimes usually referred to as "kala-azar" (KA) is the deadliest form of leishmaniasis and is the causal reason for about 20,000 to 40,000 deaths worldwide, with total reported VL cases between 200,000 to 400,000 [1]. Visceral leishmaniasis (VL) has been targeted by the WHO for elimination as it is fatal, if left untreated. According to WHO, most of the cases of VL have been reported in the Indian subcontinent, Sudan, Ethiopia, and Brazil [2].
Apart from the heterogeneity in the clinical manifestations, the problem is further complicated by the presence of "asymptomatic" infections where the patients do not display any symptoms of the disease [3] and hence, their detection poses a challenge. In addition, some VL treated patients (6 months to several years after the treatment regimen) show a macular, maculopapular, and nodular rash that contain dormant parasites [2,4,5]. These individuals are themselves recovered, but serve as an active source of new infection when exposed to vectors. Such individuals are referred to as post-kala-azar dermal leishmaniasis (PKDL) infected.
The severity of this disease at such different scales thus, demands for efficient implementation of disease intervention strategies that can limit the spread of the infection among human populations. Recent studies on VL intervention strategies can be broadly classified under four groups; studies relating to-animal reservoir control, vector population control, human reservoir control and finally a group that includes studies where multiple interventions were conducted concurrently [6]. Animal reservoirs are eliminated from infection spread by either culling them [7], or use of canine vaccines [8] and insecticide impregnated collars [9]. The sandfly vector populations are controlled by spray of insecticides [10,11] and use of treated bednets [12]. As part of human reservoir control, drug-based treatment of infected individuals is one of the large scale programmes undertaken by WHO to reduce the incidence of VL cases in highly concentrated areas [13].
Although these methods are the most commonly used, there are problems and issues with the use of each aforementioned strategy. Lack of recording the actual number of reservoirs in a given area poses a major challenge for animal intervention strategies [14,15]. An important problem of insecticide sprays or insecticide treated bednets is the increasing resistance of sandflies to insecticides like DDT and deltamethrin [16,17]. Similarly, continued treatment of infected individuals with antimonials and miltefosine lead to drug resistance within the parasites, where patients stop responding to drug treatments [18,19]. The issues of implementation of strategies, and the associated problems with each of these strategies, further emphasizes the requirement of novel intervention strategies or the use of effective combinations of the existing strategies to target the elimination of the disease.
Mathematical models provide a detailed framework to study, analyze the dynamics of VL disease transmission and further contribute towards optimal choice of intervention strategies to control VL spread [3]. Historically, epidemiological models have long been proposed to understand leishmaniasis disease transmission and its control. The inter-epidemic periods between 1875 and 1950 in Assam, India was studied by Dye [20] using a deterministic model to describe the dynamics of VL. This model was extended to canine VL in Malta [21] to explain the efficacy of various control methods [22]. After this pioneering work, many mathematical models have been employed to understand the VL transmission dynamics [1,[23][24][25][26][27][28][29][30] but only few articles attempt to describe the detailed VL disease transmission dynamics [1,[25][26][27][28][29][30].
Transmission of VL can also affect animals. Leishmania donovani has been found in animals in East Africa; in Brazil, high levels of infection occur in dog populations (canine or CVL) [3]. Through anthroponotic medium or the zoonotic medium; visceral leishmaniasis can be transmitted between human to human or between animal and human respectively [5] via sandflies [31]. Burattini et al. [32], developed an SEIR type model between sandfly, animal and human populations for zoonotic transmission of visceral leishmaniasis. Following up, many other models [1, 25, 28-30, 33, 34] consider zoonotic transmission along with potential PKDL progression rate in humans via the addition of another infective stage. It is important to note here that most of these studies were dedicated to understand the visceral leishmaniasis transmission and do not directly focus on control of disease spread. The model by Stauch et al. [28] predicts the role of long-term intervention strategies, combined with active case detection and include efficacious treatment. Later ELmojtaba et al. [34] also, formulated a compartmental model and applied two optimal controls, namely treatment and vaccination within the model to investigate optimal strategies for controlling the spread of the disease. A recent mathematical model [35] attempted to understand the effect of specific optimal strategies for controlling anthroponotic cutaneous leishmaniasis in human populations; although this model only focuses on individual control strategies which are ineffective in controlling disease spread and is also not validated with real population data.
With this background, we modify the model of ELmojtaba et al. [25] and propose a detailed compartment-based mathematical model to explain the anthroponotic VL transmission dynamics in three distinct populations-the human and the animal as hosts, and sandfly as the vector, with the detailed analysis of the human infected population into its clinically distinct classes namely, asymptomatic, symptomatic, and PKDL infected. Adding to the complex nature of transmission, the disease also can have a seasonal fluctuation depending on the vector population [31]. To capture the presence of possible seasonal pattern observed within the data of reported VL cases, we include a periodicity factor in our model. We analyze this model to derive conditions for positive invariance, global stability of the unique disease free equilibrium, expression for the basic reproduction number in a periodic environment, and the conditions for the existence and permanence of at least one positive endemic periodic solution. We attempt to fit the model outcomes to the number of new VL cases occurred in South Sudan to estimate model parameters. We also perform parameter sensitivity analysis to identify the most sensitive parameters in our model.
Although the aforementioned studies do propose intervention techniques for disease control; none of the studies previously investigate the effects of the interventions in the disease dynamics while accounting for periodic seasonal fluctuations, their efficacy and costeffectiveness which may sometimes be limited by availability of resources. We further propose combinations of known preventive measures, namely, (i) treated bednets, (ii) treatment of infective humans and (iii) spray of insecticides and compare them with each other to understand the short and long term effects of the interventions on infected vector and human populations. We also use optimal control theory to study the efficacy of different strategic combinations of VL interventions in disease elimination as well as their cost-effectiveness. We explore the effects of proposed combinations and calculate the Infection Averted Ratio (IAR) with the Incremental Cost-Effectiveness Ratio (ICER) to investigate the efficacy of the combinations to eliminate disease from the population and their cost-effectiveness. Specifically, we emphasize that, by carrying out such a comparative analysis; predicting the involved costs and the corresponding outcomes of alternative control strategies can be useful to decision makers, who are often faced by the challenge of resource allocation. An optimal balance among different types of interventions may significantly reduce the number of VL cases and deaths at a minimal cost. This can lead to a well-coordinated effort and an effective implementation of strategies for disease control.

Model construction
The transmission dynamics of VL in the Indian subcontinent was modeled by a system of ordinary differential equations as per [1,28,29]. Following this framework, we consider the basic SIR type model with respect to history of infection within the human population. The model comprises of the human, reservoir and sandfly populations with seasonally forced biting rates on the sandfly population. Similarly, let the reservoir host population be divided into two categories, susceptible reservoir, S R (t), and infected reservoir, I R (t), such that The total vector (sandfly) population, denoted by N V is subdivided into susceptible sandflies S V (t), and infected sandflies I V (t), such that All the humans initially remain susceptible to infection and are assumed to grow in number with a constant birth Λ H and death μ h rate. After being bitten by an infectious sandfly, susceptible humans (S H ) are considered to become asymptomatically infected (I A ) with force of infection ab I V N H where a is the mean rate of bites per sandfly and b is the sandfly to human (reservoir) transmission probability. The per-capita biting rate of sandflies a is equal to the number of bites received per human from sandfly due to conservation of bites mechanism. Asymptomatic stages can include those humans with sub-symptomatic and non symptomatic early infection. Here, the asymptomatic stage (I A ) describes the subset of all those humans capable of contributing towards disease transmission. If alive, they remain asymptomatic for 1 g h days. A fraction of these individuals (ρ 1 ) develop symptomatic KA (I H ), some (ρ 2 ) become PKDL-infected (P H ) and the remaining ρ 3 = 1 − ρ 1 − ρ 2 recover from the asymptomatic infected stage (I A ). Symptomatic humans (I H ) are eligible for treatment. These individuals die due to VL at an average rate δ, or get treated at an average rate α 1 . A proportion of patients σ, from I H stage successfully combat the parasites and recover from the infection (R H ). The remaining proportion (1 − σ) of patients putatively enter into the dormant stage (T H ); followed by a PKDL (P H ) infection if on an average they survive upto 1 d p days [3]. Humans with PKDL get treated at an average rate α 2 , or recover naturally at an average rate β. Cellular immunity remains for 1 p r days, after which recovered humans (R H ) again become susceptible (S H ). Susceptible reservoirs are recruited from the population at a constant rate Λ R , and acquire infection VL following contacts with infected sandflies at a rate ab I V N H where a and b as described above. We assume that the transmission probability per bite is the same for human and reservoir because sandflies do not distinguish between humans and reservoirs. It is also assumed that reservoirs do not die due to the disease, but are limited by a per capita natural mortality rate μ r . Susceptible sandflies are recruited at a constant rate Λ V , and acquire VL infection following contacts with human infected with visceral leishmaniasis(I A and I H ) or human having PKDL (P H ) or reservoir infected (I R ) with visceral leishmaniasis at an average rate equal to where a is the per-capita biting rate, and c is the transmission probability for sandfly infection. The infection probabilities of sandfly depend on the stage of infection [28]. We assume that μ 1 and μ 2 are the respective infection probabilities of sandfly for biting humans and reservoir in the stages I A , I H , P H , I R . Sandflies suffer natural mortality at a per-capita rate μ v regardless of their infection status.
Abubakar et. al. [36] described that, the monthly distribution of VL cases reflected the general pattern observed in South Sudan, with less cases during the transmission season (April to June) and a peak during the dry season, beginning continuously in September. Hence, we assumed a periodic form for the biting rate as follows: aðtÞ ¼ a 0 ð1 þ d r sin 2pt 12 Þ. The biting rate a(t) of the sandfly population varies periodically with different temperatures which is assumed to be time periodic for a period of 12 months. a 0 denote the average biting rate and δ r denotes the amplitude of seasonality [37][38][39][40].
With this assumption and the description of the terms, we get the following system of differential equations:

Model properties
Let C denote all continuous functions on the real line. Given f 2 C, and if f is ω-periodic, then the average value of f on a time interval [0, ω] can be defined as: The maximum and minimum values of f on a time interval [0, ω] is denoted as f M and f m , respectively and defined as and All parameters of the model (1) are assumed to be nonnegative. Furthermore since the above model monitors living populations, it is assumed that all the state variables are nonnegative at time t = 0. It is noted that in the absence of the disease (δ = 0), the total human population size, This shows that the biologically-feasible region: Here, m is taken as a constant because it is well known that a vector takes a fixed number of blood meals per unit time independent of the population density of the host. Similarly, we let n ¼ N V N R be the female sandfly vector reservoir ratio defined as the number of female sandflies per reservoir host.
Disease free equilibrium calculation and its stability analysis, mathematical description of basic reproduction number (R 0 ), existence and permanence of endemic periodic solution and its global stability has been discussed as separate sections in Section A in S1 File.

Model calibration (Case study)
To estimate model parameters, we use the monthly VL infected case records reported in South Sudan for the year 2012 [36]. The estimated parameters from the data are - The initial human demographic parameters S H (0), I H (0), P H (0), T H (0), R H (0) as well as initial infected reservoir population I R (0), initial infected sandfly population I V (0) were also estimated.
We assume that The carrying capacity (N V ) of the sandfly population is taken to be a multiple of the total human population at the beginning, i.e. N V (0) = k 2 × N H (0), where k 2 is the total number of sandfly per human. Similarly, initially N R (0) = k 1 × N H (0), where k 1 is the total number of reservoir per human. We estimated k 1 and k 2 from the given data of Visceral Leishmaniasis. Initial susceptible sandfly population S V (0) = N V (0) − I V (0) and initial susceptible reservoir According to the [36], 28,328 new cases of VL were reported from September 2009 until December 2012. Relapses represented 8.6% cases in 2012 and PKDL was noted in 4.6% of patients in 2012. We added a compartment I C to our model (1) to calculate the cumulative number of new notified VL infections. The number of new VL cases I C (symptomatic infected I H +PKDL infected P H +Transient T H ) from the model (1) has the following form which represents the rate of cumulative new VL cases. Here, While simulating our model (1) with above mentioned assumptions, the numerical solutions for the above equation gives the predicted monthly cumulative VL incidence. For performing numerical simulations, ode15s and ode45 MATLAB ODE solvers were used. Here, I C (0) = number of new notified cases at the first time point of the data (C(0)).
We minimize the sum of the squared error between the model and data, which is given by where Y i is the cumulative VL data, and Y i ¼ I C ðt i ; b yÞ þ , * N(0, Iσ 2 ) where t i = 0, 31, 60, 91. . . days [39], and be the error of fit, which follows an independent Gaussian distribution having unknown variance σ 2 . The prior distribution can be viewed as representing the current state of knowledge, or uncertainty, about the model parameters prior to data being observed. From previous literature, we observe that all our unknown parameters are non negative and bounded. It is more realistic to assume that the prior distributions are "bell curve" shape rather than flat one. So, an independent Gaussian prior specification is assumed for the unknown parameters b y [41] of the model (1) i.e. θ j * N(ν j , ψ j ), where j = 1,2,3,. . .. . .N. Realizations of Gaussian processes with a proper covariance function can provide nearly all functions we can encounter in "real life". Also, they are convenient and provide exact inference and marginal distribution.
We also assume that the inverse of the error variance follows a gamma distribution as prior with the following form: where S 2 0 , n 0 are the prior mean and prior accuracy of σ 2 , respectively.
Using conditional conjugacy property of Gamma distribution, the conditional distribution is rðs À 2 jY; b yÞ also a Gamma distribution with This conditional conjugacy property makes it possible to sample and update within each Metropolis-Hastings simulation step for the other parameters. Since, we assume independent Gaussian prior specification for b y, therefore we can now calculate the prior sum-of-squares for the given b y as: Then, for a fixed value of σ 2 , the posterior distribution of b y is as follows: and the posterior ratio needed in the Metropolis-Hastings acceptance probability can be written as: MCMC tool box in MATLAB version R2011b was used to estimate the unknown b y for the model [41,42]. Geweke's Z-scores [43] were used to ensure the chain convergence (Table D in S1 File).
The estimated model parameters, including human, reservoir and sandfly demographic parameters, for South Sudan in 2012 are given in Table A and B of S1 File. Further, plots for the posterior distributions of the estimated parameters, including human, reservoir and sandfly demographic parameters of the model (1) are given in S1.A Fig. The basic reproductive number (R 0 ), the expected number of secondary cases produced by a single infection in a completely susceptible population, was also calculated for each parameter set so as to get a posterior distribution of R 0 (S1. B Fig). Also, trace plots of parameters (S1. C Fig) were generated to identify whether the Markov chain has converged or not.

Model sensitivity analysis
Using the standard combination of Latin Hypercube Sampling (LHS) and Partial Rank Correlation Coefficient (PRCC) multivariate analysis, we performed the sensitivity analysis for the model (1). LHS is a stratified Monte Carlo sampling method, where the random parameter distributions are divided into N equal probability intervals and samples are taken from each [44], where N is the sample size. PRCC is an efficient method for measuring the nonlinear, monotonic relationship between inputs and the model outcome of interest.
The inputs are the estimated parameters as given in the "Model calibration" section and as described "Parameter Estimation" of Section C in S1 File. The sensitivity analysis was performed for 2000 random samples of each parameter in the model for ranges given in Table A of S1 File. The details of sensitivity analysis are discussed in the "Sensitivity Analysis-Results" section of Section C in S1 File. Whereas, the model outcome is the cumulative proportion of infectious individuals, which is the solution of dI C dt ¼ r 1 g h I A þ ð1 À sÞa 1 I H þ d p T H þ r 2 g h I A and R 0 .
The optimal control problem Construction of the problem. We incorporated three different popularly used VL intervention strategies in the model (1) namely, the use of treated bed nets, treatment of infective individuals using antibiotics and spray of insecticides [6]. The control functions, u 1 , u 2 , u 3 and u 4 ; represent time dependent efforts of use of treated bednets, treatment of symptomatic KA patients, treatment of PKDL patients and spray of insecticides respectively. The controls are practiced on a time interval [0; t f ], where t f is the final time.
System of non-linear differential equations representing the effect of different interventions on our basic model (1) is given as follows:  Table A and Table B of S1 File. The control functions u 1 (t), u 2 (t), u 3 (t) and u 4 (t) are bounded and Lebesgue integrable functions.
Our control problem involves that in which the number of infected individuals with visceral leishmaniasis and the cost of applying different controls are minimized subject to the system (5). The objective function to be minimized is defined as subject to the state Eq (5) and the total cost is given by where t f is the final time, σ 1 is the discount rate applied to future years and A 1 , A and A 2 are weight constants of the I A , I H , P H group, respectively, whereas, B, C, D and E are weight constants for treated bed nets, treatment (for I H , P H ) and spray of insecticide efforts respectively which regularize the optimal control. The discounting procedure reflects inherent uncertainty about the future. It is assumed that, there is no linear relationship between the coverage of these interventions and their corresponding costs [45]. We seek an optimal control u Ã 1 ðtÞ, u Ã 2 ðtÞ, u Ã 3 ðtÞ and u Ã 4 ðtÞ such that where Analysis of the optimal control problem. The necessary conditions that an optimal control must satisfy, come from the Pontryagin's Maximum Principle [46]. The Hamiltonian H, with respect to u 1 , u 2 , u 3 and u 4 can be written as: where g i is the right hand side of the differential equation of the i th state variable and λ i are the adjoint variables. By applying Pontryagin's Maximum Principle [46,47] we get the following result: Proposition 1 There exists an optimal control u Ã 1 , u Ã 2 , u Ã 3 and u Ã 4 and corresponding solution, where Details of adjoint functions @l i @t is given in the "Optimal Control" Section A of S1 File. Numerical simulations. The optimality system is a two-point boundary problem, because of the initial condition on the state system (5), and the terminal condition on the adjoint system (10). Using the values for the parameters given in Tables A and B in S1 File, and the initial conditions in Table C in S1 File, we first solve the initial valued state system (5) forward in time, using an initial guess for the above defined control functions (11)(12)(13)(14). Then, using the same initial guess for the control functions, we solve the adjoint system (10) with terminal conditions backward in time. The controls are updated at each iteration using the optimality conditions (11)(12)(13)(14) The iterations continue until the system converges.

Cost-effectiveness analysis
Although the choice of optimal combinations of control strategies and need for a proper implementation is important, it is also equally important to choose a combination that is also cost-effective while implementing on a large scale. Cost-effectiveness analysis is a method to evaluate the benefits associated with interventions in infections (treatment of infected individuals, use of treated bednets and spray of insecticides)with respect to the strategy's involved cost [45]. To calculate the cost-effectiveness of our strategies, we used two methods, namely the infection averted ratio (IAR) and the incremental cost-effectiveness ratio (ICER).
IAR. The infection averted ratio (IAR) can be defined as the ratio of number of infections averted to the number of recovered. The number of infection averted is given as the difference between the total number of infective individuals present without control and total infective individuals present with control.
ICER. The incremental cost-effectiveness ratio (ICER) can be defined as the additional cost per additional health outcome. ICER is an incremental ratio of the difference in total cost between one strategy and the next best alternative to the difference in total number of averted infections through each strategy.

Model fitting and validation
Monthly VL infected cases reported in South Sudan for the year 2012 [36] was chosen to estimate our model parameters. The data shows a peak during months of Jan-Feb, after which the reported VL cases decline to lower numbers ( Fig 1A). This data was chosen such that the model predictions using the estimated model parameters almost accurately fit with the cumulative number of VL cases reported for the year 2012. The cumulative number of cases are almost accurately predicted by the model for the year 2012 ( Fig 1B). Further, using the fitted model, predictive simulations for the next year Jan 2013-Dec 2013 was also performed to signify the use of the model to predict future VL cases. Our model predicted that cumulatively 4,394 and 6,587 persons are diagnosed as new KA cases at the end of December 2012 and December 2013 respectively. Thus, approximately 6,587−4,394 = 2,193 new KA cases were predicted for the year 2013. WHO [48] reported that around 2,364 new VL cases were reported in South Sudan for the year 2013. The initial values for the simulations were taken according to the reported yearly census data of South Sudan. The total human populations (N H ) obtained after simulations for the 2012 and 2013 cases are indicated in S6 Fig and is suggestive of around 4.3% yearly growth of population in South Sudan, which is also comparable with the reported population growth (approximately 4% yearly). The cumulative number of predicted VL cases from our model were comparable to the actual reported VL cases in South Sudan for 2013 further substantiating the predictive power of the model. Note that the above mentioned simulations were performed assuming a periodic sandfly biting rate due to the periodic nature of the data, which was further verified using seasonality test. Performing the seasonality test, the null hypothesis that "there is no effect of seasonality on disease dynamics" was rejected at a significance level of p<0.05, thus, supporting the notion of seasonality influence on VL disease dynamics, modeled as a periodic sandfly biting rate, albeit with low amplitude of seasonality (see Section B in S1 File for details).

Choice of optimal control strategy combinations
The most popularly used VL intervention strategies include the treatment of infected, use of treated bed nets and spray of insecticides for vector control [6]. Initially, each control strategy and its effect on VL disease transmission was investigated individually. S3 S4 and S5 Figs display the model-based predictions of effects of different intervention strategies on a yearly basis and its comparison with a scenario where no optimal control is introduced in the population. Simulations were performed on the previously standardized model with estimated parameters (Tables A and B in S1 File) and initial values (Table C in S1 File) to compare the normal and optimally treated scenarios. Model predictions suggest that the treatment policy is very useful to optimally control KA infected individuals (S3. D Fig). As the number of humans harboring infection reduce due to medical treatment, susceptible sandfly vectors are also deprived of infective human hosts for their blood meal and hence, do not become infected. The amount of infected vector populations thus, reduce exponentially, although the rate at which they decrease is low (S3. E Fig). The optimal use of treated bednets on the other hand performs even worse as sandfy populations reduce only intermittently, after which the rate of infected vectors increase(S4. E Fig). Similarly, optimal use of insecticides are relatively weaker in controlling infected populations as compared to both the aforementioned policies (S5 Fig). It is important to note here that none of these policies are by themselves effective enough in reducing VL infection spread and hence, require better strategies to combat the VL infection. We propose from our model that instead of using each of these strategies separately, optimal combinations of the aforementioned control strategies would be better alternatives to reduce VL disease spread. So as to test this contention, we introduced different optimal combination strategies in our model and numerically compare their effects on infected populations.
We devise and test for the following combinations: 1. Strategy A: combination of use of treated bednets, treatment of infective individuals and spray of insecticides.

Strategy B: combination of use of treated bednets and treatment of infective individuals.
3. Strategy C: combination of use of treated bednets and spray of insecticides.

Strategy D: combination of treatment of infective individuals and spray of insecticides.
Strategy A. Under this strategy, we use all the four controls u 1 , u 2 , u 3 and u 4 to optimize the objective function J. Comparing Strategy A with a situation where no control strategy was used, it was observed that the susceptible and recovered human populations (S H + R H ) increase in number (Fig 2A), asymptomatic KA populations (I A ) marginally reduce (Fig 2B and 2C), symptomatic KA and PKDL (I H + P H ) reduce at an exponential rate (Fig 2D), and infected sandfly populations (I V ) (Fig 2E) decreases significantly at a near exponential rate where sandfly populations reduce below 1000 within 15 days. At t = 100 days, comparing Strategy A with no strategy, there is an increase in S H + R H by 9975 individuals, decrease in I A by 2409 individuals and I V by 9694 individuals (Fig 2) respectively. With this strategy, I H + PKDL will be eliminated from the system within t = 50 days in humans. The control profile shown in Fig 2F, shows that the control u 1 is at 10% initially, after that it drops slowly to the lower bound while control u 2 and u 3 decreases from the maximum of 100% to the lower bound in 98 days and the control u 4 is at 36% for beginning time before dropping slowly to the lower bound in the 92 th day. This suggests that a low effort is required on treated bednet and insecticide spray under this strategy.
Strategy B. In this strategy, the treated bednet control u 1 and the treatment control u 2 , u 3 are used to optimize the objective function J while we set the spray of insecticides control u 4 , to zero. Comparing Strategy B with a situation where no control strategy was used, it was observed that the susceptible and recovered population (S H + R H ) increase in number (Fig 3A), asymptomatic KA populations (I A ) marginally reduce (Fig 3B and 3C) suggesting that asymptomatics become symptomatic after which they are cured, and symptomatic KA, PKDL (I H + P H ) cases reduce at an exponential rate (Fig 3D) very similar to Strategy A. Infected sandfly populations (I V ) (Fig 3E) decreases significantly after introduction of Strategy B, albeit at a  rate slower than Strategy A (Fig 3E). At t = 100 days, comparing Strategy B with no strategy, there is an increase in the S H + R H by 8991 individuals and decrease in I A by 1429 individuals, I H + P H by 3596 individuals, and I V by 5430 individuals respectively (Fig 3). In Fig 3F, the control u 2 , u 3 is at the upper bound of 100% and drops gradually until reaching the lower bound, while control on treated bednets u 1 is at the maximum of 32% initially before dropping gradually to the lower bound in the 98th day. This suggests that there is a low effort for the use of treated bednets and a higher effort required for medical treatment of individuals under this strategy.
Strategy C. Here under this strategy treated bednet control u 1 and the spray of insecticide control u 4 are used to optimize the objective function J while we set treatment control u 2 , u 3 = 0. Comparing Strategy C with a situation where no control strategy was used, it was observed that the susceptible and recovered population (S H + R H ) (Fig 4A and 4B) increase and the numbers of infected humans I A (Fig 4C and 4D), I H + PKDL (Fig 4E and 4F) with optimal strategy show only a marginal decrease as compared to the numbers in the case without control. The number of infected humans are reduced by a low amount in this strategy when compared with Strategy A and B. Infected sandflies I V (Fig 4G) reduce to considerably low values which is better than any other combination strategy. At t = 100 days, comparing Strategy C with no Optimal control of visceral leishmaniasis disease transmission strategy, S H + R H increase by 2511, I A decrease by 2426, I H + P H decrease by 59 and I V decrease by 9553 individuals respectively (Fig 4). The control profile is shown in Fig 4H; we see that both the control u 1 and u 4 start from their upper bounds steadily decline and converge to its lower bound within 96 days. This suggests that there is a low effort involved for the use of both treated bednets and spray of insecticides when applied in combination.
Strategy D. Under this strategy, we optimize the objective function J using the treatment control u 2 , u 3 and the spray of insecticides controls u 4 while the treated bed net control u 1 = 0. Comparing Strategy D with the no control situation, it was observed that the susceptible and recovered population (S H + R H ) significantly increased Fig 5A, the asymptomatic population I A (Fig 5B and 5C) also decreased, I H + PKDL (Fig 5D) and infected sandflies I V (Fig 5E) also reduce tremendously. Strategy D decreases the infected human and infected sandfly population drastically while increase the number of susceptible and recovered human population (Fig 5). At t = 100 days, comparing Strategy D with no strategy, S H + R H increase by 9974 individuals, I A decrease by 2405 individuals and I V decrease by 9696 individuals respectively. With this strategy, I H + PKDL will be eliminated from the system within t = 50 days in humans. Fig 5F highlights the optimal control profile for strategy D. It was observed that the treatment control u 2 , u 3 starting from its upper bound drop at a very slow rate and reach its lower bound taking around 96 days while u 4 is at the maximum of 40% before dropping to the lower bound in the 95 days. This again suggests that the spray of insecticide controls require less effort for implementation when used in combination with treatment policies.

Effect on R 0
The above applied optimal control strategies applied in our model are unable to predict the long term dynamics of the disease as a whole. As and when the intervention strategies are stopped, a few remaining infectious people/vectors can initiate a new outbreak of the disease [28]. In the context of epidemiology, the basic reproduction number (R 0 ), that describes the number of secondary infections arising from a single individual during period of infection, is Optimal control of visceral leishmaniasis disease transmission an effective measure for understanding long term endemicity [49]. Hence, the effect of applied strategies on R 0 was investigated. Fig 6 describes the numerical simulation results of R 0 under different control strategy combinations. Assuming that optimal control combinations are implemented in the beginning of the year; in the early stages, strategy C performs well and suppresses the spread of disease, strategy B performs almost similarly throughout and do not affect R 0 indicating long term persistence of disease for these two combinations, followed by strategies A and D. Eventually after completing 100 days, as the provided optimal control has reached its lower bound, R 0 increases displaying disease persistence (Fig 6B). This analysis provides us with the fact that implementation of intervention strategies at different time points leads to different disease dynamics. Hence, it is important to determine the time point at which optimal control needs to be implemented. Also, it further points to the fact that complete elimination of the disease by any control strategy combination can be achieved only if the control strategies continue for long periods of time (determined by the upper bounds of the controls u 1 , u 2 , u 3 and u 4 ).

Cost-effectiveness of the suggested strategies
In order to understand the cost-effectiveness of each of the optimal combinations, ICER and IAR were calculated for each strategy (see "Cost-effective analysis" in Methods section). The cost involved for the treatment strategies on a per person basis are given in Table B of S1 File. The comparisons of the calculated ICER and IAR for the different scenarios are given in Table 1. The lowest ICER and highest IAR for strategy D indicates that strategy D outperforms all other intervention strategies with respect to cost-effectiveness and efficiency. Strategy B is the next best strategy as it has the second lowest ICER and second highest IAR. Strategy A is also effective in eliminating the disease but has highest involved cost for implementation (highest ICER), and hence, is the least cost-effective strategy. Strategy C performs worst with

Conclusion
Visceral leishmaniasis, being a deadly tropical disease requires community efforts and largescale elimination programmes for effective control of disease spread among populations. Although efforts are being taken, the choice of the best strategy or combination of strategies, that largely depends on the regional, seasonal and temperature variations still remains a formidable task at hand [6]. In this paper, we introduce a general non-autonomous anthroponotic visceral leishmaniasis model that considers the human (infected compartments divided into symptomatic, asymptomatic, PKDL-infected classes) and sandfly populations and probe further to understand the effect of optimal control strategy combinations on the infected populations considered within the model. To capture the realistic situations, we considered the effect of seasonal variations on the biting rate of sandfly as an integrable periodic function in the model. Further, we mathematically analyse the model to derive the basic reproduction number R 0 for the model and show that the disease-free equilibrium of the proposed model is globally asymptotically stable if R 0 < 1. If R 0 > 1, then the proposed system has at least one positive periodic solution, and the solution is uniformly persistent. We have also proven that if R 0 > 1, then the positive periodic solution is globally asymptotically stable. We also numerically analyze the model by estimating parameters fitted for the VL cases reported in South Sudan for the year 2012 [36]. The estimated value of basic reproduction number (R 0 ) in periodic environment in South Sudan for the year 2012 was 2.67 (with 95% CI). We further use this standardized model to predict the number of VL infected cases reported for the year 2013 in South Sudan [48]. The model predictions were highly comparable with the number of actual reported cases at the end of 2013 substantiating the predictive capability of the model. Time dependent intervention strategies can be implemented to curtail a vector-borne disease on a finite time interval. Using optimal control analysis, we further investigated the effects of popular intervention strategies and their optimal combinations on infected human and vector populations in the model under periodic seasonal fluctuations, thereby depicting control in realistic situations. In this article, we consider three types of control, i.e. use of treated bednets, treatment of infectives and spray of insecticides. Tracking the time-dependent changes, it was observed that the combination of spray of insecticides & drug-based treatment of infected individuals (Strategy D) and the combination of treated bednet, spray of insecticides & drugbased treatment (Strategy A) performs well for the time period of intervention. To observe long term effect of control on disease spread, the effects of control strategies were observed on the basic reproduction number (R 0 ) for entire period of intervention. This analysis indicated significant changes in the number of possible secondary infections from an infected individual, with respect to implementation of intervention strategies at different time points of the infection in population. Hence, it is important to determine the time point at which optimal control needs to be implemented. Also, complete elimination of the disease by any control strategy combination can be achieved only if the control strategies continue for long periods of time. These results pose a realistic view of visceral leishmaniasis disease spread, its control and its resurrection, once the control strategies are removed.
As opposed to the previous model on cutaneous leishmaniasis [35], we have considered the model specific to anthroponotic visceral leishmaniasis by incorporating the symptomatic, asymptomatic and transient infectious classes separately and implemented combinatory optimal control strategies and cost-effective analysis to propose an efficient elimination programme. Further, our predictions for visceral leishmaniasis are completely opposite to the model-based control strategies of cutaneous leishmaniasis [35], which predict complete elimination of the disease contradicting the endemicity observed in the real disease scenario after removal of a intervention strategy. Moreover, we have validated our study with real disease incidences and have shown the predictive capability of our model. From our study, the effect of the different strategies on R 0 further indicated that the combination of treated bednets & spray of insecticides (Strategy C) performs the best to control disease, followed by strategies A & D which have comparable performance. Whereas, cost-effectiveness analysis using IAR and ICER, indicate that the combination of drug-based treatment of infective individuals and spray of insecticides (Strategy D) is the most optimal, cost-effective, and efficacious strategy followed by the combination of treated bednets and drug-based treatment control (Strategy B) to control disease dynamics. Strategy A also displays a high efficacy in eliminating the disease comparable to Strategy D, but is accompanied by a higher involved cost. Thus, we conclude that for cases where the severity of disease is low, it is favourable to choose both an efficacious and cost-effective strategy (Strategies B and D) whereas, when the intensity of the disease is high and the priority is to control the disease spread, strategies A or C, which are less-cost effective but immediately efficacious to control disease spread within a short span of time, can be applied. Thus, our model can be useful for decision makers, who are often faced by the challenge of resource allocation, to choose different control strategies with respect to the severity of the disease.
Supporting information S1 File. Appendix. This file contains the mathematical analysis and calculations for the disease free equilibrium and its stability analysis, mathematical description of basic reproduction number (R 0 ), existence and permanence of endemic periodic solution and its global stability from the model (Section A). Seasonality testing (Section B), parameter estimation and model sensitivity analysis (Section C). The file also contains the estimated model parameters and their description (Tables A and B) and initial values (Table C) of the variables considered in the model and Geweke's Z-score for each parameter (Table D)