A Novel Multi-scale Immuno-epidemiological Model of Visceral Leishmaniasis in Dogs

Leishmaniasis is a neglected and emerging disease prevalent in Mediterranean and tropical climates. As such, the study and development of new models are of increasing importance. We introduce a new immuno-epidemiological model of visceral leishmaniasis in dogs. The within-host system is based on previously collected and published data, showing the movement and proliferation of the parasite in the skin and the bone-marrow, as well as the IgG response. The betweenhost system structures the infected individuals in time-sinceinfection and is of vector-host type. The within-host system has a parasite-free equilibrium and at least one endemic equilibrium, consistent with the fact that infected dogs do not recover without treatment. We compute the basic reproduction number R0 of the immunoepidemiological model and provide the existence and stability results of the population-level disease-free equilibrium. Additionally, we prove existence of an unique endemic equilibrium when R0 > 1, and evidence of backward bifurcation and existence of multiple endemic equilibria when R0 < 1.


I. INTRODUCTION
The leishmaniases are a group of diseases found in over 90 countries around the world, spread by over 30 species of the phlebotomine sand flies and infecting a variety of hosts including humans and dogs.While cutaneous leishmaniasis is more common, visceral leishmaniasis (VL) is lethal if untreated.We focus on zoonotic visceral leishmaniasis (ZVL), which has symptoms including enlarged spleen and liver and non-specific symptoms such as fever, weight loss, and anemia [1].The non-specificity makes diagnosis challenging, particularly in the case of dogs [15].The leishmaniases are classified as a Neglected Tropical Disease (NTD), with an estimated 0.2-0.4 million new human cases per year [14] and hundreds of millions at risk of new infection [19].The WHO has stated that it is one of the most significant tropical diseases in the world [19].As such, it is imperative to continue the study of VL, including its epidemiology, immunology, control measures, and identification.
Visceral leishmaniasis is usually caused by the L. donovani and L. infantum protozoa.The L. donovani-induced VL is more common in Africa and Asia, while the L. infantum-induced VL is more common in the Americas and the Mediterranean [10].While dogs and other mammals are infected by the L. donovani-induced VL, it has been shown that dogs are a primary reservoir only for the disease caused by the L. infantum species [2].In fact, previous work has postulated that dogs are the main contributor to the spread of VL, claiming that 20% of the infected dogs cause 80% of transmission [2].For these reasons, we focus primarily on the role of dogs in the persistence of VL.
To date, there have been few mathematical models for the between-host dynamics of VL [15].Some models have chosen to compartmentalize asymptomatic or latently infected, resistant, and/or recovered classes, but almost all models have been ODEs.Some models, including the first model of VL in dogs by Dye includes resistance from birth [3].Another model was developed by Ribas et al [4] which studies the population level dynamics between dogs and humans.The model was used to argue that treatment in dogs does not reduce significantly human illness.Shimozako et al considered a model of leishmaniasis in dogs and humans [16] and concluded that latent dogs contribute more to the illness than clinically ill dogs.Seva et al [12] investigated an outbreak of VL in Spain through a model that involved multiple host classesrabbits, hares, dogs, and humans.All of the VL epidemic models studied to date are single-scale ODE models.
Even fewer models have been made for the within-host dynamics.Länger et al [8] developed a model examining the effect of IgG1, IgG2a, and lymphocytes on the parasite load, concluding that their model could be used in identifying biologically significant parameters [8].Siewe et al [17] modeled macrophages, parasite loads, dendritic cells, T cells, and cytokines, and simulated the effects of various control measures.As a result, they stated that an increase in IFN-γ production should lead to a decrease in parasite load; an implication for potential therapy.
We develop, for the first time, an immuno-epidemiological model for VL.While still in the early stages of development, we intend for this model to display the effect that the within-host dynamics may have on infection of the vector.Additionally, we hope to assess the infectivity of a vector based on parasite reproduction and, eventually, control measure efficacy.Similarly, we plan to examine the parasite reproduction inside hosts with respect to the efficacy towards a vector, immune response, and treatment.Since these processes occur at a drastically different rate, we adopt the multi-scale approach.
Our within-host system was designed to fit data from Courtenay et al [2].Our between-host system is of vector-host type, with the infected host class structured by time-since-infection.The betweenhost system contains susceptible, infected, and recovered host classes and susceptible, carrier, and infectious vector classes.Courtenay et al [2] conclude that, while samples were taken from both the skin tissue and bone marrow to record parasite loads, it is the parasite load in the skin tissue that is the most reliable indicator of VL infection [2].We use this fact in the linking of the within-and between-host systems.
Following the introduction of the model in Section 2, we present a parasite-free equilibrium of the within-host system and prove it to be unstable in Section 3.1.Then, in Section 3.2, we introduce the basic reproduction number R 0 , and prove the disease-free equilibrium of the immunoepidemiological model to be locally asymptotically stable when R 0 < 1.We show the existence of an endemic equilibrium when R 0 > 1, and the existence of multiple endemic equilibria and the presence of backward bifurcation, when R 0 < 1.
In Section 4 we discuss our conclusions.

