ANALYSIS AND SIMULATIONS WITH A MULTI-SCALE MODEL OF CANINE VISCERAL LEISHMANIASIS

Visceral leishmaniasis in dogs is believed to have an impact on the prevalence of the disease in human populations. Here, we continue the analysis of the nested immuno-epidemiological model of visceral leishmaniasis in dogs, including a proof of well-posedness using functional analytical methods. Once well-posedness is established, we continue stability analysis of the endemic equilibria and provide necessary and sufficient conditions for the presence of backward bifurcation, and prove the instability of the lower endemic equilibrium in the presence of backward bifurcation. Lastly, we provide a number of simulations of the model using a number of control strategies. Control measures currently in use attempt to reduce the parasite load in the host, reduce the vector population, reduce the vector biting rate, and remove infected hosts. We examine various combinations of these strategies and conclude that a strategy combining culling infected dogs and removing vectors from the population by means such as insecticide will be the most effective. Mathematics Subject Classification. 92D30. Received January 12, 2020. Accepted June 29, 2020.

differences within an individual compartment. Nested models consist of an ODE within-host system using timesince-infection and a between-host system using time-since-infection and chronological time. This method was extended to vector-host models by Rock, Wood, and Keeling [18].
The decision to use a multi-scale approach is informed by the work of Williams and Dye, noting that the vector lifespan is significantly shorter than that of the host [24]. Additionally, we can use a nested model to examine the between-host dynamics using a limited amount of within-host data [10].
Canine visceral leishmaniasis has been modeled for over thirty years with a variety of models, beginning with a standard SIR model. The first model was introduced by Dye [6], and was comprised entirely of ordinary differential equations. Since then, most models of visceral leishmaniasis have been ODE systems as well [17]. Some models, such as [14], included an exposed but not infectious class, while others, such as [3,12,19] included asymptomatic hosts. Previous work has posed that asymptomatic hosts may be one of the driving forces of the disease [4].
New models require new analysis. In biomathematical models, the well-posedness is often overlooked, instead displaying traditional results. In the classical sense, a system is well-posed if a unique solution which depends continuously on the input exists for each initial value. Additionally, in biological systems, we hope to have realistic and biologically relevant solutions, e.g., nonnegative inputs yielding nonnegative outputs.
In this paper, we prove well-posedness of the multi-scale immuno-epidemiological model found in [22]. For context, we will first introduce and give a brief motivation of the immuno-epidemiological model of zoonotic visceral leishmaniasis established in [22]. Some of the value of ODE systems is the relative simplicity in its analysis. However, the introduction of time-since-infection in an epidemiological model necessitates the use of partial differential equations, which complicates analysis. In systems with PDEs, showing well-posedness is a nontrivial task that has often relied on the tools of semi-group theory, an example of which can be seen in Browne [1]. We introduce, instead, a method to display well-posedness of PDE systems through functionalanalytical methods, yielding a, perhaps, more tangible method. The foundation of this method is accessible in [11]. In particular, we provide a proof that the between-host subsystem, structured by time-since-infection, is well-posed.
Following the proof of well-posedness, we continue the stability analysis of equilibria started in [22]. This includes the derivation of the necessary and sufficient conditions for the presence of backward bifurcation, as well as a proof of instability of the lower of the two equilibria.
Finally, we present a number of simulations of the full model. We begin by showing initial values exist that are consistent with previous analysis, as well as the occurrence of bistability in the presence of backward bifurcation. We then examine the effects of drug therapy on the model and conclude with the results of applying multiple control measures simultaneously.

Model formulation
We begin with a brief introduction of the model. While the state variables and parameters will be described, some of the specifics will be omitted, and can be examined further in previous work by Welker and Martcheva [22]. A full list of state variables can be found in Table 1. The within-host system examines parasite loads and immune response, while the between-host system is of vector-host type structured by time-since-infection.

The within-host system
The within-host system consists of three first order ordinary differential equations and focuses on parasite loads in the skin tissue and bone marrow of dogs, and the immunoglobulin G (IgG) immune response. In the system below, P S (x) is the parasite load in the skin, P O (x) is that of the bone marrow, and G(x) is the IgG concentration triggered by the presence of the parasite, where x is the time-since-infection. A list of the parameters can be found in Table 2. The subscripts S and O will indicate skin tissue and bone marrow, We assume reproduction of the parasite is limited by carrying capacities K S and K O . Recruitment is modeled logistically with rates r S and r O . The parasite is deposited in the skin, and can travel to the bone marrow at rate k S . The parasite in the bone marrow can travel inside a blood cell back to the skin at rate k O . The parameter ρ is a density conversion coefficient.
The natural death rate of the parasite is included in the logistic terms, but the clearance rate due to immune response is separate. The parameters ε S and ε O are the clearance rates of the parasite induced by IgG. Parasiteinduced IgG production occurs at rates a S and a O , and d is the IgG clearance rate.

