Examining the reservoir potential of animal species for Leishmania infantum\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varvec{\textit{Leishmania infantum}}$$\end{document} infection

Leishamaniasis is a parasitic disease transmitted by sandfly Phlebotomus spp. and is seen in tropical and subtropical countries in which an estimated 12 million persons are infected. Among various types of leishmaniasis, zoonotic visceral leishmaniasis (ZVL) caused by Leishmania infantum is an important amphixenosis shared by human and other animals. Although identifying the natural reservoir host would help better understand the transmission dynamics of Leishmania spp., little effort has been made to quantitatively clarify the dynamics involving the reservoir host of ZVL. The present study investigated the reservoir potential of four wild animals in maintaining ZVL, using prevalence data from Latin American countries in Amazons and examining the role of crab-eating fox, spiny rat, common opossum and black rat in maintaining the transmission. Reflecting frequent reinfections, a susceptible-infected-susceptible model was employed, enabling us to estimate model parameters from endemic prevalence data. The next generation matrix of the multi-host system was computed, permitting us to theoretically examine the reservoir potential of each animal species. Our estimates indicated that there is no unique reservoir host consisting of single animal species. Crab eating fox was considered to play an important role in maintaining L. infantum transmission, but this was the case only in combination with other hosts. The present study indicates that animal species other than canine play important roles in maintaining transmission of Leishmania infantum, which is different from conventional wisdom that centered on the importance of canine only. Greater sample size with additional entomological and genetic insights into inter-specific contact would be required to implement more explicit assessments.


