Model Formulation and Statistical Inference.

We explore four models of rabies dynamics within single bat colonies, representing alternative hypotheses based on field and experimental studies for the biology of rabies infection. Each model is described through modifications of the classic susceptible–exposed–infected–recovered paradigm (22, 23), and assumes that following exposure to VBRV, bats enter an exposed state (E) and subsequently either enter a noninfectious state with temporary immunity (T) with probability 1 − α, or enter a lethal infectious state with probability α (see SI Appendix for further discussion on immunizing exposures and infectious states). Bats that develop a lethal infection initially enter an infectious, but clinically silent state, I N , where bats have normal social and feeding behaviors before progressing to a rabid infectious state, I R , where they are aggressive or withdrawn (24). Model I made no further assumptions. Additional models captured recovery from an infectious state (model II), immune-boosting upon reexposure (model III), and lifelong immunity (model IV). Schematic representations for each model are provided in Fig. 2A and full details provided in the SI Appendix.

Fig. 2. (A) Schematic diagram of each model where the force of infection . Demographic processes are omitted for clarity. Filled states have seroconverted, and white states are seronegative. Model IV (lifelong immunity) is a base model that includes all black lines and states. Model I (temporary immunity and lethal infection) additionally includes the blue text, model II (nonlethal infection) includes red and blue text where 1 − ρ is the probability of recovery from infection, and model III (immune boosting) includes yellow and blue where c is the strength of immune boosting that depends on the force of infection. (B) Here 95% confidence intervals (CIs) for lethal infection probability versus the effect of immigration on the force of infection. As in A, model I is colored blue, model II is colored red, model III is colored yellow, and model IV is represented by black. Models I–III arrive at the same MLE for α and ϕ (blue star), and the MLE for model IV is indicated by a black star. (C) Change in likelihood from the MLE for each model across all R 0 values. Points above the horizontal gray line fall within the 95% CI (found using the likelihood ratio test), and the vertical line indicates . The MLE of R 0 and associated 95% CIs are also shown. (D) Log likelihood (log ℒ) and ΔAIC score (with zero associated with the best-fitting model) corresponding to the MLE for each model. The ΔAIC score accounts for differences in the number of estimated parameters.

In each model, transmission to susceptible individuals S occurs by bites from infectious bats (I N or I R ). Transmission is assumed to be frequency dependent because seroprevalence in the field was not significantly associated with the size of these bat colonies (19). Immigration of infectious bats also contributes to the force of infection λ at a constant rate ϕ so that where N is the total population size and β N and β R are the transmission rates from bats in the I N and I R states, respectively. Here, ϕ is a term that arises externally and represents infectious bats entering a colony and exposing susceptible bats before either leaving or dying. It can also capture interactions with infectious noncolony members during foraging.

Although many model parameters were inferred from challenge studies and knowledge of vampire bat life history (specifically, parameters describing the time bats spend in each state and all mortality parameters; SI Appendix, Table S2), the probability of developing lethal infection (α) versus temporary immunity following exposure and the impact of movement or interactions between colonies ϕ are unknown, but likely crucial determinants of transmission dynamics. Therefore, we used likelihood-based statistical inference methods (25) to confront the seroprevalence data from Peru with our transmission models to obtain maximum likelihood estimates (MLEs) and associated confidence bounds for α and ϕ when optimized over the transmission parameters, β N and β R (Materials and Methods).

We compared parameter estimates for each model (I–IV) by overlaying the 95% confidence intervals for α versus ϕ (Fig. 2B). Although model I (temporary immunity and lethal infection) is the best fitting model, the associated ΔAkaike Information Criterion (ΔAIC) indicates that models II and III (recovery from an infectious state and immune boosting) also create plausible representations of the data. It is important to note that the confidence intervals and the associated MLE for models I–III consistently point to a low probability of developing lethal infection—implying that most exposures are subclinical and immunizing—and frequent immigration of infectious bats. The role of immigration in VBRV dynamics is highlighted by the intrinsic basic reproductive ratio R 0 . In the absence of imports, the R 0 values corresponding to the MLE are all less than one (0.61, 0.34, and 0.51 for models I, II, and III, respectively), and confidence intervals associated with each model confirm that enzootic persistence relies on spatial structure and immigration of infectious individuals (Fig. 2C).