The between-host system
The between-host system is of vector-host type. The subscripts H and V will indicate the host and vector, respectively. The host system consists of two first order ordinary differential equations for the susceptible  (S H (t)) and recovered (R H (t)) hosts, as well as a first order partial differential equation for the infected host class (I H (t)), which is structured by time-since-infection, i.e., In the above equation, i(t, x) represents the density of infected hosts at time t with time-since-infection x. The total number of hosts is given by N H (t) = S H (t) + I H (t) + R H (t). The vector system consists of three first order ordinary differential equations for susceptible (S V (t)), carrier (C V (t)), and infectious (I V (t)) classes. The between-host system is given below, and a list of the parameters can be found in Table 3. The time unit used will be t unless x is specified. The flow chart of the model is given by Figures 1 and 2.
The flow-chart of the between-host system. (2.10) Reproduction rates are given by Λ H and Λ V . Natural death rates are given by m H and m V . The average biting rate is given by a. Hosts become infected at rate β H after being bitten by an infectious vector, and vectors become carriers at rate β V after biting an infectious host. Infected hosts, in the presence of treatment, will move to the recovered class at rate σ. Infected hosts will die at rate µ. Recovered hosts will lose immunity at rate γ and become susceptible again. Carrier vectors will become infectious at rate τ . Lastly, k u is a conversion constant between the time scales, i.e., x = k u t.

Linking the within-and between-host systems
The time scale of the vector occurs at a much faster rate than that of the hosts. Thus, we replace certain parameters in the between-host system with functions dependent on time-since-infection. The full list of parameters used in the linking functions is given in Table 4. Courtenay et al. [4] concluded that the best indicator of infectiousness to a susceptible vector was the parasite load in the skin of the host. Hence P S (x) is used to determine the transmission rate β V : where ξ is the rate of decay, c V is the maximal transmission coefficient, and ν = −2/ ln(10) [13]. The recovery rate σ, in the presence of treatment, is also dependent on time-since-infection. It is given by where δ 0 , θ S , θ O , and κ are constants [21]. Lastly, disease-induced mortality µ is also dependent on time-sinceinfection and is given by where ψ S , ψ O , and η are constants [9].

Well-posedness
In this section we show that the immuno-epidemiological model of Leishmaniasis introduced in the previous section is well-posed, that is for every set of non-negative initial conditions there exists a unique solution which is also nonnegative. We first prove the well-posedness of the within-host model.
Then, for every y • ∈ R n + , there exists a unique solution of y = F (t, y), y(0) = y • , with values in R n + , which is defined on some To these ends, let y = (P S , P O , G), F 1 (t, y) = P S (t), F 2 (t, y) = P O (t), and F 3 (t, y) = G (t). Then, for T > 0, y ∈ C[0, T ] × C[0, T ] × C[0, T ] =: Ω 1 , and we note that Ω 1 is a Banach space. Since F (t, y) is continuously differentiable, it is locally Lipschitz, as shown in the following Lemma, which makes use of the Contraction Mapping Theorem.
Proof. Let f be continuously differentiable on X. Then, for all p ∈ X, there exists a convex ball B r (p) such that all partial derivatives are bounded by some K > 0. Hence, for ||x − p|| ≤ r, there exists a constant C(r) ≥ 0 such that ||Df (x)|| ≤ C(r), where Df is the Jacobian of f . Let x, y ∈ B r (p). Define Hence for all x, y ∈ B r (p). Therefore f is locally Lipschitz with constant C(r).
Next, assume y ∈ R 3 + . Suppose y 1 = P S = 0. Then Next, consider the modified system Let H S (P S ) and H O (P O ) be the first and second terms respectively in (3.4). Note that H S (P S ) is quadratic with a maximum at P S = K S /2. Hence there exist non-unique ζ S , ω S > 0 such that ζ S − ω S ρP S ≥ H S (P S ) for all P S (see the Lemma below).
for all x.
Proof. Let Then Hence f has an extremum at which is a minimum since f (x) = 2r/K > 0 for all x. Let ζ = (ω + r)K. Then for ω < 3r. Therefore ζ = (ω + r)K and 0 < ω < 3r satisfy the inequality for all x.
By a similar argument, there exist non-unique Let P = ρP S + P O and H(P ) = ζ − ωP . Note that H(P ) = 0 for P = ζ/ω. Integrating P ≤ ζ − ωP , we have Then Since P S , P O , and G are bounded, we have that b = ∞, and the solutions exist and are nonnegative for all t.
These results lead to the following theorem.
Lastly, suppose (P S (0), P O (0), G(0)) / ∈ S . Then there exists a t * such that Thus we conclude that the within-host system is well-posed.