Introduction
Leishamaniasis is a parasitic amphixenosis caused by Leishamania spp. that belongs to Trypanosomatidae. The disease is transmitted by sandfly Phlebotomus spp. and classified as visceral and cutaneous leishmaniases [1][2][3][4][5][6][7][8][9][10]. Leishmaniasis is prevalent in tropical and subtropical 88 countries, and an estimated 12 million persons are infected. The population at risk is as large as 350 million. Visceral leishmaniasis takes several years before the infected individual develops symptoms. Among the cases of visceral leshmaniasis, hypertrophy of spleen and liver is observed. The patient may often naturally recover, but if it is left untreated, severe cases could result in death within two years of illness onset [1][2][3][4][5][6][7][8][9][10]. Among various types of leishmaniasis, zoonotic visceral leishmaniasis (ZVL) caused by Leishmania infantum is an impotant amphixenosis that is shared by both human and other animal species. Historically, ZVL was observed in Asia, Middle East and Europe. It was then geographically spread to Latin America, and in the present day, it is mainly seen in India, Nepal, Bangladesh, Sudan and Brazil [1].
The reservoir host of an infectious disease is the animal species that can maintain transmission, and thus, allow persistence of the disease. Such animal species can act as the source of disease for humans. If the reservoir host is identified, public health experts could consider potential countermeasures against the infection. Moreover, one can better understand the mechanisms of transmission dynamics [10,13]. In mathematical sense, one can identify the reservoir host using various approaches such as by analyzing epidemiological datasets [10,[12][13][14]. For instance, Nishiura et al. [10] have defined the reservoir host by examining the next generation matrix of a multi-host system. In a population with four different animal species, a 4-by-4 next generation matrix K = R i j is defined, i.e., characterizing secondary transmissions through interactions within and between four animal species, where R i j gives the average number of secondary cases in host i caused by a single primary case in host j in a fully susceptible population. The basic reproduction number of this system is given by the dominant eigenvalue of K [11].
To understand the role of host 1, one has to consider two hypothetical settings, i.e., when the population consists of host 1 alone, and when the host 1 does not exist in the population. With regard to the former condition, the secondary transmission is described by R 11 only, and if R 11 > 1, it indicates that the host type 1 can maintain transmission on its own. Such host is referred to as the maintenance host. As for the latter condition, the population consists of remaining 3 species, and the next generation matrix of such hypothetical population K reads If the dominant eigenvalue of K is greater than 1, it indicates that the transmission could be maintained even in the absence of host type 1. If the dominant eigenvalue of K is smaller than 1, the host type 1 is regarded as an essential host for maintaining transmission. We regard the host species that can act as both maintenance and essential host as the reservoir host [10]. The reservoir host of ZVL has yet to be explicitly identified. Without understanding the reservoir potential of various animal species, it is difficult to consider possible countermeasures against animals for the prevention of human infection. The purpose of the present study is to quantitatively assess the reservoir dynamics of L. infantum infection in wildlife through the analysis of empirical data. Since publications of ZVL in wildlife setting are commonly and consistently seen in Amazons, Brazil and its surrounding areas, the present study focuses on published evidence from northwestern part of Latin America.

Materials and methods
While visceral leishmaniasis has been seen in various countries in tropical and subtropical zones, the present study focused on Amazon area from which particularly high incidence in humans has been reported. Since the infectious disease spreads across borders, not only Brazil but also Colombia and Venezuela are also included as the subject of our study. From these countries, wild animals have been sampled and researchers isolated L. infantum. Figure 1 shows the observed prevalence data (i.e. proportion positive based on cross sectional surveys) of four animal species, namely, crab-eating fox (Cerdocyon thous), spiny rat (Trichomys apereoides), common opossum (Didelphis marsupialis) and black rat (Rattus rattus) [1,15]. Along with the sample positives, the confidence intervals are also given, reflecting sample size for each species. The present study focuses on prevalence surveys that used polymerase chain reaction (PCR) method, because the prevalence of four species was commonly examined using PCR and the datasets based on other laboratory methods (e.g. parasitemia) were not consistently available across all species.
In Fig. 1, the confidence intervals are wide, because the prevalence surveys rested on small number of samples. Nevertheless, at first sight of the figure, one can observe that crab-eating fox yields the highest prevalence. A quick thought in relation to the reservoir dynamics may be that crab-eating fox plays an important role in maintain-

Fig. 1
Observed prevalence data of leishmaniasis by different animal species. From an animal species with the highest prevalence, we label the species in ascending order. The whiskers represent confidence intervals derived from binomial distribution ing transmission of L. infantum. However, it is more fruitful to explicitly identify the reservoir host based on objective analysis and thoroughly examine the reservoir potential of other animal species. To achieve an explicit and objective analysis, the present study employs a simple epidemiological model.
From experimental studies, it has been known that these wild animals do not often reveal symptoms and experience frequent re-infections [1,16,17]. Such evidence helps us to assume that all individuals never die of the infection and acquire immunity, and thus, to adopt an susceptible-infectious-susceptible (SIS) model [18][19][20][21][22][23][24]. Let x i and y i be proportions susceptible and infectious of species i, respectively, we have where γ i is the recovery rate of species i, μ i the natural death rate of i, and λ i is the force of infection of species i described by where β i j is the transmission rate from species j to i. The natural death rate is assumed as identical to the natural birth rate, and both x and y equally contribute to the natural birth event. Accordingly, in the first subequation of system (3), the demographic term The mean infectious period (or recovery rate) and the natural death rate are assumed known and can be extracted from literature. Assumed infectious periods and life expectancies are shown in Table 1 [25][26][27][28][29][30][31].  Five different contact matrixes for leishmaniasis in wild animals. Matrix B 1 assumes that the contact rates between different types of animals are smaller than contacts within the same species. The within-group contact is weighted by θ that can be interpreted as the proportion of contacts that are spent for within-group mixing. B 3 assumes that half of contacts are spent for between-group mixing and B 2 is intended to be between B 1 and B 3 . B 4 assumes that the contact rates of animal types that yielded higher prevalence than others would be more influential than others. B 5 assumes that withing-group mixing tends to be much higher than between-group mixing and that the between-group mixing is approximated by the parameter for type 4 Doing so, unknown parameters are only β i j and the estimates would help determine the next generation matrix of the multi-host system for L. infantum infection. Nevertheless, while we have 4 data inputs from empirical observation (Fig. 1) and assumes a stationary state [19,21,23], the matrix of β i j is 4 × 4 and the degree of freedom is insufficient. For this reason, we construct plausible contact matrices using only 4 parameters. Constructing multiple contact matrices, we address the uncertainty with respect to unobserved contact, as was also similarly practiced elsewhere [10,12,33]. Figure 2 shows the assumed five patterns of contact matrices. In principle, it is natural to assume that the transmission rate within the same species is greater than those occurring between species [10][11][12]32,33]. Matrices B 1 , B 2 and B 3 are referred to as the formulation that rests on assortative mixing assumption. The parameter θ is interpreted as the proportion of contacts spent for within group mixing. For B 1 , B 2 and B 3 , we assume θ at 0.9, 0.7 and 0.5, respectively. n i stands for the relative population size among all four species. Since we never have an access to such data, n i = 0.25 for any i for mathematical convenience. Matrices B 4 and B 5 qualitatively intends to capture the underlying transmission mechanism by allocating four parameters in sixteen entries [and such matrices are referred to as the matrices that describe who acquire infection from whom (WAIFW)]. Matrix B 4 assumes that the contact rates of animal types that yielded lower prevalence would be influential on others. Matrix B 5 assumes that the contacts between different animal species are very infrequent and non-diagonal elements are compensated by the contact rate of species with the lowest prevalence.
Maximum likelihood method was employed to infer unknown parameters of contact matrix. From the SIS model (4), it is evident that the prevalence at stationary state is given by Let k i be the sample size of species i among which m i were positive. The likelihood function to estimate β i j is Minimizing the negative logarithm of (6), parameters β i j are estimated, and subsequently, the next generation matrix K = R i j is quantified as R i j = β i j /(γ j + μ j ) [11,34,35].
Let P i be the projection matrix on type i, i.e., p ii = 1, and p i j = 0 for all other entries. The host-specific reproduction number U i [10] is given by the dominant eigenvalue of the matrix that includes only the type(s) of interest, i.e., where ρ(.) denotes the largest eigenvalue. It should be noted that the species i in this context can be either single animal species or a combination of multiple animal species. Whereas Funk et al. [13] discussed the next generation matrix K by explicitly separating the vector species from others, the present study simplifies the model and considers the R i j as the average number of secondary cases in animal species i caused by a single infected animal in j 'through sandfly bites'. This simplification was conducted because there aren't known multiple vector species for leishmaniasis. Let I be identity matrix. The host-excluded reproduction number Q i is given by the dominant eigenvalue of the matrix that excludes only the type(s) of interest, i.e., For a species or a combination of species, those satisfying U i > 1 is referred to as the maintenance host as it indicates that the presence of species can allow transmission to be maintained. Similarly, those satisfying Q i < 1 is referred to as the essential host, because the transmission cannot be continued in the absence of that (those) host(s). The minimum set of hosts that satisfy both U i > 1 and Q i < 1 is defined as the reservoir community [10,13]. If the condition is satisfied by a single species, such host is referred to as the unique reservoir host. Since the sampling distributions of U i and Q i are unclear, the 95 percent confidence intervals (CI) were derived from bootstrapping method. As mentioned above, we assumed that the natural death rate and recovery rate are known and fixed (Table 1). Whereas the natural death rate may not have a large impact on quantitative fate of the reservoir identification exercise, it is fruitful to examine the sensitivity of U i and Q i to the recovery rate γ i . Thus, we varied the recovery rate by multiplying 0.50, 0.75, 1.25 and 1.50 to the baseline value and examined how U i and Q i changes. For each set of assigned value, we implemented the maximum likelihood estimation and calculated the dominant eigenvalue of the next generation matrix.

