Enhancing reservoir control in the co-dynamics of HIV-VL: from mathematical modeling perspective

HIV patients are vulnerable to developing active visceral leishmaniasis (VL). To understand this complication, we studied a mathematical model for HIV and visceral leishmaniasis coinfection. In this approach, we reckoned two distinct equilibria: the disease-free and the endemic equilibria. The local and global stability of the disease-free equilibrium were thoroughly investigated. To further support the qualitative findings, we performed simulations to quantify the changes of the dynamical behavior of the full model for variation of relevant parameters. Increasing the rate of VL recovery (φ1), the recovery rate for VL–HIV Co-infection (φ2), removing reservoirs (c1), minimizing the contact rate (βh) are important in controlling the transmission of individual and co-infection disease of VL and HIV. In conclusion, possible measures should be implemented to reduce the number of infected individuals. Therefore, we recommend that policy makers and stakeholders incorporate these measures during planing and implementation phases to control the transmission of VL–HIV co-infection.


Introduction
Visceral leishmaniasis (VL) also known as 'Kala-azar' is a vector borne, zoonotic disease caused by Leishmania donovani species [1,2]. There are more than 20 species of leishmania that can cause human infection. The infection is transmitted following a successful bite and inoculation by the infected phlebotomine female sand flies [3,4]. The World Health Organization(WHO) considers leishmaniasis as the sixth most important endemic disease in the world [5], and it is distributed around the world in 90 countries [6], most of which are developing countries associated with malnutrition, population displacement, poor housing, a weak immune system, and lack of financial resources [6]. WHO estimated that from about 900,000 to 1.3 million new cases of leishmaniasis are reported per year [1,2,6,7]; of these, approximately 0.2-0.4 million are of visceral leishmaniasis (VL) [6,7]. The spread of the disease is linked to environmental changes such as deforestation, building of dams, irrigation schemes, and urbanization [7].
Human immunodeficiency virus (HIV) is the etiological agent responsible for the acquired immunodeficiency syndrome (AIDS) [8]. There are multiple modes of HIV transmission including sexual intercourse, sharing needles with HIV-infected persons, or via HIV-contaminated blood transfusions and others [8]. It is estimated that 36.7 million people worldwide are living with HIV [7] and 2.0 million new infections are reported per year [7,9]. Also approximately one million died of HIV-related causes globally [10].
It is expected that there is an overlap between the transmission areas of HIV and leishmaniasis. Due to this fact there have been an increasing number of cases of VL-HIV coinfection, which has spread throughout the world [4,7]. The co-infection has been reported in 35 countries [11,12]. So that VL-HIV co-infection is an emerging new threat to global public health and development [4,11,12].
Mathematical modeling has a great role in describing the dynamics of infectious diseases in a community [13,14]. Several scholars have developed different models for HIV, VL, and their co-infection with other diseases to study their transmission dynamics. For example, many mathematical models have been developed to understand the transmission nature of HIV [9,[15][16][17][18] and VL [19][20][21][22][23][24][25] explicitly. In addition, some co-infection models for HIV and malaria [26], HIV-TB [8] and VL and malaria [5] were proposed and analyzed. Hussaini et al. [27] recently presented a mathematical model to study the transmission dynamics of HIV and VL co-infection. In developing their model, they took into account the human and sand fly populations. However, they did not consider the reservoir population for VL and co-infection transmission dynamics in their assumption. Reservoirs are important for maintaining the life cycle of many leishmania species and hence are important for transmission of zoonotic and rural/sylvatic infections. Reservoirs acquire infection with leishmaniasis following contact with infected sand flies, and they are infected to their life time. Treating humans while removing reservoirs from the system reduces the fraction of infected sand flies, which gives a good control strategy against the disease [3,21].
In this study, we formulated and analyzed mathematical model for VL-HIV coinfection, which incorporates the key epidemiological and biological features of each of the two diseases by considering the reservoir population for VL transmission.
The paper is organized as follows: In Sect. 2, the mathematical model of HIV-VL is presented together with a set of definitions and basic underlying hypotheses. The formulated HIV-VL mathematical model is analyzed in Sect. 3. In Sect. 4 numerical simulations to give a better interpretation of the analytical results were reported. The last part, Sect. 5, is devoted to the discussion and conclusions.

Baseline model formulation
In this work, the population that has been considered is classified into three subpopulations, namely the human population, the sand fly population, and the reservoir population; each of them are again divided into different states. The human population N h (t) that are sexually active in a certain community are sub-divided into susceptible individuals S h (t), individuals infected with visceral leishmaniasis (VL) only I hl (t), HIV only infected individuals without clinical symptoms of AIDS I hh (t), individuals co-infected with VL and HIV with no clinical symptoms of AIDS I hhl (t), HIV infected individuals with AIDS symptoms A h (t), and co-infected individuals both with clinical symptoms of AIDS and VL A hl (t).
The number of susceptible individuals can be increased by a constant rate of recruitment h and recovery from only VL infected class with a rate φ 1 and diminished by natural death A susceptible individual can acquire VL infection following the effective contact with VL infected single sand fly at a rate of λ l , and he/she also acquires HIV at a rate of λ h if the individual is in effective contact with an HIV infected individual. Once the susceptible individuals are infected with either of infections, they will transfer to the respective infection classes (either I hh (t) or I hl (t)). In this co-infection model the natural death rate μ h is assumed equal for all human population classes. For a more detailed description of the parameters, refer to Table 1. Thus, the total human population is given by The second sub-population of the model, the sand fly population N s (t), has two sub-classes denoted by S s (t) and I s (t) at time t representing susceptible and infected sand flies, respectively. The susceptible sand flies can be recruited at a rate of s through birth to the sand fly population. The natural death rate of the sand fly is denoted by μ s , it can contribute to the reduction of the sand fly susceptible class. This sums the total sand fly population to The last sub-population considered in this disease dynamics study is the reservoir population N r (t) at time t, which is also categorized into two groups; the susceptible class S r (t) and the infected reservoir I r (t). Thus, the total reservoir population is given by The forces of infection associated with HIV/AIDS, VL, sand flies and reservoir population are denoted by λ h , λ l , λ s , and λ r , respectively, and are given as follows: Considering the formulations and assumptions above with the compartmental diagram in Fig. 1, the HIV/AIDS-VL co-infection model can be described by the following system of ordinary differential equations: The rates of change for the total human population, the sand fly population, and the reservoir population with change in time are as follows:

Invariant region
The system in (1) describes the epidemiological dynamics of the human population, the sand fly population, and the reservoir population. If disease specific induced death rates are assumed negligible ( Thus, the feasible region is All the model parameters and variables are nonnegative. If nonnegative initial values of the state variables are taken from this region, then the solution of system (1) will remain in . Therefore, this region is a positively invariant set, and hence studying the dynamics of model (1) in the biologically meaningful and mathematically well-posed region is sufficient.

Analysis of the HIV/AIDS-VL co-infection model
Model analysis is an important part in modeling the epidemiological phenomenon to find the qualitative and theoretical results. Before proceeding to further analysis, nondimensionalizing the sub-populations as seen below helps to study the proportion of each of the state variables.
Here, it can be defined that p = N s N h and q = N s N r are the female sand fly-human ratio and female sand fly-reservoir ratio [5], respectively, and they can be assumed as constants [28]. Then the new deterministic model for the above re-scaled classes is as follows: with the biologically feasible region of Region is a positively invariant domain, and it is mathematically well-posed under system (2). Now it is possible to proceed to the detailed analysis of epidemiological model (2).

Disease-free equilibrium and basic reproduction number
Studying the possibility of the population to be invaded has paramount importance. Here the disease-free state that shows the absence of the diseases is investigated. The DFE for system (2) is given by Using the next generation approach (see [29,30]), F and V stand for new infections and transition of infection rates, respectively. The partial derivatives of F and V with respect to each non-susceptible class are denoted by F and V , respectively, where Then the basic reproduction number associated with system (2) is defined by the spectral radius of the next generation matrix FV -1 and obtained as Furthermore, it can also be given as is the basic reproduction number for an HIV only model, i.e., at the absence of VL. And R 0h = a l c l (μ r a l b l p+a l b l q(φ 1 +δ l +μ h )) μ s μ r (φ 1 +δ l +μ h ) represents the reproduction number for a VL only model at the absence of HIV. Proof To prove Theorem 3.1, it is apparent to start with linearizing model (2) around E 0 . Thus, the Jacobian matrix is given as follows: Hence the eigenvalues of the matrix J E 0 are λ 1 = -μ h , λ 4 = -(φ 2 + ε 1 κ + τ δ l + μ h ), λ 6 = -(φ 3 + δ hl + μ h ), λ 7 = -μ s , λ 9 = -μ r , and the rest can be obtained from the following characteristic polynomials: where +μ h ), B 2 = μ r +μ s +μ h +φ 2 +δ l , B 1 = (μ h +φ 2 +δ l )(μ r +μ s )+μ r μ s -a l c l a l b l q-a l c l pa l b l , B 0 = (μ h +φ 2 +δ l )(μ r μ s -a l c l a l b l q)a l b 1 pa l c l μ r .
For a quadratic polynomial to have negative roots, the Routh-Hurwitz stability criterion states that both A 1 and A 0 must be greater than zero. Thus, from the quadratic equation of system (3) if R 0 < 1, A 1 and A 0 yield And for cubic polynomial the criterion is B 2 > 0, B 0 > 0 and B 2 B 1 -B 0 > 0. Thus, from the cubic polynomial of equation (3), B 2 is apparently positive; as can be seen it is the sum of positive parameters. It is given for B 0 > 0 that (μ h + φ 2 + δ l )(μ r μ sa l c l a l b l q)a l b 1 pa l c l μ r > 0, which simplified to R 0 < 1. And for B 2 B 1 -B 0 > 0, (μ r + μ s + μ h + φ 2 + δ l ) (μ h + φ 2 + δ l )(μ r + μ s ) + μ r μ sa l c l a l b l qa l c l pa l b l -(μ h + φ 2 + δ l )(μ r μ sa l c l a l b l q) + a l b 1 pa l c l μ r > 0 also holds. Therefore, whenever R 0 < 1, the DFE becomes locally asymptotically stable.

Global stability of DFE
Applying the Castillo-Chavez et al. [31] (Appendix A), let us rewrite the model in system (2) as The column vector X contains uninfected classes, while the components of vector Z are the infected individuals. Thus, for the transformed system, model (2) is From the first condition, dX dt = F(X, 0), X * is globally asymptotically stable, the reduced system is given by Therefore, for the DFE of system (2), say E 0 = (X * , 0), to be globally stable, the following two conditions must be satisfied [31]. From the second condition, G(X, Z) = AZ -Ĝ(X, Z) withĜ(X, Z) ≥ 0, where A = D z G(X * , 0) = ( ∂G i ∂z j (X * , 0)) is the linearization matrix of system (2) evaluated at E 0 . Hence, the matrix A is where only v 22 is changed to β hκμ h , the rest are as defined previously. Now the column vectorĜ(X, Z) is given bŷ This implied that the third and the fifth entries ofĜ(X, Z) are negative if all the parameters and the state variables in those entries are taken strictly positive. HenceĜ(X, Z) is not greater than zero. Thus, the second condition (H 2 ) stated above is not fulfilled. This concludes that the DFE may not be globally asymptotically stable. Apart from the following special cases, this may be a guarantee for bifurcation analysis. For two special cases, we have the following sub-results.
be the bifurcation parameter solved at R 0l < R 0h = 1 (R 0 = 1). Then the Jacobian of system (5) at DFE is given by where the shorthand representations are as given previously, but β h is given by the bifurcation parameter β * . There are two eigenvectors associated with J β * at its simple (zero) eigenvalue. Thus, the right eigenvector denoted by (w = (w 1 , w 2 , w 3 , w 4 , w 5 , w 6 , w 7 , w 8 , w 9 , w 10 ) T ) has the entries as given below.
And the left eigenvector represented by , v 10 ) T ) is given with the following values of entries: The center manifold approach described in [31] is used to study the direction of bifurcation analysis with computation of a and b, where The first, the seventh, and the ninth entries of the left eigenvector are zero, thus the partial derivatives of the corresponding drift functions (f 1 , f 7 , and f 9 ) are not important in computing a and b than wasting time. Thus, from system (5), the only crucial ones and nonzero derivative of the drift evaluated at (E 0 , β * ) are given below.
= a l b l p, = -ε 2 a l b l p, ∂ 2 f 6 ∂x 5 ∂x 8 = ∂ 2 f 6 ∂x 8 ∂x 5 = ε 2 a l b l p, Additionally, Thus, using equations (6), (7), (8), and (9), it can be found that and where It is clear from equation 11 that the bifurcation coefficient, b, is automatically positive. Thus, it follows from Theorem 4.1 in [31] that the the HIV-VL co-infection model will undergo backward bifurcation if the bifurcation coefficient, a, given by equation 10, is positive.