Well-posedness of the between-host system
We have the following theorem: The proof of this theorem can be found in Appendix B.

Analysis of equilibria
In [22] we defined the basic reproduction number as: The between-host system has a disease-free equilibrium which is locally asymptotically stable when R 0 < 1 and unstable when R 0 > 1 [22].

Endemic equilibria
Recall that 10) Substituting (4.8) into (4.2), we get which yields (4.14) Next, note that Substituting (4.11), (4.14), and (4.15) into the initial condition yields Dividing both sides of the above by i 0 and m H , we see that We rearrange (4.16) as follows:

Necessary and sufficient condition for backward bifurcation
To obtain the condition for backward bifurcation, we differentiate (4.18) with respect to β H , noting that i 0 is a function of β H . The derivative of the left-hand side is The derivative of the right-hand side is Next, we note that at the critical value where R 0 = 1, we have that i 0 = 0. Evaluating both the left-and right-hand sides of the derivative at i 0 = 0, we obtain Denote the bracketed term above and the right-hand side as G and H , respectively. Then Since H is always positive, G must be negative, i.e., Next, we solve the equation R 0 = 1 for β H and obtain which must be negative. Thus, we can multiply the above by m H /Λ H to obtain which yields the bifurcation condition Denote the right hand side of (4.17) as F(i 0 ). Then which may be positive or negative. Then with stability if the left hand side is greater than the right hand side.

Linearization and stability of endemic equilibria
In this section we show the lack of stability of the lower of the two endemic equilibria when backward bifurcation occurs. Consider the following between-host system linearized around a generic equilibrium: (At this point, we will abandon the * notation for convenience.) Then and thus Furthermore, n V = s V + c V + y V , and From the above equation, we have either λ = −m V or n V = 0, which would imply s V + c V + y V = 0. Suppose the last equation holds. Then Hence Recall that From (4.29), The above equation substituted into (4.37) yields From (4.28), we have We can rewrite equation (4.25) as Substituting (4.36), (4.38), and (4.40) into (4.27), we see that or, dividing by y(0), Next, recall We first work with the first term on the right hand side of (4.42). First, note that Then, since we have that Thus Next, we note that the second and third terms of (4.42) can be written as respectively. Substituting the terms into (4.42), we have that (4.43) Using (4.24) in the above, we see that (4.44) Thus for the lower endemic equilibrium F (i 0 ) < 0 and therefore χ(0) > 1. Since for λ real χ(λ) → 0 as λ → ∞, then there exists λ * real and positive such that χ(λ * ) = 1. Hence, the lower equilibrium is unstable.

Simulations
We now present a number of simulations beginning with standard results and concluding with the implementation of control strategies. For each scenario, the figure displaying the number of infected hosts (I H (t)) is shown, but other results may be discussed as well. Parameter values can be found in Table A.1 unless otherwise noted.

Disease-free equilibrium
First, we present the result of convergence to the disease-free equilibrium E 0 when R 0 < 1. In Figure 3, we observe that it takes less than 1000 days to eradicate the disease.
This simulation includes no treatment, but the population of the vector is significantly smaller than that of the hosts (N V (0) = 500 compared to N H (0) > 10 5 ). This may lead to the susceptible vector being less likely to bite an infected host, thereby lowering the likelihood of the disease spreading to the vector population.

Endemic equilibria
5.2.1. Endemic equilibrium with R 0 > 1 Next, we examine the first case to result in convergence to an endemic equilibrium. Figure 4 presents the result in which R 0 > 1. We observe that the system converges to an endemic equilibrium in less than 1000 days, following an outbreak that occurs at approximately day 500.
The major difference between this system and the one presented in Section 5.1 is the size of the initial populations. Here, we have a significantly larger vector population than host population. This results in a much greater chance that an infected host will be bitten by a susceptible vector. Additionally, we observe that the final host population size is less than one-tenth of the initial size -a significant reduction.