A. The Within-Host System
Our within-host model is motivated by time series data in Courtenay et al [2] pertaining to the parasite loads in the skin tissue and bone marrow of dogs, as well as the immunoglobulin G (IgG) concentration.We derive a system coupling the A full list of parameter meanings can be found in Table I and the variable meanings in Table IV: For the parasite loads, we assume reproduction is limited by a carrying capacities K S and K O .Hence the equations for P S and P O use logistic terms to model recruitment, with rates r S and r O .
We note that to infect a host, a sand fly must deposit parasites into the skin, when it takes a blood meal.Similarly, for a sand fly to become infected, it must take a blood meal from an infected host's skin.As there are parasites in the bone marrow, we can assume mobility of the parasite between the skin tissue and bone marrow.Hence our model contains "travel terms," with rates k O and k S and density conversion coefficient ρ.Thus, for every time within-host unit, a fraction of the parasites in the skin move to the bone marrow, and vice versa.
While not much is known about the parasite, we assume that the life span is short, and we include the natural death as part of the logistic term.However, we include the clearance of the parasite induced by the immune response as a separate term.Hence ε S and ε O are the IgG induced clearance rates of the parasite.
Lastly, we assume that the IgG response is caused by the presence of the foreign parasite, i.e., the basal level of IgG present before the introduction of the parasite is not counted towards G.As such, a S and a O are the IgG production rates caused by the parasite loads.We assume that the clearance rate is the only way IgG leave the system, which occurs at rate d.[2] obtained data for the parasite loads (PL) and IgG concentration.We present simulations of our within-host system.The behaviors  of these curves are similar to that of the data provided.The result of our simulation is given by the curves in Figures 2(A)-2(C).The parameter values for the simulation can be found in Table II.

B. The Between-Host System
Few models exist for the between-host dynamics of VL, and even fewer for zoonotic VL in dogs.Previous models have largely been ODE in structure.Some models included an asymptomatic or latently infected class [15].However, obtaning data on infectious dogs is obstructed by the fact that identifying infectious dogs can be very difficult [13].Since an infectious host is considered only infectious to the vector, the most reliable way to test for VL is through xenodiagnosis, a process in which a susceptible vector population bites a possibly infected host, and is then tested for the presence of the parasite [2].
However, xenodiagnosis is not always feasible or practical [2].Since the symptoms of VL are non-specific, it is difficult to separate the latently infected dogs from the infectious dogs.In our multi-scale model, we structure the infectious hosts by time-since-infection, with the assumption that hosts are less infectious and less likely to display symptoms closer to when they first contract the parasite.The time-since-infection structure provides flexibility; allowing for fitting data given in different time units in the two different scales -the within-host and the between-host.We first introduce the system for the between-host dynamics of VL.Definitions of the parameters used can be found in Table III and definitions of the variables used in Table IV: where x is the time-since-infection and the total number of infected hosts, I H (t), is The host system consists of susceptible S H (t), infected i(t, x), and recovered/resistant R H (t) classes.The constant k u in ( 5) and ( 6) accounts for the difference in the rates that the time t and timesince-infection x occur, i.e., x = k u t.For the sake of initial analysis, we let k u = 1.Susceptible hosts are born at rate Λ H , and move to the infected class with standard incidence β H aS H I V /N .Recovery takes place at rate σ.The integral term in (7) is the total number of recovered individuals per unit time.Hosts exit the system either through natural death, m H , or disease-induced mortality, µ.The existence of relapse and reinfection shows the necessity of waning immunity at rate γ (see Table III).
The vector system consists of susceptible vectors S V (t), carrier vectors C V (t), who are infected      While there are many unknowns about the sand flies and leishmania, it is known that the parasite must make its way through the sand fly after a blood meal before it can be deposited in a host.The time elapsed for the parasite to potentially infect a new host, called extrinsic incubation period, is comparable to the life span of the sand fly.This requires the carrier class for the vector.
Vectors are born at rate Λ V , and exit the system only through natural death rate m V .Vectors move from the carrier class to the infectious class at rate τ .The integral terms in ( 9) and ( 10) represent the force of infection of humans to susceptible vectors.Since the rate of infection (ROI) is assumed to be dependent on x, the rate of recovery of hosts, σ, the disease-induced death rate of hosts, µ, and the rate of an infected host infecting a susceptible vector at the time of the blood meal, β V , are also dependent on x.

