Analysis of a Mathematical Model of Zoonotic Visceral Leishmaniasis (ZVL) Disease

Abstract

This research paper attempts to describe the transmission dynamic of zoonotic visceral leishmaniasis with the aid of a mathematical model by considering the asymptomatic stages in humans and animals. The disease is endemic in several countries. Data used in the research are obtained from the literature while some are assumed based on the disease dynamic. The consideration of both asymptomatic and the symptomatic infected individuals is incorporated in both humans and animals (reservoir), as well as lines of treatment for the human population. It is found that the model has two fixed points; the VL-free fixed point and the VL-endemic fixed point. Stability analysis of the fixed points shows that the VL-free fixed point is globally asymptotically stable whenever the basic reproduction number is less than one and the VL-endemic fixed point is globally asymptotically stable whenever the basic reproduction number is greater than one. Sensitivity analysis is conducted for the parameters in the basic reproduction number, and the profile of each state variable is also depicted using the data obtained from the literature and those assumed. The transmission probability from infected sandflies to animals, transmission probability from infected animals to sandflies, per capita biting rate of sandflies of animals, and rate of transfer from symptomatic infected animals to the recovered class are among the most sensitive parameters that have the greatest influence on the basic reproduction number. Moreover, the value of the basic reproduction number is obtained to be 0.98951, which may require further study, as the margin between potential disease control and outbreak is thin.

1. Introduction

Visceral leishmaniasis (VL) is one of the Neglected Tropical Diseases (NTDs). The protozoan parasites that cause leishmaniasis are spread by the bite of infected female phlebotomine sandflies. The disease has a connection to poverty and is linked to malnutrition, substandard housing, armed conflict, population displacement, illiteracy, gender discrimination, a weakened immune system, and a lack of resources [1,2].

Over 88 countries in the world have some areas where leishmaniasis is present. In these regions, there are about 350 million residents. In the tropics and subtropics are where most of the countries impacted by leishmaniasis can be found. Leishmaniasis can be found in rain forests in Central and South America as well as in deserts in West Asia. India, Bangladesh, Nepal, Sudan, and Brazil account for more than 90% of the global cases of visceral leishmaniasis. Mexico, Central America, and South America, from northern Argentina to Texas (but not in Uruguay, Chile, or Canada); Asia (but not Southeast Asia); the Middle East; and Africa (mainly East and North Africa, with some cases elsewhere) are the regions where leishmaniasis is found [3].

The primary reservoir of the parasite in these places is the domestic dog, which is spread by the phlebotomine sandfly, Lutzomyia longipalpis [4,5]. Upon introduction into a community, the parasite sustains a dog–insect–dog peri-domestic transmission cycle [6], wherein infected flies sporadically bite humans, resulting in zoonotic visceral leishmaniasis (ZVL). Because of this, the incidence and prevalence of ZVL are crucial epidemiological factors for limiting transmission [7], and the estimates of these parameters rely on the accurate identification of dogs that are affected [8]. Dogs that are affected are therefore frequently eliminated as part of ZVL control operations.

In the Mediterranean, ZVL is a veterinary and medical issue, where dogs are valued as “valuable” animals [9]. Human cases of ZVL are common in Brazil, where 90% of cases reported are in humans, dogs are viewed as “less valuable”, and ZVL is primarily regarded as a medical issue. Because of this, 850,000 dogs in Brazil undergo screening each year, and 20,000–25,000 [9,10] of those canines are put to death after receiving a positive test. This official position that dogs are undervalued is held by some governments, such as the Brazilian Health Authorities, which prohibit treating infected dogs with drugs intended for human use and mandate the mandatory elimination of dogs that test positive for drugs. Thousands of dogs are euthanized each year without a counterproof diagnosis because health officials threaten to fine dog owners severely if they refuse to have their animals’ serabiological analyses performed. Therefore, rather than having their dogs submitted to serological screening, owners either forbid health agents from entering their houses or move their pets to locations without serological screening. Some proprietors even avoid taking their dogs to veterinary doctors, fearing an examination could result in a death sentence for their animals [11].

The two parts of ethics that need to be stressed as regards leishmanisis are those pertaining to animals and human relations, where human health is at stake and animals and their link with human society are valued more. It is important that these ethic issues are addressed so that they do not impede the effort towards the control of leishmaniasis.

During World War II, leishmaniasis and sandfly fever were highly prevalent among soldiers stationed in the Persian Gulf region. This region saw the deployment of some 697,000 US troops during the Gulf War (1990–1991). In this cohort, only 12 cases of visceral and 19 cases of cutaneous leishmaniasis were found to be diagnosed. The application of insecticides and repellents, reduced summertime transmission rates, and increased urbanization all contributed to the improvement [5,6]. It is estimated that in 2003, approximately 150 American soldiers fighting in Iraq were diagnosed with leishmaniasis, and additional cases are anticipated. Preliminary information on 22 instances of leishmaniasis treated at Walter Reed Army Medical Center that were contracted by US troops in Afghanistan, Kuwait, and Iraq between August 2002 and September 2003 was recently released.

In order to effectively combat visceral leishmaniasis, early detection and proper treatment are essential. Visceral leishmaniasis symptoms and signs are non-specific; thus, a diagnosis is only made by combining clinical symptoms and laboratory testing that are specific to Leishmania. The level of the health system determines the diagnostic approach for medical services in endemic areas. For application in the field in the majority of endemic locations, two serological tests—the Direct Agglutination Test (DAT) and the rK39 antigen-based immunochromatographic tests—have been developed [12]. A Recombinant DNA Technology (RDT) for visceral leishmaniasis is an easy test that may be used both peripherally and centrally. It identifies antibodies. According to a scientific evaluation of published data, the sensitivity of RDTs varies with the eco-epidemiological regions, particularly in East Africa, where it is low. The drawbacks of all serological tests include their inability to reliably diagnose relapse and the fact that a sizable portion of healthy individuals living in endemic regions without a history of visceral leishmaniasis are positive for antileishmanial antibodies as a result of asymptomatic infections. These limitations are caused by the persistence of antibodies for prolonged periods after cure. Because of this, it is necessary to diagnose visceral leishmaniasis by combining antibody-based testing with a consistent clinical case definition.

Pentavalent antimonials have been the first-line treatment for VL for more than 70 years, although there are currently four medications available for the condition. This treatment takes 20 to 30 days, is toxic (3–5% deaths due to treatment), and is accompanied by increasing failure rates as noted in [13]. The first oral medication against VL is miltefosine, which has been recommended as a first-line medication in the VL elimination initiative; however, because of its long half-life, it is teratogenic and may quickly lead to resistance. Amphotericin B is used in two formulations: “Conventional” amphotericin B and “Liposomal” amphotericin B (AmBisome), and finally, paromomycin (PMM), which was registered in India in 2006, and is currently being tested in a phase IV trial [8,14].

An in vitro point-of-care test is necessary to confirm or rule out active cases for early diagnosis since untreated cases operate as reservoirs for infection and endanger the community to contracting leishmaniasis. Comparably, to ascertain whether therapy for visceral leishmaniasis has been successful, a laboratory test is necessary [15].

Despite the fact that model projections are predicated on precise parameter estimations and scenarios, actual conditions vary greatly, contributing a variable complexity that is difficult to measure. Numerous features of VL disease remain poorly understood, and its elimination has been identified as a public health concern [16]. Specifically, the progression of the disease through many stages, including PKDL and asymptomatic infection, remains unclear, particularly with respect to infectivity and diagnosis. Since recent efforts to manage VL appear to be working, the 2020 deadline for ending VL as a public health concern was rescheduled to 2017 [4,16]. However, modeling and studies of gearbox have indicated a likely scenario. However, modeling and transmission studies have raised the conceivable potential that those who are asymptomatic could obstruct the removal of the infection or hide the true number of leishmaniasis infections. Leishmaniasis infections are a major worry for any group in societies worldwide. Many studies have presented models to understand the disease, as well as the importance of eliminating VL infection independent of population-level asymptomatic infection classifications.

Using ordinary differential equations, a compartment-based mathematical model of zoonotic visceral leishmaniasis transmission and its control across three distinct populations—human, animal, and sandfly—was created in [17]. Asymptomatic, symptomatic, post-kala-azar cutaneous leishmaniasis, and transiently infected individuals made up the human population. The study examined the effects of the asymptomatic stages, but it did not address the treatment phases for the infected human population.

The treatment approach for VL is centered on managing patients who have manifested symptoms of the disease. This implies that since laboratory testing is the only method available to identify those who are infected but have not yet displayed any symptoms, asymptomatic individuals will continue to transmit the disease unchecked. This indicates that the current treatment approach might not be enough on its own to bring the condition under control and sufficiently lower its occurrence. The ratio of symptomatic transmission to sandflies that receive medical attention is the main focus of the argument on disease transmission. If affected people are the main agent transmitting the disease, then treatment will have a major effect on additional transmission. Primary indirect evidence for this disease comes from the long-term cycles seen in many regions and the spatial clustering of leishmaniasis cases; both of these are best explained if Kala-Azar cases are the source of transmission. On the other hand, if asymptomatic infections account for most of the transmission, this means that early treatment of leishmaniasis may not significantly reduce transmission. Currently, there is less evidence to support this scenario, but it is not implausible if multiple asymptomatic infections propagate slowly to sandflies, increasing the rate of transmission to humans and animals. This emphasizes the necessity of conducting additional study and gathering precise data on the VL history and infectiousness relativity of various infection stages in order to enhance control strategies.

This work therefore will contribute to the methodology of vector control in the management of VL disease. Through this research, it is hoped that maximal benefits and new strategies will be developed, especially given that our model considers the critical states of infectiousness of the disease.

Furthermore, the study will lessen the lack of models and data available to address VL disease. It will also be very helpful in understanding how the population of asymptomatic animals, which typically impedes management of the disease, interacts with sandfly populations, which appear to be crucial to the eradication effort. Furthermore, it will advance our knowledge of the biology of sandflies and how transmission might alter over the course of a year. Future studies could make use of the fundamental frameworks established by those who have already studied the seasonality of sandflies in VL and the broader body of knowledge on the seasonal dynamics of other disease vectors.

In this paper, therefore, a comprehensive compartment-based mathematical model that examines the transmission process of VL in three different populations (humans and animals acting as hosts, and sandflies acting as the vector) is considered, taking into account the precise categorization of the human infected population into three distinct clinical classes (symptomatic, asymptomatic, and PKDL infected) and subdividing asymptomatic classes into early and late asymptomatic classes, and considering classes of lines of treatment. The model includes a total of 20 compartments (represented as variables) with regard to various demographic categories, and it represents the rate of disease transmission as a flow from one compartment to another. For this big system, the model is examined both analytically and numerically in order to determine the potential control mechanism and show the impact of the disease transmission amongst the different populations.

The paper is organized in this order: In Section 2, we formulate the model together with the description of the parameters defined in the model. In Section 3, we discuss the basic properties of the model. In Section 4, we obtain the fixed points, compute the basic reproduction number, analyze the local and global stability of the fixed points obtained, and conduct sensitivity analysis on some key parameters. In Section 5, we present the numerical analysis vis-à-vis sensitivity analysis and numerical simulation. Finally, in Section 6, we discuss the results and conclusion.

2. Model Formulation

In this research paper, a model of zoonotic visceral leishmaniasis disease incorporating lines of treatment for the infected human population is considered. Let $N_H\left(t\right)$ stand for the total human population, $N_A\left(t\right)$ for the total animal population (reservoir), and $N_F\left(t\right)$ for the total sandfly population (vectors). In order to represent a biologically realistic complex scenario among populations, the total human population is further divided into twelve compartments: susceptible $S_H\left(t\right)$, exposed $E_H\left(t\right)$, early asymptomatic infected $I_H^E\left(t\right)$, late asymptomatic infected $I_H^L\left(t\right)$, early recovered stage who are DAT-positive and not yet LST-positive $R_H^E\left(t\right)$, late recovered stage who are DAT-negative but still LST-positive $R_H^L\left(t\right)$, symptomatic infected $I_H^S\left(t\right)$, infected individuals who are receiving first-line treatment $I_H^1\left(t\right)$, infected individuals who are receiving second-line treatment $I_H^2\left(t\right)$, PKDL-infected $I_H^P\left(t\right)$, recovered humans who have cleared the parasite $R_H^1\left(t\right)$, and putative recovered human $R_H^2\left(t\right)$. The first and second lines of treatment, which aim to reduce the parasite population within each host, are the only ones that helped the people in this class recover. Similarly, the animal or reservoir populations are divided into five compartments: susceptible animal $S_A\left(t\right)$, exposed animal $E_A\left(t\right)$, asymptomatic infected animal $I_A^A\left(t\right)$, symptomatic infected animal $I_A^S\left(t\right)$, and recovered asymptomatic infected animal $R_A\left(t\right)$. The sandfly (vector) population is divided into three compartments: susceptible sandflies $S_F\left(t\right)$, exposed sandflies $E_F\left(t\right)$, and infected sandflies $I_F\left(t\right)$. The regular SEIR model is used to simulate the flow of infection between these compartments in the human and animal populations, whereas the SEI carrier-type model is used to simulate the flow of infection through the vector. Each arrow between the compartments in the model Figure 1 below indicates the passage of infection from one compartment to the next and illustrates the model structure incorporating the numerous components considered in the model. The arrows entering each susceptible population represent the birth rates of humans, animals, and sandflies $({\mathsf{\Pi }}_H,{\mathsf{\Pi }}_A,{\mathsf{\Pi }}_F)$, respectively, while the arrows moving out from each compartment that do not enter any compartment represent the death rates $({\mu }_H,{\mu }_A,{\mu }_F)$ of each individual compartment considered in the model. Using a standard incidence function as defined by [18] the interactions of susceptible humans and animals with the sandfly populations are modeled, considering the likelihood that an infected sandfly may spread infection to a susceptible human or animal host $({\lambda }_F,{\eta }_F)$ and also that the sandfly populations can contract infection by an infected animal with incidence functions $({\rho }_A,{\xi }_A)$, respectively. A fraction of exposed humans $\left({\theta }_1\right)$ become asymptomatic first at an early stage before late stage at a rate of ${\theta }_2$. After infection, exposed humans enter the early asymptomatic stage $I_H^E\left(t\right)$ and become PCR-positive in peripheral blood. As the specific antibodies develop in the blood serum due to the infection, the fraction of humans in the early asymptomatic stage move to the late asymptomatic stage $I_H^L\left(t\right)$, which is characterized by the onset of DAT-positivity. If not dying, humans remain in this asymptomatic stage. Infection with Leishmania parasites proceeds in most cases asymptomatically, with only a minor fraction of cases subsequently developing KA. A major fraction ${\tau }_1{\psi }_1$ does not develop disease, becomes PCR-negative in the early recovered stage $R_H^E\left(t\right)$, and develops LST positivity $R_H^L\left(t\right)$ at a rate of $\epsilon$. A remaining fraction of ${\tau }_2{\psi }_1$ asymptomatic infections become symptomatic $I_H^S\left(t\right)$, and a tiny fraction ${\tau }_3{\psi }_1$ of putatively recovering humans develop a state of PCR-negativity in peripheral blood, while still harboring a non-detectable number of parasites; this is denoted as $R_H^2\left(t\right)$, from where relapse to PKDL ($I_H^P\left(t\right)$) follows. A fraction ${\tau }_2{\psi }_1$ develops symptomatic KA ($I_H^S$) whilst staying PCR-positive. Those developing symptoms are eligible for treatment. If not dying, these patients receive first-line treatment ($I_H^1$). A proportion ${\omega }_1$ of patients clear parasites under first-line treatment, recover ($R_H^1$), and finally become LST-positive ($R_H^L$) like those with asymptomatic infections. A fraction $\phi$ of the LST-positive ($R_H^L$) patients become susceptible. The remaining proportion of patients represent the treatment failures that are split into a PCR-positive proportion of ${\omega }_3$ patients receiving second-line treatment ($I_H^2$) and a proportion of ${\omega }_2$ patients putatively recovering into a state of PCR-negativity. The second proportion still harbors a non-detectable number of parasites ($R_H^2$, from where relapse to PKDL will follow). For second-line treatment (in state $I_H^2$), as with first-line treatment, a proportion ${\gamma }_2$ of patients under second-line treatment recover ($R_H^1$) and become LST-positive ($R_H^L$). The remaining proportion of ${\gamma }_1$ patients recover putatively into a state of PCR-negativity ($R_H^2$) from which, again, relapse to PCR-positivity and PKDL ($I_H^P$) follows for those who do not die. All PKDL patients ($I_H^P$) are treated until full recovery ($R_H^1$).