Endemic equilibrium with
In [22], the presence of backward bifurcation was shown to occur under certain circumstances. Here, we observe such a case. In Figure 5, we observe convergence to an endemic equilibrium when R 0 < 1. It is worth noting that the condition shown to indicate the presence of backward bifurcation in (4.2) has been met.
While the system does converge to an endemic equilibrium, it does so in a drastically different manner than the system presented in (5.2.1). Over the same timespan, it would appear as though no outbreak occurs, and  the number of infected hosts equilibrates. However, we observe an outbreak that occurs approximately 10000 days later, followed by convergence to a higher endemic equilibrium.

Bistability
In the Figure 6, we show convergence to separate equilibria based solely on a change in the initial amount of susceptible hosts. This occurs when R 0 < 1 and backward bifurcation is present.  It is worth noting when, with everything else held constant, I H (t) switches from converging to the endemic equilibrium to the disease-free equilibrium. With the given values, we see that if S H (0) ≤ 10750, then I H (t) converges to the endemic equilibrium. Otherwise, I H (t) converges to the disease-free equilibrium. This can be observed in Figure 6. Hence, according to this simulation, we can infer that a larger susceptible population can force the disease to die out. This could be due to factors such as host choice. The infected vector has a larger susceptible host population to infect, but the susceptible vector has a lower probability of biting an infected host.

Introducing treatment by increasing κ
Up to this point, no simulations shown have involved treatment. In this section, we present a simulation involving the treatment of visceral leishmaniasis. In the Figure 7, we introduce treatment simply by changing the value of κ, the coefficient of the linking function for recovery, while holding all other parameters constant. This does not have a real world parallel, but instead is used to show the effects of introducing some hypothetical method of treatment on the prevalence of the disease. In this case, we also set γ, the waning immunity rate, to 1.
The reproduction number corresponding to each value of κ, all of which are less than 1, can be found in Table 5. In this simulation, each solution converges to the disease-free equilibrium, but the behaviors of the solutions are worth noting.
First, it is clear that the initial minor change in κ, from zero to 10 −10 , greatly reduces the peak of the initial outbreak. In fact, this change results in a maximum number of infected hosts approximately one-third the size of an untreated population. Further increases in κ decrease the maximum value of I H (t) even more. Obviously, the value of decreasing the number of infected hosts cannot be overstated. Fewer hosts having to deal with symptoms and treatment is important, but the decrease in the total number of infected hosts is vital to lowering the chance a vector can be infected when biting a randomly chosen host.
The second observation relates to when the peak of the outbreak occurs. Differences between consecutive values of κ are not as obvious, but when κ = 0, the peak occurs after 1000 days, while the significantly smaller peak when κ = 5 × 10 −9 occurs very close to the beginning of the simulation.
Lastly, we observe the reduction in time it takes for the solutions to converge to the disease-free equilibrium. When κ = 0, we see I H (t) is approximately zero at around 3900 days. The initial increase of κ to 10 −10 reduces Table 5. Effect of κ on R 0 . Increasing κ seems to decrease the reproduction number.
κ R 0 0 0.362185 10 −10 0.358694 1 × 10 −9 0.328863 2 × 10 −9 0.298840 3 × 10 −9 0.271773 4 × 10 −9 0.247359 5 × 10 −9 0.225325 this to approximately 3000 days. Increasing κ further clearly decreases the amount of time it takes for I H (t) to be close to zero, with κ = 5 × 10 −9 reducing the time to approximately 1000 days. The two previous observations combined imply that increasing treatment efficacy should, sometimes greatly, reduce the amount of time a population must work towards disease eradication. This is particularly important in populations with higher poverty levels that may not be prepared for the disease financially.

