Study area and data collection

We developed simulations within a 64-km² ORV management zone in urban Hamilton, Ontario, Canada (Fig. 1). The city of Hamilton is located in southern Ontario at the western end of Lake Ontario and has a population of 536,917. Hamilton is part of a contiguous urban area surrounding the city of Toronto, which has a population of over 7.8 million people (Statistics Canada 2016). The study area was located just east of downtown Hamilton and included roughly 25% nonresidential and 75% residential (Fig. 1). Reincursion of raccoon rabies virus into Ontario was first detected in December 2015 in Hamilton, possibly caused by a long-distance translocation from southeastern New York State (Nadin-Davis et al. 2020).

Figure1

Study area in Hamilton, Ontario, Canada, for which simulations were run to assess the efficiency of local rabies vaccination strategies for raccoons (Procyon lotor) in an urban setting. Residential and nonresidential areas are indicated in dark blue and light yellow, respectively. The inset 8 × 8-km (64-km²) zone is where the simulations were conducted and where shading indicates residential and nonresidential areas.

In response to the detection of RRV, the Ontario Ministry of Natural Resources and Forestry (MNRF) started an ORV program using ONRAB^® baits in December 2015 to control and eliminate RRV in Hamilton (Lobo et al. 2018). Since 2015, ORV baits have been distributed by hand, bait stations, helicopter, and fixed-wing aircraft using various target densities and combinations of strategies (see Vaccination Strategy Design in  Supplementary Material Appendix 1 (JWD-D-23-00059R2_Bastille-Rousseau_REGULAR_Suppl_Final.docx)). Additionally, TVR was implemented starting in 2017 in targeted areas. Fixed-wing aerial baiting is limited to rural areas located outside of our study area. In the study area, helicopter baiting and bait stations were used to target nonresidential urban areas (e.g., large public parks, buffered areas along highways and ravines) whereas hand baiting and TVR were generally limited to residential urban areas. Raccoon density in the area was assumed to be at 18/km².

The Southern Ontario Land Resource Information System (SOLRIS) was used for land-cover classification. Land-cover categories “Transportation”, “Built-up Area Pervious,” and “Built-up Area Impervious” were reclassified as residential whereas other categories were reclassified as nonresidential (mostly forested areas). Pixels were therefore either residential or nonresidential. Individual raccoon movement data were not available for Hamilton, so auxiliary raccoon movement data were obtained from five raccoons (three female, two male) fitted with GPS collars in Toronto, Ontario, Canada. Collars were active during July–September 2010 and locations were subsampled every 30 minutes. Raccoons tracked in Toronto experienced a similar range of urban conditions to the Hamilton study area, representing residential areas and nonresidential areas. The raccoon project was approved by the York University Animal Care Committee (approval number 2010-25W).

Simulation structure

Our framework followed the approaches developed in McClure et al. (2022). Code for this framework is available as supporting information in McClure et al. (2022). Briefly, we analyzed raccoon GPS data to extract raccoon movement speed (step length) by fitting a gamma distribution and estimated resource selection coefficients to parameterize the movement simulations. These coefficients were estimated using a third-order resource selection function (RSF; Johnson 1980; Boyce et al. 2002) in which GPS locations were compared with 10,000 generated locations within each individual home range (estimated using a kernel density estimate). We assumed that the five tracked raccoons were representing a random sample of the overall population and used the mean and standard deviation in the observed parameters (coefficients from RSF and shape and scale parameters of the gamma distribution) to generate realistic movement parameters for 1,000 individuals. We then developed a spatially explicit model of raccoon movement informed by the SOLRIS land cover and our estimates of movement and habitat selection. We simulated raccoon movement for a 1-mo period (August) within our 64-km² study area and evaluated how different vaccination strategies might lead to different vaccine-based seroprevalences. We validated simulated raccoon movement by comparing home-range size of simulated raccoons with that of the five GPS-collared raccoons. Compared with McClure et al. (2022), we simulated only one type of raccoon movement, based on observed movement patterns based on the GPS-tracked individuals and observed density in the study areas (18/km² or 1,152 raccoons).

