Evaluation of the Worldwide Occurrence of Rabies in Dogs and Cats Using a Simple and Homogenous Framework for Quantitative Risk Assessments of Rabies Reintroduction in Disease-Free Areas through Pet Movements

Dog and cat rabies cases imported from rabies enzootic countries represent a major threat for areas that have acquired rabies-free status and quantitative risk analyses (QRAs) are developed in order to assess this risk of rabies reintroduction through dog and cat movements. Herein we describe a framework to evaluate dog and cat rabies incidence levels in exporting countries along with the associated uncertainty for such QRAs. For enzootic dog rabies areas (EDRAs), we extended and adapted a previously published method to specify the relationship between dog rabies vaccination coverage and canine rabies incidence; the relationship between dog and cat rabies incidences; and then to predict annual dog and cat rabies incidences. In non-enzootic dog rabies areas (nEDRAs), we provided annual incidence based on declared dog and cat rabies cases. For EDRAs, we predicted an annual incidence potentially greater than 1.5% in dogs and about ten times lower in cats with a high burden in Africa and Asia but much lower in Latin America. In nEDRAs, the occurrence of rabies was lower and of similar magnitude in dogs and cats. However, wildlife could still potentially infect dogs and cats through spillover events. This framework can directly be incorporated in QRAs of rabies reintroduction.