Drug treatment
According to the WHO [23], current drug treatments include amphotericin B deoxycholate (AB), paromomycin (PM), pentavalent antimonials (Sb v ), and miltefosine (MF). Each treatment has complications related to efficacy, toxicity, price, and schedule. The only drug treated orally is MF.
We first observe simulations below in which Sb v is administered via intramuscular injection for the required treatment period of 30 days. Much is unknown regarding the pharmacodynamics of Sb v , but a broad generalization of the believed mechanism is the limiting of available cellular energy [23], thereby reducing the ability of the parasite to reproduce within the host. All of the remaining simulations use κ = 10 −10 and γ = 1/24.
After Sb v , we look at simulations using a different drug, miltefosine (MF). Not only is MF the only VL treatment administered orally, but it seems to be more effective than Sb v , which may have a declining effect compared to when it was originally used [23].

Successful pentavalent antimony treatment (Treatment with Sb v )
In Figure 8, we observe the treatment resulting in the removal of the disease from the population. The within-host system is modified to where represents the effect of Sb v on the reproduction of the parasite. The between-host system is left unchanged. During the 30 days of treatment, the reproduction rates r O and r S are multiplied by 0 < s 0 < 1. When treatment ceases, the system above reduces to the original within-host system for the remainder of the simulation. The simulations of the within-host system seem as expected. The general shapes of each graph are the same, perhaps with a slightly lower peak. However, the impact of treatment on the between-host system is drastic. Not only is maximum prevalence approximately halved, but we observe treatment forcing convergence to the disease-free equilibrium. Therefore, we should assume that treatment with pentavalent antimony, if successful to a certain point, can remove the disease from a population.

Unsuccessful pentavalent antimony treatment
Unfortunately, a rather high efficacy of treatment was considered above. We consider an alternate situation. It is possible that asymptomatic infected dogs contribute to the persistence of the disease. Is it possible that treated dogs, who may not exhibit any signs of the disease, may negatively impact the host population as a whole? In Figure 9, we observe the case in which everything is the same as the previous situation except for a decreased efficacy of the Sb v treatment.
As in the previous case, the within-host simulations are as expected. Unfortunately, we see a negative effect of treatment regarding the infected host population. Not only do we have persistence of the disease, but the endemic equilibrium of the treatment system is in fact higher than that of the untreated system. It is hard to say if this is a realistic situation. However, this simulation should inspire caution with ineffective treatment.

Delayed pentavalent antimony treatment
A strong, unrealistic assumption of this treatment model is that treatment begins upon initial infection. This can be remedied by altering S(t) to  for some t 0 > 0. Figure 10 illustrates the impact of a 5 day delay (t 0 = 5) on the successful Sb v treatment simulation shown in Figure 8.
This result shows the necessity of quickly implementing Sb v treatment. Whereas immediate treatment would lead to the removal of the disease from the population, a delay (in this case, of 5 days) shows convergence to a higher endemic equilibrium, similar to that of the previous simulation.

Successful miltefosine treatment -Case 1
Like Sb v , the actual mechanisms of MF are not well understood. However, it is possible that MF directly affects the parasites themselves, in contrast with the possible mechanism of Sb v . For this reason, we observe two different effects on the within-host system, both of which result in elimination of the disease.
In the first case, we simply adjust the last terms in the P S and P O equations, i.e., The result of the treatment by this method can be seen in Figure 11. An odd consequence of MF treatment in this model is a higher maximum of the skin parasite load that occurs later than that of the untreated model. However, we see that treatment causes prevalence to peak at less than half the value of the untreated system, followed by convergence to the disease-free equilibrium.

Successful miltefosine treatment -Case 2
On the other hand, we could augment the original within-host system by adding a new term, i.e., The result of treatment by this method can be found in Figure 12.
The behavior shown in Figure 12 is similar to the first case of MF treatment shown in Figure 11. However, we will continue with the assumption that this formulation of MF treatment is accurate.