The sandfly population is considered in the susceptible ($S_F$), exposed ($E_F$), or infectious ($I_F$),) stage. Sandflies can become infected by blood meals taken from an infectious animal through a certain interaction. The infection rate of sandflies is determined by the following; the biting rate, the infection probabilities of sandflies dependent on the infection status of the hosts, and the number of infectious hosts. We assume that each blood meal of a susceptible sandfly leads to a sandfly infection if taken from either an asymptomatic or symptomatic infected animal.

The following assumptions are considered in the formulation of the model:

i.

Entire human infected class includes PKDL-infected people as well as early asymptomatic infected, late asymptomatic infected, and symptomatic infected people.

ii.

According to each person’s immunogenic potential, both early and late asymptomatic humans may eventually show symptoms and join the symptomatic infected group, and the symptomatic infected group may eventually become PKDL-infected, or they may gradually recover [19].

iii.

The case where humans from the early asymptomatic class recover is not considered in this model.

iv.

The capacity of the population of human, animals, and sandflies to fend off infection also plays a role in how smoothly the population moves from one compartment to another.

v.

Humans acquire the disease but do not transmit.

vi.

Transmission between animals and sandflies is assumed to be indirect.

vii.

A fraction of late asymptomatic infected humans that are developing symptoms will receive first treatment.

Through the schematic diagram depicted in Figure 1, a system of non-linear differential equations is obtained and presented below:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{d{S}_H}{dt}}& =& {\mathsf{\Pi }}_H-{\lambda }_F{S}_H-{\mu }_H{S}_H,\hfill \\ \hfill {\displaystyle \frac{d{E}_H}{dt}}& =& {\lambda }_F{S}_H-({\theta }_1+{\mu }_H)E_H,\hfill \\ \hfill {\displaystyle \frac{d{I}_H^E}{dt}}& =& {\theta }_1{E}_H-({\theta }_2+{\mu }_H)I_H^E,\hfill \\ \hfill {\displaystyle \frac{d{I}_H^L}{dt}}& =& {\theta }_2{I}_H^E-({\psi }_1+{\mu }_H)I_H^L,\hfill \\ \hfill {\displaystyle \frac{d{R}_H^E}{dt}}& =& {\tau }_1{\psi }_1{I}_H^L-(\epsilon +{\mu }_H)R_H^E,\hfill \\ \hfill {\displaystyle \frac{d{R}_H^L}{dt}}& =& \epsilon R_H^E+{\varrho }_1{\psi }_4{R}_H^1-{\mu }_H{R}_H^L,\hfill \\ \hfill {\displaystyle \frac{d{I}_H^S}{dt}}& =& {\tau }_2{\psi }_1{I}_H^L-(\sigma +{\mu }_H)I_H^S,\hfill \\ \hfill {\displaystyle \frac{d{I}_H^1}{dt}}& =& \sigma I_H^S-({\psi }_2+{\mu }_H)I_H^1,\hfill \\ \hfill {\displaystyle \frac{d{I}_H^2}{dt}}& =& {\omega }_3{\psi }_2{I}_H^1-({\psi }_3+{\mu }_H)I_H^2,\hfill \\ \hfill {\displaystyle \frac{d{I}_H^P}{dt}}& =& \delta R_H^2-(\varphi +{\mu }_H)I_H^P,\hfill \\ \hfill {\displaystyle \frac{d{R}_H^1}{dt}}& =& {\omega }_1{\psi }_2{I}_H^1+\varphi I_H^P+{\gamma }_2{\psi }_3{I}_H^2-({\varrho }_1{\psi }_4+{\mu }_H)R_H^1,\hfill \\ \hfill {\displaystyle \frac{d{R}_H^2}{dt}}& =& {\tau }_3{\psi }_1{I}_H^L+{\omega }_2{\psi }_2{I}_H^1+{\gamma }_1{\psi }_3{I}_H^2-(\delta +{\mu }_H)R_H^2,\hfill \\ \hfill {\displaystyle \frac{d{S}_A}{dt}}& =& {\mathsf{\Pi }}_A-{\eta }_F{S}_A-{\mu }_A{S}_A,\hfill \\ \hfill {\displaystyle \frac{d{E}_A}{dt}}& =& {\eta }_F{S}_A-({\zeta }_1+{\mu }_A)E_A,\hfill \\ \hfill {\displaystyle \frac{d{I}_A^A}{dt}}& =& {\varpi }_1{\zeta }_1{E}_A-({\zeta }_2+{\mu }_A)I_A^A,\hfill \\ \hfill {\displaystyle \frac{d{I}_A^S}{dt}}& =& {\varpi }_2{\zeta }_1{E}_A+{\nu }_2{\zeta }_2{I}_A^A-(\kappa +{\mu }_A)I_A^S,\hfill \\ \hfill {\displaystyle \frac{d{R}_A}{dt}}& =& {\nu }_1{\zeta }_2{I}_A^A+\kappa I_A^S-{\mu }_A{R}_A,\hfill \\ \hfill {\displaystyle \frac{d{S}_F}{dt}}& =& {\mathsf{\Pi }}_F-({\rho }_A+{\xi }_A+{\chi }_A)S_F-{\mu }_F{S}_F,\hfill \\ \hfill {\displaystyle \frac{d{E}_F}{dt}}& =& ({\rho }_A+{\xi }_A+{\chi }_A)S_F-(\pi +{\mu }_F)E_F,\hfill \\ \hfill {\displaystyle \frac{d{I}_F}{dt}}& =& \pi E_F-{\mu }_F{I}_F,\hfill \end{array}$$

(1)

subject to the following initial conditions:

$$\begin{array}{c}\hfill \begin{array}{c}{S}_H\left(0\right)=S_{H_0}>0,E_H\left(0\right)=E_{H_0}\ge 0,I_H^E\left(0\right)=I_{H_0}^E\ge 0,I_H^L\left(0\right)=I_{H_0}^L\ge 0,\hfill \\ R_H^E\left(0\right)=R_{H_0}^E\ge 0,R_H^L\left(0\right)=R_{H_0}^L\ge 0,I_H^S\left(0\right)=I_{H_0}^S\ge 0,I_H^1\left(0\right)=I_{H_0}^1\ge 0,\hfill \\ I_H^2\left(0\right)=I_{H_0}^2\ge 0,I_H^P\left(0\right)=I_{H_0}^P\ge 0,R_H^1\left(0\right)=R_{H_0}^1\ge 0,R_H^2\left(0\right)=R_{H_0}^2\ge 0,\hfill \\ S_A\left(0\right)=S_{A_0}>0,E_A\left(0\right)=E_{A_0}\ge 0,I_A^A\left(0\right)=I_{A_0}^A\ge 0,I_A^S\left(0\right)=I_{A_0}^S\ge 0,\hfill \\ R_A\left(0\right)=R_{A_0}\ge 0,S_F\left(0\right)=S_{F_0}>0,E_F\left(0\right)=E_{F_0}\ge 0,I_F\left(0\right)=I_{F_0}\ge 0.\hfill \end{array}\end{array}$$

where:

$$\begin{array}{c}\hfill {\lambda }_F={\displaystyle \frac{{\beta }_H{g}_H{I}_F}{N_H+N_A}},{\eta }_F={\displaystyle \frac{{\beta }_A{g}_A{I}_F}{N_H+N_A}},{\rho }_A={\displaystyle \frac{{\beta }_F{g}_A{I}_A^A}{N_H+N_A}},{\xi }_A={\displaystyle \frac{{\beta }_F{g}_A{I}_A^S}{N_H+N_A}},{\chi }_A={\displaystyle \frac{{\beta }_F{g}_A{h}_A{E}_A}{N_H+N_A}}\end{array}$$

$$\begin{array}{c}\hfill {\tau }_1+{\tau }_2+{\tau }_3=1,{\omega }_1+{\omega }_2+{\omega }_3=1,{\gamma }_1+{\gamma }_2=1,{\varpi }_1+{\phi }_2=1,{\nu }_1+{\nu }_2=1.\end{array}$$

3. Basic Properties of the Integer Model

Here, we explore the basic properties of the model. All state variables are dependent on time, t. The total populations of humans ($N_H$), sandflies ($N_F$), and animals ($N_A$) are given by

$$\begin{array}{cc}\hfill {\displaystyle \frac{N_H}{dt}}& ={\mathsf{\Pi }}_H-{\mu }_H{N}_H,\hfill \\ \hfill {\displaystyle \frac{N_A}{dt}}& ={\mathsf{\Pi }}_A-{\mu }_A{N}_A,\hfill \\ \hfill {\displaystyle \frac{N_F}{dt}}& ={\mathsf{\Pi }}_F-{\mu }_F{N}_F,\hfill \end{array}$$

(2)

with the initial conditions

$$\begin{array}{c}\hfill \begin{array}{c}{S}_H\left(0\right)=S_{H_0}>0,E_H\left(0\right)=E_{H_0}\ge 0,I_H^E\left(0\right)=I_{H_0}^E\ge 0,I_H^L\left(0\right)=I_{H_0}^L\ge 0,\hfill \\ R_H^E\left(0\right)=R_{H_0}^E\ge 0,R_H^L\left(0\right)=R_{H_0}^L\ge 0,I_H^S\left(0\right)=I_{H_0}^S\ge 0,I_H^1\left(0\right)=I_{H_0}^1\ge 0,\hfill \\ I_H^2\left(0\right)=I_{H_0}^2\ge 0,I_H^P\left(0\right)=I_{H_0}^P\ge 0,R_H^1\left(0\right)=R_{H_0}^1\ge 0,R_H^2\left(0\right)=R_{H_0}^2\ge 0,\hfill \\ S_A\left(0\right)=S_{A_0}>0,E_A\left(0\right)=E_{A_0}\ge 0,I_A^A\left(0\right)=I_{A_0}^A\ge 0,I_A^S\left(0\right)=I_{A_0}^S>0,\hfill \\ R_A\left(0\right)=R_{A_0}\ge 0,S_F\left(0\right)=S_{F_0}\ge 0,E_F\left(0\right)=E_{F_0}\ge 0,I_F\left(0\right)=I_{F_0}\ge 0.\hfill \end{array}\end{array}$$

Therefore, in order for the solutions to the model (1) with non-negative initial data to have biological significance, it is necessary to demonstrate that they will always remain non-negative $t>0$.

Positivity of Solution

Since the model (1) describes the human population, here, we show that all the state variables that are non-negative for all $t\ge 0$ solutions are all bounded. Before analyzing the model, it is important to show that the state variables of the model remain non-negative for all non-negative initial conditions. We claim the following result.

Theorem 1.

Let $S_H>0,E_H\ge 0,...,I_F\ge 0$. The solutions $S_H,E_H,...,I_F$ of the model system (1) for $t\ge 0$ are positive. For the model system (1), the region Ω is positively invariant and all solutions starting in Ω approach, enter, or stay in Ω , where $\mathsf{\Omega }={\mathsf{\Omega }}_H\cup {\mathsf{\Omega }}_A\cup {\mathsf{\Omega }}_F\in {\mathbb{R}}_{+}^{12}\times {\mathbb{R}}_{+}^5\times {\mathbb{R}}_{F+}^3$:

$$\begin{array}{c}\hfill \begin{array}{c}{\mathsf{\Omega }}_H=(S_H,E_H,I_H^E,I_H^L,I_H^S,R_H^E,R_H^L,I_H^1,I_H^2,I_H^P,R_H^1,R_H^2)\in {\mathbb{R}}_{H+}^{12},\hfill \\ {\mathsf{\Omega }}_A=(S_A,E_A,I_A^A,I_A^S,E_A,R_A)\in {\mathbb{R}}_{A+}^5,\hfill \\ {\mathsf{\Omega }}_F=(S_F,E_F,I_F)\in {\mathbb{R}}_{F+}^3.\hfill \end{array}\end{array}$$

Proof. 

We use the method of contradiction as in [20] to prove Theorem (1). Under the given initial conditions, it is straightforward to prove that the components of the solutions of the model system (1) are always positive; otherwise, we assume a contradiction. We claim there exists a first time, $t_1$

$t_1$: $S_H\left(t_1\right)=0,{\displaystyle \frac{d{S}_H}{d{t}_1}}<0,E_H\left(t\right)>0,I_H^E\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_2$ or claim ∃

$t_2$: $E_H\left(t_2\right)=0,{\displaystyle \frac{d{E}_H}{d{t}_2}}<0,S_H\left(t\right)>0,I_H^E\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_3$ or claim ∃