We simulated 560 vaccination strategy designs that combined the different ORV baiting and TVR vaccination strategies used by MNRF since 2015 (Table 1). Hand baiting and TVR were simulated only in residential areas, whereas bait stations and helicopter baiting were limited to nonresidential areas. No hand baiting was applied if TVR coverage was 100%.

Table 1

Summary of vaccination strategies and equivalencies among units (percentage, actual number, actual density) used in simulations run to assess the efficiency of local rabies vaccination strategies for raccoons (Procyon lotor) in Hamilton, Ontario, Canada. For each strategy, the reference level is presented first and indicates the intensity generally used. The actual number of baits used in the simulation is specified as well as the resulting actual density (baits/km²). This density is reflective of where the specific strategy is being deployed. Spatial vaccination data specific to the study area were used to estimate actual number and actual density and any mismatch between actual and target density. The area column indicates whether the strategy is applied to a residential area (R) or nonresidential area (NR).

We then overlaid simulated vaccination strategies on raccoon movement and estimated the proportion of raccoons that encountered and consumed a vaccine bait or that were trapped and vaccinated when in the same grid cell, resulting in seroconversion (rabies antibody response; Sobey et al. 2010; Brown et al. 2012; Gilbert et al. 2018). Oral bait consumption and seroconversion involved three processes: 1) colocation of a raccoon and bait within the same grid cell (30 × 30 m), 2) bait encounter and consumption given colocation that accounted for nontarget and conspecific bait competition and reduced bait availability over time, and 3) seroconversion given bait consumption. In the simulation, TVR was applied in a subsequent step in which a specific number of raccoons were randomly selected as trapped, vaccinated, and released back into the study area. As there was no disease in our simulated raccoon populations, we estimated the vaccine-induced population seroprevalence as the number of seroconverted individuals divided by the raccoon population across the study area. For each management strategy, we estimated the population seroprevalence achieved through vaccination and the associated monetary cost of the strategy. More details regarding each step of the simulations are provided in  Supplementary Material Appendix 1 (JWD-D-23-00059R2_Bastille-Rousseau_REGULAR_Suppl_Final.docx).

Statistical analysis

We used a linear mixed-effects model with a logit link to analyze the relative impacts of baiting strategy features on vaccine-based seroprevalence (McClure et al. 2022). All covariates were considered as categorical (i.e., “dummy” variables). Strategy was added as a random intercept because each strategy was repeated 250 times. Parameters (hand baiting, TVR, helicopter baiting, and bait-station baiting) were included as main effects. Observed vaccine-based seroprevalence was also summarized graphically to highlight the potential interactive effects of hand baiting and TVR in residential areas and helicopter baiting and bait stations in nonresidential areas. We evaluated the magnitude of the correlation between predicted vaccine-based seroprevalence and cost of the strategy (cost calculations for each strategy are presented in  Supplementary Material Appendix 2 (JWD-D-23-00059R2_Bastille-Rousseau_REGULAR_Suppl_Final.docx)). Using a separate linear mixed-effects model, we examined the effects of vaccination strategies on efficiency, where the efficiency of a strategy s was defined as

and relative efficiency was calculated as

where min(eff) and max(eff) represent the minimum or maximum values across all strategies. Independent variables included each vaccination strategy (hand baiting, TVR, helicopter baiting, and bait station baiting) as categorical variables and with strategy as a random intercept. We identified the three most efficient strategies across a range of target levels of vaccine-based seroprevalence, including 20%, 30%, 40%, 50%, and 60% seroprevalence. All simulations and statistical analyses were run on a workstation using R 4.0.3 (R Core Team 2022).

Raccoon movement and model parameterization

