www.biomathforum.org/biomath/index.php/biomath ORIGINAL ARTICLE A novel multi-scale immuno-epidemiological model of visceral leishmaniasis in dogs J. Shane Welker, Maia Martcheva Department of Mathematics, University of Florida, FL 32611, USA spane@ufl.edu, maia@ufl.edu Received: 8 August 2018, accepted: 2 January 2019, published: 23 January 2019 Abstract—Leishmaniasis is a neglected and emerging disease prevalent in Mediterranean and tropical climates. As such, the study and develop- ment 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 pub- lished data, showing the movement and proliferation of the parasite in the skin and the bone-marrow, as well as the IgG response. The between-host system structures the infected individuals in time-since- infection 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 immuno-epidemiological model and provide the existence and stability results of the population-level disease-free equilibrium. Additionally, we prove ex- istence of an unique endemic equilibrium when R0 > 1, and evidence of backward bifurcation and existence of multiple endemic equilibria when R0 < 1. Keywords-leishmaniasis in dogs, backward bifur- cation, immuno-epidemiological model, stability, pa- rameter estimation, immune dynamics AMS SUBJECT CLASSIFICATION: 92D30 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 leish- maniasis (ZVL), which has symptoms including enlarged spleen and liver and non-specific symp- toms such as fever, weight loss, and anemia [1]. The non-specificity makes diagnosis challenging, particularly in the case of dogs [15]. The leish- maniases 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 Copyright: © 2019 Welker et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Citation: J. Shane Welker, Maia Martcheva, A novel multi-scale immuno-epidemiological model of visceral leishmaniasis in dogs, Biomath 8 (2019), 1901026, http://dx.doi.org/10.11145/j.biomath.2019.01.026 Page 1 of 12 http://www.biomathforum.org/biomath/index.php/biomath https://creativecommons.org/licenses/by/4.0/ https://creativecommons.org/licenses/by/4.0/ http://dx.doi.org/10.11145/j.biomath.2019.01.026 J. Shane Welker, Maia Martcheva, A novel multi-scale immuno-epidemiological model of visceral ... and Asia, while the L. infantum-induced VL is more common in the Americas and the Mediter- ranean [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 signif- icantly 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 classes – rabbits, 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 biologi- cally 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 vec- tor, 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 between- host 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 R0, and prove the disease-free equilibrium of the immuno- epidemiological model to be locally asymptoti- cally stable when R0 < 1. We show the existence of an endemic equilibrium when R0 > 1, and the existence of multiple endemic equilibria and the presence of backward bifurcation, when R0 < 1. In Section 4 we discuss our conclusions. II. THE MODEL 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 Biomath 8 (2019), 1901026, http://dx.doi.org/10.11145/j.biomath.2019.01.026 Page 2 of 12 http://dx.doi.org/10.11145/j.biomath.2019.01.026 J. Shane Welker, Maia Martcheva, A novel multi-scale immuno-epidemiological model of visceral ... Vectors SV (t) CV (t) IV (t) de ath de ath de ath bite τbirth Hosts SH(t) ∫ i(t,x)dx RH(t) birth de ath de ath de ath transmission recovery waning immun ity Parasite load in skin (PS(x)) Parasite load in bone marrow (PO(x)) IgG level (G(x)) Within-Host Dynamics Fig. 1: The flow-chart of the immuno-epidemiological model of VL in dogs dynamics of the parasite load in the skin tissue (PS(x)), the parasite load in the bone marrow (PO(x)), and the IgG concentration (G(x)), where x is the time-since-infection. The subscript S will indicate skin tissue, while the subscript O will indicate the bone marrow. We first introduce the model below and then we explain and motivate it. A full list of parameter meanings can be found in Table I and the variable meanings in Table IV: P ′S =rSPS ( 1 − PS KS ) + 1 ρ kOPO −kSPS −εSPSG (1) P ′O =rOPO ( 1 − PO KO ) −kOPO + kSρPS −εOPOG (2) G′ =aSρPS + aOPO −dG (3) For the parasite loads, we assume reproduction is limited by a carrying capacities KS and KO. Hence the equations for P ′S and P ′ O use logistic terms to model recruitment, with rates rS and rO. We note that to infect a host, a sand fly must de- posit 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 kO and kS 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, aS and aO 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. Courtenay et al [2] obtained data for the parasite loads (PL) and IgG concentration. We present sim- ulations of our within-host system. The behaviors Biomath 8 (2019), 1901026, http://dx.doi.org/10.11145/j.biomath.2019.01.026 Page 3 of 12 http://dx.doi.org/10.11145/j.biomath.2019.01.026 J. Shane Welker, Maia Martcheva, A novel multi-scale immuno-epidemiological model of visceral ... TABLE I: Parameters for the within-host system. Parameter Description Unit aO PO induced IgG production rate IgG AU / [parasites · time x] aS PS induced IgG production rate IgG AU / [parasites · time x] d Natural clearance rate of IgG 1 / [time x] εO IgG induced PO clearance 1 / [time x · IgG AU/mL] εS IgG induced PS clearance 1 / [time x · IgG AU/mL] KO Parasite carrying capacity in bone marrow parasites/mL KS Parasite carrying capacity in skin tissue parasites/g kO Parasite travel rate from bone marrow to skin tissue 1 / [time x] kS Parasite travel rate from skin tissue to bone marrow 1 / [time x] rO Parasite reproduction rate in bone marrow 1 / [time x] rS Parasite reproduction rate in skin tissue 1 / [time x] ρ Conversion Coefficient [g skin] / [mL bone marrow] 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 dynam- ics 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 dif- ficult [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 con- tract the parasite. The time-since-infection struc- ture 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: S′H = ΛH − βHaSHIV N −mHSH + γRH, (4) it + kuix = −(σ(x) + µ(x) + mH)i(t,x), (5) kui(t, 0) = βHaSHIV N , (6) R′H = ∫ ∞ 0 σ(x)i(t,x)dx− (γ + mH)RH, (7) where x is the time-since-infection and the total number of infected hosts, IH(t), is IH(t) = ∫ ∞ 0 i(t,x)dx. (8) The host system consists of susceptible SH(t), infected i(t,x), and recovered/resistant RH(t) classes. The constant ku in (5) and (6) accounts for the difference in the rates that the time t and time- since-infection x occur, i.e., x = kut. For the sake of initial analysis, we let ku = 1. Susceptible hosts are born at rate ΛH, and move to the infected class with standard incidence βHaSHIV /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, mH, 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 vec- tors SV (t), carrier vectors CV (t), who are infected Biomath 8 (2019), 1901026, http://dx.doi.org/10.11145/j.biomath.2019.01.026 Page 4 of 12 http://dx.doi.org/10.11145/j.biomath.2019.01.026 J. Shane Welker, Maia Martcheva, A novel multi-scale immuno-epidemiological model of visceral ... (a) (b) (c) Fig. 2: Figures (a) and (b) show simulations of the parasites in the skin and bone marrow, respectively, in log10 values. Figure (c) shows simulations of log10 IgG antibody units. TABLE II: Parameter values used. Parameter Value Unit aO 4.550060 × 10−1 IgG AU / [parasites · day] aS 1.739817 × 10−5 IgG AU / [parasites · day] d 4.136157 × 10−3 parasites/day εO 8.927579 × 10−7 1 / [day · IgG AU/mL] εS 6.968101 × 10−7 1 / [day · IgG AU/mL] KO 1.034155 × 106 parasites/mL KS 1.007303 × 108 parasites/g kO 5.723339 × 10−9 1 / day kS 5.228469 × 10−4 1 / day rO 2.614700 × 10−2 1 / day rS 2.965272 × 10−2 1 / day ρ 1 [g skin] / [mL bone marrow] TABLE III: Parameters for the between-host system. Parameter Description Units a Average biting rate bites / [time t · vectors] βH Rate of host becoming infected after bite 1 / [bites/hosts] βV Rate of vector becoming infected after bite 1 / [bites/vectors] γ Rate of becoming susceptible after recovery 1 / [time t] ku Time scaling constant [time x] / [time t] ΛH Birth rate of hosts hosts / [time t] ΛV Birth rate of vectors vectors / [time t] mH Natural death rate of hosts 1 / [time t] mV Natural death rate of vectors 1 / [time t] µ(x) Disease induced death rate of hosts 1 / [time t] σ(x) Rate of recovery of hosts 1 / [time t] τ Rate of moving from carrier to infectious 1 / [time t] Biomath 8 (2019), 1901026, http://dx.doi.org/10.11145/j.biomath.2019.01.026 Page 5 of 12 http://dx.doi.org/10.11145/j.biomath.2019.01.026 J. Shane Welker, Maia Martcheva, A novel multi-scale immuno-epidemiological model of visceral ... TABLE IV: Definitions of dependent variables. Variable Description Unit SH(t) Susceptible hosts at time t hosts IH(t) Infected hosts at time t hosts i(t,x) Density of hosts infected x time units ago at time t hosts / [time x] RH(t) Recovered hosts at time t hosts SV (t) Susceptible vectors at time t vectors CV (t) Carrier vectors at time t vectors N(t) Host population hosts IV (t) Infectious vectors at time t vectors PO(x) Parasite load of bone marrow at time x parasites/mL PS(x) Parasite load of skin tissue at time x parasites/g G(x) IgG concentration at time x IgG AU/mL TABLE V: Parameters for linking. Parameter Description Unit ξ Rate of exponential decay 1 / [parasites/g] cV Maximal transmission coefficient 1 / [bites/vector] δ0 Constant 1 / mL κ Constant 1 / [IgG AU/mL · time t] η Constant 1 / [IgG AU/mL · time t] ν Constant (unitless) ψO Constant 1 / [parasites/mL · time t] ψS Constant 1 / [parasites/g · time t] θO Constant 1 / parasites θS Constant 1 / parasites but not infectious yet, and infectious vectors IV (t). 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 in- fect 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. S′V = ΛV −mV SV − aSV N ∫ ∞ 0 βV (x)i(t,x)dx (9) C′V = aSV N ∫ ∞ 0 βV (x)i(t,x)dx−(τ +mV )CV (10) I′V = τCV −mV IV , (11) Vectors are born at rate ΛV , and exit the system only through natural death rate mV . 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 dif- ferent 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 PS(x) to determine the rate of transmission βV (x): βV (x) = cV e −ξPνS (x), (12) where ξ is the rate of exponential decay, cV is the maximal transmission coefficient, and ν = Biomath 8 (2019), 1901026, http://dx.doi.org/10.11145/j.biomath.2019.01.026 Page 6 of 12 http://dx.doi.org/10.11145/j.biomath.2019.01.026 J. Shane Welker, Maia Martcheva, A novel multi-scale immuno-epidemiological model of visceral ... −2/ ln(10). This relationship was derived by Li et al [9], using data for dengue. Assuming treatment, the recovery rate also de- pends on the time-since-infection. We assume re- covery occurs when the within-host parasite load becomes zero. Thus, we let σ(x) = κG δ0 + θSρPS + θOPO , (13) where δ0 is a small constant, and θS, θO, and κ are constants [18]. Note that when PS and PO approach 0, σ becomes large due to δ0. To link the disease-induced mortality, µ, we let µ(x) = ψSρPS + ψOPO + ηG, (14) where ψS, ψO, and η are constants, [6]. III. ANALYSIS A. Analysis of the Within-Host System The within-host system (1)-(3) has a parasite- free equilibrium E0 = (0, 0, 0). To determine its stability, we consider the Jacobian at the parasite- free equilibrium. We have the following result. Theorem 1. The parasite-free equilibrium E0 is always unstable. Proof: Let k̂O = 1 ρ kO, k̂S = kSρ, and âS = aSρ. The Jacobian of the within-host system evaluated at E0 is J0 :=  rS −kS k̂O 0k̂S rO −kO 0 âS aO −d   , which has the eigenvalue λ1 = −d. The remaining eigenvalues are eigenvalues of the submatrix J1 := ( rS −kS k̂O k̂S rO −kO ) . Note that k̂Ok̂S = kOkS. Then E0 is stable if and only if det(J1) = (rS−kS)(rO−kO)−kSkO > 0 and Tr(J1) < 0. Suppose that det(J1) > 0. Note that det(J1) = rSrO −rSkO −rOkS = rS(rO −kO) −rOkS. So rS(rO−kO)−rOkS > 0 if and only if rS(rO− kO) > rOkS. Then rO > kO. Similarly, we can get rS > kS. Hence Tr(J1) > 0 when det(J1) > 0. Therefore E0 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) al- ways has at least one parasite equilibrium E∗. This equilibrium is unique if rS < kS. If rS > kS, the equilibrium is unique if k̂O ( εS d âS + rS KS ) > (rS −kS) εS d aO (15) Proof: To show existence, we set the within- host system equal to zero and reduce the system to 0 =rSPS ( 1 − PS KS ) + k̂OPO −kSPS − εS d PS(âSPS + aOPO) (16) 0 =rOPO ( 1 − PO KO ) −kOPO + k̂SPS − εO d PO(âSPS + aOPO). (17) Solving (16) for PO, we obtain PO = PSf1(PS), where f1(PS) = 1 ΦO(PS) [ kS + εS d âSPS−rS ( 1− PS KS )] and ΦO(PS) = k̂O − εS d aOPS. Substituting PO in (17), we obtain the following equation for PS: F(PS) := f1(PS)G1(PS) = 0 with G1(PS) =rO ( 1 − PSf1(PS) KS ) −kO + k̂S f1(PS) − εO d (âSPS + aOPSf1(PS)). Biomath 8 (2019), 1901026, http://dx.doi.org/10.11145/j.biomath.2019.01.026 Page 7 of 12 http://dx.doi.org/10.11145/j.biomath.2019.01.026 J. Shane Welker, Maia Martcheva, A novel multi-scale immuno-epidemiological model of visceral ... Denote by P̂S the value of PS such that ΦO(P̂S) = 0. Further, denote by P̄S the value of PS such that f1(P̄S) = 0. We have P̄S = rS −kS εSâS d + rS KS , P̂S = k̂O εS d aO . We consider the following cases Case 1:] We have rS < kS or rO < kO. Then, f1(PS) > 0 and ΦO(PS) > 0 iff PS ∈ (0, P̂S). We have that f1(PS) is an incrasing function of PS. Further G1(PS) is a decreasing function of PS. As f1(PS) > 0 on PS ∈ (0, P̂S) the roots of F(PS) = 0 are the same as the roots of G1(PS) = 0. Since G1 is decreasing, if a root exists, it must be unique. Since G1(0) = rO−kO + kSkO kS −rS = rO + kOrS kS −rS > 0 if kS > rS, as in this case. On the other hand lim PS→P̂S − G1(PS) = −∞ Since G1(PS) is continuous on (0, P̂S), then there is at least one solution of G1(PS) = 0. Hence, there exists a unique P∗S ∈ (0, P̂S) such that F(P∗S) = 0 and P ∗ O = P ∗ Sf1(P ∗ S) > 0. In this case, it is easy to see that G∗ = âSP ∗ S + aOP ∗ O d is positive as well. Case 2: We have rS > kS. Since f1(0) < 0, the solution, if it exists, lies in a different interval. Note f1(PS) > 0 iff PS ∈ (P̄S, P̂S). However, in this case we don’t know whether P̄S < P̂S or vice versa. So we have to consider two subcases. Case 2A: Assume inequality (15), that is, as- sume P̄S < P̂S. Then f1(PS) is an increasing function of PS with f1(P̄S) = 0, lim PS→P̂S − f1(PS) = ∞. It is easy to see that in this case we also have lim PS→P̄S + G1(PS) =∞, lim PS→P̂S − G1(PS) = −∞. Further, G1(PS) is also monotone and con- tinuous as in Case 1. Hence, there exists a unique P∗S ∈ (P̄S, P̂S) such that F(P ∗ S) = 0 and P∗O = P ∗ Sf1(P ∗ S) > 0. In this case it is easy to see that G∗ is positive as well. We also note that F(P̄S) = k̂S > 0. Thus, P̄S is not a solution. Case 2B: Assume P̄S > P̂S. Then f1(PS) is a decreasing function of PS. f1(P̄S) = 0, lim PS→P̂S − f1(PS) = ∞ It is easy to see that in this case we also have lim PS→P̄S + G1(PS) =∞, lim PS→P̂S − G1(PS) = −∞. Since G1(PS) is also continuous as in Case 1, there exists at least one solution in (P̂S, P̄S). However, in this case G1(PS) may not be monotone and the equilibrium may not be unique. Hence, there exists at least one P∗S ∈ (P̂S, P̄S) such that F(P∗S) = 0 and P ∗ O = P∗Sf1(P ∗ S) > 0. In this case it is easy to see that G∗ is positive as well, and P̄S is not a solution to F(PS) = 0. B. Analysis of the Full Immuno-Epidemiological Model The immuno-epidemiological model has a disease-free equilibrium E0 = ( ΛH mH , 0, 0, ΛV mV , 0, 0 ) . (18) We linearize model (4)-(11) around E0. Looking for exponential solutions of the form ~z(t) = zeλt, Biomath 8 (2019), 1901026, http://dx.doi.org/10.11145/j.biomath.2019.01.026 Page 8 of 12 http://dx.doi.org/10.11145/j.biomath.2019.01.026 J. Shane Welker, Maia Martcheva, A novel multi-scale immuno-epidemiological model of visceral ... we obtain the characteristic equation F(λ) = 1, where λ ∈ C and F(λ) = τ mV + λ · am τ + mV + λ βHa·∫ ∞ 0 βV (x)e −λxe− ∫ x 0 (σ(ξ)+µ(ξ)+mH)dξdx. (19) Then the basic reproduction number is F(0), or R0 = % of vectors that become infectious︷ ︸︸ ︷ τ τ + mV · aβH mV︸ ︷︷ ︸ RV · RH︷ ︸︸ ︷ am ∫ ∞ 0 βV (x)e − ∫ x 0 (σ(ξ)+µ(ξ)+mH)dξdx, (20) where RV is the basic reproduction number of the vectors and RH is the basic reproduction number of the hosts [11]. The parameter m denotes the ratio of the vector to hosts and is defined as m = ΛV mV mH ΛH . It should be noted that R0 is dependent on the within-host system. Theorem 3. If R0 < 1, then the disease-free equilibrium E0 is locally asymptotically stable. If R0 > 1, then E0 is unstable. Proof: Suppose R0 < 1. Then F(λ) = 1 has a unique negative solution λ∗ ∈ R. For a complex λ, let λ = α1 + iα2, and assume α1 ≥ 0. Then |F(λ)| ≤ a2βHmτ |mV + λ| |τ + mV + λ| · ∫ ∞ 0 βV (x) ∣∣∣e−λx∣∣∣π(x)dx ≤ a2βHmτ (mV + α1)(τ + mV + α1) · ∫ ∞ 0 βV (x)e −α1xπ(x)dx =F(α1) ≤ F(0) = R0 < 1. Since 1 = |F(λ)| ≤ R0 < 1, a contradiction, we must have α1 < 0. Hence every complex solution to F(λ) = 1 must have a negative real part. Therefore, E0 is locally asymptotically stable when R0 < 1. Now suppose that R0 > 1. Then for positive λ ∈ R, F ′(λ) = (mV + λ)(τ + mV + λ)(a 2βHmτ) (mV + λ)2(τ + mV + λ)2 · ∫ ∞ 0 (−x)βV (x)e−λxπ(x)dx − a2βHmτ(τ + 2mV + 2λ) (mV + λ)2(τ + mV + λ)2 · ∫ ∞ 0 βV (x)e −λxπ(x)dx = − a2βHmτ (mV + λ)2(τ + mV + λ)2 · [ (mV + λ)(τ + mV + λ) · ∫ ∞ 0 xβV (x)e −λxπ(x)dx + (τ +2mV +2λ) ∫ ∞ 0 βV (x)e −λxπ(x)dx ] < 0, since the bracketed expression is always positive for non-negative βV (x) 6≡ 0. Since R0 > 1, then F(0) > 1. Since lim λ→∞ F(λ) = 0 and F is decreasing, F(λ) = 1 has a unique positive solution λ∗ ∈ R. Thus E0 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 0 = ΛH − βHaSHIV N −mHSH + γ γ + mH · βHaSHIV N Σ, (21) fSV (IV )SH N2 = Q, (22) N2 − ΛH mH N = − βHaSHIV mH M, (23) Biomath 8 (2019), 1901026, http://dx.doi.org/10.11145/j.biomath.2019.01.026 Page 9 of 12 http://dx.doi.org/10.11145/j.biomath.2019.01.026 J. Shane Welker, Maia Martcheva, A novel multi-scale immuno-epidemiological model of visceral ... where fSV (IV ) := SV = ΛV mV − τ + mV τ IV , Q = (τ + mV )mV τ · 1 βHa2BV , M = ∫ ∞ 0 µ(x)π(x)dx, Σ = ∫ ∞ 0 σ(x)π(x)dx, BV = ∫ ∞ 0 βV (x)π(x)dx, π(x) = exp ( − ∫ x 0 (σ(ξ) + µ(ξ) + mH)dξ ) , N = SH + ∫ ∞ 0 i(t,x)dx + RH = ΛH mH − βHaSHIV N M. Solving (21) for SH gives fSH (IV /N). Substitut- ing that into (23) yields N = ΛH mH − βHa [ fSH ( IV N )] IV N mH M =: fN (IV /N). We let X = βHaIV /N and p = 1−γΣ/(γ+mH). We redefine fSH and fN as functions of X, and obtain fSH (X) = ΛH mH + pX , fN (X) = ΛH mH [ 1 − MX mH + pX ] . We have that QN = SV SH/N = fSV (X)fSH (x)/fN (X). Expanding and rearranging, we obtain a0X 2 + b0X + c0 = 0, (24) where a0 = (p−M)2 mH m R0 + τ + mV τβHa (p−M), (25) b0 = 2m R0 (p−M) −mp + τ + mV τβHa mH, (26) c0 = ( 1 R0 − 1 ) mHm, (27) and m = S0V /N 0. Theorem 4. When R0 > 1, there exists a unique positive endemic equilibrium. Proof: We have that a0 > 0, since p−M > 0. If R0 > 1, then c0 is necessarily negative. Hence, the equation (24) has exactly one positive solution. It is not hard to see that in this case SH = fSH (X) > 0 and N = fN (X) > 0. This also implies that SV > 0. Hence, a unique positive equilibrium exists. On the other hand, if R0 < 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: p + aBV mH < 2M, noting that M is the disease-induced mortality. Thus, backward bifurcation in this model can occur only if M > 0. Theorem 5. If R0 < 1 and b0 is negative, then backward bifurcation occurs and two endemic equilibria exist. If R0 < 1 and b0 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 R0 is plotted on the x axis. The parameter varied in the figure is the mortality rate of the vector mV . It may be important to note the decrease in the level of the upper equilibrium upon increasing mV . 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 immuno- epidemiological model for zoonotic visceral leish- maniasis (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. Biomath 8 (2019), 1901026, http://dx.doi.org/10.11145/j.biomath.2019.01.026 Page 10 of 12 http://dx.doi.org/10.11145/j.biomath.2019.01.026 J. Shane Welker, Maia Martcheva, A novel multi-scale immuno-epidemiological model of visceral ... Fig. 3: Backward bifurcation of (24) when R0 < 1 and b0 < 0. Various values of mV are shown. 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 equi- librium’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 di- rectly affecting the vector. The basic reproduction number for the immuno- epidemiological model R0 was introduced, and the disease-free equilibrium of the between-host system was shown to be locally asymptotically stable when R0 < 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 := βHaIV /N was used to establish the ex- istence and characterize the endemic equilibria of the immuno-epidemiological model. When R0 > 1, the model was shown to have a unique positive endemic equilibrium. However, when R0 < 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 justifica- tion 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 R0 below one with temporary control measures, such as mosquito spraying, to put the disease on elimination path [7]. ACKNOWLEDGEMENTS The authors acknowledge partial support from NSF grant DMS-1515661, and the contributions of George Pu in the fitting of the within-host model. We also thank the referees for their comments and suggestions. REFERENCES [1] About Leishmaniasis - DNDi. 2016. URL: http : / / www . dndi . org / diseases - projects / leishmaniasis/. [2] Orin Courtenay et al. “Heterogeneities in Leishmania infantum infection: Using skin parasite burdens to identify highly infec- tious dogs”. In: PLoS Neglected Tropi- cal Diseases 8.1 (2014), pp. 1–9. ISSN: 19352735. DOI: 10 . 1371 / journal . pntd . 0002583. [3] Christopher Dye. “The epidemiology of canine visceral leishmaniasis in southern France: classical theory offers another ex- planation of the data”. In: Parasitology 96.1 (1988), pp. 19–24. DOI: 10 . 1017 / S0031182000081622. Biomath 8 (2019), 1901026, http://dx.doi.org/10.11145/j.biomath.2019.01.026 Page 11 of 12 http://www.dndi.org/diseases-projects/leishmaniasis/ http://www.dndi.org/diseases-projects/leishmaniasis/ http://dx.doi.org/10.1371/journal.pntd.0002583 http://dx.doi.org/10.1371/journal.pntd.0002583 http://dx.doi.org/10.1017/S0031182000081622 http://dx.doi.org/10.1017/S0031182000081622 http://dx.doi.org/10.11145/j.biomath.2019.01.026 J. Shane Welker, Maia Martcheva, A novel multi-scale immuno-epidemiological model of visceral ... [4] “Estimating the Optimal Control of Zoonotic Visceral Leishmaniasis by the Use of a Mathematical Model”. In: The Scientific World Journal 2013.810380 (2013). [5] SM Garba, AB Gumel, and MRA Bakar. “Backward bifurcations in dengue transmis- sion dynamics”. In: Math. Biosci. 215(1) (2008), pp. 11–25. [6] Michael A. Gilchrist and Akira Sasaki. “Modeling host-parasite coevolution: A nested approach based on mechanistic mod- els”. In: Journal of Theoretical Biology 218.3 (2002), pp. 289–308. ISSN: 0022- 5193. DOI: https : / / doi . org / 10 . 1006 / jtbi . 2002 . 3076. URL: http : / / www . sciencedirect . com / science / article / pii / S0022519302930766. [7] H. Gulbudak and M. Martcheva. “Forward hysteresis and backward bifurcation caused by culling in avian influenza model”. In: Math. Biosci. 246(1).9 (2013), pp. 202–212. [8] Bettina M Länger et al. “Modeling of leish- maniasis infection dynamics: novel applica- tion to the design of effective therapies.” In: BMC systems biology 6.1 (2012), p. 1. ISSN: 1752-0509. DOI: 10.1186/1752- 0509- 6- 1. URL: http://www.biomedcentral.com/1752- 0509/6/1. [9] Xue-Zhi Li, JunYuan Yang, and Maia Martcheva. “Age Structured Epidemic Mod- eling”. In: (expected 2019). [10] J. Lukes et al. “Evolutionary and geo- graphical history of the Leishmania dono- vani complex with a revision of current taxonomy”. In: Proceedings of the Na- tional Academy of Sciences 104.22 (2007), pp. 9375–9380. ISSN: 0027-8424. DOI: 10. 1073/pnas.0703678104. URL: http://www. pnas.org/cgi/doi/10.1073/pnas.0703678104. [11] Maia Martcheva. An Introduction to Math- ematical Epidemiology. 1st ed. Vol. 61. Springer US, 2015. ISBN: 978-1-4899- 7611-6. [12] A da Paixo Sev et al. “Efficacies of preven- tion and control measures applied during an outbreak in Southwest Madrid, Spain”. In: PloS one 12 (). [13] CA Petersen and SC Barr. “Canine leishma- niasis in North America: Emerging or newly recognized?” In: 76.October 2009 (2012), pp. 211–220. ISSN: 00092665. DOI: 10 . 1007 / s11103 - 011 - 9767 - z . Plastid. arXiv: NIHMS150003. [14] CDC Centers for Disease Control and Pre- vention. CDC - Leishmaniasis. URL: https: //www.cdc.gov/parasites/leishmaniasis/. [15] K.S. Rock et al. “Progress in the Math- ematical Modelling of Visceral Leishma- niasis”. In: Elsevier Ltd., 2016. ISBN: 9780128050606. [16] Helio Junji Shimozako, Jianhong Wu, and Eduardo Massad. “Mathematical modelling for Zoonotic Visceral Leishmaniasis dynamics: A new analysis considering updated parameters and notified human Brazilian data”. In: KeAi Infectious Disease Modelling 2.2 (2017), pp. 143–160. DOI: 10 . 1016 / j . idm . 2017 . 03 . 002. URL: http://ac.els-cdn.com/S2468042716300173/ 1 - s2 . 0 - S2468042716300173 - main . pdf ? {\ } tid = c8ed3798 - 8db2 - 11e7 - b1dc - 00000aacb360{\ & }acdnat = 1504118647{\ } 14d164bb702c07064ef2e80eeff2e974. [17] Nourridine Siewe et al. “Immune response to infection by Leishmania : A mathemat- ical model”. In: Math. Biosci. 276 (2016), pp. 28–43. ISSN: 00255564. DOI: 10.1016/ j . mbs . 2016 . 02 . 015. URL: http : / / linkinghub . elsevier . com / retrieve / pii / S0025556416000468. [18] N. Tuncer et al. “Structural and prac- tical identifiability issues of immuno- epidemiological vector-host models with application to Rift Valley Fever”. In: Bull. Math. Biol. 78.9 (2016), pp. 1796–1827. [19] WHO — Leishmaniasis. 2017. URL: http : //www.who.int/leishmaniasis/en/. Biomath 8 (2019), 1901026, http://dx.doi.org/10.11145/j.biomath.2019.01.026 Page 12 of 12 http://dx.doi.org/https://doi.org/10.1006/jtbi.2002.3076 http://dx.doi.org/https://doi.org/10.1006/jtbi.2002.3076 http://www.sciencedirect.com/science/article/pii/S0022519302930766 http://www.sciencedirect.com/science/article/pii/S0022519302930766 http://www.sciencedirect.com/science/article/pii/S0022519302930766 http://dx.doi.org/10.1186/1752-0509-6-1 http://www.biomedcentral.com/1752-0509/6/1 http://www.biomedcentral.com/1752-0509/6/1 http://dx.doi.org/10.1073/pnas.0703678104 http://dx.doi.org/10.1073/pnas.0703678104 http://www.pnas.org/cgi/doi/10.1073/pnas.0703678104 http://www.pnas.org/cgi/doi/10.1073/pnas.0703678104 http://dx.doi.org/10.1007/s11103-011-9767-z.Plastid http://dx.doi.org/10.1007/s11103-011-9767-z.Plastid http://arxiv.org/abs/NIHMS150003 https://www.cdc.gov/parasites/leishmaniasis/ https://www.cdc.gov/parasites/leishmaniasis/ http://dx.doi.org/10.1016/j.idm.2017.03.002 http://ac.els-cdn.com/S2468042716300173/1-s2.0-S2468042716300173-main.pdf?{\_}tid=c8ed3798-8db2-11e7-b1dc-00000aacb360{\&}acdnat=1504118647{\_}14d164bb702c07064ef2e80eeff2e974 http://ac.els-cdn.com/S2468042716300173/1-s2.0-S2468042716300173-main.pdf?{\_}tid=c8ed3798-8db2-11e7-b1dc-00000aacb360{\&}acdnat=1504118647{\_}14d164bb702c07064ef2e80eeff2e974 http://ac.els-cdn.com/S2468042716300173/1-s2.0-S2468042716300173-main.pdf?{\_}tid=c8ed3798-8db2-11e7-b1dc-00000aacb360{\&}acdnat=1504118647{\_}14d164bb702c07064ef2e80eeff2e974 http://ac.els-cdn.com/S2468042716300173/1-s2.0-S2468042716300173-main.pdf?{\_}tid=c8ed3798-8db2-11e7-b1dc-00000aacb360{\&}acdnat=1504118647{\_}14d164bb702c07064ef2e80eeff2e974 http://ac.els-cdn.com/S2468042716300173/1-s2.0-S2468042716300173-main.pdf?{\_}tid=c8ed3798-8db2-11e7-b1dc-00000aacb360{\&}acdnat=1504118647{\_}14d164bb702c07064ef2e80eeff2e974 http://ac.els-cdn.com/S2468042716300173/1-s2.0-S2468042716300173-main.pdf?{\_}tid=c8ed3798-8db2-11e7-b1dc-00000aacb360{\&}acdnat=1504118647{\_}14d164bb702c07064ef2e80eeff2e974 http://dx.doi.org/10.1016/j.mbs.2016.02.015 http://dx.doi.org/10.1016/j.mbs.2016.02.015 http://linkinghub.elsevier.com/retrieve/pii/S0025556416000468 http://linkinghub.elsevier.com/retrieve/pii/S0025556416000468 http://linkinghub.elsevier.com/retrieve/pii/S0025556416000468 http://www.who.int/leishmaniasis/en/ http://www.who.int/leishmaniasis/en/ http://dx.doi.org/10.11145/j.biomath.2019.01.026 Introduction The Model The Within-Host System The Between-Host System Linking the Within- and Between-Host Systems Analysis Analysis of the Within-Host System Analysis of the Full Immuno-Epidemiological Model Discussion and Conclusion