$t_3$: $I_H^E\left(t_3\right)=0,{\displaystyle \frac{d{I}_H^E}{d{t}_3}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_1$ or claim ∃

$t_4$: $I_H^L\left(t_4\right)=0,{\displaystyle \frac{d{I}_H^L}{d{t}_4}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_4$ or claim ∃

$t_5$: $R_H^E\left(t_5\right)=0,{\displaystyle \frac{d{R}_H^E}{d{t}_5}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_5$ or claim ∃

$t_6$: $R_H^L\left(t_6\right)=0,{\displaystyle \frac{d{R}_H^L}{d{t}_6}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_6$ or claim ∃

$t_7$: $I_H^S\left(t_7\right)=0,{\displaystyle \frac{d{I}_H^S}{d{t}_7}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_7$ or claim ∃

$t_8$: $I_H^1\left(t_8\right)=0,{\displaystyle \frac{d{I}_H^1}{d{t}_1}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_8$ or claim ∃

$t_9$: $I_H^2\left(t_9\right)=0,{\displaystyle \frac{d{I}_H^2}{d{t}_8}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_9$ or or claim ∃

$t_{10}$: $I_H^P\left(t_{10}\right)=0,{\displaystyle \frac{d{I}_H^P}{d{t}_{10}}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_{10}$ or claim ∃

$t_{11}$: $R_H^1\left(t_{11}\right)=0,{\displaystyle \frac{d{R}_H^1}{d{t}_{11}}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_{11}$ or claim ∃

$t_{12}$: $R_H^2\left(t_{12}\right)=0,{\displaystyle \frac{d{R}_H^2}{d{t}_{12}}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_{12}$ or claim ∃

$t_{13}$: $S_A\left(t_{13}\right)=0,{\displaystyle \frac{d{S}_A}{d{t}_{13}}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_{13}$ or claim ∃

$t_{14}$: $E_A\left(t_{14}\right)=0,{\displaystyle \frac{d{E}_A}{d{t}_{14}}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_{14}$ or claim ∃

$t_{15}$: $I_A^A\left(t_{15}\right)=0,{\displaystyle \frac{d{I}_A^A}{d{t}_{15}}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_{15}$ or claim ∃

$t_{16}$: $I_A^S\left(t_{16}\right)=0,{\displaystyle \frac{d{I}_A^S}{d{t}_{16}}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_{16}$ or claim ∃

$t_{17}$: $R_A\left(t_{17}\right)=0,{\displaystyle \frac{d{R}_A}{d{t}_{17}}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_{17}$ or claim ∃

$t_{18}$: $S_F\left(t_{18}\right)=0,{\displaystyle \frac{d{S}_F}{d{t}_{18}}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_{18}$ or claim ∃

$t_{19}$: $E_F\left(t_{19}\right)=0,{\displaystyle \frac{d{E}_A}{d{t}_{19}}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,I_F\left(t\right)>0$ for $0<t<t_{19}$ or claim ∃

$t_{20}$: $I_F\left(t_{20}\right)=0,{\displaystyle \frac{d{I}_A}{d{t}_{20}}}<0,S_H\left(t\right)>0,E_H\left(t\right)>0,...,E_F\left(t\right)>0$ for $0<t<t_{20}$.

In the first claim, we have

$${\displaystyle \frac{S_H}{d{t}_1}}={\mathsf{\Pi }}_H>0,$$

(3)

which is a contradiction to the assumption. That is, $S_H$ remains positive.

In the second claim, we have

$${\displaystyle \frac{d{E}_H}{d{t}_2}}={\lambda }_F{S}_H>0,$$

(4)

which is a contradiction to the assumption. That is, $E_H$ remains positive.

In the third claim, we have

$${\displaystyle \frac{d{I}_H^E}{d{t}_3}}={\theta }_1{E}_H>0,$$

(5)

which is a contradiction to the assumption. That is, $I_H^E$ remains positive.

In the fourth claim, we have

$${\displaystyle \frac{d{I}_H^L}{d{t}_4}}={\theta }_2{I}_H^E>0,$$

(6)

which is a contradiction to the assumption. That is, $I_H^L$ remains positive.

In the fifth claim, we have

$${\displaystyle \frac{d{R}_H^E}{d{t}_5}}={\tau }_1{\psi }_1{I}_H^L>0,$$

(7)

which is a contradiction to the assumption. That is, $R_H^E$ remains positive.

In the sixth claim, we have

$${\displaystyle \frac{d{R}_H^L}{d{t}_6}}=\epsilon R_H^E+{\varrho }_1{\psi }_4{R}_H^1>0,$$

(8)

which is a contradiction to the assumption. That is, $R_H^L$ remains positive.

In the seventh claim, we have

$${\displaystyle \frac{d{I}_H^S}{d{t}_7}}={\tau }_2{\psi }_1{I}_H^L>0,$$

(9)

which is a contradiction to the assumption. That is, $I_H^S$ remains positive.

In the eighth claim, we have

$${\displaystyle \frac{d{I}_H^1}{d{t}_8}}=\sigma I_H^S>0,$$

(10)

which is a contradiction to the assumption. That is, $I_H^1$ remains positive.

In the ninth claim, we have

$${\displaystyle \frac{d{I}_H^2}{d{t}_9}}={\omega }_3{\psi }_2{I}_H^1>0,$$

(11)

which is a contradiction to the assumption. That is, $I_H^2$ remains positive.

In the tenth claim, we have

$${\displaystyle \frac{d{I}_H^P}{d{t}_{10}}}=\delta R_H^2>0,$$

(12)

which is a contradiction to the assumption. That is, $I_H^P$ remains positive.

In the eleventh claim, we have

$${\displaystyle \frac{d{R}_H^1}{d{t}_{11}}}={\omega }_1{\psi }_2{I}_H^1+\varphi I_H^P+{\gamma }_2{\psi }_3{I}_H^2>0,$$

(13)

which is a contradiction to the assumption. That is, $R_H^1$ remains positive.

In the twelfth claim, we have

$${\displaystyle \frac{d{R}_H^2}{d{t}_{12}}}={\tau }_3{\psi }_1{I}_H^L+{\omega }_2{\psi }_2{I}_H^1+{\gamma }_1{\psi }_3{I}_H^2>0,$$

(14)

which is a contradiction to the assumption. That is, $R_H^2$ remains positive.

In the thirteenth claim, we have

$${\displaystyle \frac{d{S}_A}{d{t}_{13}}}={\mathsf{\Pi }}_A>0,$$

(15)

which is a contradiction to the assumption. That is, $S_A$ remains positive.

In the fourteenth claim, we have

$${\displaystyle \frac{d{E}_A}{d{t}_{14}}}={\eta }_F{S}_A>0,$$

(16)

which is a contradiction to the assumption. That is, $E_A$ remains positive.

In the fifteenth claim, we have

$${\displaystyle \frac{d{I}_A^A}{d{t}_{15}}}={\varpi }_1{\zeta }_1{E}_A>0,$$

(17)

which is a contradiction to the assumption. That is, $I_A^A$ remains positive.

In the sixteenth claim, we have

$${\displaystyle \frac{d{I}_A^S}{d{t}_{16}}}={\varpi }_2{\zeta }_1{E}_A+{\nu }_2{\zeta }_2{I}_A^A>0,$$

(18)

which is a contradiction to the assumption. That is, $I_A^S$ remains positive.

In the seventeenth claim, we have

$${\displaystyle \frac{d{R}_A}{d{t}_{17}}}={\nu }_1{\zeta }_2{I}_A^A+\kappa I_A^S>0,$$

(19)

which is a contradiction to the assumption. That is, $R_A$ remains positive.

In the eighteenth claim, we have

$${\displaystyle \frac{d{S}_F}{d{t}_{18}}}={\mathsf{\Pi }}_F,$$

(20)

which is a contradiction to the assumption. That is, $S_F$ remains positive.

In the nineteenth claim, we have

$${\displaystyle \frac{d{E}_F}{d{t}_{19}}}=({\rho }_A+{\xi }_A+{\chi }_A)S_F>0,$$

(21)

which is a contradiction to the assumption. That is, $E_F$ remains positive.

In the twentieth claim, we have

$${\displaystyle \frac{d{I}_F}{d{t}_{20}}}=\pi E_F>0,$$

(22)

which is a contradiction to the assumption. That is, $I_F$ remains positive.

Therefore, all solutions of the model (1) remain positive for all non-negative initial conditions as required. □

Also, since ${\displaystyle \frac{d{N}_H}{dt}}={\mathsf{\Pi }}_H-{\mu }_H{N}_H,{\displaystyle \frac{d{N}_A}{dt}}={\mathsf{\Pi }}_A-{\mu }_A{N}_A,{\displaystyle \frac{d{N}_F}{dt}}={\mathsf{\Pi }}_F-{\mu }_F{N}_F$, this means that $N_H⟶{\displaystyle \frac{{\mathsf{\Pi }}_H}{{\mu }_H}},N_A⟶{\displaystyle \frac{{\mathsf{\Pi }}_A}{{\mu }_A}}$ and $N_F⟶{\displaystyle \frac{{\mathsf{\Pi }}_F}{{\mu }_F}}$ are bounded. Based on biological considerations, the model system (1) is studied in the following feasible region:

$$\begin{array}{c}\hfill \mathsf{\Omega }=\left\{\right(S_H,E_H,I_H^E,I_H^L,R_H^E,R_H^L,I_H^S,I_H^1,I_H^2,I_H^P,R_H^1,R_H^2,S_A,E_A,\\ \hfill I_A^A,I_A^S,R_A,S_F,E_F,I_F)\\ \hfill \in {\mathbb{R}}_{+}^{20}={\mathbb{R}}_{H+}^{12}\times {\mathbb{R}}_{A+}^5\times {\mathbb{R}}_{F+}^3:N_H\left(t\right)\le {\displaystyle \frac{{\mathsf{\Pi }}_H}{{\mu }_H}},N_A\left(t\right)\le {\displaystyle \frac{{\mathsf{\Pi }}_A}{{\mu }_A}},N_F\left(t\right)\le {\displaystyle \frac{{\mathsf{\Pi }}_F}{{\mu }_F}}\}.\end{array}$$

(23)

This region is positively invariant with respect to the model system (1). This also mean that, all solutions of the model (1), with initial conditions in $\mathsf{\Omega }$, will remain in $\mathsf{\Omega }$ for all $t\ge 0$. We can conveniently consider the solutions of the model in $\mathsf{\Omega }$. (see also, [21]). Combining this result and Theorem (1), we have the following lemma:

Lemma 1.

The region Ω is positively-invariant for the model 1 with initial conditions in ${\mathbb{R}}_{+}^{20}={\mathbb{R}}_{H+}^{12}\times {\mathbb{R}}_{A+}^5\times {\mathbb{R}}_{F+}^3$.

Additionally, it is easy to see that each of the differential equations of the model system (1) is Lipschitz continuous with the given initial conditions and has solutions. Moreover, the solution is unique and since $\mathsf{\Omega }$ is a positively invariant region, the solution exist for any time $t\ge 0$ (see, [22]).

4. Fixed Points of the Model and Their Stability Analysis

In this section, we explore the existence and stability of fixed points of the model (1).

4.1. VL-Free Fixed Point

Let $x=x^{*}$ be the VL-free fixed point for the system (1). Then, $f\left(x^{*}\right)=0$, ${\displaystyle \frac{dx}{dt}}=f\left(x\right)$, $x={(x_1,x_2,\dots ,x_{20})}^T$, $x_1=S_H,x_2=E_H,\dots ,x_{20}=I_F$ in order of the model variables. This implies that ${\displaystyle \frac{d{x}^{*}}{dt}}=f\left(x^{*}\right)$, $x^{*}={(x_1^{*},x_2^{*},\dots ,x_{20}^{*})}^T$, $x_1^{*}=S_H^{*},x_2^{*}=E_H^{*},\dots ,x_{20}^{*}=I_F^{*}$.

Thus, at fixed point $f\left(x^{*}\right)=0$, implies

$$\begin{array}{c}\hfill {\mathsf{\Pi }}_H-{\lambda }_F{S}_H^{*}-{\mu }_H{S}_H^{*}=0,\end{array}$$

(24)

$$\begin{array}{c}\hfill {\lambda }_F{S}_H^{*}-({\theta }_1+{\mu }_H)E_H^{*}=0,\end{array}$$

(25)

$$\begin{array}{c}\hfill {\theta }_1{E}_H^{*}-({\theta }_2+{\mu }_H)I_H^{E*}=0,\end{array}$$

(26)

$$\begin{array}{c}\hfill {\theta }_2{I}_H^{E*}-({\psi }_1+{\mu }_H)I_H^{L*}=0,\end{array}$$

(27)

$$\begin{array}{c}\hfill {\tau }_1{\psi }_1{I}_H^{L*}-(\epsilon +{\mu }_H)R_H^{E*}=0,\end{array}$$

(28)

$$\begin{array}{c}\hfill \epsilon R_H^{E*}+{\varrho }_1{\psi }_4{R}_H^{1*}-{\mu }_H{R}_H^{L*}=0,\end{array}$$

(29)

$$\begin{array}{c}\hfill {\tau }_2{\psi }_1{I}_H^{L*}-(\sigma +{\mu }_H)I_H^{S*}=0,\end{array}$$

(30)

$$\begin{array}{c}\hfill \sigma I_H^{S*}-({\psi }_2+{\mu }_H)I_H^{1*}=0,\end{array}$$

(31)

$$\begin{array}{c}\hfill {\omega }_3{\psi }_2{I}_H^{1*}-({\psi }_3+{\mu }_H)I_H^{2*}=0,\end{array}$$

(32)

$$\begin{array}{c}\hfill \delta R_H^{2*}-(\varphi +{\mu }_H)I_H^{P*}=0,\end{array}$$

(33)

$$\begin{array}{c}\hfill {\omega }_1{\psi }_2{I}_H^{1*}+\varphi I_H^{P*}+{\gamma }_2{\psi }_3{I}_H^{2*}-({\varrho }_1{\psi }_4+{\mu }_H)R_H^{1*}=0,\end{array}$$

(34)

$$\begin{array}{c}\hfill {\tau }_3{\psi }_1{I}_H^{L*}+{\omega }_2{\psi }_2{I}_H^{1*}+{\gamma }_1{\psi }_3{I}_H^{2*}-(\delta +{\mu }_H)R_H^{2*}=0,\end{array}$$

(35)

$$\begin{array}{c}\hfill {\mathsf{\Pi }}_A-{\eta }_F{S}_A^{*}-{\mu }_A{S}_A^{*}=0,\end{array}$$

(36)

$$\begin{array}{c}\hfill {\eta }_F{S}_A^{*}-({\zeta }_1+{\mu }_A)E_A^{*}=0,\end{array}$$

(37)

$$\begin{array}{c}\hfill {\varpi }_1{\zeta }_1{E}_A^{*}-({\zeta }_2+{\mu }_A)I_A^{A*}=0,\end{array}$$

(38)

$$\begin{array}{c}\hfill {\varpi }_2{\zeta }_1{E}_A^{*}+{\nu }_2{\zeta }_2{I}_A^{A*}-(\kappa +{\mu }_A)I_A^{S*}=0,\end{array}$$

(39)

$$\begin{array}{c}\hfill {\nu }_1{\zeta }_2{I}_A^{A*}+\kappa I_A^{S*}-{\mu }_A{R}_A^{*}=0,\end{array}$$

(40)

$$\begin{array}{c}\hfill {\mathsf{\Pi }}_F-({\rho }_A+{\xi }_A+{\chi }_A)S_F^{*}-{\mu }_F{S}_F^{*}=0,\end{array}$$

(41)

$$\begin{array}{c}\hfill ({\rho }_A+{\xi }_A+{\chi }_A)S_F^{*}-(\pi +{\mu }_F)E_F^{*}=0,\end{array}$$

(42)

$$\begin{array}{c}\hfill \pi E_F^{*}-{\mu }_F{I}_F^{*}=0.\end{array}$$

(43)

At the VL-free fixed point, it means there is no presence of the disease, meaning $E_H=0$, $I_H^E=0$, $I_H^L=0,R_H^E=0$, $R_H^L=0,I_H^S=0$, $I_H^1=0,I_H^2=0$, $I_H^P=0,R_H^1=0,R_H^2=0$, $E_A=0,I_A^A=0,I_A^S=0,R_A=0,E_F=0,I_F=0$.

Thus, we have

$${\mathsf{\Pi }}_H-{\mu }_H{S}_H^{*}=0⟺S_H^{*}={\displaystyle \frac{{\mathsf{\Pi }}_H}{{\mu }_H}},$$

$${\mathsf{\Pi }}_A-{\mu }_A{S}_A^{*}=0⟺S_A^{*}={\displaystyle \frac{{\mathsf{\Pi }}_A}{{\mu }_A}},$$

$${\mathsf{\Pi }}_F-{\mu }_F{S}_F^{*}=0⟺S_F^{*}={\displaystyle \frac{{\mathsf{\Pi }}_F}{{\mu }_F}}.$$

Therefore, the VL-free fixed point is given by

$$\begin{array}{ccc}\hfill x^{*}& =& {(x_1^{*},x_2^{*},\dots ,x_{20}^{*})}^T,\hfill \\ & =& {(S_H^{*},E_H^{*},\dots ,I_F^{*})}^T.\hfill \end{array}$$

Hence, the VL-free fixed point, denoted by $E^0$, is given by

$$E^0=\left({\displaystyle \frac{{\mathsf{\Pi }}_H}{{\mu }_H}},0,0,0,0,0,0,0,0,0,0,0,{\displaystyle \frac{{\mathsf{\Pi }}_A}{{\mu }_A}},0,0,0,0,{\displaystyle \frac{{\mathsf{\Pi }}_F}{{\mu }_F}},0,0\right).$$

(44)

4.2. Basic Reproduction Number $R_0$

The number of secondary infections caused by a single infectious individual over the course of their whole infectious period is known as the basic reproduction number. The basic reproduction number can be expressed mathematically as a spectral radius. The number of new infections caused by a single infected person in a community that is totally susceptible is defined by the spectral radius $R_0$, a threshold parameter for disease control [23]. The number of VL infections caused by an active VL case is what we refer to as the basic reproduction number, or $R_0$, in this instance. To ascertain the model system’s basic reproductive number, we employ the method of next generation matrix described in [23] to determine the basic reproductive number of the system (1). The matrices $\mathcal{F}$ and $\mathcal{V}$, for the new infection terms and the remaining transfer terms, are, respectively, given by

$$\begin{array}{c}\hfill \mathcal{F}=\left(\begin{array}{c}{\lambda }_F{S}_H\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\\ {\eta }_F{S}_A\\ 0\\ 0\\ ({\rho }_A+{\xi }_A+{\chi }_A)S_F\\ 0\end{array}\right),\mathcal{V}=\left(\begin{array}{c}({\theta }_1+{\mu }_H)E_H\\ -{\theta }_1{E}_H+({\theta }_2+{\mu }_H)I_H^E\\ -{\theta }_2{I}_H^E+({\psi }_1+{\mu }_H)I_H^L\\ -{\tau }_2{\psi }_1{I}_H^L+(\sigma +{\mu }_H)I_H^S\\ -\sigma I_H^S+({\psi }_2+{\mu }_H)I_H^1\\ -{\omega }_3{\psi }_2{I}_H^1+({\psi }_3+{\mu }_H)I_H^2\\ -\delta R_H^2+(\varphi +{\mu }_H)I_H^P\\ ({\zeta }_1+{\mu }_A)E_A\\ -{\varpi }_1{\zeta }_1{E}_A+({\zeta }_2+{\mu }_A)I_A^A,\\ -{\varpi }_2{\zeta }_1{E}_A-{\nu }_2{\zeta }_2{I}_A^A+(\kappa +{\mu }_A)I_A^S\\ (\pi +{\mu }_F)E_F\\ -\pi E_F+{\mu }_F{I}_F\end{array}\right).\end{array}$$

We then obtain the matrices $\mathit{F}$ and $\mathit{V}$ at the VL-free fixed point, $E^0$. Thus, we have

$$\mathit{F}=\left(\begin{array}{cccccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& b_1\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& b_2\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& b_3& b_4& b_4& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right),$$

(45)

$b_1={\displaystyle \frac{{\beta }_H{g}_H{\mathsf{\Pi }}_H{\mu }_A}{{\mathsf{\Pi }}_H{\mu }_A+{\mathsf{\Pi }}_A{\mu }_H}},b_2={\displaystyle \frac{{\beta }_A{g}_A{\mathsf{\Pi }}_A{\mu }_H}{{\mathsf{\Pi }}_H{\mu }_A+{\mathsf{\Pi }}_A{\mu }_H}},b_3={\displaystyle \frac{{\beta }_F{g}_A{\mathsf{\Pi }}_F{\mu }_A{\mu }_H}{{\mu }_F({\mathsf{\Pi }}_H{\mu }_A+{\mathsf{\Pi }}_A{\mu }_H)}},b_4={\displaystyle \frac{{\beta }_F{g}_A{h}_A{\mathsf{\Pi }}_F{\mu }_A{\mu }_H}{{\mu }_F({\mathsf{\Pi }}_H{\mu }_A+{\mathsf{\Pi }}_A{\mu }_H)}}$

$$\mathit{V}=\left(\begin{array}{cccccccccccc}{a}_{11}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ -{\theta }_1& a_{12}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& -{\theta }_2& a_{13}& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& -a_{22}& a_{14}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& -\sigma & a_{15}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& -a_{23}& a_{16}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& a_{17}& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& a_{18}& 0& 0& 0& \\ 0& 0& 0& 0& 0& 0& 0& -a_{24}& a_{19}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& -a_{25}& -a_{26}& a_{20}& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& a_{21}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& -\pi & {\mu }_F\end{array}\right),$$

(46)

where

$$\begin{array}{c}\hfill a_{11}={\theta }_1+{\mu }_H,a_{12}={\theta }_2+{\mu }_H,a_{13}={\psi }_1+{\mu }_H,a_{14}=\sigma +{\mu }_H,a_{15}={\psi }_2+{\mu }_H,\\ \hfill a_{16}={\psi }_3+{\mu }_H,a_{17}=\varphi +{\mu }_H,a_{18}={\zeta }_1+{\mu }_A,a_{19}={\zeta }_2+{\mu }_A,a_{20}=\kappa +{\mu }_A,\\ \hfill a_{21}=\pi +{\mu }_A,a_{22}={\tau }_2{\psi }_1,a_{23}={\omega }_2{\psi }_2,a_{24}={\varpi }_1{\zeta }_1,a_{25}={\varpi }_2{\zeta }_1,a_{26}={\nu }_2{\zeta }_2.\end{array}$$

the determinant of $\mathit{V}$ is obtained to be

$$det\left(\mathit{V}\right)=a_{11}{a}_{12}{a}_{13}{a}_{14}{a}_{15}{a}_{16}{a}_{17}{a}_{18}{a}_{19}{a}_{20}{a}_{21}{\mu }_F$$

$$V^{-1}=\left(\begin{array}{cccccccccccc}{\displaystyle \frac{1}{a_{11}}}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ A_1{\displaystyle \frac{1}{a_{12}}}0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ A_2{A}_3{\displaystyle \frac{1}{a_{13}}}0& 0& 0& 0& 0& 0& 0& 0& 0\\ A_4{A}_5{A}_6{\displaystyle \frac{1}{a_{14}}}0& 0& 0& 0& 0& 0& 0& 0\\ A_4{A}_5{A}_6{\displaystyle \frac{1}{a_{14}}}0& 0& 0& 0& 0& 0& 0& 0\\ A_7{A}_8{A}_9{A}_{10}-a_{15}& 0& 0& 0& 0& 0& 0& 0\\ A_{11}{A}_{12}{A}_{13}{A}_{14}{A}_{15}{\displaystyle \frac{1}{a_{16}}}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& {\displaystyle \frac{1}{a_{17}}}0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& {\displaystyle \frac{1}{a_{18}}}0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& a_{16}& {\displaystyle \frac{1}{a_{19}}}0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& a_{17}& a_{18}& {\displaystyle \frac{1}{a_{20}}}0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& {\displaystyle \frac{1}{a_{21}}}0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& a_{19}& {\displaystyle \frac{1}{{\mu }_F}}\end{array}\right)$$

(47)

$$\begin{array}{c}\hfill \begin{array}{c}{A}_1={\displaystyle \frac{{\theta }_1}{a_{11}{a}_{12}}},A_2={\displaystyle \frac{{\theta }_1{\theta }_2}{a_{11}{a}_{12}{a}_{13}}},A_3={\displaystyle \frac{{\theta }_2}{a_{12}{a}_{13}}},A_4={\displaystyle \frac{{\theta }_1{\theta }_2{a}_{22}}{a_{11}{a}_{12}{a}_{13}{a}_{14}}},A_5={\displaystyle \frac{{\theta }_2{a}_{22}}{a_{12}{a}_{13}{a}_{14}}},\hfill \\ A_6={\displaystyle \frac{a_{22}}{a_{13}{a}_{14}}},A_7={\displaystyle \frac{{\theta }_1{\theta }_2\sigma a_{22}}{a_{11}{a}_{12}{a}_{13}{a}_{14}{a}_{15}}},A_8={\displaystyle \frac{{\theta }_2\sigma a_{22}}{a_{12}{a}_{13}{a}_{14}}},A_9={\displaystyle \frac{\sigma a_{22}}{a_{13}{a}_{14}}},A_{10}={\displaystyle \frac{\sigma }{a_{14}}},\hfill \\ A_{11}={\displaystyle \frac{{\theta }_1{\theta }_2\sigma a_{22}{a}_{22}}{a_{11}{a}_{12}{a}_{13}{a}_{14}{a}_{16}}},A_{12}={\displaystyle \frac{{\theta }_2\sigma a_{22}{a}_{29}}{a_{12}{a}_{13}{a}_{14}{a}_{16}}},A_{13}={\displaystyle \frac{\sigma a_{22}{a}_{23}}{a_{13}{a}_{14}{a}_{16}}},A_{14}={\displaystyle \frac{\sigma a_{23}}{a_{14}{a}_{16}}},A_{15}={\displaystyle \frac{a_{29}}{a_{16}}},\hfill \\ A_{16}={\displaystyle \frac{a_{24}}{a_{18}{a}_{25}}},A_{17}={\displaystyle \frac{a_{19}{a}_{25}+a_{24}{a}_{26}}{a_{18}{a}_{15}{a}_{20}}},A_{18}={\displaystyle \frac{a_{26}}{a_{25}{a}_{20}}},A_{19}={\displaystyle \frac{\pi }{{\mu }_F{a}_{21}}}.\hfill \end{array}\end{array}$$

Thus, we have

$${\mathit{FV}}^{-\mathit{1}}=\left(\begin{array}{cccccccccccc}0& 0& 0& 0& 0& 0& 0& 0& 0& 0& {\displaystyle \frac{\pi b_1}{{\mu }_F{a}_{21}}}& {\displaystyle \frac{b_1}{{\mu }_F}}\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& {\displaystyle \frac{\pi b_2}{{\mu }_F{a}_{21}}}& {\displaystyle \frac{b_2}{{\mu }_F}}\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& & 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& B_1& B_2& {\displaystyle \frac{b_4}{a_{20}}}& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right),$$

(48)

where

$$\begin{array}{c}\hfill B_1={\displaystyle \frac{b_4{a}_{24}(a_{20}+a_{26})+a_{25}{b}_4(a_{20}+a_{25})}{a_{18}{a}_{25}{a}_{20}}},B_2={\displaystyle \frac{a_{20}+a_{26}}{a_{25}{a}_{20}}}.\end{array}$$

The spectral radius $\rho \left(F{V}^{-1}\right)=R_0$, the basic reproduction number to the system, which is the dominant eigenvalue, after simplifying, is given as

$$R_0=\sqrt{{\displaystyle \frac{\pi (a_{20}{b}_3+b_2{b}_4{a}_{24}{a}_{25})+\pi b_4{a}_{24}(a_{26}+a_{20}{b}_2)}{{\mu }_F{a}_{18}{a}_{20}{a}_{21}{a}_{25}}}}.$$

(49)

4.3. Local Stability Analysis of VL-Free Fixed Point

To analyze the stability of the VL-free fixed point, we linearize the non-linear system (1) by taking a small perturbation about the fixed points.

Since the decomposed matrices $\mathcal{F}$ and $\mathcal{V}$ exist, which satisfy conditions (A1)–(A5) in [23], it follows from Theorem 2 of the same [23] that we have the following.

Lemma 2.

The VL-free $E^0$ of the model given by (1) is locally asymptotically stable (LAS) whenever $R_0<1$, and unstable when $R_0>1$.

The epidemiological implication of Lemma 2 is that a small influx of infectious individuals will not generate large outbreaks in the population if $R_0<1$. Basically, $R_0$ calculates the typical number of secondary cases in a population that is fully susceptible that are caused by a single infectious individual. According to Lemma 2, if $R_0<1$, then a modest excess of infectious individuals will not cause significant outbreaks in the community. For the VL-free $E^0$ of the model when $R_0<1$, a global asymptotic stability (GAS) property needs to be provided for the disease elimination to be independent of the initial sizes of the populations of the model. We examine now the global asymptotic stability (GAS) of the VL-free fixed point, $E^0$.

4.4. Global Stability Analysis of the VL-Free Fixed Point

The basic reproduction number of the model is the threshold quantity $R_0$ as obtained in (49). It measures how many new infected cells, on average, are produced [21,23]. It is implied by Lemma 1 that the disease may die out when $R_0<1$ and the initial sizes of the model’s subpopulations are within $\mathsf{\Omega }$. It is crucial to demonstrate that the VL-free fixed point is globally asymptotically stable to guarantee the independence of disease elimination when $R_0<1$.

We explore the global asymptotic stability of the VL-free fixed point $E^0$ of the model using a Lyapunov functional.

Theorem 2.

Suppose that ${\mu }_F{a}_{18}{a}_{19}{a}_{20}{a}_{21}\ge \pi {\beta }_A{\beta }_F{g}_A^2[a_{24}(a_{20}+a_{26})+a_{19}(a_{20}{h}_A+a_{25})]$ and $\pi a_{24}{b}_2{b}_4>{\mu }_F{a}_{18}{a}_{20}{a}_{21}$, the VL-free fixed point ($E^0$) is globally asymptotically stable in the region Ω whenever $R_0<1$.

Proof. 

Consider the Lyapunov function

$$\mathit{V}=g_1{E}_H+I_H^E+I_H^L+I_H^1+I_H^2+I_H^S+g_2{I}_H^P+g_3{E}_A+g_4{I}_A^A+g_5{I}_A^S+g_6{E}_F+g_7{I}_F,$$

(50)

where

$$\begin{array}{c}\hfill \begin{array}{c}{g}_1={\displaystyle \frac{{\mu }_F{a}_{18}{a}_{19}{a}_{20}{a}_{21}-\pi {\beta }_A{\beta }_F{g}_A^2[a_{24}(a_{20}+a_{26})+a_{19}(a_{20}{h}_A+a_{25})]}{a_{18}\pi {\beta }_H{g}_H}},\hfill \\ g_2=\pi a_{20}{b}_3+\pi a_{24}{a}_{26}{b}_4+\pi a_{20}{a}_{24}{b}_2{b}_4+a_{2}5(\pi a_{24}{b}_2{b}_4-{\mu }_F{a}_{18}{a}_{20}{a}_{21}),\hfill \\ g_3={\displaystyle \frac{{\beta }_F{g}_A[a_{24}(a_{20}+a_{26})+a_{19}(a_{20}{h}_A+a_{25})]}{a_{18}}},\hfill \\ g_4={\beta }_F{g}_A(a_{20}+a_{26}),\hfill \\ g_5=a_{19}{\beta }_F{g}_A,\hfill \\ g_6=a_{19}{a}_{20},\hfill \\ g_7={\displaystyle \frac{a_{19}{a}_{20}{a}_{21}}{\pi }},\hfill \end{array}\end{array}$$

with the Lyapunov derivative along the solution curve

$$\dot{\mathit{V}}=g_1\dot{E_H}+\dot{I_H^E}+\dot{I_H^L}+\dot{I_H^1}+\dot{I_H^2}+\dot{I_H^S}+g_2\dot{I_H^P}+g_3\dot{E_A}+g_4\dot{I_A^A}+g_5\dot{I_A^S}+g_6\dot{E_F}+g_7\dot{I_F},$$

$$\begin{array}{c}\hfill \begin{array}{c}{\dot{E}}_H={\displaystyle \frac{d{E}_H}{dt}},{\dot{I}}_H^E={\displaystyle \frac{d{I}_H^E}{dt}},{\dot{I}}_H^L={\displaystyle \frac{d{I}_H^L}{dt}}{\dot{I}}_H^1={\displaystyle \frac{d{I}_H^1}{dt}},{\dot{I}}_H^2={\displaystyle \frac{d{I}_H^2}{dt}},{\dot{I}}_H^S={\displaystyle \frac{d{I}_H^S}{dt}},{\dot{I}}_H^P={\displaystyle \frac{d{I}_H^P}{dt}},\hfill \\ {\dot{E}}_A={\displaystyle \frac{d{E}_A}{dt}},{\dot{I}}_A^A={\displaystyle \frac{d{I}_A^A}{dt}},{\dot{I}}_A^S={\displaystyle \frac{d{I}_A^S}{dt}},{\dot{E}}_F={\displaystyle \frac{d{E}_F}{dt}}{\dot{I}}_F={\displaystyle \frac{d{I}_F}{dt}}.\hfill \end{array}\end{array}$$

Substituting the derivatives of each from the system (1) in $\dot{\mathit{V}}$, we have

$$\begin{array}{ccc}\hfill \dot{\mathit{V}}& =& g_1[{\lambda }_F{S}_H-({\theta }_1+{\mu }_H)E_H]+[{\theta }_1{E}_H-({\theta }_2+{\mu }_H)I_H^E]+[{\theta }_2{I}_H^E-({\psi }_1+{\mu }_H)I_H^L]\hfill \\ & +& [\sigma I_H^S-({\psi }_2+{\mu }_H)I_H^1]+[{\omega }_3{\psi }_2{I}_H^1-({\psi }_3+{\mu }_H)I_H^2]+[{\tau }_2{\psi }_1{I}_H^L-(\sigma +{\mu }_H)I_H^S]\hfill \\ & +& g_2[\delta R_H^2-(\varphi +{\mu }_H)I_H^P]+g_3[{\eta }_F{S}_A-({\zeta }_1+{\mu }_A)E_A]+g_4[{\varpi }_1{\zeta }_1{E}_A-({\zeta }_2+{\mu }_A)I_A^A]\hfill \\ & +& g_5[{\varpi }_2{\zeta }_1{E}_A+{\nu }_2{\zeta }_2{I}_A^A-(\kappa +{\mu }_A)I_A^S]+g_6[({\rho }_A+{\xi }_A+{\chi }_A)S_F-(\pi +{\mu }_F)E_F]\hfill \\ & +& g_7[\pi E_F-{\mu }_F{I}_F],\hfill \\ & =& g_1{\lambda }_F{S}_H+g_3{\eta }_F{S}_A+g_6({\rho }_A+{\xi }_A+{\chi }_A)S_F+[{\theta }_1-({\theta }_1+{\mu }_H)]E_H+[{\theta }_2\hfill \\ & -& ({\theta }_2+{\mu }_H)]I_H^E+[{\tau }_2{\psi }_1-({\psi }_1+{\mu }_H)]I_H^L+[\sigma -(\sigma +{\mu }_H)]I_H^S+[{\omega }_3{\psi }_2\hfill \\ & -& ({\psi }_2+{\mu }_H)]I_H^1-({\psi }_3+{\mu }_H)I_H^2-g_2(\varphi +{\mu }_H)I_H^P+g_2\delta R_H^2+[g_4{\varpi }_1{\zeta }_1\hfill \\ \hfill & -& g_3({\zeta }_1+{\mu }_A)+g_5{\varpi }_2{\zeta }_1]E_A+[g_5{\nu }_2{\zeta }_2-g_4({\zeta }_2+{\mu }_A)]I_A^A-g_5(\kappa +{\mu }_A)I_A^S\hfill \\ & +& [g_7\pi -g_6(\pi +{\mu }_F)]E_F-g_7{\mu }_F{I}_F.\hfill \end{array}$$

(51)

Now,

$$\begin{array}{c}{g}_1{\lambda }_F{S}_H={\displaystyle \frac{g_1{\beta }_H{g}_H{I}_F}{N_H+N_A}}{S}_H<g_1{\beta }_H{g}_H{I}_F,\hfill \\ g_3{\eta }_F{S}_A={\displaystyle \frac{g_3{\beta }_A{g}_A{I}_F}{N_H+N_A}}{S}_A<g_3{\beta }_A{g}_A{I}_F,\hfill \\ g_6({\rho }_A+{\xi }_A+{\chi }_A)S_F={\displaystyle \frac{g_6({\beta }_F{g}_A{I}_A^A+{\beta }_F{g}_A{I}_A^S+{\beta }_F{g}_A{h}_A{E}_A)}{N_H+N_A}}{S}_F<g_6({\beta }_F{g}_A{I}_A^A+{\beta }_F{g}_A{I}_A^S+{\beta }_F{g}_A{h}_A{E}_A).\hfill \end{array}$$

Also, using the fact that if $a,b,c>0$, then $a+b-c<a+b$, dropping $g_2(\varphi +{\mu }_H)I_H^P$, and simplifying the term according to the variable, (51) becomes

$$\begin{array}{ccc}\hfill & \le & [{\theta }_1-({\theta }_1+{\mu }_H)]E_H+[{\theta }_2-({\theta }_2+{\mu }_H)]I_H^E+[{\tau }_2{\psi }_1-({\psi }_1+{\mu }_H)]I_H^L+[\sigma -(\sigma +{\mu }_H)]I_H^S\hfill \\ \hfill & +& [{\omega }_3{\psi }_2-(ps{i}_2+{\mu }_H)]I_H^1-({\psi }_3+{\mu }_H)I_H^2+g_2\delta R_H^2+[g_4{\varpi }_1{\zeta }_1-g_3({\zeta }_1+{\mu }_A)+g_5{\varpi }_2{\zeta }_1\hfill \\ \hfill & +& g_6{\beta }_F{g}_A{h}_A]E_A+[g_5{\nu }_2{\zeta }_2-g_4({\zeta }_2+{\mu }_A)+g_6{\beta }_F{g}_A]I_A^A+[g_6{\beta }_F{g}_A-g_5(\kappa +{\mu }_A)I_A^S\hfill \\ \hfill & +& [g_7\pi -g_6(\pi +{\mu }_F)]E_F+[g_1{\beta }_H{g}_H+g_3{\beta }_A{g}_A-g_7{\mu }_F]I_F,\hfill \\ & =& -{\mu }_H{E}_H-{\mu }_H{I}_H^E-({\tau }_1{\psi }_1+{\tau }_3{\psi }_1+{\mu }_H)I_H^L-{\mu }_H{I}_H^S-({\omega }_1{\psi }_2+{\omega }_2{\psi }_2+{\mu }_H)I_H^1\hfill \\ \hfill & -& ({\psi }_3+{\mu }_H)I_H^2+g_2\delta R_H^2+[g_4{\varpi }_1{\zeta }_1-g_3({\zeta }_1+{\mu }_A)+g_5{\varpi }_2{\zeta }_1+g_6{\beta }_F{g}_A{h}_A]E_A\hfill \\ \hfill & +& [g_5{\nu }_2{\zeta }_2-g_4({\zeta }_2+{\mu }_A)+g_6{\beta }_F{g}_A]I_A^A+[g_6{\beta }_F{g}_A-g_5(\kappa +{\mu }_A)I_A^S+[g_7\pi -g_6(\pi +{\mu }_F)]E_F\hfill \\ & +& [g_1{\beta }_H{g}_H+g_3{\beta }_A{g}_A-g_7{\mu }_F]I_F.\hfill \end{array}$$

(52)

Substituting $g_1,g_2,g_3,g_4,g_5,g_6,g_7$ in (52), we have

$$[g_4{\varpi }_1{\zeta }_1-g_3({\zeta }_1+{\mu }_A)+g_5{\varpi }_2{\zeta }_1+g_6{\beta }_F{g}_A{h}_A]E_A=0,$$

$$[g_5{\nu }_2{\zeta }_2-g_4({\zeta }_2+{\mu }_A)+g_6{\beta }_F{g}_A]I_A^A=0,$$

$$[g_6{\beta }_F{g}_A-g_5(\kappa +{\mu }_A)I_A^S=0,$$

$$[g_7\pi -g_6(\pi +{\mu }_F)]E_F=0,$$

$[g_1{\beta }_H{g}_H+g_3{\beta }_A{g}_A-g_7{\mu }_F]I_F=0$ and (52) becomes

$$\begin{array}{ccc}& =& -{\mu }_H{E}_H-{\mu }_H{I}_H^E-(a_{27}+a_{32}+{\mu }_H)I_H^L-{\mu }_H{I}_H^S-({\omega }_1{\psi }_2+a_{23}+{\mu }_H)I_H^1\hfill \\ & -& ({\psi }_3+{\mu }_H)I_H^2+[\pi a_{20}{b}_3+\pi a_{24}{a}_{26}{b}_4+\pi a_{20}{a}_{24}{b}_2{b}_4+a_{25}(\pi a_{24}{b}_2{b}_4\hfill \\ & -& {\mu }_F{a}_{18}{a}_{20}{a}_{21}\left)\right]\delta R_H^2,\hfill \\ & =& -{\mu }_H{E}_H-{\mu }_H{I}_H^E-(a_{27}+a_{32}+{\mu }_H)I_H^L-{\mu }_H{I}_H^S-({\omega }_1{\psi }_2+a_{23}+{\mu }_H)I_H^1\hfill \\ & -& ({\psi }_3+{\mu }_H)I_H^2+[\pi (a_{20}{b}_3+a_{24}{a}_{25})+\pi a_{24}{b}_4(a_{26}+a_{20}{b}_2)-{\mu }_F{a}_{18}{a}_{20}{a}_{21}{a}_{25})]\delta R_H^2,\hfill \\ & =& -{\mu }_H{E}_H-{\mu }_H{I}_H^E-(a_{27}+a_{32}+{\mu }_H)I_H^L-{\mu }_H{I}_H^S-({\omega }_1{\psi }_2+a_{23}+{\mu }_H)I_H^1\hfill \\ & -& ({\psi }_3+{\mu }_H)I_H^2+\left({\displaystyle \frac{\pi (a_{20}{b}_3+a_{24}{a}_{25})+\pi a_{24}{b}_4(a_{26}+a_{20}{b}_2)}{{\mu }_F{a}_{18}{a}_{20}{a}_{21}{a}_{25}}}-1\right)\delta {\mu }_F{a}_{18}{a}_{20}{a}_{21}{a}_{25}{R}_H^2,\hfill \\ & =& -{\mu }_H{E}_H-{\mu }_H{I}_H^E-(a_{27}+a_{32}+{\mu }_H)I_H^L-{\mu }_H{I}_H^S-({\omega }_1{\psi }_2+a_{23}+{\mu }_H)I_H^1\hfill \\ & -& ({\psi }_3+{\mu }_H)I_H^2+(R_0^2-1)\delta {\mu }_F{a}_{18}{a}_{20}{a}_{21}{a}_{25}{R}_H^2,\hfill \\ & <& 0.\hfill \end{array}$$

Clearly, $\dot{V}\le 0$ for $R_0<1$ with $V=0$ only at the VL-free fixed point, $E^0$. Thus, it follows that by LaSalle’s invariance principle (as stated in [22]), $E^0$ is globally asymptotically stable. □

Theorem 2 holds biological significance, as it states that VL can be eliminated from an infected individual whenever $R_0<1$. Specifically, for the system (1), $R_0<1$ is both necessary and sufficient for VL elimination. Additionally, in a large population, the disease dies out.

4.5. Existence of VL Endemic Fixed Point

In order to determine the presence of a non-zero (endemic) model fixed point, we let

$$\begin{array}{ccc}\hfill E^{**}& =& (S_H^{**},E_H^{**},I_H^{E**},I_H^{L**},R_H^{E**},R_H^{L**},I_H^{S**},I_H^{1**},I_H^{2**},I_H^{P**},R_H^{1**},R_H^{2**},S_A^{**},\hfill \\ & & E_A^{**},I_A^{A**},I_A^{S**},R_A^{**},S_F^{**},E_F^{**},I_F^{**}).\hfill \end{array}$$

(53)

denote any arbitrary fixed point of the model. We then solve the equations of the model in terms of the associated forces of infection (${\lambda }_F,{\eta }_F,{\rho }_A,{\zeta }_A,{\xi }_A$). It follows from (24)–(43), that the system has a unique endemic fixed point, $E^{**}$, where

$$\begin{array}{cc}\hfill S_H^{**}& ={\displaystyle \frac{{\mathsf{\Pi }}_H}{{\lambda }_F^{**}+{\mu }_H}},\hfill \\ \hfill E_H^{**}& ={\displaystyle \frac{{\mathsf{\Pi }}_H{\lambda }_F^{**}}{a_{11}({\lambda }_F^{**}+{\mu }_H)}},\hfill \\ \hfill I_H^{E**}& ={\displaystyle \frac{{\theta }_1{\mathsf{\Pi }}_H{\lambda }_F^{**}}{a_{11}{a}_{12}({\lambda }_F^{**}+{\mu }_H)}},\hfill \\ \hfill I_H^{L**}& ={\displaystyle \frac{{\theta }_1{\theta }_2{\mathsf{\Pi }}_H{\lambda }_F^{**}}{a_{11}{a}_{12}{a}_{13}({\lambda }_F^{**}+{\mu }_H)}},\hfill \\ \hfill R_H^{E**}& ={\displaystyle \frac{a_{27}{\theta }_1{\theta }_2{\mathsf{\Pi }}_H{\lambda }_F^{**}}{a_{11}{a}_{12}{a}_{13}{a}_{35}({\lambda }_F^{**}+{\mu }_H)}},\hfill \\ \hfill I_H^{1**}& ={\displaystyle \frac{a_{22}{\theta }_1{\theta }_2\sigma {\mathsf{\Pi }}_H{\lambda }_F^{**}}{a_{11}{a}_{12}{a}_{13}{a}_{14}{a}_{15}({\lambda }_F^{**}+{\mu }_H)}},\hfill \\ \hfill I_H^{2**}& ={\displaystyle \frac{a_{22}{a}_{23}{\theta }_1{\theta }_2\sigma {\mathsf{\Pi }}_H{\lambda }_F^{**}}{a_{11}{a}_{12}{a}_{13}{a}_{14}{a}_{15}{a}_{16}({\lambda }_F^{**}+{\mu }_H)}},\hfill \\ \hfill R_H^{2**}& ={\displaystyle \frac{c_1{\psi }_1{\theta }_1{\theta }_2{\mathsf{\Pi }}_H{\lambda }_F^{**}}{a_{11}{a}_{12}{a}_{13}{a}_{14}{a}_{15}{a}_{16}{a}_{38}({\lambda }_F^{**}+{\mu }_H)}},\hfill \\ \hfill I_H^{S**}& ={\displaystyle \frac{a_{22}{\theta }_1{\theta }_2{\mathsf{\Pi }}_H{\lambda }_F^{**}}{a_{11}{a}_{12}{a}_{13}{a}_{14}({\lambda }_F^{**}+{\mu }_H)}},\hfill \\ \hfill I_H^{P**}& ={\displaystyle \frac{c_1\delta \psi {\theta }_1{\theta }_2{\mathsf{\Pi }}_H{\lambda }_F}{a_{11}{a}_{12}{a}_{13}{a}_{14}{a}_{15}{a}_{16}{a}_{17}{a}_{38}({\lambda }_F^{**}+{\mu }_H)}},\hfill \\ \hfill R_H^{1**}& ={\displaystyle \frac{(c_2+c_1\delta \varphi +c_3){\psi }_1{\theta }_1{\theta }_2{\mathsf{\Pi }}_H{\lambda }_F^{**}}{a_{11}{a}_{12}{a}_{13}{a}_{14}{a}_{15}{a}_{16}{a}_{17}{a}_{37}{a}_{38}({\lambda }_F^{**}+{\mu }_H)}},\hfill \\ \hfill R_H^{L**}& ={\displaystyle \frac{c_4{\lambda }_F^{**}+c_5{\lambda }_F^{**}(c_2+c_1\delta \varphi +c_3)}{{\mu }_H{a}_{11}{a}_{12}{a}_{13}{a}_{14}{a}_{15}{a}_{16}{a}_{17}{a}_{37}{a}_{38}({\lambda }_F^{**}+{\mu }_H)}},\hfill \\ \hfill S_A^{**}& ={\displaystyle \frac{{\mathsf{\Pi }}_A}{{\eta }_F^{**}+{\mu }_A}}\hfill \\ \hfill E_A^{**}& ={\displaystyle \frac{{\mathsf{\Pi }}_A{\eta }_F^{**}}{a_{18}({\eta }_F^{**}+{\mu }_A)}},\hfill \\ \hfill I_A^{A**}& ={\displaystyle \frac{b_{14}{\mathsf{\Pi }}_A{\eta }_F^{**}}{a_{18}{a}_{19}({\eta }_F^{**}+{\mu }_A)}},\hfill \\ \hfill I_A^{S**}& ={\displaystyle \frac{(a_{19}{b}_{15}+b_{14}{b}_{16}){\mathsf{\Pi }}_A{\eta }_F^{**}}{a_{18}{a}_{19}{a}_{20}({\eta }_F^{**}+{\mu }_A)}},\hfill \\ \hfill R_A^{**}& ={\displaystyle \frac{[a_{20}{a}_{24}{b}_{19}+\kappa (a_{19}{b}_{15}+b_{14}{b}_{16})]{\mathsf{\Pi }}_A{\eta }_F^{**}}{a_{18}{a}_{19}{a}_{20}{\mu }_A({\eta }_F^{**}+{\mu }_A)}},\hfill \\ \hfill S_F^{**}& ={\displaystyle \frac{{\mathsf{\Pi }}_F}{{\rho }_A^{**}+{\zeta }_A^{**}+{\xi }_A^{**}-{\mu }_F}},\hfill \\ \hfill E_F^{**}& ={\displaystyle \frac{({\rho }_A^{**}+{\zeta }_A^{**}+{\xi }_A^{**}){\mathsf{\Pi }}_F}{a_{21}({\rho }_A^{**}+{\zeta }_A^{**}+{\xi }_A^{**}-{\mu }_F)}},\hfill \\ \hfill I_F^{**}& ={\displaystyle \frac{({\rho }_A^{**}+{\zeta }_A^{**}+{\xi }_A^{**})\pi {\mathsf{\Pi }}_F}{{\mu }_F{a}_{21}({\rho }_A^{**}+{\zeta }_A^{**}+{\xi }_A^{**}-{\mu }_F)}},\hfill \end{array}$$

(54)

where,

$$\begin{array}{c}\hfill \begin{array}{c}{a}_{27}={\tau }_1{\psi }_1,a_{28}={\varrho }_1{\psi }_4,a_{29}={\omega }_3{\psi }_2,a_{30}={\omega }_1{\psi }_1,a_{31}={\gamma }_2{\psi }_3,\hfill \\ a_{32}={\tau }_3{\psi }_1,a_{33}={\gamma }_1{\psi }_3,a_{34}={\nu }_1{\zeta }_1,a_{35}=ϵ+{\mu }_H,a_{36}=\varphi +{\mu }_H,\hfill \\ a_{37}={\varrho }_1+{\mu }_H,a_{38}=\delta +{\mu }_H\hfill \\ c_1={\tau }_3{a}_{14}{a}_{15}{a}_{16}+a_{23}{\tau }_2\sigma +a_{23}{a}_{33}{\tau }_2\sigma ,\hfill \\ c_2={\psi }_2\sigma a_{16}{a}_{17}{a}_{22}{a}_{38},\hfill \\ c_3=a_{23}{a}_{31}{\tau }_2\sigma .\hfill \\ c_4=ϵ{\psi }_1{\theta }_1{\theta }_2{a}_{14}{a}_{15}{a}_{16}{a}_{17}{a}_{37}{a}_{38}{\mathsf{\Pi }}_H,\hfill \\ c_5=a_{28}{\psi }_1{\theta }_1{\theta }_2{\mathsf{\Pi }}_H.\hfill \end{array}\end{array}$$

4.6. Local Stability Analysis of VL Endemic Fixed Point

We rely on the linearization approach to establish the local asymptotic stability of the VL endemic fixed point, that is, finding the eigenvalues of the linearized system around the fixed point.

Theorem 3.

If $R_0>1$, then the VL endemic fixed point $E^{**}$ of system (1) is locally asymptotically stable (LAS).

Proof. 

Linearize the model system around the VL endemic fixed point. Also, to contain the elements of the matrix of the Jacobian, we redefine the parameter representations. Thus, we have the Jacobian matrix at the endemic fixed point $E^{**}$ denoted by $J_{E^{**}}$ as

$$J_{E^{**}}=\left(\begin{array}{cccccccccccccccccccc}{e}_1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& e_{11}\\ {\lambda }_F^{**}& e_2& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& m_4\\ 0& {\theta }_1& e_3& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& {\theta }_2& e_4& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& a_{27}& e_5& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& ϵ& e_6& 0& 0& 0& a_{28}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& a_{22}& 0& 0& e_7& 0& 0& 0& a_{28}& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& \sigma & e_8& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& a_{29}& e_9& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& e_{10}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& a_{30}& a_{31}& \varphi & e_{12}& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& a_{32}& 0& 0& 0& a_{23}& a_{32}& 0& 0& e_{13}& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& e_{14}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& {\eta }_F^{**}& e_{15}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& a_{24}& e_{16}& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& a_{33}& a_{26}& e_{17}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& a_{34}& \kappa & e_{18}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& e_{22}& e_{23}& e_{23}& 0& e_{19}& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& b_3& b_2& b_2& 0& 0& e_{20}& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& \pi & e_{21}\end{array}\right),$$

(55)

where

$$\begin{array}{c}\hfill \begin{array}{c}{e}_1=-m_1,e_2=-a_{11},e_3=-a_{12},e_4=-a_{13},e_5=-a_{35},e_6=-{\mu }_H,e_7=-a_{14},\hfill \\ e_8=-a_{15},e_9=-a_{16},e_{10}=-a_{17},e_{11}=-m_4,e_{12}=-a_{37},e_{13}=-a_{38},\hfill \\ e_{14}=-m_2,e_{15}=-a_{18},e_{16}=-a_{19},e_{17}=-a_{20},e_{18}=-{\mu }_A,e_{19}=-m_3,\hfill \\ e_{20}=-a_{21},e_{21}=-{\mu }_F,e_{22}=-b_3,e_{23}=-b_2.\hfill \end{array}\end{array}$$

The characteristic equation associated with the model system (1) is given by $|J_{E^{**}}-\lambda I|=0$, where $\lambda (={\lambda }_i,i=1,2,3,\dots ,20)$ denotes the eigenvalues of the $J_{E^{*}*}$, $I=I_{20\times 20}$ identity matrix.

This implies

$$|J_{E^{*}*}-\lambda I|=0.$$

(56)

Thus,

$$\begin{array}{c}\hfill \begin{array}{c}(e_1-{\lambda }_1)(e_2-{\lambda }_2)(e_3-{\lambda }_3)(e_4-{\lambda }_4)(e_5-{\lambda }_5)(e_6-{\lambda }_6)(e_7-{\lambda }_7)(e_8-{\lambda }_8)(e_9-{\lambda }_9)(e_{10}-{\lambda }_{10})\hfill \\ (e_{12}-{\lambda }_{11})(e_{13}-{\lambda }_{12})(e_{14}-{\lambda }_{13})(e_{15}-{\lambda }_{14})(e_{16}-{\lambda }_{15})(e_{17}-{\lambda }_{16})(e_{18}-{\lambda }_{17})(e_{19}-{\lambda }_{18})\hfill \\ (e_{20}-{\lambda }_{19})(e_{21}-{\lambda }_{20})=0.\hfill \end{array}\end{array}$$

This implies

${\lambda }_1=e_1=-m_1=-({\lambda }_F^{**}+{\mu }_H)<0\hspace{1em}\hspace{1em}{\lambda }_{11}=e_{12}=-a_{37}=-({\varrho }_1{\psi }_4+{\mu }_H)<0$,

${\lambda }_2=e_2=-a_{11}=-({\theta }_1+{\mu }_H)<0\hspace{1em}\hspace{1em}\hspace{1em}{\lambda }_{12}=e_{13}=-a_{38}=-(\delta +{\mu }_H)<0$,

${\lambda }_3=e_3=-a_{12}=-({\theta }_2+{\mu }_H)<0\hspace{1em}\hspace{1em}\hspace{1em}{\lambda }_{13}=e_{14}=-m_2=-({\eta }_F^{**}+{\mu }_H)<0$,

${\lambda }_4=e_4=-a_{13}=-({\psi }_1+{\mu }_H)<0\hspace{1em}\hspace{1em}\hspace{1em}{\lambda }_{14}=e_{15}=-a_{18}=-({\zeta }_1+{\mu }_A)<0$,

${\lambda }_5=e_5=-a_{35}=-(ϵ+{\mu }_H)<0\hspace{1em}\hspace{1em}\hspace{1em}{\lambda }_{15}=e_{16}=-a_{19}=-({\zeta }_2+{\mu }_A)<0$,

${\lambda }_6=e_6=-{\mu }_H<0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\lambda }_{16}=e_{17}=-a_{20}=-(\kappa +{\mu }_A)<0$,

${\lambda }_7=e_7=-a_{14}=-(\sigma +{\mu }_H)<0\hspace{1em}\hspace{1em}\hspace{1em}{\lambda }_{17}=e_{18}=-{\mu }_A<0$,

${\lambda }_8=e_8=-a_{15}=-({\psi }_2+{\mu }_H)<0\hspace{1em}\hspace{1em}\hspace{1em}{\lambda }_{18}=e_{19}=-m_3<0$,

${\lambda }_9=e_9=-a_{16}=-({\psi }_3+{\mu }_H)<0\hspace{1em}\hspace{1em}\hspace{1em}{\lambda }_{19}=e_{20}=-a_{21}=-(\pi +{\mu }_A)<0$,

${\lambda }_{10}=e_{10}=-a_{17}=-(\varphi +{\mu }_H)<0\hspace{1em}\hspace{1em}\hspace{1em}{\lambda }_{20}=e_{21}=-{\mu }_F<0$,

${\lambda }_i,i=1,2,3,\dots ,20$ are all negative.

We therefore conclude that the VL endemic fixed point $E^{**}$ is locally asymptotically stable. Hence, the proof is complete. □

This result indicates that the VL disease is likely to spread quickly and the number of infected individuals may rise exponentially if $R_0$ remains greater than unity for some time $t>0$.

4.7. Global Stability Analysis of VL Endemic Fixed Point

Here, we explore the global stability of the VL endemic fixed point $E^{**}$ to help understand the spread of the VL disease in a population.

Theorem 4.

There is no periodic orbit for the system (1).

Proof. 

Applying the Dulac’s criterion as used in [24], let

$$\begin{array}{c}\hfill X=(S_H,E_H,I_H^E,I_H^L,R_H^E,R_H^L,I_H^S,I_H^1,I_H^2,I_H^P,R_H^1,R_H^2,S_A,E_A,I_A^A,I_A^S,R_A,S_F,E_F,I_F).\end{array}$$

Taking the Dulac’s function as

$$G={\displaystyle \frac{1}{I_F}},$$

from the system (1) equations, we have

$$\begin{array}{cc}\hfill G{\displaystyle \frac{d{S}_H}{dt}}& ={\displaystyle \frac{{\mathsf{\Pi }}_H}{I_F}}-{\lambda }_F{\displaystyle \frac{S_H}{I_F}}-{\mu }_H{\displaystyle \frac{S_H}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{E}_H}{dt}}& ={\lambda }_F{\displaystyle \frac{S_H}{I_F}}-({\theta }_1+{\mu }_H){\displaystyle \frac{E_H}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{I}_H^E}{dt}}& ={\theta }_1{\displaystyle \frac{E_H}{I_F}}-({\theta }_2+{\mu }_H){\displaystyle \frac{I_H^E}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{I}_H^L}{dt}}& ={\theta }_2{\displaystyle \frac{I_H^E}{I_F}}-({\psi }_1+{\mu }_H){\displaystyle \frac{I_H^L}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{R}_H^E}{dt}}& ={\tau }_1{\psi }_1{\displaystyle \frac{I_H^L}{I_F}}-(\epsilon +{\mu }_H){\displaystyle \frac{R_H^E}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{R}_H^L}{dt}}& =\epsilon {\displaystyle \frac{R_H^E}{I_F}}+{\varrho }_1{\psi }_4{\displaystyle \frac{R_H^1}{I_F}}-{\mu }_H{\displaystyle \frac{R_H^L}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{I}_H^S}{dt}}& ={\tau }_2{\psi }_1{\displaystyle \frac{I_H^L}{I_F}}-(\sigma +{\mu }_H){\displaystyle \frac{I_H^S}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{I}_H^1}{dt}}& =\sigma {\displaystyle \frac{I_H^S}{I_F}}-({\psi }_2+{\mu }_H){\displaystyle \frac{I_H^1}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{I}_H^2}{dt}}& ={\omega }_3{\psi }_2{\displaystyle \frac{I_H^1}{I_F}}-({\psi }_3+{\mu }_H){\displaystyle \frac{I_H^2}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{I}_H^P}{dt}}& =\delta {\displaystyle \frac{R_H^2}{I_F}}-(\varphi +{\mu }_H){\displaystyle \frac{I_H^P}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{R}_H^1}{dt}}& ={\omega }_1{\psi }_2{\displaystyle \frac{I_H^1}{I_F}}+\varphi {\displaystyle \frac{I_H^P}{I_F}}+{\gamma }_2{\psi }_3{\displaystyle \frac{I_H^2}{I_F}}-({\varrho }_1{\psi }_4+{\mu }_H){\displaystyle \frac{R_H^1}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{R}_H^2}{dt}}& ={\tau }_3{\psi }_1{\displaystyle \frac{I_H^L}{I_F}}+{\omega }_2{\psi }_2{\displaystyle \frac{I_H^1}{I_F}}+{\gamma }_1{\psi }_3{\displaystyle \frac{I_H^2}{I_F}}-(\delta +{\mu }_H){\displaystyle \frac{R_H^2}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{S}_A}{dt}}& ={\displaystyle \frac{{\mathsf{\Pi }}_A}{I_F}}-{\eta }_F{\displaystyle \frac{S_A}{I_F}}-{\mu }_A{\displaystyle \frac{S_A}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{E}_A}{dt}}& ={\eta }_F{\displaystyle \frac{S_A}{I_F}}-({\zeta }_1+{\mu }_A){\displaystyle \frac{E_A}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{I}_A^A}{dt}}& ={\varpi }_1{\zeta }_1{\displaystyle \frac{E_A}{I_F}}-({\zeta }_2+{\mu }_A){\displaystyle \frac{I_A^A}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{I}_A^S}{dt}}& ={\varpi }_2{\zeta }_1{\displaystyle \frac{E_A}{I_F}}+{\nu }_2{\zeta }_2{\displaystyle \frac{I_A^A}{I_F}}-(\kappa +{\mu }_A){\displaystyle \frac{I_A^S}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{R}_A}{dt}}& ={\nu }_1{\zeta }_2{\displaystyle \frac{I_A^A}{I_F}}+\kappa {\displaystyle \frac{I_A^S}{I_F}}-{\mu }_A{\displaystyle \frac{R_A}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{S}_F}{dt}}& ={\displaystyle \frac{{\mathsf{\Pi }}_F}{I_F}}-({\rho }_A+{\xi }_A+{\chi }_A){\displaystyle \frac{S_F}{I_F}}-{\mu }_F{\displaystyle \frac{S_F}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{E}_F}{dt}}& =({\rho }_A+{\xi }_A+{\chi }_A){\displaystyle \frac{S_F}{I_F}}-(\pi +{\mu }_F){\displaystyle \frac{E_F}{I_F}},\hfill \\ \hfill G{\displaystyle \frac{d{I}_F}{dt}}& =\pi {\displaystyle \frac{E_F}{I_F}}-{\mu }_F.\hfill \end{array}$$

Thus,

$$\begin{array}{cc}\hfill {\displaystyle \frac{dGX}{dt}}& ={\displaystyle \frac{\partial }{d{S}_H}}\left(G{\displaystyle \frac{d{S}_H}{dt}}\right)+{\displaystyle \frac{\partial }{d{E}_H}}\left(G{\displaystyle \frac{d{E}_H}{dt}}\right)+{\displaystyle \frac{\partial }{d{I}_H^E}}\left(G{\displaystyle \frac{d{I}_H^E}{dt}}\right)+{\displaystyle \frac{\partial }{dd{I}_H^L}}\left(G{\displaystyle \frac{d{I}_H^L}{dt}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{R}_H^E}}\left(G{\displaystyle \frac{d{R}_H^E}{dt}}\right)+{\displaystyle \frac{\partial }{d{R}_H^L}}\left(G{\displaystyle \frac{d{R}_H^L}{dt}}\right)+{\displaystyle \frac{\partial }{d{I}_H^S}}\left(G{\displaystyle \frac{d{I}_H^S}{dt}}\right)+{\displaystyle \frac{\partial }{d{I}_H^1}}\left(G{\displaystyle \frac{d{I}_H^1}{dt}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{I}_H^2}}\left(G{\displaystyle \frac{d{I}_H^2}{dt}}\right)+{\displaystyle \frac{\partial }{d{I}_H^P}}\left(G{\displaystyle \frac{d{I}_H^P}{dt}}\right)+{\displaystyle \frac{\partial }{d{R}_H^1}}\left(G{\displaystyle \frac{d{R}_H^1}{dt}}\right)+{\displaystyle \frac{\partial }{d{R}_H^2}}\left(G{\displaystyle \frac{d{R}_H^2}{dt}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{S}_A}}\left(G{\displaystyle \frac{d{S}_A}{dt}}\right)+{\displaystyle \frac{\partial }{d{E}_A}}\left(G{\displaystyle \frac{d{E}_A}{dt}}\right)+{\displaystyle \frac{\partial }{d{I}_A^A}}\left(G{\displaystyle \frac{d{I}_A^A}{dt}}\right)+{\displaystyle \frac{\partial }{d{I}_A^S}}\left(G{\displaystyle \frac{d{I}_A^S}{dt}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{R}_A}}\left(G{\displaystyle \frac{d{R}_A}{dt}}\right)+{\displaystyle \frac{\partial }{d{S}_F}}\left(G{\displaystyle \frac{d{S}_F}{dt}}\right)+{\displaystyle \frac{\partial }{d{E}_F}}\left(G{\displaystyle \frac{d{E}_F}{dt}}\right)+{\displaystyle \frac{\partial }{d{I}_F}}\left(G{\displaystyle \frac{d{I}_F}{dt}}\right),\hfill \\ \hfill & ={\displaystyle \frac{\partial }{d{S}_H}}\left({\displaystyle \frac{{\mathsf{\Pi }}_H}{I_F}}-{\lambda }_F{\displaystyle \frac{S_H}{I_F}}-{\mu }_H{\displaystyle \frac{S_H}{I_F}}\right)+{\displaystyle \frac{\partial }{d{E}_H}}\left({\lambda }_F{\displaystyle \frac{S_H}{I_F}}-({\theta }_1+{\mu }_H){\displaystyle \frac{E_H}{I_F}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{I}_H^E}}\left({\theta }_1{\displaystyle \frac{E_H}{I_F}}-({\theta }_2+{\mu }_H){\displaystyle \frac{I_H^E}{I_F}}\right)+{\displaystyle \frac{\partial }{d{I}_H^L}}\left({\theta }_2{\displaystyle \frac{I_H^E}{I_F}}-({\psi }_1+{\mu }_H){\displaystyle \frac{I_H^L}{I_F}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{R}_H^E}}\left({\tau }_1{\psi }_1{\displaystyle \frac{I_H^L}{I_F}}-(\epsilon +{\mu }_H){\displaystyle \frac{R_H^E}{I_F}}\right)+{\displaystyle \frac{\partial }{d{R}_H^L}}\left(\epsilon {\displaystyle \frac{R_H^E}{I_F}}+{\varrho }_1{\psi }_4{\displaystyle \frac{R_H^1}{I_F}}-{\mu }_H{\displaystyle \frac{R_H^L}{I_F}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{I}_H^S}}\left({\tau }_2{\psi }_1{\displaystyle \frac{I_H^L}{I_F}}-(\sigma +{\mu }_H){\displaystyle \frac{I_H^S}{I_F}}\right)+{\displaystyle \frac{\partial }{d{I}_H^1}}\left(\sigma {\displaystyle \frac{I_H^S}{I_F}}-({\psi }_2+{\mu }_H){\displaystyle \frac{I_H^1}{I_F}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{I}_H^2}}\left({\omega }_3{\psi }_2{\displaystyle \frac{I_H^1}{I_F}}-({\psi }_3+{\mu }_H){\displaystyle \frac{I_H^2}{I_F}}\right)+{\displaystyle \frac{\partial }{d{I}_H^P}}\left(\delta {\displaystyle \frac{R_H^2}{I_F}}-(\varphi +{\mu }_H){\displaystyle \frac{I_H^P}{I_F}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{R}_H^1}}\left({\omega }_1{\psi }_2{\displaystyle \frac{I_H^1}{I_F}}+\varphi {\displaystyle \frac{I_H^P}{I_F}}+{\gamma }_2{\psi }_3{\displaystyle \frac{I_H^2}{I_F}}-({\varrho }_1{\psi }_4+{\mu }_H){\displaystyle \frac{R_H^1}{I_F}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{R}_H^2}}\left({\tau }_3{\psi }_1{\displaystyle \frac{I_H^L}{I_F}}+{\omega }_2{\psi }_2{\displaystyle \frac{I_H^1}{I_F}}+{\gamma }_1{\psi }_3{\displaystyle \frac{I_H^2}{I_F}}-(\delta +{\mu }_H){\displaystyle \frac{R_H^2}{I_F}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{S}_A}}\left({\displaystyle \frac{{\mathsf{\Pi }}_A}{I_F}}-{\eta }_F{\displaystyle \frac{S_A}{I_F}}-{\mu }_A{\displaystyle \frac{S_A}{I_F}}\right)+{\displaystyle \frac{\partial }{d{E}_A}}\left({\eta }_F{\displaystyle \frac{S_A}{I_F}}-({\zeta }_1+{\mu }_A){\displaystyle \frac{E_A}{I_F}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{I}_A^A}}\left({\varpi }_1{\zeta }_1{\displaystyle \frac{E_A}{I_F}}-({\zeta }_2+{\mu }_A){\displaystyle \frac{I_A^A}{I_F}}\right)+{\displaystyle \frac{\partial }{d{I}_A^S}}\left({\varpi }_2{\zeta }_1{\displaystyle \frac{E_A}{I_F}}+{\nu }_2{\zeta }_2{\displaystyle \frac{I_A^A}{I_F}}-(\kappa +{\mu }_A){\displaystyle \frac{I_A^S}{I_F}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{R}_A}}\left({\nu }_1{\zeta }_2{\displaystyle \frac{I_A^A}{I_F}}+\kappa {\displaystyle \frac{I_A^S}{I_F}}-{\mu }_A{\displaystyle \frac{R_A}{I_F}}\right)+{\displaystyle \frac{\partial }{d{S}_F}}\left({\displaystyle \frac{{\mathsf{\Pi }}_F}{I_F}}-({\rho }_A+{\xi }_A+{\chi }_A){\displaystyle \frac{S_F}{I_F}}-{\mu }_F{\displaystyle \frac{S_F}{I_F}}\right)\hfill \\ & +{\displaystyle \frac{\partial }{d{E}_F}}\left(({\rho }_A+{\xi }_A+{\chi }_A){\displaystyle \frac{S_F}{I_F}}-(\pi +{\mu }_F){\displaystyle \frac{E_F}{I_F}}\right)+{\displaystyle \frac{\partial }{d{I}_F}}\left(\pi {\displaystyle \frac{E_F}{I_F}}-{\mu }_F\right),\hfill \\ \hfill & =\left(-{\lambda }_F{\displaystyle \frac{1}{I_F}}-{\mu }_H{\displaystyle \frac{1}{I_F}}\right)+\left(-({\theta }_1+{\mu }_H){\displaystyle \frac{1}{I_F}}\right)+\left(-({\theta }_2+{\mu }_H){\displaystyle \frac{1}{I_F}}\right)\hfill \\ & +\left(-({\psi }_1+{\mu }_H){\displaystyle \frac{1}{I_F}}\right)+\left(-(\epsilon +{\mu }_H){\displaystyle \frac{1}{I_F}}\right)+\left(-{\mu }_H{\displaystyle \frac{1}{I_F}}\right)+\left(-(\sigma +{\mu }_H){\displaystyle \frac{1}{I_F}}\right)\hfill \\ & +\left(-({\psi }_2+{\mu }_H){\displaystyle \frac{1}{I_F}}\right)+\left(-({\psi }_3+{\mu }_H){\displaystyle \frac{1}{I_F}}\right)+\left(-(\varphi +{\mu }_H){\displaystyle \frac{1}{I_F}}\right)\hfill \\ & +\left(-({\varrho }_1{\psi }_4+{\mu }_H){\displaystyle \frac{1}{I_F}}\right)+\left(-(\delta +{\mu }_H){\displaystyle \frac{1}{I_F}}\right)+\left(-{\eta }_F{\displaystyle \frac{1}{I_F}}-{\mu }_A{\displaystyle \frac{1}{I_F}}\right)\hfill \\ & +\left(-({\zeta }_1+{\mu }_A){\displaystyle \frac{1}{I_F}}\right)+\left(-({\zeta }_2+{\mu }_A){\displaystyle \frac{1}{I_F}}\right)+\left(-(\kappa +{\mu }_A){\displaystyle \frac{1}{I_F}}\right)\hfill \\ & +\left(-{\mu }_A{\displaystyle \frac{1}{I_F}}\right)+\left(-({\rho }_A+{\xi }_A+{\chi }_A){\displaystyle \frac{1}{I_F}}-{\mu }_F{\displaystyle \frac{1}{I_F}}\right)+\left(-(\pi +{\mu }_F){\displaystyle \frac{1}{I_F}}\right)\hfill \\ & +\left(-\pi E_F{I}_F^{-2}\right),\hfill \\ \hfill & =-({\lambda }_F+{\eta }_F+{\rho }_A+{\xi }_A+{\chi }_A+\delta +\kappa +12{\mu }_H+5{\mu }_A+2{\mu }_F+{\theta }_1+{\theta }_2\hfill \\ & +\pi +{\psi }_1+{\psi }_2+{\psi }_3+{\rho }_1{\psi }_4+\epsilon +\sigma +\varphi +{\zeta }_1+{\zeta }_2){\displaystyle \frac{1}{I_F}}-\pi {\displaystyle \frac{E_F}{I_F^2}},\hfill \\ \hfill & <0.\hfill \end{array}$$

Hence, there exists no periodic solution for the system (1). □

With $\mathsf{\Omega }$ being positively invariant, all solutions to the system (1) originate and remain in $\mathsf{\Omega }$ for all t. This is supported by the Poincaré–Bendixson theorem. Consequently, the theorem follows.

Theorem 5.

The VL endemic fixed point $E^{**}$ for the system (1) is globally asymptotically stable whenever $R_0>1$.

This result indicates that the disease may not die out whenever the threshold $R_0>1$.

4.8. Sensitivity Analysis

This subsection presents the analytical sensitivity analysis of some of the important parameters in $R_0$ (i.e., ${\beta }_A,{\beta }_F,g_A$, and $\kappa$) using the differential approach as applied by [25] to quantify their impact on the model dynamic. Moreover, the $R_0$ important parameters’ numerical and analytical values are produced using precise assumptions based on parameter values. Some insight on tracking the model’s onset at variant places is given by the analytical expressions produced. By reducing the value to less than unity, the threshold value $R_0$ is acknowledged as the primary method of stopping the disease’s progress. Positively signed parameters are considered highly and proportionally sensitive to $R_0$, whereas negatively signed parameters are not declining $R_0$. This is obtained using partial derivatives of $R_0$ with respect to the key parameter. The following are valid based on the expression for $R_0$ provided $D_8{D}_{10}<D_9{D}_{11}$:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{\partial R_0}{\partial {\beta }_A}}& =& {\displaystyle \frac{D_1\sqrt[]{({\beta }_A{D}_1+D_2)D_3{D}_4}}{2({\beta }_A{D}_1+D_2)D_3{D}_4}}>0,\hfill \\ \hfill {\displaystyle \frac{\partial R_0}{\partial {\beta }_F}}& =& {\displaystyle \frac{D_5\sqrt[]{{\beta }_F{D}_3{D}_4{D}_5}}{2{\beta }_F{D}_3{D}_4}}>0,\hfill \\ \hfill {\displaystyle \frac{\partial R_0}{\partial g_A}}& =& {\displaystyle \frac{(D_6+2g_A{D}_7)\sqrt[]{(g_A{D}_6+g_A^2{D}_7)D_3{D}_4}}{2D_3{D}_4(g_A{D}_6+g_A^2{D}_7)}}>0,\hfill \\ \hfill {\displaystyle \frac{\partial R_0}{\partial \kappa }}& =& {\displaystyle \frac{k(D_8{D}_{10}-D_9{D}_{11})\sqrt[]{D_3(\kappa D_9+D_{10})(\kappa D_8+D_{11})}}{(\kappa D_8+D_{11})(\kappa D_9+D_{10})}}<0.\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill D_1& ={\beta }_F{g}_A^2{\mu }_A{\mu }_H^2{\mathsf{\Pi }}_{A}P{i}_F,\hfill \\ \hfill D_2& ={\beta }_F{g}_A{h}_A\kappa {\mu }_A{\mathsf{\Pi }}_F({\zeta }_1{\mu }_F{\mathsf{\Pi }}_H+{\zeta }_1{\mu }_H{\mathsf{\Pi }}_A+m{u}_A^2{\mathsf{\Pi }}_H+{\mu }_A{\mu }_H{\mathsf{\Pi }}_A)\hfill \\ & +{\beta }_F{g}_A{h}_A{\mu }_A^3{\mu }_F{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H{\zeta }_1+{\beta }_F{g}_A{h}_A{\mu }_A^3{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F{\zeta }_1+{\beta }_F{g}_A{h}_A{\mu }_A^4{\mu }_F{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H\hfill \\ & +{\beta }_F{g}_A{h}_A{\mu }_A^3{\mu }_F{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F{\zeta }_1+{\beta }_F{g}_A{\mu }_A^2{\mu }_H^2{\nu }_2{\mathsf{\Pi }}_F{\varrho }_1{\zeta }_1{\zeta }_2,\hfill \\ \hfill D_3& ={\mu }_A^2{\mu }_F{\mathsf{\Pi }}_H^2+2{\mu }_A{\mu }_F{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_H+{\mu }_F{\mu }_H,\hfill \\ \hfill D_4& =\kappa {\mu }_F\pi {\zeta }_1{\zeta }_2+\kappa {\mu }_A{\mu }_F{\zeta }_1{\zeta }_2+\kappa {\mu }_A{\mu }_F\pi {\zeta }_1+\kappa {\mu }_A^2{\mu }_F{\zeta }_1+\kappa {\mu }_A\pi {\zeta }_2+\kappa {\mu }_A^2{\mu }_F{\zeta }_2\hfill \\ & +\kappa {\mu }_A^2{\mu }_F\pi +\kappa {\mu }_A^3{\mu }_F+{\mu }_A{\mu }_F\pi {\zeta }_1{\zeta }_2+{\mu }_A^2{\mu }_F{\zeta }_1{\zeta }_2+{\mu }_A^2{\mu }_F\pi {\zeta }_1+{\mu }_A^3{\mu }_F{\zeta }_1{\zeta }_2\hfill \\ & +{\mu }_A^2{\mu }_F\pi {\zeta }_2+{\mu }_A^3{\mu }_F{\zeta }_2+{\mu }_A^3{\mu }_F\pi +{\mu }_A^4{\mu }_F,\hfill \\ \hfill D_5& =\kappa g_A{h}_A{\mu }_A{\mu }_F{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H{\zeta }_1+\kappa g_A{h}_A{\mu }_A{\mu }_F{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F{\zeta }_1+{\beta }_F\kappa g_H{h}_A{\mu }_A^3{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H\hfill \\ & +\kappa g_A{h}_A{\mu }_A^2{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F+g_A{h}_A{\mu }_A^3{\mu }_F{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H{\zeta }_1+{\beta }_F{g}_A{h}_A{\mu }_A^3{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F{\zeta }_1\hfill \\ & +g_A{h}_A{\mu }_A^4{\mu }_F{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H+g_A{h}_A{\mu }_A^3{\mu }_F{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F+{\beta }_A{g}_A^2{\mu }_A{\mu }_H^2{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F\hfill \\ & +g_A{\mu }_A^2{\mu }_H^2{\nu }_2{\mathsf{\Pi }}_F{\varrho }_1{\zeta }_1{\zeta }_2+g_A{\mu }_A{\mu }_H^2{\nu }_2{\mathsf{\Pi }}_F{\varrho }_1{\zeta }_1{\zeta }_2,\hfill \\ \hfill D_6& ={\beta }_F\kappa h_A{\mu }_A{\mu }_F{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H{\zeta }_1+{\beta }_F\kappa h_A{\mu }_A{\mu }_F{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F{\zeta }_1+{\beta }_F\kappa h_A{\mu }_A^2{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F\hfill \\ & +{\beta }_F\kappa h_A{\mu }_A^3{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H+{\beta }_F{h}_A{\mu }_A^4{\mu }_F{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H+{\beta }_F{h}_A{\mu }_A^3{\mu }_F{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F\hfill \\ & +{\beta }_F{h}_A{\mu }_A^3{\mu }_F{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F{\zeta }_1+{\beta }_F{h}_A{\mu }_A^2{\mu }_F{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F{\zeta }_1\hfill \\ & +{\beta }_F{\mu }_A^2{\mu }_H^2{\nu }_2{\mathsf{\Pi }}_F{\varrho }_1{\zeta }_1{\zeta }_2+{\beta }_F{\mu }_A{\mu }_H^2{\nu }_2{\mathsf{\Pi }}_F{\varrho }_1{\zeta }_1{\zeta }_2,\hfill \\ \hfill D_7& ={\beta }_A{\beta }_F{\mu }_A{\mu }_H^2{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F,\hfill \\ \hfill D_8& ={\beta }_F{g}_A{h}_A{\mu }_A{\mu }_F{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H{\zeta }_1+{\beta }_F{g}_A{h}_A{\mu }_A{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F{\zeta }_1+{\beta }_F{g}_A{h}_A{\mu }_A^3{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H\hfill \\ & +{\beta }_F{g}_A{h}_A{\mu }_A^2{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_H,\hfill \\ \hfill D_9& ={\mu }_F\pi {\zeta }_1{\zeta }_2+{\mu }_A{\mu }_F{\zeta }_1{\zeta }_2+{\mu }_A{\mu }_F\pi {\zeta }_1+{\mu }_A^2{\mu }_F{\zeta }_1+{\mu }_A\pi {\zeta }_2+{\mu }_A^2{\mu }_F{\zeta }_2\hfill \\ & +{\mu }_A^2{\mu }_F\pi +{\mu }_A^3{\mu }_F,\hfill \\ \hfill D_{10}& ={\mu }_A{\mu }_F\pi {\zeta }_1{\zeta }_2+{\mu }_A^2{\mu }_F{\zeta }_1{\zeta }_2+{\mu }_A^2{\mu }_F{\zeta }_1+{\mu }_A^3{\mu }_F{\zeta }_1{\zeta }_2+{\mu }_A^2{\mu }_F\pi {\zeta }_2+{\mu }_A^3{\mu }_F{\zeta }_2\hfill \\ & +{\mu }_A^3{\mu }_F\pi +{\mu }_A^4{\mu }_F,\hfill \\ \hfill D_{11}& ={\beta }_F{g}_A{h}_A{\mu }_A^3{\mu }_F{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H{\zeta }_1+{\beta }_F{g}_A{h}_A{\mu }_A^3{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F{\zeta }_1+{\beta }_F{g}_A{h}_A{\mu }_A^4{\mu }_F{\mathsf{\Pi }}_F{\mathsf{\Pi }}_H\hfill \\ & +{\beta }_F{g}_A{h}_A{\mu }_A^4{\mu }_H{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F+{\beta }_A{\beta }_F{g}_A^2{\mu }_A{\mu }_H^2{\mathsf{\Pi }}_A{\mathsf{\Pi }}_F+{\beta }_F{g}_A{\mu }_A^2{\mu }_H^2{\nu }_2{\mathsf{\Pi }}_F{\varrho }_1{\zeta }_1{\zeta }_2\hfill \\ & +{\beta }_F{g}_A{\mu }_A{\mu }_H^2{\nu }_2{\mathsf{\Pi }}_F{\varrho }_1{\zeta }_1{\zeta }_2.\hfill \end{array}$$

The equations in (57) show that decreasing the transmission probability from infected sandflies to susceptible animals (${\beta }_A$), the transmission probability from infected animals to susceptible sandflies (${\beta }_F$), and the per capita biting rate of sandflies of animals ($g_A$) will help reduce the transmission of VL infection. Again, from (57), increasing the fraction of rate of transfer of symptomatic infected animals to the recovered class ($\kappa$) will also reduce the transmission of VL disease in a population.

The elasticity indices for all the $R_0$ parameters which is given by

$${\gamma }_q^{R_0}={\displaystyle \frac{\partial R_0}{\partial q}}.{\displaystyle \frac{q}{R_0}}.$$

(58)

where $R_0$ denotes the basic reproduction ratio and q is a parameter are analyzed numerically using Variable Precision Arithmetic (VPA) in MATLAB in the section for numerical simulation.

5. Numerical Simulation

In this section, the numerical outcomes of the various compartments of the model are observed using parameter values from the literature and some assumptions. All values of the parameters in

Table 1. The description of parameters in the ZVL model.

Parameter	Description
${\mathsf{\Pi }}_H$	Recruitment rate of humans
${\mathsf{\Pi }}_A$	Recruitment rate of animals
${\mathsf{\Pi }}_F$	Recruitment rate of sandflies
${\theta }_1$	Transfer rate of exposed humans to early VL infection stage
${\theta }_2$	Transfer rate of early asymptomatic VL-infected humans to late VL infection stage
${\tau }_1$	Proportion of late asymptomatic VL-infected humans moving to $R_H^E$
${\tau }_2$	Proportion of late asymptomatic VL-infected humans moving to $I_H^S$
${\tau }_3$	Proportion of late asymptomatic VL-infected humans moving to $R_H^2$ classes
${\psi }_1$	Rate of transfer of late asymptomatic infected humans to symptomatic infected class
${\psi }_2$	Rate of transfer of infected humans receiving first-line treatment to second-line treatment, to recovered humans who have cleared the parasite and to putative recovered human classes
${\psi }_3$	Rate of transfer of infected humans receiving second-line treatment to recovered humans who have cleared the parasite and to putative recovered human classes
${\omega }_1$	Proportions of infected humans receiving first-line treatment moving to recovered class who have cleared the parasite
${\omega }_2$	Proportions of infected humans receiving first-line treatment moving to putative recovered class
${\omega }_3$	Proportions of infected humans receiving first-line treatment moving to infected class receiving second-line treatment after first-line treatment failure
${\gamma }_1$	Proportion of infected humans receiving second-line treatment moving to putative recovered class
${\gamma }_2$	Proportion of infected humans receiving second-line treatment moving to recovered class who have cleared the parasite
${\varrho }_1$	Proportion of recovered humans who have cleared the parasite moving to recovered class who are DAT-positive and LST-positive
${\varrho }_2$	Proportion of recovered humans who have cleared the parasite moving to putative recovered class
${\mu }_H$	Natural death rates of humans
${\mu }_A$	Natural death rates of animals
${\mu }_F$	Natural death rates of sandflies
$\delta$	Rate of transfer from putative recovered humans to PKDL-infected class
$\epsilon$	Rate of transfer from recovered humans who are DAT-positive and not yet
	LST-positive to recovered humans who are DAT-positive and still LST-positive
$\varphi$	Rate of transfer of PKDL-infected humans to putative recovered class
$\phi$	Rate of transfer from recovered humans who are DAT-positive and but still
	LST-positive susceptible humans
$\sigma$	Rate of transfer of humans from symptomatic infected class to infected humans receiving first-line treatment class
${\varpi }_1$	Proportion of exposed animals moving to asymptomatic infected class
${\varpi }_2$	Proportion of exposed animals moving to symptomatic infected class
${\zeta }_1$	Rate of transfer of exposed animals to asymptomatic infected class and symptomatic infected class
${\zeta }_2$	Rate of transfer of asymptomatic infected animals to symptomatic infected class and recovered class
${\nu }_1$	Proportion of asymptomatic infected animals moving to symptomatic infected class
${\nu }_2$	Proportion of asymptomatic infected animals moving to recovered class
$\kappa$	Rate of transfer from symptomatic infected animals to recovered class
$\pi$	Rate of transfer of exposed sandflies to infected sandflies
${\beta }_H$	Transmission probability from infected sandflies to susceptible humans
${\beta }_A$	Transmission probability from infected sandflies to susceptible animals
${\beta }_F$	Transmission probability from infected animal to susceptible sandflies
	(kept constant for both asymptomatic and symptomatic infected class)
$g_H$	Per capita biting rate of sandflies of humans
$g_A$	Per capita biting rate of sandflies of animals (kept constant for both asymptomatic and symptomatic infected classes)
$h_A$	Modification parameter for the relative infectiousness of an animal

are given in

Table 2. Parameter values used for the simulation result.

Parameter	Value	Unit	Source
${\mathsf{\Pi }}_H$	19.5	${\mathrm{day}}^{-1}$	[17]
${\mathsf{\Pi }}_A$	8.33	${\mathrm{day}}^{-1}$	[26]
${\mathsf{\Pi }}_F$	210.62	${\mathrm{day}}^{-1}$	Assumed
${\theta }_1$	0.0111	${\mathrm{day}}^{-1}$	[27]
${\theta }_2$	0.01667	${\mathrm{day}}^{-1}$	[19]
${\tau }_1$	0.001	${\mathrm{day}}^{-1}$	Assumed
${\tau }_2$	0.699	${\mathrm{day}}^{-1}$	$1-({\tau }_1+{\tau }_3)$
${\tau }_3$	0.3	${\mathrm{day}}^{-1}$	[28,29]
${\psi }_1$	0.083	${\mathrm{day}}^{-1}$	[19]
${\psi }_2$	0.033	${\mathrm{day}}^{-1}$	[19]
${\psi }_3$	0.333	${\mathrm{day}}^{-1}$	[19]
${\omega }_1$	0.92	${\mathrm{day}}^{-1}$	$1-({\omega }_2+{\omega }_3)$
${\omega }_2$	0.03	${\mathrm{day}}^{-1}$	[19]
${\omega }_3$	0.05	${\mathrm{day}}^{-1}$	[19]
${\gamma }_1$	0.03	${\mathrm{day}}^{-1}$	[19]
${\gamma }_2$	0.97	${\mathrm{day}}^{-1}$	$1-{\gamma }_2$
${\varrho }_1$	0.0135	${\mathrm{day}}^{-1}$	[19]
${\mu }_H$	0.00795	${\mathrm{day}}^{-1}$	Assumed
${\mu }_A$	0.0169	${\mathrm{day}}^{-1}$	[17]
${\mu }_F$	0.14	${\mathrm{day}}^{-1}$	Assumed
$\delta$	0.00397	${\mathrm{day}}^{-1}$	[19]
$\epsilon$	0.0135	${\mathrm{day}}^{-1}$	[19]
$\varphi$	0.00556	${\mathrm{day}}^{-1}$	[19]
$\sigma$	1	${\mathrm{day}}^{-1}$	[19]
${\varpi }_1$	0.79	${\mathrm{day}}^{-1}$	Assumed
${\varpi }_2$	0.21	${\mathrm{day}}^{-1}$	Assumed
${\zeta }_1$	0.69	${\mathrm{day}}^{-1}$	Assumed
${\zeta }_2$	0.211	${\mathrm{day}}^{-1}$	Assumed
${\nu }_1$	0.97	${\mathrm{day}}^{-1}$	Assumed
${\nu }_2$	0.03	${\mathrm{day}}^{-1}$	$1-{\nu }_1$
$\kappa$	0.115	${\mathrm{day}}^{-1}$	[17]
$\pi$	0.2	${\mathrm{day}}^{-1}$	[17]
${\beta }_H$	1	${\mathrm{day}}^{-1}$	[19]
${\beta }_A$	1	${\mathrm{day}}^{-1}$	[19]
${\beta }_F$	1	${\mathrm{day}}^{-1}$	[19]
$g_H$	0.02856	${\mathrm{day}}^{-1}$	[30,31]
$g_A$	0.2856	${\mathrm{day}}^{-1}$	[17]
$h_A$	1.39	${\mathrm{day}}^{-1}$	[18]

. This is to gain insight into the behavior of the numerical solutions of the model.

The initial states considered for the model variables given in

Table 3. The description of variables in the VL model.

Variable	Description
$S_H$	Population of susceptible humans
$E_H$	Population of exposed humans
$I_H^E$	Population of infected humans at early asymptomatic stage
$I_H^L$	Population of infected humans at late asymptomatic stage
$I_H^S$	Population of symptomatic infected humans
$R_H^E$	Population of recovered humans who are DAT-positive and not yet LST-positive
$R_H^L$	Population of recovered humans who are DAT-positive and but still LST-positive
$I_H^P$	Population of infected humans at the PKDL stage
$I_H^1$	Population of infected humans who are receiving first-line treatment
$I_H^2$	Population of infected humans who are receiving second-line treatment
$R_H^1$	Population of recovered humans who have cleared the parasite
$R_H^2$	Population of putative recovered humans
$S_A$	Population of susceptible animals
$E_A$	Population of exposed animals
$I_A^A$	Population of asymptomatic infected animals
$I_A^S$	Population of symptomatic infected animals
$R_A$	Population of recovered animals
$S_F$	Population of susceptible sandflies (vector)
$E_F$	Population of exposed sandflies
$I_F$	Population of infected sandflies

are $S_H=250,000$, $E_H=3000I_H^E=2000I_H^L=1000I_H^S=900,R_H^E=200,R_H^L=150,I_H^1=150I_H^2=200$, $I_H^P=120,R_H^1=100,R_H^2=110,S_A=2500,E_A=700,I_A^A=200,I_A^S=190$, $R_A=180$, $S_F=10,000E_F=3500,I_F=2000$. The elasticity indices of the parameters of $R_0$ are obtained numerically in MATLAB2014a using Variable Precision Arithmetic (VPA) for precision. By employing high-precision arithmetic to reduce errors in the crucial computations, this will ensure precise assessment of the basic reproduction number’s sensitivity to changes in a particular parameter. It is adjudged that VPA makes the computations more precise, which is important when little adjustments to the parameters have a big impact on $R_0$. The elasticity index is computed more precisely using the VPA than using the standard floating-point arithmetic [32,33,34]. This is required because normal floating-point arithmetic may not be accurate enough to handle the model’s numerous parameters and intricate relationships, which could result in errors that may happen during computation. Additionally, the Variable Precision Arithmetic (VPA) used for the sensitivity analysis is a global sensitivity analysis method, as it considers the entire input space, making it more robust for complex, multi-dimensional models like in this study.

Table 4. Elasticity indices of $R_0=0.98951$ to the parameters of the model.

Parameter	Value	Elasticity Index
${\mathsf{\Pi }}_H$	19.5	−0.58423
${\mathsf{\Pi }}_A$	8.33	−0.11367
${\mathsf{\Pi }}_F$	210.62	0.6979
${\mu }_H$	0.00795	0.58423
${\mu }_A$	0.0169	0.03997
${\mu }_F$	0.14	−1.3958
${\varpi }_1$	0.79	0.029476
${\varpi }_2$	0.21	−0.69521
${\zeta }_1$	0.69	−1.3472
${\zeta }_2$	0.211	0.024342
${\nu }_2$	0.03	0.024342
$\kappa$	0.115	−0.023604
$\pi$	0.2	−0.39048
${\beta }_H$	1	0
${\beta }_A$	1	0.0051338
${\beta }_F$	1	0.6979
$g_H$	0.02856	0
$g_A$	0.02856	0.70304
$h_A$	1.39	0.029476

provides the numerical values indicating $R_0$’s relative significance. It is discovered that some parameters have a negative elasticity index, while others are positive, and a few zero. When the parameters show a positive relationship, it means that raising the values of that parameter will significantly impact how frequently the disease spreads. Conversely, if there is a negative relationship, raising these parameters would aid in slowing the disease’s spread.

The numerical values from Table 4, provides a comprehensive overview of which parameters are critical for controlling the disease; that is the impact of parameter changes on $R_0$. Understanding the elasticity indices helps in strategic planning for disease control, allowing public health officials to focus on the most influential parameters. This targeted approach can optimize resource allocation and improve the effectiveness of intervention strategies, ultimately contributing to better disease management and control efforts.

To gain more insight, the contribution of some of the parameters in $R_0$, some experiments were conducted by varying the values of the key parameters (using the basic reproduction number $R_0$ as response function). The transmission probability from infected sandflies to animals (${\beta }_A$), transmission probability from infected animals to sandflies (${\beta }_F$), per capita biting rate of sandflies of animals ($g_A$), and rate of transfer from symptomatic infected animals to the recovered class ($\kappa$) are among the parameters that have greatest influence on $R_0$. In particular, the range of $R_0\in (0.98951,1.1138)$ as the values of the critical parameters (${\beta }_A,{\beta }_F,g_A$ and $\kappa$) are varied from $(0.8,0.8,0.01856,0.115)$ to $(1.3,1.3,0.02856,0.415)$, respectively. These numerical results corroborate the analytical results obtained for the critical parameters. It is worth mentioning that attention should be given to the most sensitive parameters identified so that the value of $R_0$ remains less than one, to avoid disease outbreak.

Thus, this study shows that effective disease control entails a multi-faceted approach based on minimizing the transmission probability from infected sandflies to animals, transmission probability from infected animals to sandflies, and per capita biting rate of sandflies of animals; increasing the rate of transfer from symptomatic infected animals to the recovered class; and the early diagnosis of ZVL cases in animals (reservoirs), among others.

The simulations also show that ZVL modeling studies in communities may require the implementation of an insecticide-based treatment strategy of infected animals (reservoirs) to reduce disease transmission and burden, and increase sandfly mortality, which is the main source of ZVL transmission.

The profiles of the state variables are given in Figure 2a–t:

6. Conclusions

In this paper, a mathematical model of the transmission dynamic of ZVL incorporating several additional complexities was considered. These complexities include lines of treatment (first and second lines of treatment) to reflect the healthcare interventions for visceral leishmaniasis, depending on the progression of the disease and the response to initial treatments, and asymptomatic and symptomatic stages in animals to differentiate between asymptomatic and symptomatic stages in animal hosts, which is crucial because animals, particularly canines, play a significant role in the transmission cycle of zoonotic visceral leishmaniasis. Asymptomatic animals can still be infectious, which complicates control efforts.

Providing precise instructions for disease control, the stability analysis demonstrates that the stability of the VL-endemic and VL-free fixed points globally will make the development of strategies for the eradication of the disease easy. Public health actions can be made simpler by emphasizing the necessity of maintaining $R_0$ below one in order to eradicate the disease.

Furthermore, the research’s findings are consistent with the previous findings in the literature, particularly regarding the significance of an integrated strategy that addresses several elements of transmission, prioritizes vector (sandfly) control, and regulates reservoirs (animals). The research’s other distinctive features include an emphasis on the targeted increase in sandfly mortality as determined by the analysis, the recommendation for insecticide-based treatment specifically for infected animals, and a clear understanding of the various stages of infection in both humans and animals. These elements reflect new approaches to resource optimization, even though they are in line with more general ZVL management techniques.

The model can be used for forecasting the disease outcomes and directing management methods because of its analytical and numerical consistency, which increases its credibility. This is further improved by the sensitivity analysis, which identifies crucial variables that have a major impact on ZVL dynamics and helps to come up with more successful intervention plans.