Sensitivity analysis
In this section, sensitivity indices of R 0 with respect to the parameters are calculated, as shown in Table 2, using the formula, where y is the model parameter, following the technique described in [36,37]. These indices show how important each parameter is to the disease transmission. Since R 0 = max{R 0l , R 0h }, we obtained the sensitivity indices of R 0l and R 0h separately (see Appendix B).    Fig. 2 and Table 2, we have the following statements. If the sensitivity index is positive, then the reproduction number increases along with increasing value of the parameter. That means that positive index parameters have a power of expanding the disease if their value increases. On the other hand, if the sensitivity index is negative, then reproduction decreases with the increasing value of the parameter. This also mean that negative index parameters have a power of reducing the burden of the disease in the community as their value increases. From this policy makers and stakeholders are expected to act accordingly in combating VL infection, HIV infection, and their co-infection from the community.

Numerical simulation
In this section, we use numerical simulations to support the analytical results previously established. We used Maple 18 to show the effect of some parameters in the expansion as well as for the control of HIV only, VL only, and co-infection of HIV and VL. We used parameter values in Table 3 for simulation purpose.

Effect of human recovery rate (φ 1 ) on VL infectious population
In Fig. 3, we investigated the effect of φ 1 in reducing VL-only infectious population by maintaining the other parameters constant. The figure reflects that when the values of φ 1 increase, the number of VL only infectious population is diminished. Therefore, we should  focus on improving recovery rates either by treating infected populations or by increasing the levels of individuals' immunity to VL disease in the population. Government policy makers should take this into account as a mitigation strategy.