Raccoons varied in their responses to nonresidential areas, with three individuals showing strong avoidance of nonresidential areas (relative to residential areas) and two individuals weakly selecting for nonresidential areas ( Supplementary Material Table S1 and Appendix 1 (JWD-D-23-00059R2_Bastille-Rousseau_REGULAR_Suppl_Final.docx)). GPS-tracked raccoon home-range sizes ranged from 0.11 to 0.31 km². Simulated raccoon movement had similar properties regarding home-range size but with slightly larger variance, ranging from 0.05 to 1 km² ( Supplementary Material Fig. S1 and Appendix 1 (JWD-D-23-00059R2_Bastille-Rousseau_REGULAR_Suppl_Final.docx)).

Seroprevalence achieved by vaccination strategies

Across the 560 vaccination strategies, the predicted vaccine-based seroprevalence in raccoons averaged 33% (range 0–75%). Increasing the number of baits for hand baiting from a target density of 50% to 200% (Table 1) increased predicted seroprevalence from around 10% to 13% (Fig. 2). Trends were highly similar for helicopter baiting, with predicted increases in seroprevalence ranging from 10% to 12% for strategies with 50% and 200% of target bait densities. The number of bait stations used in the simulations had a smaller impact on seroprevalence than hand and helicopter baiting, and the number of baits per station very minimally impacted seroprevalence (Fig. 2). Combining TVR with other vaccination strategies led to an increase in seroprevalence across all strategies. The highest seroprevalence was observed when TVR was used to cover 100% of residential areas (Fig. 3). In nonresidential areas, combining helicopter baiting and bait stations led to higher seroprevalence, with an increased intensity of both strategies (Fig. 4). However, using these two techniques in nonresidential areas led to a predicted seroprevalence of only 16% in raccoons (Fig. 4).

Figure2

Effects of vaccination strategy on predicted seroprevalence for raccoons (Procyon lotor) in the urban setting of Hamilton, Ontario, Canada. Plots show predicted seroprevalence for different baiting designs relative to the absence of vaccination. For bait stations, one panel shows the effect of the number of stations while holding the number of baits per stations at 300 (100%) and the other panel is showing the effect of the number of baits per station while holding the number of bait stations at 30 (100%; the standard rates used by Ministry of Natural Resources and Forests in urban settings). For trap-vaccinate-release (TVR), the fraction of the study area covered (cvg) and the percentage of target animals vaccinated (vac) are provided. Density and number of baits associated with each strategy is presented in Table 1. Note the different y-axis scale for TVR.

Figure3

Interactive effects of vaccination strategies in residential areas on simulated seroprevalence for raccoons (Procyon lotor) in the urban setting of Hamilton, Ontario, Canada. Hand baiting was done at four different densities, whereas trap-vaccinate-release (TVR) was done at two different spatial coverages (cvg) and two different percentages of raccoons vaccinated (vac). No hand baiting was performed when TVR was done over the whole study area (100% coverage).

Figure4

Interactive effects of vaccination strategies in nonresidential areas on simulated seroprevalence for raccoons (Procyon lotor) in the urban setting of Hamilton, Ontario, Canada. Helicopter baiting was done at three different efforts, and the number of bait stations varied from 15 (50%) to 30 (100%) to 60 (200%). Given the very weak effect of the number of baits per station (Figure 2), the number of baits was set at 300 (100%) per station for the figure. Table 1 presents the equivalencies between percentage and actual numbers.

Efficiency of vaccination strategies