Introduction
Rabies is a major and widespread zoonosis with a case-fatality rate of 100%, that causes approximately 60,000 human deaths each year [1]. In some areas such as Western Europe, Oceania or Japan the enzootic circulation of rabies (associated with Rabies Virus RABV) in domestic and wild animal populations has been halted, thus preventing human exposures. Nonetheless, rabies risk persists at low levels in these areas mostly because of (re)importations of rabies-infected animals, especially dogs and cats, from rabies enzootic areas [2]. In this context, it is crucial to assess the probability of rabies reintroduction through dog and cat movements in order to provide a deeper understanding of the processes responsible for the risk persistence. Moreover, models built to assess the probability of rabies reintroduction can be used to evaluate efficacy of risk mitigation measures or to test alternative Vet. Sci. 2020, 7, 207; doi:10.3390/vetsci7040207 www.mdpi.com/journal/vetsci scenarios of risk management. As a consequence, such risk assessments may be increasingly important in a context of globalisation, with more owners likely travelling with their pets [3][4][5]. An increasing number of countries with rabies-free status may also be interested in such analysis to define rabies risk management policies. Both quantitative and semi-quantitative risk assessment analyses have been implemented for rabies reintroduction in rabies-free areas [6][7][8][9][10][11][12][13][14][15][16][17].
In these quantitative risk assessments (QRAs), one crucial step has been and still is to provide rabies incidence in the animal population of interest (i.e., dogs and/or cats) in the exporting area(s). Incidence corresponds to the number of new cases of the disease of interest during a certain period of time (e.g., one year for annual incidence) divided by the number of individuals in the population [18]. The uncertainty is also key to depict this parameter accurately and can be specified by the variance of the parameter's distribution [19]. The incidence will then define the risk level represented by one given animal species in one given area (or travelling to this area) [18], before considering implementation of risk mitigation measures (e.g., vaccination, serologic testing, border control). Providing accurate rabies incidence distributions can be challenging due to the scarcity, heterogeneity and variable reliability of available data to produce parameter distributions for QRA models. Various methodologies have been used to quantify dog and cat rabies incidences for the different QRAs of rabies reintroduction. Some rely on reported dog and cat rabies case counts by countries to official institutions (e.g., World Animal Health Organisation (OIE)) and rabies incidence is then modelled with a Gamma distribution to represent uncertainty of a mean number of events, occurring as a Poisson process, per unit time [8,[11][12][13][14]16]. This method is convenient especially when dealing with worldwide dog and cat movements since these data are available for most of the countries. It is probably a reliable assumption for areas where rabies cases occur sporadically and where surveillance systems perform well as stated by the Great Britain Advisory Group on Quarantine [20]. Nonetheless, in some areas, mostly in enzootic dog rabies areas (EDRAs), there is an underreporting of animal rabies cases which highly depends on the considered country or region and its rabies surveillance system [21]. For example, an active surveillance programme in Kenya detected >70 times more rabid dogs than the existing passive surveillance system [22]. It has thus been proposed to use the maximum annual incidence over a several years period (e.g., 3 or 4 years) to maximise the risk [8,13]. However, the use of a maximum value in this context may not compensate for the level of underreporting. Others proposed to multiply these surveillance incidence data by a coefficient according to discrepancies observed between active and passive surveillance [15]. The use of such method to account for this underreporting bias is possible if animals included in the risk assessment of rabies reintroduction only come from a single country or area since this type of coefficient is adapted to a given surveillance system. Another drawback of using rabies case counts is the need to estimate dog and cat population sizes to obtain incidence. Dog and cat population sizes are difficult to estimate especially for areas with little information and can thus add more uncertainty in QRA models. Thus, despite being common in QRAs of rabies introduction, the use of declared rabies cases and animal population sizes could led to biased rabies incidences. It has also been proposed to use other data sources such as rabies incidences (i.e., number of new rabies cases divided by the population size of the investigated area) reported in scientific literature and as the result of specific investigations, often in a context of reinforced surveillance [10,20]. These data sources can then be used to produce non-parametric distributions (e.g., Uniform or Triangular distributions) directly from available rabies incidences [10] or by working from these reported incidence values completed with subject matter expert opinions [20]. It is noteworthy that such data sources avoid the need to use population size estimates.
In order to correctly evaluate rabies incidences, it is important to recall that rabies in dog and cat populations can have different epidemiological profiles, with a simple distinction that can be made between EDRAs and non-enzootic dog rabies areas (nEDRAs). In EDRAs, mainly in Africa and Asia, dogs are considered as the primary reservoir of rabies with high incidence levels. These countries have also the highest burden of human rabies, since dogs account for more than 95% of human rabies transmissions [1]. In these areas, research on the epidemiology of rabies in cats is less frequent, probably because of less significant public health implications since cats do not substantially contribute to human rabies exposures [1]. In EDRAs, cat rabies is nonetheless suspected to be closely related to dog rabies and the occurrence of cat rabies cases is assumed to be driven by the exposition to the dog reservoir [23][24][25]. In nEDRAs, such as North America, Europe, Oceania, some Asian Countries (e.g., Japan, Korea) or Southern Cone (e.g., Argentina, Chile, Uruguay and Paraguay), dogs do not act as the primary reservoir of rabies (as the result of control measures such as mass vaccination campaigns). However, sporadic cases of dog and cat rabies can still occur through spillover from wildlife (mainly from mesocarnivores or bats), if a wildlife rabies reservoir is present (e.g., North America, Eastern Europe, Southern Cone), or following the importation of dogs and cats from EDRAs as stated before [2,[26][27][28][29][30][31].
Our objective was to develop a framework to provide annual dog and cat rabies incidences and their associated uncertainty in the different parts of the world. This framework also aimed at including the distinct rabies epidemiological profile of each area, classified as EDRAs and nEDRAs. This work was primarily designed to provide dog and cat rabies incidence distributions at the country or group of countries level for QRA models of rabies reintroduction in disease-free areas through pet movements. It could thus strengthen these models with the support of a framework taking into account the available state of knowledge about dog and cat rabies incidences worldwide and avoiding the use of poorly reliable sources of data when possible.

Framework Overview and Data Sources
To provide annual dog and cat rabies incidence values with their uncertainty for different areas of the world we distinguished EDRAs and nEDRAS given the two different underlying epidemiological contexts presented above. In EDRAs we developed prediction models based on a methodology initially proposed by Hampson et al. in 2015 [1] that was adapted and extended (especially for cat rabies incidence). To fit the two prediction models presented below, we conducted an extensive literature search to identify records providing dog and cat rabies incidence values. We searched the PubMed database using key words "rabies" in the title or the abstract AND "dog(s)" OR "canine" in the title or abstract AND "incidence" OR "prevalence". Scientific articles published before November 2019 were investigated. The search results and the selection process are presented in Figure 1. probably because of less significant public health implications since cats do not substantially contribute to human rabies exposures [1]. In EDRAs, cat rabies is nonetheless suspected to be closely related to dog rabies and the occurrence of cat rabies cases is assumed to be driven by the exposition to the dog reservoir [23][24][25]. In nEDRAs, such as North America, Europe, Oceania, some Asian Countries (e.g., Japan, Korea) or Southern Cone (e.g., Argentina, Chile, Uruguay and Paraguay), dogs do not act as the primary reservoir of rabies (as the result of control measures such as mass vaccination campaigns). However, sporadic cases of dog and cat rabies can still occur through spillover from wildlife (mainly from mesocarnivores or bats), if a wildlife rabies reservoir is present (e.g., North America, Eastern Europe, Southern Cone), or following the importation of dogs and cats from EDRAs as stated before [2,[26][27][28][29][30][31].
Our objective was to develop a framework to provide annual dog and cat rabies incidences and their associated uncertainty in the different parts of the world. This framework also aimed at including the distinct rabies epidemiological profile of each area, classified as EDRAs and nEDRAs. This work was primarily designed to provide dog and cat rabies incidence distributions at the country or group of countries level for QRA models of rabies reintroduction in disease-free areas through pet movements. It could thus strengthen these models with the support of a framework taking into account the available state of knowledge about dog and cat rabies incidences worldwide and avoiding the use of poorly reliable sources of data when possible.

Framework Overview and Data Sources
To provide annual dog and cat rabies incidence values with their uncertainty for different areas of the world we distinguished EDRAs and nEDRAS given the two different underlying epidemiological contexts presented above. In EDRAs we developed prediction models based on a methodology initially proposed by Hampson et al. in 2015 [1] that was adapted and extended (especially for cat rabies incidence). To fit the two prediction models presented below, we conducted an extensive literature search to identify records providing dog and cat rabies incidence values. We searched the PubMed database using key words "rabies" in the title or the abstract AND "dog(s)" OR "canine" in the title or abstract AND "incidence" OR "prevalence". Scientific articles published before November 2019 were investigated. The search results and the selection process are presented in Figure 1. Literature search and records' selection process to fit dog and cat rabies incidence prediction models in enzootic dog rabies areas. Model (1) and model (2) correspond to the annual dog rabies incidence and annual cat rabies incidence predictions models respectively. Literature search and records' selection process to fit dog and cat rabies incidence prediction models in enzootic dog rabies areas. Model (1) and model (2) correspond to the annual dog rabies incidence and annual cat rabies incidence predictions models respectively. We identified 18 publications (14 additional publications for the dog rabies incidence model [22,[32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50], three publications for the cat rabies incidence model and one publication for both models [34,[51][52][53][54][55]) that were used to complete the initial dataset provided by Hampson et al. in 2015 [1].
For nEDRAs, data-sources included declared dog and cat rabies case counts to official institutions (national or supra-national; e.g., OIE, PANAFTOSA-PAHO/WHO, Rabies-Bulletin-Europe) or publications reporting and compiling these data. Dog and cat population size estimates, necessary to compute incidence values when working with rabies case counts, were based on human:dog and human:cat ratios and were searched in the scientific and grey literature. Finally, human population sizes were extracted from the 2018 World Bank dataset [56].

•
Model for dog rabies incidence prediction In order to predict annual dog rabies incidence in different parts of the world, we used the functional relationship proposed by Hampson et al. in 2015 [1] which assumes that, in a given dog population, annual dog rabies incidence Idog depends on dog rabies vaccination coverage VC: (1) We estimated parameters Idog max , the annual dog rabies incidence with non-existent vaccination coverage, and S1, a shape parameter, using annual dog rabies incidence time series available for Latin America and the associated dog rabies vaccination coverage previously used by Hampson et al. in 2015 [1] (48 observations). We supplemented it with data gathered from the 15 records identified in our literature search and reported annual dog rabies incidence and associated dog vaccination coverage in a context of reinforced surveillance in limited areas. Using these records, we added 30 observations to the initial dataset. This extended dataset reflected the various epidemiological contexts of rabies, most specifically with the data from studies carried out in Asia or Africa (see Table A1 for the extension of the initial dataset and for a list of references).

•
Model for cat rabies incidence prediction In EDRAs, we assumed that the annual incidence of cat rabies depends on the annual incidence of dog rabies, since dogs act as the primary reservoir [23][24][25]. Similarly to the annual dog rabies incidence model and according to the observed linear relationship between log-transformed raw data on annual dog and cat rabies incidences, we defined the functional relationship between cat and dog rabies incidences, Icat and Idog respectively, as follows: Parameters F, a multiplicative factor, and S2, a shape parameter, were estimated using reports that concurrently provided data on annual dog and cat rabies incidences for the same area identified in our literature search (detailed data and references are presented in Table A2). This model assumes equivalent rabies surveillance pressure on dogs and cats (i.e., equal probability of detection of rabies cases in the two species) since both cats and dogs are companion animals and no data in favour of different probabilities of detection were available.
The parameters of models (1) and (2) were estimated by maximum likelihood using the Broyden-Fletcher-Goldfarb-Shanno algorithm [57]. As previously stated by Hampson et al. [1], we assumed that residuals were Gamma-distributed. In other words, for each value of dog vaccination coverage in the model (1) and for each value of annual dog rabies incidence in the model (2), uncertainty was modelled by a Gamma distribution with specific shape and scale parameters: Γ(k 1 , (1) and Γ(k 2 , Icat k 2 ) for model (2); with k 1 and k 2 , the shape parameters, being estimated along with other parameters by maximum likelihood.

•
Prediction of dog and cat rabies incidences Annual dog rabies incidence was predicted using model (1) and dog vaccination coverage data from Hampson et al. in 2015 [1] provided for country or cluster of countries and obtained by expert elicitations. A "cluster" refers to countries grouped together due to their proximity, similar socioeconomic status and epidemiological rabies situation as previously described [1], adapted to factor in recent data and expert opinions. The 23 clusters (of which 16 are EDRAs) used here are presented in Figure 2. Annual cat rabies incidence was then predicted using dog rabies incidence predictions and model (2) both at country and cluster levels.
Vet. Sci. 2020, 7, x FOR PEER REVIEW 5 of 25 Annual dog rabies incidence was predicted using model (1) and dog vaccination coverage data from Hampson et al. in 2015 [1] provided for country or cluster of countries and obtained by expert elicitations. A "cluster" refers to countries grouped together due to their proximity, similar socioeconomic status and epidemiological rabies situation as previously described [1], adapted to factor in recent data and expert opinions. The 23 clusters (of which 16 are EDRAs) used here are presented in Figure 2. Annual cat rabies incidence was then predicted using dog rabies incidence predictions and model (2) both at country and cluster levels. Uncertainty about the incidence predictions (i.e., distributions of and ) was assessed for each cluster with the following iterative process which is also summarised in Figure 3: (1) Random selection of the dog rabies vaccination coverage using a PERT distribution (defined using minimum, modal and maximum values provided in [1]) for the chosen cluster; (2) Determination of the Gamma distribution parameters for dog rabies incidence using the selected level of vaccination coverage. The shape parameter was provided by the fitted model (1) and the scale parameter by the value of (considered as the distribution mean) divided by the shape parameter; (3) Random selection and storage of a value in the generated annual dog rabies incidence distribution; (4) Determination of the Gamma distribution parameters for cat rabies incidence through the same process as in (2) using the selected level of annual dog rabies incidence but using model (2); (5) Random selection and storage of a value in the generated annual cat rabies incidence distribution. This iterative process was repeated 10,000 times. Uncertainty about the incidence predictions (i.e., distributions of Idog and Icat) was assessed for each cluster with the following iterative process which is also summarised in Figure 3: (1) Random selection of the dog rabies vaccination coverage using a PERT distribution (defined using minimum, modal and maximum values provided in [1]) for the chosen cluster; (2) Determination of the Gamma distribution parameters for dog rabies incidence using the selected level of vaccination coverage. The shape parameter k 1 was provided by the fitted model (1) and the scale parameter by the value of Idog (considered as the distribution mean) divided by the shape parameter; (3) Random selection and storage of a value in the generated annual dog rabies incidence distribution; (4) Determination of the Gamma distribution parameters for cat rabies incidence through the same process as in (2) using the selected level of annual dog rabies incidence but using model (2); (5) Random selection and storage of a value in the generated annual cat rabies incidence distribution. This iterative process was repeated 10,000 times.

nEDRAs
In nEDRAs (corresponding to seven out 23 clusters, see Figure 2) the assumption that dog rabies incidence depends on dog rabies vaccination coverage is not sustainable since the enzootic cycle is assumed to be disrupted and since rabies cases appear sporadically through exposure to another reservoir (i.e., wildlife) or through importations of infected pets from EDRAs. Similarly, cat rabies was supposed to be independent on dog rabies. In this specific context we used surveillance data, assuming no under-reporting. New dog and cat rabies cases were postulated to occur under a Poisson process and as previously used [8,11,13,16], we defined a Gamma distribution to model annual rabies incidence (a mean number of events per unit time, i.e., year) as follows [19]: Annual rabies incidence ~ Γ(Number of rabies cases during , 1 ) Human population size Human:animal ratio where stands for the period (in years) during which rabies cases are enumerated. Uncertainty regarding human:animal ratios was taken into account. We used a PERT distribution if three or more ratios were available for the cluster using the mean of the ratios as modal value for this nonparametric distribution, the lowest value as minimum and the highest value as maximum. If only two ratios were available we defined a Uniform distribution using the values as extrema. Then 10,000 values were drawn from this distribution for each cluster of countries in nEDRAs. Analyses were performed in the R environment (version R-3.6.1) [58] using R studio software [59]. Packages "bbmle" [57] for maximum likelihood estimates and "mc2d" [60] for PERT distributions and simulations were used.

nEDRAs
In nEDRAs (corresponding to seven out 23 clusters, see Figure 2) the assumption that dog rabies incidence depends on dog rabies vaccination coverage is not sustainable since the enzootic cycle is assumed to be disrupted and since rabies cases appear sporadically through exposure to another reservoir (i.e., wildlife) or through importations of infected pets from EDRAs. Similarly, cat rabies was supposed to be independent on dog rabies. In this specific context we used surveillance data, assuming no under-reporting. New dog and cat rabies cases were postulated to occur under a Poisson process and as previously used [8,11,13,16], we defined a Gamma distribution to model annual rabies incidence (a mean number of events per unit time, i.e., year) as follows [19]: Annual rabies incidence ∼ Γ(Number of rabies cases during t , 1 t ) Human population size Human:animal ratio (3) where t stands for the period (in years) during which rabies cases are enumerated. Uncertainty regarding human:animal ratios was taken into account. We used a PERT distribution if three or more ratios were available for the cluster using the mean of the ratios as modal value for this non-parametric distribution, the lowest value as minimum and the highest value as maximum. If only two ratios were available we defined a Uniform distribution using the values as extrema. Then 10,000 values were drawn from this distribution for each cluster of countries in nEDRAs. Analyses were performed in the R environment (version R-3.6.1) [58] using R studio software [59]. Packages "bbmle" [57] for maximum likelihood estimates and "mc2d" [60] for PERT distributions and simulations were used.
Vet. Sci. 2020, 7, x FOR PEER REVIEW 7 of 25 Figure 4. Model best fits for annual dog and cat rabies incidences in enzootic dog rabies areas. Model best fit for incidence of annual dog rabies is presented in (a) and that for annual cat rabies incidence in (b). The points represent the observations and are coloured by continents in which studies were conducted. Black lines represent best fits. 95%CI: 95% confidence interval of incidences (2.5th-97.5th percentile interval of the Gamma distribution of the residuals).
Using these parameter estimates, mean annual dog and cat rabies incidences were then predicted for each country ( Figure 5 and see Table A3 for detailed results).

Figure 5.
Mean predicted annual dog and cat rabies incidences by country in enzootic dog rabies areas. Map (a) for mean annual dog rabies incidence was produced using model (1) and map (b) for mean annual cat rabies incidence was produced using model (2). If dog rabies vaccination coverage was not available for a specific country, a global cluster estimate was used to predict canine rabies incidence. . Model best fits for annual dog and cat rabies incidences in enzootic dog rabies areas. Model best fit for incidence of annual dog rabies is presented in (a) and that for annual cat rabies incidence in (b). The points represent the observations and are coloured by continents in which studies were conducted. Black lines represent best fits. 95%CI: 95% confidence interval of incidences (2.5th-97.5th percentile interval of the Gamma distribution of the residuals).
Using these parameter estimates, mean annual dog and cat rabies incidences were then predicted for each country ( Figure 5 and see Table A3 for detailed results).
Vet. Sci. 2020, 7, x FOR PEER REVIEW 7 of 25 Figure 4. Model best fits for annual dog and cat rabies incidences in enzootic dog rabies areas. Model best fit for incidence of annual dog rabies is presented in (a) and that for annual cat rabies incidence in (b). The points represent the observations and are coloured by continents in which studies were conducted. Black lines represent best fits. 95%CI: 95% confidence interval of incidences (2.5th-97.5th percentile interval of the Gamma distribution of the residuals).
Using these parameter estimates, mean annual dog and cat rabies incidences were then predicted for each country ( Figure 5 and see Table A3 for detailed results).

Figure 5.
Mean predicted annual dog and cat rabies incidences by country in enzootic dog rabies areas. Map (a) for mean annual dog rabies incidence was produced using model (1) and map (b) for mean annual cat rabies incidence was produced using model (2). If dog rabies vaccination coverage was not available for a specific country, a global cluster estimate was used to predict canine rabies incidence.

Figure 5.
Mean predicted annual dog and cat rabies incidences by country in enzootic dog rabies areas. Map (a) for mean annual dog rabies incidence was produced using model (1) and map (b) for mean annual cat rabies incidence was produced using model (2). If dog rabies vaccination coverage was not available for a specific country, a global cluster estimate was used to predict canine rabies incidence.
Annual dog and cat rabies incidences were then simulated at cluster level as shown in Figure 6a. The annual dog rabies incidence in EDRAs ranged from~0 up to >1.5% (Figure 6a). The highest dog rabies incidences were predicted in Africa and Asia, where inter-country variability was greater due to heterogeneity in dog rabies vaccination coverage (e.g., mean predicted annual dog rabies incidence: 152/100,000 in Thailand versus 571/100,000 in Nepal). The predicted annual cat rabies incidence in EDRAs followed, in the proposed model, the same trends as annual canine rabies incidence but were tenfold lower than for dogs, ranging from~0 up to~1% in some African and Asian countries ( Figures 5  and 6a). As a result of high vaccination coverage, dog and cat rabies occurrence was predicted to be low in Latin America.
Vet. Sci. 2020, 7, x FOR PEER REVIEW 8 of 25 Annual dog and cat rabies incidences were then simulated at cluster level as shown in Figure  6a. The annual dog rabies incidence in EDRAs ranged from ~0 up to >1.5% (Figure 6a). The highest dog rabies incidences were predicted in Africa and Asia, where inter-country variability was greater due to heterogeneity in dog rabies vaccination coverage (e.g., mean predicted annual dog rabies incidence: 152/100,000 in Thailand versus 571/100,000 in Nepal). The predicted annual cat rabies incidence in EDRAs followed, in the proposed model, the same trends as annual canine rabies incidence but were tenfold lower than for dogs, ranging from ~0 up to ~1‰ in some African and Asian countries (Figures 5 and 6a). As a result of high vaccination coverage, dog and cat rabies occurrence was predicted to be low in Latin America. Figure 6. Annual dog and cat rabies incidences by cluster of countries. Annual incidence in enzootic dog rabies areas (a) and in non-enzootic dog rabies areas (b). Boxes represent the interquartile range (IQR); vertical bars represent the median; upper whiskers extend to the 75th percentile and the highest value that is within the 75th percentile + 1.5 × IQR; lower whiskers extend to the 25th percentile and the lowest value that is within the 25th percentile − 1.5 × IQR. Dots represent the mean annual incidence of the clusters (completed with mean values in (b)). *: no reported cases. Figure 6. Annual dog and cat rabies incidences by cluster of countries. Annual incidence in enzootic dog rabies areas (a) and in non-enzootic dog rabies areas (b). Boxes represent the interquartile range (IQR); vertical bars represent the median; upper whiskers extend to the 75th percentile and the highest value that is within the 75th percentile + 1.5 × IQR; lower whiskers extend to the 25th percentile and the lowest value that is within the 25th percentile − 1.5 × IQR. Dots represent the mean annual incidence of the clusters (completed with mean values in (b)). *: no reported cases.

nEDRAs
Data gathered to compute incidences by using Equation (3) (i.e., rabies case counts and human:animal ratios) are summarised at the cluster level for the seven clusters considered as nEDRAs (Table 1). Declared dog and cat rabies case counts were collected from national or supranational organisms (OIE, PANAFTOSA-PAHO/WHO, Rabies-Bulletin-Europe, Centers for Disease Control and Prevention, Canadian Food Inspection Agency) for the 2013-2017 period. References used to define human:animal ratio distributions (PERT or Uniform) are also provided.
Results obtained for annual incidences in the seven clusters are reported in Figure 6b. Mean annual incidences are also presented by country using count data at the country level (Figure 7, see Table A4 for detailed results). Annual dog and cat rabies incidence values appeared to be low (<1/100,000 except for the Eastern Europe cluster where annual dog rabies incidence reached 3.51/100,000) even if wildlife can contribute to a number of spillover events in some clusters where a wildlife rabies reservoir exists (North America, Dog Rabies-Free Latin America, Eastern Europe). Rabies-Free Asia, Rabies Free African Islands and Rabies-Free Oceania clusters reported no imported cases contrary to Rabies-Free European cluster for the considered period, indicating heterogeneous risk levels associated with dog and cat translocations.

Discussion
We have proposed here a framework to define dogs and cats rabies incidence distributions as inputs for QRAs of rabies reintroduction in rabies-free areas through dog and cat movements. Provided parameters can directly be used to define QRA models. Moreover, such framework can allow updates if more data become available to fit the two prediction models for EDRAs, if other dog

Discussion
We have proposed here a framework to define dogs and cats rabies incidence distributions as inputs for QRAs of rabies reintroduction in rabies-free areas through dog and cat movements. Provided parameters can directly be used to define QRA models. Moreover, such framework can allow updates if more data become available to fit the two prediction models for EDRAs, if other dog rabies vaccination coverage are applied or if other case counts (e.g., over other periods) are used to define Gamma distributions for nEDRAs. It is also important to highlight that the distinction between EDRAs and nEDRAs in this framework is a simplification of a more complex reality. For example in some nEDRAs rabies transmission chains between dogs could occur without persisting over time (e.g., in Eastern Europe where there are high levels of spillover from wildlife). Conversely, it is also possible that in EDRAs rabies transmission chains start to be disrupted because of high vaccination coverage, making the classification complex (e.g., In Latin America countries still reporting dog rabies, the incidence is low and transmission in dog populations is believed to be stopped soon [76]).

Model and Results for EDRAs
The method used for EDRAs is well suited for homogeneously evaluating annual dog and cat rabies incidences as it ignores discrepancies between countries in animal rabies surveillance and case reporting when using national surveillance reports [21]. The strength of our approach also lies in its ability to predict uncertainty ( Figure 6), which is critical when performing quantitative risk analyses [19]. "Uncertainty" here is used in its broadest meaning since it accounts for both biological variability and limited knowledge and was particularly high in this context when incidence levels were high. However, bringing together studies carried out in different regions of the world illustrated various epidemiological contexts. Limited knowledge appeared to have an important impact when evaluating cat rabies occurrence for high dog rabies incidences since little information was available, leading to high-leverage data (Figure 4). It is also important to consider that this method, for the purpose of rabies reintroduction QRAs, aims to provide the global occurrence of rabies in dogs and cats at country or cluster level, but does not account for the spatial heterogeneity of rabies occurrence within a country or incidence time variations [77][78][79]. Spatial variability (e.g., as observed in Latin America with sometimes very focal areas where dog rabies is present [31]) and time variations are probably captured by the high level of uncertainty at low dog vaccination coverage for annual dog rabies incidence and, as a consequence, at high annual dog rabies incidences for annual cat rabies incidences. Also, animal subpopulation specificities (e.g., indoor pets versus free-roaming animals) could not be included in this model, whose purpose was to provide global incidence values at the country or cluster of countries level. Finally, in both models (1) and (2), the annual incidence of dog and cat rabies depends on only one variable, dog rabies vaccination coverage. To account for the effect of other covariates, more complex models need to be developed. Dog density and population turnover have been reported to be major predictors [77,79,80], but they could not be included in this framework since they were not available in most references. Such limitations were already identified by Hampson et al. [1]. Moreover, vaccination coverage in dog populations are sometimes difficult to estimate. To take this into account, we included uncertainty for this parameter when predicting incidences (by using a PERT distribution). However an overestimation of vaccination coverage would induce a lower rabies incidences in dogs and cats, and vice versa. Nonetheless, this method suits well the needs for quantitative risk analysis providing a homogenous way to provide annual rabies incidences for dogs and cats in exporting areas considered as EDRAs. Moreover, uncertainty associated with these values does not rely on non-parametric distributions arbitrarily chosen but aims at providing the actual level of knowledge and the biological variability observed in the field in an objective and standardised way. When compared to values used in risk analyses, we found similar results to QRAs also using similar data sources (i.e., incidence data obtained through investigations in a limited area with reinforced surveillance). For example, Hudson et al. in 2017 [10] used an Uniform distribution of 0.0001-0.05 (mean = 0.025) to model annual dog rabies incidence (i.e., annual probability of infection) in Indonesia which cover range of values similar to our results for Indonesia ( Figure 6-mean: 3.38 × 10 −3 ; 5th-95th percentiles: 2.07 × 10 −5 -1.27 × 10 −2 ) but with a lower mean. However, another report using similar sources of data, produced non-parametric Triangular distributions (minimum, mode, maximum) with lower variance and with a lower range of values: for example they used a Triangular distribution (1 × 10 −5 ; 1 × 10 −4 ; 1 × 10 −3 ) for Western Pacific while we obtained a mean dog rabies incidence of 5.66 × 10 −3 (5th-95th percentiles: 4.09 × 10 −5 -2.06 × 10 −2 ); or a Triangular distribution (1 × 10 −5 ; 1 × 10 −4 ; 1 × 10 −4 ) for Sub Saharan Africa while we obtained a mean of 5.68 × 10 −3 (5th-95th percentiles: 3.93 × 10 −5 -1.99 × 10 −2 ) for Southern Africa. One explanation is that the authors decided to lower incidence level for the purpose of the QRAs assuming that imported dogs came from lower-risk areas within these EDRAs [20]. For QRAs using declared rabies counts and declared animal populations sizes in EDRAs, direct comparison is difficult. Indeed, data sources are of a different kind and incidence values are rarely directly provided since it is an intermediate step in the risk calculation. Nonetheless one QRA presented instantaneous prevalence (IP) (i.e., annual incidence multiplied by mean incubation period (≈35 days)/365 days) and the ranking based on these IPs between EDRAs seems different from our results (e.g., Caribbean had the second highest IP in this other study) [13]. The reason for such differences is not traceable since source of variation could arise from rabies case counts and/or dog and cat population size estimates.

Model and Results for nEDRAs
In nEDRAs, where cat rabies is not considered to be dependent on dog rabies, since cat rabies occurs through exposure to infected wildlife or the importation of infected pets (and not through exposure to infected dogs), the incidence in cats may be higher than the incidence in dogs. This implies that, in such contexts, the occurrence of dog and cat rabies is mainly related to the level of exposure to wildlife and that cats, especially in countries with a high Human Development Index (HDI) (e.g., North America), certainly roam freely more often than dogs [81]. According to this hypothesis, it was consistent to observe a similar annual dog and cat rabies incidence in countries with a lower HDI (e.g., Eastern Europe), where there are probably more free-roaming dogs [82]. For these incidence levels to be valid, we assumed that dog and cat rabies cases were well detected and reported in this context. This is a reasonable assumption since in nEDRAs, which often have a high HDI, dog and cat rabies occurs only sporadically because most dogs and cats are owned and surveillance systems are efficient. However, in case of animal rabies underreporting, we expect incidences to be underestimated with our framework. It is also worth noting that biased population sizes would affect the validity of incidence values. In order to limit the impact of this bias, we included uncertainty on human:animal ratios. The most common bias in the estimation of animal population sizes occurs when some categories of animals are ignored, e.g., non-owned dogs and cats, and leads to an underestimation of population sizes. Such underestimation would conduct to overestimated incidences. Nonetheless, this bias would have minor impact since incidence values were already very low in nEDRAs. Moreover, the use of this method to compute rabies incidences concerned a minority of clusters of countries (seven out of 23). When considering rabies-free areas with no wildlife rabies reservoir, very low residual rabies incidence (Rabies Free Europe) or absence of rabies occurrence (Rabies Free Asia, Rabies Free African Islands and Rabies Free Oceania) were observed. This discrepancy could arise from the implementation of different mitigation measures and the associated level of compliance or from different kinds of pet flows in terms of volume and origin to these rabies-free clusters. For rabies incidence in nEDRAs we found similar values and similar level of uncertainty to the other reports mentioning IPs (with respect to coefficient to obtain IP from incidence values) [13]. This is consistent with the fact that the same methodology was used in the specific case of nEDRAs, with a Gamma distribution to model a mean number of events per unit time (assuming new rabies cases follow a Poisson process); even if we used human:animal ratios to estimate dog and cat population sizes (allowing to follow human population growth) instead of raw counts. Non parametric distributions using similar data sources for nEDRAs produced same range of values but with more uncertainty (i.e., greater distribution variances). For example, the Great Britain Advisory Group on Quarantine [20] used for dog and cat rabies incidence in North America a Triangular distribution (0; 1 × 10 −6 ; 1 × 10 −5 ) which includes the range of values that we found when combining dog and cat rabies cases: mean incidence of 2.13 × 10 −6 (5th-95th percentiles: 2.00 × 10 −6 -2.27 × 10 −6 ).

Conclusions
The framework developed was used to evaluate annual dog and cat rabies incidences, along with their corresponding uncertainties, in the different parts of the world for the purpose of QRAs of rabies reintroduction in disease-free areas through dog and cat movements. In EDRAs (especially in African and Asian countries), we predicted an annual incidence of up to 1.5% in dogs and 1% in cats with high uncertainty levels. In nEDRAs, residual dog and cat rabies occurrence was observed due to infected dog and cat importations and, more importantly, to exposure to wildlife in areas where it is infected (with values >2/100,000 in Eastern Europe for example). In these areas incidence levels appeared to be similar in dogs and cats, or higher in cats compared to dogs. This framework can directly be implemented in QRAs of rabies reintroduction and avoid the use of non-parametric distributions to model dog and cat rabies incidences. We thus provide an objective approach to define parameters' distributions for such models that take into account the current level of knowledge. More generally, such framework (and associated results) could also interest anyone wanting to obtain dog and cat rabies incidences at the country of or group of countries level.

Conflicts of Interest:
The authors declare no conflict of interest.   [50] Annual incidence values correspond to an annual number of new rabies cases divided by a population size. † : Reference specifically added to provide dog vaccination coverage. ‡ : Reference specifically added to provide dog population size in order to calculate incidence.  Annual incidence values correspond to an annual number of new rabies cases divided by a population size. † : Reference specifically added to provide animal population sizes in order to calculate incidence.    Annual incidence values correspond to an annual number of new rabies cases divided by a population size.