Delayed miltefosine treatment
Since MF treatment is believed to be more effective than other drug treatments currently, we examine multiple cases of delayed MF treatment. The treatment period still lasts for 28 days, but the start time is varied. Thus If t 0 = 0, i.e., treatment starts immediately, the model reduces to that presented in Figure 12. This model converges to the disease-free equilibrium at approximately 5900 days. Figure 13 shows the effect of delaying MF treatment. With just a ten day delay in beginning treatment, convergence to the disease-free equilibrium takes approximately 4000 more days than the immediate treatment model.
A delay in just one more day causes the disease to persist. Relative to the base model, an 11 day delay will cause prevalence to peak at a much later date and a slightly higher level. Additionally, an 11 day delay resulted in the highest endemic equilibrium presented. Increasing the delay to 50 days yielded a result much closer to the untreated model, but still converged to a higher endemic equilibrium.
Increasing the delay further yielded results moving closer to the untreated model. A 200 day delay resulted in an almost identical I H (t) curve to that of the untreated model. However, while most of the figures showed similarities between the untreated and 200 day delay scenarios, such a large delay resulted in a massive increase in infectious vectors. Furthermore, a 200 day delay resulted in interesting behavior in the skin and bone marrow parasite loads, including an obvious loss of differentiability.
Relative to Sb v , it would appear that MF treatment may be more forgiving regarding a delay in treatment. This would be consistent with the belief that MF is a more effective treatment. However, since these parameters are synthetic, such a result is not guaranteed. Since Sb v use is declining, we will move forward with the assumption that MF is, in fact, more effective than Sb v . Thus, we will assume MF is used in all within-host treatment. In order to examine the effect of multiple control measures, we augment the full model as follows: where the control parameters (in bold) are given on Table 6. In Section 5.3, we examined the effect of drug treatment for infected dogs. In addition to drug treatment (m), we now examine different control strategies including the use of culling infected dogs (c), increasing vector mortality with methods such as insecticide-treated bed nets (ITNs) and insecticide spray (n), and decreasing bite rate with methods such as repellent-impregnated collars (r). Since these parameters are arbitrarily chosen, we include the percent increase needed to switch convergence from an endemic equilibrium (EE) to the disease-free equilibrium (DFE).
Additionally, we note that authors such as Dye [7] believed culling strategies would be ineffective, while others such as Reithinger et al. [15] believed it could reduce prevalence of ZVL. Thus, we assume that culling will not be solely sufficient in eradicating the disease, but may be effective in conjunction with other control measures. As shown on Table 7, letting c = .0014 and setting all other control parameters to zero would eradicate the disease. Thus, for the simulations in which culling is utilized, we let c = .0013 (unless otherwise noted), which would not itself remove the disease from the population.

Miltefosine treatment and insecticide/bed nets
We first consider the control strategy in which infected dogs are treated with Miltefosine upon infection and vector mortality is increased by control measures such as insecticide and bed nets. We assume the Miltefosine removes 1% of the parasite per unit time x. This is not necessarily practical since ZVL is considered a veterinary problem instead of a public health problem, i.e., public funds are not typically used for interventions. As shown on Table 8, we see that increasing vector mortality by an order of 10 −4 is necessary to eradicate the disease.

Culling and repellent
Since treating ZVL may not be practical, we consider alternative control strategies. In Table 9, we observe the results of using control measures such as insect repellent to reduce the biting rate in the presence of a culling strategy. To successfully eliminate the disease from the population, the biting rate must be reduced by a factor of the order of 10 −3 . It is not predicted that repellent on its own would be enough to eradicate the disease.

Culling, repellent, and insecticide/bed nets
Authors such as Ribas et al. [16] have predicted that vector control may be an essential component of control strategies. In Table 10, we see that control measures that increase vector mortality by an order of 10 −5 are effective in the presence of culling. Compared to the results shown on Table 8, we see that culling, in lieu of drug treatment, reduces the burden of control measures such as insecticides and bed nets that increase vector mortality.
Including repellent further reduces the burden on insecticides and bed nets, as shown in Table 11. However, it is not a significant reduction from the results shown in Table 10. On the other hand, the efficacy of repellent is of the order 10 −4 -a decrease from the results shown in Table 9. Additionally, it is worth noting that a small decrease in n requires a larger increase in r to eliminate the disease, which is consistent with the notion of vector control being more effective. Lastly, we consider reducing the efficacy of culling with values of n and r comparable to those found in Table 11. In Table 12, we observe that even a slight reduction in culling efficacy results in the persistence of the disease. This suggests that culling could play an important role in an effective control strategy.