There was a strong correlation (R²=0.94) between predicted seroprevalence and the associated cost of management strategies (Fig. 5). Strategies in which 100% of the residential areas were treated with TVR had distinctly higher seroprevalence (>0.6) in raccoons, but TVR strategies also had a markedly higher associated cost. When standardizing the predicted seroprevalence in raccoons by strategy cost (relative efficiency), strategies in which lower seroprevalence was achieved were often more efficient (Fig. 6), especially for target population thresholds under 40%. Helicopter-only strategies tended to be the most efficient (Fig. 6 and Table 2), although the highest seroconversion in raccoons achieved by this strategy alone was <0.2 (Fig. 6). The TVR strategies were overall the least efficient, although high-density hand baiting had similar efficiency to TVR with 50% coverage (Fig. 6,  Supplementary Material Appendix 3 and Table S1 (JWD-D-23-00059R2_Bastille-Rousseau_REGULAR_Suppl_Final.docx)). A smaller number of bait stations was more efficient than strategies using more stations (Fig. 6,  Supplementary Material Appendix 3 and Table S1 (JWD-D-23-00059R2_Bastille-Rousseau_REGULAR_Suppl_Final.docx)). Even though TVR was inefficient when considering strategy costs, incorporating some level of TVR was necessary to achieve population seroprevalence higher than 30% (Table 2). Treating the entire residential area solely with TVR was the only method that led to seroprevalence higher than 60%. Combining an intermediate level of TVR with other vaccination methods was able to achieve seroprevalence between 30% and 60%.

Figure5

Relationship between averaged predicted rabies seroprevalence in raccoons and cost of the vaccination strategy (in Canadian dollars) based on simulations run to assess the efficiency of local rabies vaccination strategies for raccoons (Procyon lotor) in Hamilton, Ontario, Canada. The cluster of filled-in points represents strategies in which the entire residential area was treated with trap-vaccinate-release (TVR; 100% coverage) with regular effort (100% vaccination).

Figure6

Relationship between averaged predicted seroprevalence and relative efficiency for each vaccination strategy based on simulations run to assess the efficiency of local rabies vaccination strategies for raccoons (Procyon lotor) in Hamilton, Ontario, Canada. Strategies in which only one method was used are coded by symbol to highlight their overall efficiency. Each point on the graph represents one of the 560 vaccination strategies. The open circles represent mixtures of strategies that do not include trap-vaccinate-release (TVR).

Table 2

Top three lowest-cost rabies vaccination strategies to achieve seroprevalences of 20, 30, 40, 50, and 60% in a 64 km² study area of urban Hamilton, Ontario, Canada, based on simulations run to assess the efficiency of local rabies vaccination strategies for raccoons (Procyon lotor).

Broader applications and caveats

A critical question regarding our simulations is how well they could be transferred to other urban areas of RRV management outside the 64-km² focal area in Hamilton, Ontario, Canada. In terms of efficiency, we contend that our simulations should be directly applicable to other urban areas of southern Ontario and nearby areas in the US where raccoon density and behavior remain similar to those in southern Ontario. Indeed, raccoon movement was parameterized using radio-collar data from raccoons in Toronto, Ontario, Canada, and a simple land-cover classification system (residential and nonresidential). However, because different vaccination strategies were applied to residential areas (hand baiting and TVR) versus nonresidential areas (bait stations and helicopter baiting) and these strategies vary in their efficiency, it would be risky to directly extrapolate our results to areas with a markedly different ratio of residential to nonresidential areas (e.g., 75% to 25% in our case study). Our results may, however, inform what to expect if vaccinating an area with a larger fraction of nonresidential area. In such an instance, it is likely that the efficiency of the overall vaccination program would be increased, because a larger fraction of the area could be covered via helicopter baiting. This could mean that in an area where 50% of the area could be baited via helicopter, seroprevalence of up to 24% could be achieved via helicopter baiting alone. However, if higher seroprevalence is needed, TVR might also be required in the nonresidential areas.

