Notes:1Calibrated with distances in meters.2SI stands for Supporting Information.

Model calibration

The model was calibrated to the data available from the region of São Paulo acquired in 2017-2018 when vaccination of farms, but no roost control, was performed in the area. Two model parameters, the roost-to-roost transmission rate, βRR , and the roost-to-farm transmission rate, βRF , were fitted (see below); the remaining parameters were extracted from the literature (Table 1).
Model fitting was carried out using a regression-based conditional density Approximate Bayesian Computation algorithm, such as implemented in Prada et al. (2014), following Beaumont et al. (2002) and Lopes and Beaumont (2010). Briefly, summary statistics were calculated from the 2017-2018 data, and we ran a total of 7,800 simulations to calibrate the roost-to-roost and roost-to-farm transmission rates,\(\beta_{\text{RR}}\) and \(\beta_{\text{RF}}\), respectively. Due to the high number of nodes in the network, the fitting process was carried out in two phases. First, the roost-to-roost transmission rate parameter was calibrated to reach 1% prevalence across all roosts in the network by simulation of 100 of years of transmission between roosts. The spillover to farms was not considered in this phase, therefore, no intervention was allowed. Second, the roost-to-farm transmission rate parameter was fitted to generate 226 outbreaks in farms across a five-year period, as we calibrated the parameters using the data collected in 2017-2018 when only farm vaccination was performed, only this intervention was allowed in the second phase of calibration, see Model calibration in Supporting Information for more details.

Bat-rabies control scenarios