C. Linking the Within-and Between-Host Systems
While much analysis is left, we note the different time scales that will be utilized.That of the parasites and vectors will occur much faster than that of the hosts.We introduce the methods in which we initially plan to incorporate the faster time scale into the spread of the virus.
Since Courtenay et al [2] concluded that the parasite load in the dog skin tissue was the best indicator of infectiousness to the vector, we use P S (x) to determine the rate of transmission β V (x): where ξ is the rate of exponential decay, c V is the maximal transmission coefficient, and ν = −2/ ln (10).This relationship was derived by Li et al [9], using data for dengue.
Assuming treatment, the recovery rate also depends on the time-since-infection.We assume recovery occurs when the within-host parasite load becomes zero.Thus, we let where δ 0 is a small constant, and θ S , θ O , and κ are constants [18].Note that when P S and P O approach 0, σ becomes large due to δ 0 .

A. Analysis of the Within-Host System
The within-host system (1)-( 3) has a parasitefree equilibrium E 0 = (0, 0, 0).To determine its stability, we consider the Jacobian at the parasitefree equilibrium.We have the following result.
Theorem 1.The parasite-free equilibrium E 0 is always unstable.
Proof: Let kO = 1 ρ k O , kS = k S ρ, and âS = a S ρ.The Jacobian of the within-host system evaluated at E 0 is which has the eigenvalue λ 1 = −d.The remaining eigenvalues are eigenvalues of the submatrix and Tr(J 1 ) < 0. Suppose that det(J 1 ) > 0. Note that Similarly, we can get r S > k S .Hence Tr(J 1 ) > 0 when det(J 1 ) > 0. Therefore E 0 is unstable.This stability result heuristically makes sense as, without the introduction of treatment, infected hosts stay infected.
Theorem 2. The within-host system (1)-( 3) always has at least one parasite equilibrium E * .This equilibrium is unique if r S < k S .If r S > k S , the equilibrium is unique if Proof: To show existence, we set the withinhost system equal to zero and reduce the system to ) Solving (16) for P O , we obtain P O = P S f 1 (P S ), where Substituting P O in (17), we obtain the following equation for P S : Denote by PS the value of P S such that Φ O ( PS ) = 0. Further, denote by PS the value of P S such that f 1 ( PS ) = 0. We have .
We consider the following cases Case 1:] We have r S < k S or r O < k O .Then, f 1 (P S ) > 0 and Φ O (P S ) > 0 iff P S ∈ (0, PS ).
We have that f 1 (P S ) is an incrasing function of P S .Further G 1 (P S ) is a decreasing function of P S .As f 1 (P S ) > 0 on P S ∈ (0, PS ) the roots of F (P S ) = 0 are the same as the roots of G 1 (P S ) = 0. Since G 1 is decreasing, if a root exists, it must be unique.Since , as in this case.On the other hand ) is continuous on (0, PS ), then there is at least one solution of G 1 (P S ) = 0. Hence, there exists a unique P * S ∈ (0, PS ) such that F (P * S ) = 0 and P * O = P * S f 1 (P * S ) > 0. In this case, it is easy to see that is positive as well.
Case 2: We have r S > k S .Since f 1 (0) < 0, the solution, if it exists, lies in a different interval.Note f 1 (P S ) > 0 iff P S ∈ ( PS , PS ).However, in this case we don't know whether PS < PS or vice versa.So we have to consider two subcases.Case 2A: Assume inequality (15), that is, assume PS < PS .Then f 1 (P S ) is an increasing function of P S with It is easy to see that in this case we also have Further, G 1 (P S ) is also monotone and continuous as in Case 1.Hence, there exists a unique P * S ∈ ( PS , PS ) such that F (P * S ) = 0 and P * O = P * S f 1 (P * S ) > 0. In this case it is easy to see that G * is positive as well.We also note that F ( PS ) = kS > 0. Thus, PS is not a solution.Case 2B: Assume PS > PS .Then f 1 (P S ) is a decreasing function of P S .
It is easy to see that in this case we also have Since G 1 (P S ) is also continuous as in Case 1, there exists at least one solution in ( PS , PS ).However, in this case G 1 (P S ) may not be monotone and the equilibrium may not be unique.Hence, there exists at least one P * S ∈ ( PS , PS ) such that F (P * S ) = 0 and P * O = P * S f 1 (P * S ) > 0. In this case it is easy to see that G * is positive as well, and PS is not a solution to F (P S ) = 0.