Discussion and conclusion
In this paper, we continued work on the model for zoonotic visceral leishmaniasis originally presented in [22]. We first presented a proof of well-posedness of the model using functional-analytical methods. Then, we derived the necessary and sufficient condition for the presence of backward bifurcation and analyzed the stability of the lower endemic equilibria. Lastly, we presented numerous simulations, concluding with observations regarding control strategies.
The proof of well-posedness of the ODE within-host system is found using standard methods. In the proof of well-posedness of the between-host system, existence is derived as a fixed point of a system of integral operators, where the system is written as integral sequences. As we are working in a Banach space, we show that the sequence generated iteratively from starting points is Cauchy in the Banach space. We then show uniqueness by assuming two fixed points and proving they are the same using Grönwall's inequality.
We then confirm the presence of backward bifurcation stated in [22] and derive the necessary and sufficient condition for the presence of backward bifurcation in the model -a consequence of standard incidence and disease-induced mortality. This is consistent with the results found in [8]. We then use a variant of the betweenhost system linearized around a generic equilibrium to prove the lack of stability of the lower of the two endemic equilibria when backward bifurcation occurs.
Lastly, we presented simulations of the model. First, we established the existence of initial values such that the simulations agreed with theoretical results. We showed convergence to the disease-free equilibrium when R 0 < 1, convergence to an endemic equilibrium when R 0 > 1, and multiple endemic equilibria when R 0 < 1 and the condition for backward bifurcation is met. A noteworthy result is that increasing the ratio of vectors to hosts seems to lead to persistence of the disease. Consequentially, it stands to reason that vector control may be one of the strongest control measures.
We then began to apply control measures to the model. Simulations displaying successful and unsuccessful immediate drug treatment were shown. However, since diagnostics for ZVL are currently unreliable, immediate treatment is unrealistic. Thus, we considered the impact of a delayed start of treatment, and found that even a short delay of less than two weeks could result in a failed treatment and an even higher endemic equilibrium. Hence we conclude that, due to the difficulties in reliably and quickly diagnosing ZVL and the lack of public funds, drug treatment may not be a strong control measure.
The last simulations presented the impact of control strategies utilizing multiple control measures. These control measures were grouped by their effect. Culling (c) removes infected hosts from the population, control measures such as bed nets and insecticide remove vectors from the population (n), and repellents which decrease the biting rate (r). We assumed that culling would not successfully eliminate the disease on its own, and considered its effectiveness when used alongside other control measures.
Given a fixed culling rate, we observed that the order of n is smaller than the order of r necessary to eliminate the disease in the population. In other words, it will take less vector control than repellent to achieve a successful control strategy, which agrees with conclusions drawn from simulations in Section 5.4.4.

Appendix A. Parameters
Appendix B. Proof of well-posedness of the between-host system Proof. Consider the modified between-host system Take 0 < ε < 1. Define the set First, we show that F maps S into itself. Consider Then, since we have that the inequality is trivially satisfied. Finally, To show the output of the map F is also bounded, we show that for F 1 . Similar ideas can be applied to F i for i = 2, ..., 6. To these ends, we have Then It is important to note that S is a subset of a Banach space X. To complete showing that F(u(t)) maps S to itself, we note that Given u 0 = 0, define iteratively the sequence u n+1 = F (u n ). where σ = sup x σ(x). Considering the last term, we have  Thus, if we let N n+1 (t) = F 1 (u n ) − F 1 (u n−1 ) + ∞ 0 F 2 (u n ) − F 2 (u n−1 ) dx + · · · + F 6 (u n ) − F 6 (u n−1 ) , or u n+1 X ≤ ||u n || X t.
Next, define Next, we consider the sequence N n (t) = N n H (t) + N n V (t) for n = 0, 1, . . . where u n (t) = F(u n−1 (t)) with some given u 0 (t) which may be taken as zero. Then Let ||·|| C = max 0≤t≤T |·|. We have . . . N n − N n−1 C ≤ K(α − m H ) n−1 t n n! + K(α − m V ) n−1 t n n! ≤ K(K 1 ) n−1 (T ) n n! , which tends to 0 as n → ∞. In fact, we have We have the following definition: Definition B.1. Let X be a Banach space with norm ||·|| X . A Cauchy sequence (x n ) ⊂ X is a sequence such that for every ε > 0, there exists a N > 0 such that for all m, n > N , ||x m − x n || X < ε.
which goes to 0 as m → ∞. Thus, the sequence is Cauchy and it has a limit. The limit of this sequence is the fixed point of F. It is important to note that in the case of our time-since-infection structured model (as opposed to an age-structured model), the above proof only holds if the recovery rate σ(x) and the disease-induced mortality µ(x) are bounded. Lastly, to prove uniqueness, assume that F has two fixed points, i.e., u = F (u) and u = F (u). Let As before, we obtain Then Grönwall's inequality implies N (t) ≡ 0. Hence S H = S H , i = i, R H = R H , I V = I V , C V = C V , and S V = S V . Thus u = u, and therefore the fixed point is unique, and we conclude that the model is well-posed.