In contrast to models I–III, the model including lifelong immunity (model IV) suggests viral persistence across most feasible values of α with lower levels of immigration (Fig. 2B). However, the ΔAIC score of 7.38 indicates that this model provides a considerably less plausible explanation of the data (26). Further, although the presence of lifelong immunity would imply a different relationship between ϕ and α than that observed in models I–III, the conclusion that immigration is required for persistence within a colony is maintained; the estimated intrinsic R 0 remains substantially less than 1 (0.61) at the MLE.

Given the consistency of parameter estimates across models I–III, we further explored the dynamics of the most parsimonious model (model I): temporary immunity and lethal infection. To confirm the efficacy of our maximum likelihood estimation, we averaged over 1,000 simulations using the parameter values associated with the MLE across sample sites, which yields a mean global seroprevalence of ∼10.8%—nearly identical to the observed data (∼10.5%). Further, Fig. 3 displays a contour plot of the complete likelihood surface of the lethal infection probability α and the contribution of immigration to the force of infection ϕ. The MLE for the probability of lethal infection is surprisingly low (0.1) with immigration contributing notably to the force of infection. The simulated dynamics from this region of parameter space suggest that long-term rabies persistence is enabled by a high frequency of immunizing exposures, with the chain of transmission maintained by immigration events that lead to sporadic but lethal VBRV outbreaks (Fig. 3D). These dynamics are consistent with three observations of VBRV in the field. First, they capture the sporadic nature of local rabies outbreaks in livestock, as described by anecdotal reports from ranchers in Peru. Second, by again averaging over 1,000 simulations, the prevalence of active infection in vampire bats is found to be low (<1%), as has been reported in a variety of bat species in the Americas (27). Third, consistent with our field observations (Fig. 1), such dynamics allow a continuous presence of seropositive individuals with relatively stable colony sizes.

Fig. 3. Likelihood of VBRV dynamics in response to changes in α and ϕ shows that the most likely parameter combination is a low lethal infection probability (α) combined with a high effect of immigration on the force of infection (ϕ). (Center) Contour plot of the log likelihood for α versus ϕ for model 1. The reported likelihood values are maximized over β N and β R . The star indicates the MLE and the thick gray line represents the 95% CI. Above the contour plot is the profile likelihood for ϕ and to the right of the contour plot is the profile likelihood for α. The MLE corresponds to and with 95% CIs (10−3.51, 10−2.83) and (0, 0.29), respectively, noting 0 is not included in the CI for α. Parameters for sample simulations: (A) , , , ; (B) , , , ; (C) , , , ; and (D) , , , (the MLE). In A–D, shading indicates the proportion of seropositive individuals with temporary immunity (T, light red) and the proportion of seropositive individuals that are infectious ( , light blue). In contrast to D, simulations in A and C fail to consistently produce seropositive bats, while the dynamics in B drive the colony to extinction. At the MLE, ∼ 2 bats for every 10 susceptible bats present in a colony are expected to be exposed to VBRV by immigrant bats each year. Based on simulations this accounts for approximately half (53%) of all exposures, with the remaining exposures arising from bats within the colony. Department-specific parameters were consistent with the dynamics at the national scale (SI Appendix).

Thus, our estimation work indicates that the enzootic viral persistence suggested by serology in bats and recurrent spillover to livestock requires both movement of infectious bats among colonies and a high frequency of immunizing nonlethal exposures. It is important to note that these observations contrast sharply from the dynamics expected under qualitatively different and less likely regions of parameter space (Fig. 3 A–C). When both the infection probability and contribution of immigration to the force of infection (ϕ and α) are low (Fig. 3C), seropositive bats are observed only sporadically and at very low levels throughout the time series. In contrast, if α is increased with low immigration, rabies is present only for the duration of the infectious period—on average less than two weeks—and in contrast to the field data, seropositive individuals are rarely present (Fig. 3A). Under this scenario, infectious cases occasionally cause large epizootics that deplete the susceptible pool but such large-scale bat die-offs from rabies have not been reported. If both α and ϕ are high, seroprevalence falls within the range observed in the field, but the population approaches extinction owing to the high rate of lethal infections (Fig. 3B).