Although the use of animal movement data and individual-based modeling is increasing in disease ecology (Dougherty et al. 2018; Holbrook et al. 2019), the analyses remain challenging, and caution is required when interpreting and generalizing the results. Our movement model integrates many aspects of raccoon behavior, such as confined home ranges, habitat-driven movement, avoidance of areas with high density of conspecifics, and individual variation in animal movement. However, this model remains a simplification of the raccoon movement and bait uptake process and is also based on some assumptions regarding the values of specific parameters (which were generally based on local agency knowledge). A prior study from an urban area of RRV management also reported generalist tendencies of raccoons across habitats (McClure et al. 2022). More complex habitat-based modeling of raccoon movement in urban areas is computationally unrealistic and would be unlikely to lead to major differences in the results. Likewise, given that some parameters were based on expert opinion, simulating bait and vaccine distributions may not perfectly reflect how baiting is generally done on the landscape, even if existing baiting data were used to parameterize the models. To reduce these potential differences, our baiting simulation not only replicated the same number of baits but also ensured that the same spatial coverage (number of pixels receiving baits) remained similar to what is achieved by MNRF. Previous work has shown that altering the spatial coverage could drastically change vaccine-based seroprevalence (McClure et al. 2022). Even if we use the best information available to parameterize bait uptake given bait discovery on the landscape and seroconversion given bait uptake, these parameters directly impact our estimated seroprevalence. In such cases, it may be safer to compare relative seroprevalence instead of absolute seroprevalence. Studies have also demonstrated a cumulative impact of baiting annually across years for reaching seroprevalence targets and case elimination (Davis et al. 2019; Acheson et al. 2023), whereas a single-season population simulation without disease may yield more limited predictions in the context of long-term RRV control programs.

Additional analyses may be warranted to further expand our understanding of how to optimize vaccination programs to control and eliminate RRV across a broader context. For example, replicating our analyses over different urban areas with varying configuration and composition (e.g., Tardy et al. 2018) of residential and nonresidential areas could expand our findings on efficiency to a broader range of spatial contexts. Relatedly, although the study focused on urban vaccination, eliminating rabies in rural areas may present its own set of challenges. Rural baiting is often done using a fixed-wing plane flying in a straight line, but it remains unclear if the baiting distribution in this type of landscape could be improved by better matching baiting flight lines with raccoon movement behavior and habitat selection at practical scales of implementation across broad landscapes. The current simulation framework could be easily adapted to different spatial scales and levels of urbanization. Additionally, the current simulation exercise focused only on the spatial aspect of vaccination, with little consideration to the temporal aspect. Combining the current approach with a temporal disease-transmission model might provide better information regarding where and at what intensity to vaccinate. More explicit consideration of interacting disease transmission and population demographic processes would allow us to assess the cumulative effect of vaccination over multiple seasons and overall efficacy. Lastly, expanding our framework for other mesocarnivores such as skunks and other species would benefit local efforts in Hamilton and could also expand the application to other areas of the world. Efforts to carry out ORV targeting free-roaming dogs has been a main priority to supplement parenteral mass-vaccination strategies in order to reach adequate fractions of the dog population for rabies control and elimination (Yale et al. 2022). Simulation-based frameworks could aid in the design of interventions targeting free-roaming dogs in regions with endemic canine rabies, along with other bioeconomic simulation tools (Anderson et al. 2019).

Although assessing efficiency (here defined as the cost associated per percentage increase in seroprevalence) of various vaccination approaches is the first step in informing on-the-ground interventions, ultimately the goal is to scale this information to cost-effectiveness (cost associated with a desired outcome of rabies management). Even if seroprevalence level and effectiveness are tightly linked, some factors can impact the scaling from one to the other. Considering disease transmission and the potential cumulative impact of repeating vaccination over several seasons might show that efficient approaches are less cost-effective in the long term (or vice versa) as well as inform on how to distribute the effort temporally (e.g., high vaccination rate in the first years vs. a lower rate over many years). Similarly, some approaches are more scalable spatially than others. For example, our results indicated TVR as being more efficient than hand baiting in residential areas, but it could be financially or logistically infeasible to perform TVR at a much larger extent than cities (e.g., at the state or provincial level). In mixed landscapes with large rural patches, incorporating other methods such as aerial baiting could be most efficient. However, for small urban areas as in our study, TVR combined with other approaches might be the most efficient strategy for contingency actions to contain local outbreaks.