Effect of effective contact rate(β h ) on HIV infectious and co-infectious population
In this section, we examine the influence of effective contact rate β h on HIV infectious and HIV-VL co-infectious population. Figure 4 reflects that as the value of contact rate of β h is increased, the HIV infectious and HIV-VL co-infectious population is increased, which leads to the increased expansion of co-infection of HIV and VL. Consequently, to control HIV infected and co-infection of HIV and VL, decreasing the contact rate of β h is significant. Therefore, stakeholders need to focus on reducing the contact rate of HIV infectious by using an appropriate method of prevention mechanism to bring down the expansion of co-infection in the community. ulation by VL also grows up. Therefore, it is advisable to use treated bed net, chemical techniques to reduce the expansion. This is also another good control strategy to reduce VL infection.

Effect of recovery rate of HIV-VL (φ 2 ) on the co-infectious population
In this section, we investigate the effect of recovery rate of HIV and VL (φ 2 ) on the coinfectious population. As we clarified in the model description, due to treatment or other mechanisms, the co-infectious population recovers from VL only diseases with their own probability and join the recovered compartment. Therefore, Fig. 6 shows that increasing the rate of recovery of the co-infectious population has a great advantage in reducing VL diseases in the population. Figure 7 shows the effect of VL progression rate on the sand fly reservoir population. It illustrates that as the reservoir increases around the sand fly population, the infection of these reservoirs also increases. From this we can understand that it is important to remove

Results and conclusions
We developed a transmission dynamics model for VL-HIV co-infection, and the population is subdivided into ten compartments. Before starting the qualitative analysis of the model, we proved the existence of a region where the model is mathematically and epidemiologically well posed. Basic reproduction numbers, disease-free equilibrium, endemic equilibrium, and stability analysis of equilibrium points were analyzed. Numerically, we experimented on the effect of basic parameters in the expansion or control of VL only, HIV only, and their co-infection. From the result, we conclude that an increase in the rate of VL recovery (φ 1 ) contributes greatly to reducing VL infectious individuals in the community. Similarly, increasing the recovery rate for VL-HIV co-infection (φ 2 ) contributes to the reduction of co-infection in the population. Also it is important to remove reservoirs (c 1 ) from the system so that the number of infected reservoirs and infected sand flies is reduced. Reducing the contact rate β h is significant in bringing down the number of HIV and VL-HIV co-infected infectious population. The rate of recovery for co-infection also has an influence on dropping co-infectious population if its value has been improved. Therefore, stakeholders should focus on the above basic parameters by using an appropriate method of prevention mechanism to reduce the expansion of infection in VL only, HIV only, and VL-HIV co-infection in the community.