Using the calibrated model, we explored several VBR control scenarios, assessing the effect over the spread of VBR from a single introduction in a randomly selected roost. We considered three initial settings of suitability environments, depending on whether the single initial introduction was in either a high, middle, or low suitability environment (limited to roosts connected to at least five other roosts, to ensure simulations are not initiated in isolated locations). We define these three initial sets of roosts, so that the roosts with the suitability index, calculated through the ENM (Anderson, 2013; Owens et al., 2013; Soberón and Peterson, 2005), within upper decile (after excluding isolated locations) form a set of roosts in high suitability environment, roosts with the suitability index of +/-5% around the median form the middle, and within bottom decile form the low suitability environment sets of roosts (Supporting Information Figure S6).
In each initial setting, in response to VBR being detected in farms, we consider all combinations of two reactive interventions included in the model: a combination of roost control and farm vaccination, each intervention alone, or no intervention. Consequently, we simulate 12 different control scenarios (three initial settings of suitability environment with four different intervention strategies, Figure 3). This enables a comparison of the impact of environmental suitability on virus transmission and intervention effectiveness. With the model calibration, we selected the best posterior draws (107 selected), and we ran 50 simulations with each posterior, 5,350 simulations in total per control scenario.
We assessed two different outcomes: (i) the number of detected outbreaks in farms, and (ii) the distance of virus spread from initial infection in a roost to a farm in one year, for the different control scenarios. We determined whether there are any statistically significant differences between the means of the outcomes for different intervention strategies by the Welch’s one-way heteroscedastic F test, an alternative to ANOVA robust to the violation of variance homogeneity assumption, which we observed for both outcomes. The Welch’s test has one of the highest adjusted power among one-way tests for positively skewed data, which we observed for the numbers of outbreaks, and for approximately normally distributed data, that we observed for the distances (Dag et al., 2018). Since we do not confirm the hypotheses of equal means, we perform the Games-Howell post-hoc tests to recognize which pairs of intervention strategies significantly differ, assessed through the Holm-corrected p-values. Tests and visualization are performed using the ggstatsplot package of R (Patil, 2021).
Areas at persistently high risk of VBR transmission and spillover in the state of São Paulo after random introductions can be highlighted by mapping spillover events. We divided the state of São Paulo into squares of 3’ latitude times 3’ longitude (30km2). Spillover risk of farms was calculated as a proportion of (detected and undetected) infections among all simulations of a particular scenario.
Results
The average number of detected outbreaks in farms in a year, from a single introduction, is decreased significantly when an intervention strategy is implemented (being either cattle vaccination, roost control, or both), across all three suitability environments considered (Figure 4A-C). The maximal distances of virus spread from a single infection in a roost to a farm in one year, for the different intervention strategies are shown in Figure 4D-F.
The F statistics and the p -values of Welch’s F test are summarized for each comparison in Table S3 in Supporting Information, with all p -values close to zero, i.e. for both outcomes, across all three initial suitability settings, we do not confirm the hypothesis of equal means in the four intervention strategies. The Games-Howell post-hoc tests identify which pairs of intervention strategies significantly differ. The Holm-correctedp -values indicate that the outcomes for the combination of farm vaccination and roost control, versus roost control alone, are not significantly different, Figure 4A-F. Additionally, when infection starts in a low suitability environment, the number of outbreaks in farms do not significantly differ between the combination of farm vaccination and roost control, versus farm vaccination alone, Figure 4C.
The most ecologically suitable areas for bats, and thus where spread is likely to be higher, are concentrated in the east side of São Paulo state. The infection risk decreases dramatically with any intervention (whether it is farm vaccination, roost control, or both); the probability of an outbreak occurring in farms, after a single introduction, can be as high as 3.81% of the simulations ran without intervention and 1.02% of the simulations with the roost control, the most effective intervention strategy (Figure 5). High infection risk probabilities in farms were also observed in the middle and low suitability environments, which could be as high as 7.12% (Supporting Information Figure S7).
Discussion
The aim of this study was to explore the spatio-temporal dynamics of vampire-bat-driven rabies (VBR) in São Paulo, Brazil, and identify high-risk areas of spillover to cattle farms. This was achieved through the development of a novel stochastic network two-species metapopulation model. The model was used to explore the impact of current interventions, ring vaccination of farms and/or ring roost control (either bat culling or bat vaccination) around a positive farm. Our results suggest that either strategy can prevent substantial number of on-farm outbreaks, as well as significantly reduce the geographical spread of the virus. However, roost control alone or combined with farm vaccination in general leads to more significant control results than cattle vaccination alone. Interestingly, combination of both intervention did not provide a significant benefit comparable to roost control alone. We also found areas of consistently high infection risk in high roost suitability environments, and in middle and low suitability environments for bat roosting.
As possibly, the most diverse, abundant, and geographically dispersed vertebrate, bats are unique in their ability to fly, long lifespans, migratory patterns, and in hosting a diverse suite of pathogens including rabies virus (Calisher et al., 2006; Luis et al., 2013). Some of these factors contribute towards the efficacy of bats as zoonosis transmitters, but also towards the lack of data about pathogen circulation, in particular their high level of mobility and vast geographic ranges, as field data are often collected from a subset of a species geographic range over a small timescale (Benavides et al., 2016). While keeping the model relatively simple in terms of bat demography, we reproduce several important environmental drivers of disease transmission, such as elevation driving the explicit range of contact between roots and farms, male-driven transmission between bat roosts, flight distance and environmental suitability. This is key to generate useful risk maps that can support policy implementation. For example, Benavides et al. (2020) highlighted the challenge of applying bat vaccines across many roosts, which could be mitigated by focusing efforts on the areas estimated by the model to be at higher risk, which could in addition reduce cross-species exposure while reducing the impact on bat communities.
Bat culling remains a controversial approach to VBR control (Streicker et al, 2012). Alternatively, a spreadable vaccine may be administrated similar as the vampiricide. Laboratory and model results showed that the oral vaccination could be effective (Standing et al., 2017; Bakker et al., 2019). In the model we considered the implementation of a roost control, which can either represent bat culling or bat vaccination. Either way, it is modeled so that if a roost receives the intervention, all transmission ceases for at least one year. In the case of culling, the spread of the poison due to intensive grooming leads to an empty roost (Linhrat et al, 1972), which would likely be repopulated in the future, but this could potentially take a long time and has dangerous ecological implications. Furthermore, it was suggested that culling may increase recruitment of susceptible juveniles into the system, making the intervention ineffective or counterproductive, therefore, the efficacy simulated here is likely overestimated in case of culling (Streicker et al., 2012; Gentles et al., 2020). To study this in more detail, model would need to be modified to include within roost dynamics. Bat vaccination, being spread the same way, will lead to the entire roost population immune, arguably for at least one year. This type of control would not change the population structure within the roost, on the other hand, it will not reduce the impacts of bats as pests causing harm to animals by bat bites independently of rabies, including skin damage, anemia, loss of vision, loss of weight and productivity, and predisposition to other infection (Delpietro et al., 2021).
Nevertheless, cattle vaccination has also achieved considerable reduction in on-farm outbreaks and geographical spread of infection across the three initial suitability environments. Either way, a spatial mathematical model simulating the impact of these interventions, for example extending the one presented here, could be used before hand to evaluate the consequences of their introduction, and identify the most suitable locations to cover with the campaign for the successfully control, or even eradication, of the virus. As concluded by Blackwood et al. (2013), who developed several stochastic SEIR models examining viral persistence, bat population migration, and the effects of bat population culling; the mechanisms to reduce spillover via viral elimination, likely need to be spatially coordinated to be effective as we demonstrated here.
In the model we considered the minimum delay in the detection of outbreaks in the cattle farms to be 25 days, with a mean detection time around 75 days, and the most frequently observed delay (mode) of 30 days. We assumed a relatively over dispersed distribution to capture both the latency period and delay in detection. As the interventions simulated here are reactive, reducing the delay in detection could generate significant gains in reducing transmission. Alternatively, the farm or roost control could be administered in a prospective manner, for example focusing on high-risk areas. The challenge would be to justify the investment to stakeholders (whether it is the farmers paying for the cattle vaccine, or the government paying for either farm or roost control), when the risk might not be perceived.
We followed prior work made in the region (Rocha et al., 2020; Rocha and Dias, 2020), building a similar contact network as in Rocha and Dias (2020), with a consistent assumption of up to 10 km flight distance (Benavides at al., 2016), and dependence of bat foraging migration pattern on altitude (Rocha et al., 2020). How far within the 10 km distance the bats fly is determined by the number of individuals in the roost, since individuals may fly to more distant feeding sources and/or roosts to minimize competition with conspecifics (Kunz and Fenton, 2003; Rocha et al., 2020). We address these spatial interactions by utilizing a gravity model. In addition, we incorporate knowledge of favorable conditions for bats using the roost locations and an ecological niche model to capture the environmental suitability (Ecological niche model in Supporting Information). Our approach has however a number of limitations; the contact network is assumed to be time-invariant and we are examining outbreaks over a one year time-period from a single introduction. Assuming a unique infected roost at random (potentially in the middle of the region) as a starting point is unrealistic, however it allows us to better capture the expected spatial spread from a single point. The model focuses only on spillover to cattle, however, other animals are in risk (e.g. horses), and since rabies virus is zoonosis also spillover to humans occurs.
The reactive interventions depend on reports from the producers which is influenced by many socio-ecological factors, similarly adherence to intervention and thus vaccination of the animals when infection nearby is reported might be conditioned by various factors (Benavides et al., 2017). The behavior effects on intervention needs to be accounted for in model if we want more realistic predictions. Furthermore, the intervention strategies effective to reach programmatic goals needs to be evaluated in economic manner as the government and farmers financial sources are limited (Janoušková et al., 2022). For example, anemia from bat bites may reduce livestock productivity (Bakker et al., 2019), hence making a difference in bat culling compare to bat vaccination. The roost control might be more cost-efficient to the official service since a smaller number of locations should be visited and vaccine delivery (for example as a paste) is more straight-forward than cattle vaccination. The model presented here, does not evaluate the economic implications, therefore distinguish only between susceptible, exposed, and infected farms, ignoring how many number of animals and to which extent are affected by bites and/or infection. The cumulative losses due to bites if no culling is performed, and even deaths if no roost control is in place might markedly change the cost-effectiveness. Last, but not least, the trust, support and commitment of stakeholders and involved institutions is necessary to reach the expected results (Janoušková et al., 2022). For instance, vaccination of vampire bats without population reduction will be unacceptable to some stakeholders since uncontrolled bat depredation sustains exposures to non-rabies pathogens (Bakker et al., 2019). The stakeholder’s preferences have to be taken into account when assessing the sustainability of the interventions.
Conclusion
We have developed a novel stochastic network two-species metapopulation model, that captures transmission of VBR between bat roost, as well as spillover events to cattle farms. After exploring two alternative control strategies, namely reactive ring roost control (i.e. bat culling or bat vaccination) and reactive ring cattle farm vaccination, we found no large differences in their expected efficacy, however interventions in roosts were statistically significantly better in all settings considered across both outcomes (number of outbreaks and spatial spread from initial introduction). Such mathematical frameworks can prove useful to inform control interventions, particularly identifying high-risk areas where prospective vaccination, either in cattle or in bats, could take place. This will support ongoing programs, leading to more effective control. Nonetheless, to reach long-term strategies and sustainability that could move beyond control to potential local elimination and eradication, human behavior, for example, in context of interventions uptake and response to VBR infection in farms, needs to be incorporated in model to get more accurate predictions. In addition to assessing intervention strategies effectiveness and high-risk areas such as provided in this study, economic evaluation is essential before decision is made on interventions.
Acknowledgements
This work was funded by the University Global Partnership Network Research Collaboration Fund (UGPN)