Results
Quantifying the next generation matrix, the dominant eigenvalue yields the basic reproduction number R 0 , interpreted as the average number of secondary cases caused by a typical primary case, in the presence of all four species. The maximum likelihood estimates of R 0 for matrices B 1 , B 2 , B 3 , B 4 and B 5 were 1.32, 1.35, 1.43, 1.77 and 1.36, respectively, reflecting that the assumed stationary state is an endemic equilibrium. Table 2 shows maximum likelihood estimates of the host-specific reproduction number U i , while Table 3 shows maximum likelihood estimates of the host-excluded reproduction number Q i . U i for type 1 (crab-eating fox) was greater than 1 for matrices B 1 , B 2 , B 3 and B 5 , but not for B 4 . Moreover, even though type 1 was both maintenance and essential hosts using B 3 and B 5 , it was not the case for B 1 and B 2 . Thus, the reservoir potential of type 1 was not consistent across matrices, and moreover, type 1 was not regarded as the reservoir for more than half of the assumed matrices. Type 3 satisfied U i > 1 using B 1 , but otherwise U i > 1 was not observed for any other types of animal on their own.
As for combinations of host, a combination of types 1 and 2 as well as 1 and 4 was regarded as both maintenance and essential hosts using matrices B 2 , B 3 and B 5 . Another combination, types 1 and 3 was also shown to be the reservoir community except for matrix B 1 . Thus, only the combination of types 1 and 3 satisfied to be the reservoir community when matrix B 4 was used, presumably due to a mathematical reason that the third column of B 4 is calculated to be too small due to shorter infectious period of common oppossum compared with other animal species.
When B 1 (highly assortative contact matrix) was employed, all single type except for types 2 and 4 as well as all combinations of host satisfied U i > 1 due partly to the assumed extent of independence in transmission dynamics from other animal species. Nevertheless, all of them did not satisfy Q i < 1 (Table 3) indicating that the essential host does not exist for the combinations (of up to two animal species) we examined.
Combination of three species can be inspected from Tables 2 and 3 by conversely reading U i and Q i for single type of host. Namely, if we would like to know U i of a host combination of 1, 2 and 3, we should look at Q i for type 4. It is clear from the table that combinations of types 1, 2 and 3 and types 1, 3 and 4 satisfied both U i > 1 and Q i < 1 for all different types of matrices. Nevertheless, 95 % CI in Tables 2 and 3 2 1.30 (0.63, 3.70) 1.21 (0.56, 3.26) 1.12 (0.48, 2.97) 0.88 (0.08, 2.72) 1.15 (0.21, 3.70)   1 and 3 1.31 (0.65, 2.46) 1.29 (0.63, 2.45) 1.32 (0.61, 2.55) 1.53 (0.62, 2.   If Q < 1, the corresponding estimate is given in bold letters. Numbers in parentheses represent the 95 % confidence intervals were very wide, ranging from below to above the value of 1, not allowing us to fully judge the reservoir potential. Results from sensitivity analysis are shown in Fig. 3. Except for type 1 (Fig. 3a), the results of U i > 1 or U i < 1 did not greatly vary by substantially varying the infectious period γ i plus and minus 50 % of its original value. Similarly, except for type 1 (Fig. 3b), we did not observe any host with Q i < 1 by varying infectious period. Different reservoir dynamics were observed for type 1 mainly for assortative mixing, especially using B 1 , perhaps because of high sensitivity of U i and Q i to diagonal element of the next generation matrix, and thus, to the infectious period of Fig. 3 Sensitivity of the reservoir potential to infectious period of animals. The vertical axes represent the reproduction numbers that allow us to assess the reservoir potential, while horizontal axes represent relative change in the rate of recovery for each animal species. Panels a, c, e and g examines the sensitivity of the host-specific reproduction number, while b, d, f and h show the results of the host-excluded reproduction number. In each panel, the horizontal solid line represents the threshold value at which the reproduction number takes the value 1 the corresponding host type. Except for highly assortative matrices, the identification of the reservoir was not sensitively influenced by the mean infectious period.

Discussion
The present study examined the reservoir potential of four animal species, i.e., crabeating fox, spiny rat, common opossum and black rat, in allowing persistence of L. infantum transmission. While crab-eating fox yielded the highest prevalence among four species, it was not consistently regarded as the unique reservoir host across different contact matrices. All other types did not satisfy U i > 1 consistently. As for the combination of two types, i.e., types 1 and 2, 1and 3 and 1 and 4, both U i > 1 and Q i < 1 were satisfied for a part of assumed matrices, and as for the combination of three types, types 1, 2 and 3 and types 1, 3 and 4 satisfied both U i > 1 and Q i < 1 for all matrices. However, there has been no strong biological or ecological indication that type 1 (crab-eating fox) has some unique interaction with other animal species, especially in maintaining transmission of L. infantum. Moreover, Q i > 1 was consistently the case using B 1 which assumes that 90 % of contacts are spent for within-group mixing. Since Q i > 1 indicates that the remaining types of host could maintain transmission, the combinatory dynamics involving type 1 is not regarded as essential. All these results yield new insights into the reservoir dynamics of Leishmania spp.: previously, only canine was thought to act as the reservoir host of leishmaniasis as has been seen with L. donovani [1][2][3][4][5]. The present study is the first to explicitly indicate that animal species other than canine are likely to play important roles in maintaining transmission of L. infantum. Monitoring infections in animal species other than canine is deemed important to prevent human infection.
Matrix B 1 adopts highly assorattive mixing, and one could imagine that this matrix might accurately reflect the most realistic situation of multi-host transmission dynamics of L. infantum in wild life. It depends on the detailed behavior of sanflies (Phlebotomus spp.), but it is likely that different types of animals do not frequently interact from each other through common sandflies. As our analysis revealed, either single species or combination of two or more hosts allow transmission to be maintained. However, they are not regarded as essential. Strictly adhering to our definition of the reservoir host, it is likely that all types 1-4 act as the reservoir community as a single group. R 0 > 1 is in line with U i > 1 for four animals in combination, and the possible absence of transmission without four species is consistent with Q i < 1.
It should be noted that all uncertainty bounds (i.e. 95 % CI) in Tables 2 and 3 were very wide, ranging from below to above the value of 1. These did not enable us to more precisely interpret our findings, and it does indicate that the sample size should be increased in future studies.
If we have an opportunity to further examine the reservoir dynamics of Leishmania spp., there are number of points that could considerably improve our understanding and should be kept in our mind. First, an important source of uncertainty was seen in the assumed contact matrices that capture within-and between-host mixing through sandflies. One could survey the host preference of Phlebotomus spp. by exploring the biting behavior (e.g. how many are biting crab-eating fox and how many for others) and investigating their host specificities (e.g. how often inter-species transmission could occur by sharing an identical sandfly). It is anticipated that the share rate of sandfly is high within the same host type and low between different types, but we have yet to understand the plausible quantitative value. Second, not necessarily directly surveying these characteristics through entomologic investigations, but also one could examine genetic or molecular characteristics of Leishmania spp. For instance, one could implement phylogenetic analysis to accurately capture the route of transmission and evolution, thereby permitting us to track the inter-specific transmission in an explicit manner [18][19][20][21][22][23][24][36][37][38]. Such analysis could also shed light on our assumption of SIS model (e.g. if there is immune reaction through frequent re-infections and if there is an indication of evolution of the pathogen) and analysis of genetic data does not force us to adopt a stationary state assumption [18][19][20][21][22][23][24].
Enumerating limitations, we have noted that the following points are regarded as the weakness of the present study: (1) limited sample size, (2) unknown mixing pattern, (3) SIS-type assumption and (4) stationary state assumptions. In addition, it is unlikely, but there is a possibility that some important animals have not been covered by empirical studies and we have missed their contribution to the transmission dynamics. Thus, not only the number of samples and genetic data, but an ecological survey with broader research perspectives would be required in the future.
Despite these limitations, we stress out that our study has shown that animal species other than canine play important role in maintaining the transmission of L. infantum, which is different from the conventional wisdom that tended to focus on canine only. Future research studies should be conducted to clarify the extent of inter-specific interactions. Developing an indirect estimation method of such mixing pattern is one of our ongoing studies.