B. Analysis of the Full Immuno-Epidemiological Model
The immuno-epidemiological model has a disease-free equilibrium We linearize model ( 4)- (11)  we obtain the characteristic equation F (λ) = 1, where λ ∈ C and Then the basic reproduction number is F (0), or where R V is the basic reproduction number of the vectors and R H is the basic reproduction number of the hosts [11].The parameter m denotes the ratio of the vector to hosts and is defined as It should be noted that R 0 is dependent on the within-host system.
Theorem 3. If R 0 < 1, then the disease-free equilibrium E 0 is locally asymptotically stable.If R 0 > 1, then E 0 is unstable.
To study existence of endemic equilibria, we set the time derivatives in the between-host system equal to zero and reduce the system to where Solving (21) for S H gives f SH (I V /N ).Substituting that into (23) yields We let X = β H aI V /N and p = 1−γΣ/(γ +m H ).
We redefine f SH and f N as functions of X, and obtain We have that Expanding and rearranging, we obtain where and m = S 0 V /N 0 .Theorem 4. When R 0 > 1, there exists a unique positive endemic equilibrium.
Proof: We have that a 0 > 0, since p − M > 0. If R 0 > 1, then c 0 is necessarily negative.Hence, the equation (24) has exactly one positive solution.It is not hard to see that in this case S H = f SH (X) > 0 and N = f N (X) > 0. This also implies that S V > 0. Hence, a unique positive equilibrium exists.
On the other hand, if R 0 < 1, (24) may have two positive soluitons or no positive solutions.Two positive solutions are obtained if equation (24) exhibits backward bifurcation.We find a necessary and sufficient condition for the existence of two equilibria: noting that M is the disease-induced mortality.Thus, backward bifurcation in this model can occur only if M > 0.
Theorem 5.If R 0 < 1 and b 0 is negative, then backward bifurcation occurs and two endemic equilibria exist.If R 0 < 1 and b 0 is positive, there are no endemic equilibria.
The existence of the two endemic equilibria is established by the backward bifurcation shown in Figure 3 where X is plotted on the y-axis and R 0 is plotted on the x axis.The parameter varied in the figure is the mortality rate of the vector m V .It may be important to note the decrease in the level of the upper equilibrium upon increasing m V .This could imply that a stronger control measure on the vector could greatly affect the level of persistence of the disease.

IV. DISCUSSION AND CONCLUSION
In this paper, we present a new immunoepidemiological model for zoonotic visceral leishmaniasis (ZVL) in dogs, in which the within-host model simulated infectiousness based on parasite load and the between-host model was structured by time-since-infection.The within-host system examined the parasite loads in the skin and bone marrow, as well as the IgG concentration.This system agrees well with data provided by Courtenay et al [2].The within-host model was shown to have an unstable parasite-free equilibrium, in which the parasite population dies out within the host.This equilibrium's stability is in the absence of treatment, consistent with the persistence of ZVL without treatment.
While the agreement of the model solutions with the data was satisfactory, in future work we will fit the model to the data and examine the biological significance of the values found for the parameters.Upon establishing the existence of equilibria and their stability, it is of great importance to examine the effect of control measures on the parasite population, as originally presented by Dye [3].This would include existing medicinal treatments, a hypothetical vaccine, and control measures directly affecting the vector.
The basic reproduction number for the immunoepidemiological model R 0 was introduced, and the disease-free equilibrium of the between-host system was shown to be locally asymptotically stable when R 0 < 1.We then derived a quadratic equation for the equilibria of the full system, based on a reduction of the system.This equation in X := β H aI V /N was used to establish the existence and characterize the endemic equilibria of the immuno-epidemiological model.When R 0 > 1, the model was shown to have a unique positive endemic equilibrium.However, when R 0 < 1 and the coefficient of the linear term of the quadratic equation was negative, the presence of disease-induced mortality allowed for backward bifurcation to occur, consistent with the results found in ODE cases [5].This provided justification for the existence of two endemic equilibria.The presence of subthreshold equilibria generally obstructs disease eradication.Control measures in this case should be directed towards (A) removing the cause of the backward bifurcation which in this case is the disease-induced mortality in dogs or (B) coupling sustain control measures that bring R 0 below one with temporary control measures, such as mosquito spraying, to put the disease on elimination path [7].

Fig. 2 :
Fig. 2: Figures (a) and (b) show simulations of the parasites in the skin and bone marrow, respectively, in log 10 values.Figure (c) shows simulations of log 10 IgG antibody units.

TABLE I :
Parameters for the within-host system.

TABLE II :
Parameter values used.

TABLE III :
Parameters for the between-host system.

TABLE IV :
Definitions of dependent variables.

TABLE V :
Parameters for linking.
but not infectious yet, and infectious vectors I V (t).