Mechanisms for lyssavirus persistence in non-synanthropic bats in Europe: insights from a modeling study

Bats are natural reservoirs of the largest proportion of viral zoonoses among mammals, thus understanding the conditions for pathogen persistence in bats is essential to reduce human risk. Focusing on the European Bat Lyssavirus subtype 1 (EBLV-1), causing rabies disease, we develop a data-driven spatially explicit metapopulation model to investigate EBLV-1 persistence in Myotis myotis and Miniopterus schreibersii bat species in Catalonia. We find that persistence relies on host spatial structure through the migratory nature of M. schreibersii, on cross-species mixing with M. myotis, and on survival of infected animals followed by temporary immunity. The virus would not persist in the single colony of M. myotis. Our study provides for the first time epidemiological estimates for EBLV-1 progression in M. schreibersii. Our approach can be readily adapted to other zoonoses of public health concern where long-range migration and habitat sharing may play an important role.


Empirical data
Empirical data on demographic sizes of M. schreibersii and M. myotis colonies (Table S1) and on migration movements of M. schreibersii (Table S2) were obtained from a long-term survey of bat populations in natural colonies in Spain. The data we present hereafter are a continuation of the fieldwork and results already presented in Refs. 1,2 , so we refer the reader to those papers for details on data collection. In the higher resolution experimental scenario, we model explicitly the migration within different caves denoted as Summer refuges. In the following Table we provide the estimate for the associated migration flows.  (infected cells)), cellculture-adapted, EBLV-1 challenge virus (8918 FRA) was incubated with 3-fold dilutions of the sera to be titered. After incubation of the serum-virus mixtures, a suspension of BSR cells was added. Twenty-four hours later, the cell monolayer was acetone-fixed and labeled with a fluoresceinated anti-nucleocapsid antibody to detect the presence of non-neutralized virus (fluorescent foci). The optimal challenge dose (the dilution giving 80% infected cells for each virus production) is calculated. Further, titers are expressed as the arithmetic means of two independent repetitions. Samples were considered positive when the number of fluorescent foci was reduced by 50% at the 1∶27 dilution (starting dilution). This cutoff value is similar to that applied in other studies, see Refs. [2][3][4] for additional details.
While estimates are yielded from relatively small sample sizes, they are rather consistent across the 5 years of the study and with the estimates obtained from previous portions of the long-term survey 2-4 . Positive animals were found every year, except in 2014 when only one animal of the M. myotis species was sampled. Also, sampling rates are compatible with other field surveys in Europe 5-7 . Model 1 assumes the following transmission dynamics: infected bats become exposed ( ), then infectious ( ) and able to transmit the disease, and recover with rate . After an average immunity period !! , they become again susceptible. Disease progression is mathematically described by the following equations, written in deterministic continuous differential equations for sake of simplicity. Simulations are instead discrete and stochastic (see Section 3). Here we consider a frequency-dependent transmission.
where is the total bats population and the superscript p refers to the patches. In Can Palomeres we model the birth rate pulse b, and we add the cross species interaction through the factor !" !" indicates the number of infected of the other species and is the cross-species mixing interaction: The intrinsic ! ! value, the basic reproductive number in the absence of migration per each patch p, is obtained using the next generation matrix approach 8 : The force of infection for intra-species infection in a given patch p, at time t is: In Can Palomeres, where cross-species transmission between M. schreibersii and M. myotis leads to an additional contribution to the force of infection, !
In our numerical simulations, we fix a value for ! !" (and then explore a full range of values), compute the associated reproductive numbers ! ! in all other patches, and calculate the corresponding forces of infection.
The equations for the density-dependent transmission are as follows: In Can Palomeres with the cross-species mixing interaction:

Model 2
Model 2 assumes the following transmission dynamics: once infected, bats can experience a lethal infection with rate , becoming exposed ( ! ), then infectious ( ) and finally dying with rate . With a nonlethal infection (probability (1 − )), they become exposed ( ! ) without developing symptoms and then become permanently immune to the virus ( ). Here we consider a frequency-dependent transmission.
is the density dependent death rate, is the total bats population and is the carrying capacity. The carrying capacity is parameterized considering the parameterization used by Ref. 10 , and in our case it corresponds to = 7 • 10 ! for M. myotis and = 2.55 • 10 ! for M. schreibersii.
In Avenc Davì, Castanya, Summer refuges and Other caves the birth rate is = 0 so = (1 − / ). In Can Palomeres we model the birth rate pulse b, so > . As above, in Can Palomeres, the force of infection also considers the cross-species contribution. The intrinsic ! ! value, the basic reproductive number in the absence of migration per patch p, is obtained using the next generation matrix approach: The force of infection for intra-species transmission in a given patch p at time t is: As above, in Can Palomeres, the force of infection also considers the cross-species contribution,

Model 3
The transmission dynamics of model 3 is equivalent to the one described for model 2 with the addition of temporary immunity of average duration !! . Here we consider a frequency-dependent transmission.
Equations for Can Palomeres have in addition the birth rate pulse b and the cross species interaction, as for model 2. The reproductive number per patch ! ! in absence of migration is equal to Eq. (3).

Parameter values
We provide here a table listing the model parameters and their values.  *From empirical estimates, see Table S1 3. Simulation details

Stochastic and discrete integration of the disease dynamics
In each patch p for a compartment ! we extract a random variable from a multinomial distribution for each possible transition out of the compartment in the discrete time interval Δ . Such variable, ! ! , ! , determines the number of transitions from the compartment ! to a given compartment ! occurring in Δ . The change of population size of a compartment in the time interval because of disease dynamics is given by the sum over all random variables extracted, each for a specific transition: Here we consider the evolution of the susceptible compartment ! of model 1 as a concrete example. All possible transitions from this compartment are: to the exposed ! ! and to the natural death given by the demography.
The random variables for these transitions from ! are extracted from the multinomial distribution: with the transition probabilities: These three transitions cause a reduction of the size of that compartment ! . The increase is given by the transition from the recovered ! to the susceptible ! after the loss of immunity (here we do not consider Can Palomeres where an additional transition refers to the birth rate). This transition is modeled by a random number extracted from a binomial distribution: with probability ! ! →! ! = Δ , and ! being the number of recovered at time t in patch p. After extracting these numbers from the corresponding distributions, we can calculate the stochastic variation of the population size of compartment ! : Δ is equal to 1 day.

Stochastic migration of M. schreibersii
The yearly migration of M. schreibersii is defined by seasonal movements (Table S2) with an integration time scale of 1 day. The number of M. schreibersii in the compartment ! traveling from patch to patch ! in the discrete time interval Δ is an integer random variable extracted from a multinomial distribution. As a concrete example let us consider the spring migration out of Avenc Davì of bats in the !" compartment.
The two possible destinations are the compartment ! in Castanya and the compartment !" in Other caves. The random variables are extracted from: with migration probabilities: where: !" is the number of M. schreibersii in Avenc Davì in the compartment at time t; !"→!" and !"→! are the daily migration rates from Avenc Davì to Castanya and from Avenc Davì to Other caves, respectively (Table S2); !"→! and (1 − !"→! ) are the estimated percentages of bats that migrate to the two destinations, respectively (Table S2). After extracting these numbers from the corresponding distributions, we calculate the change in the !" population.

Maximum Likelihood
Likelihood was computed by assuming a binomial probability of testing positive, namely where !,! is the number of samples collected in the cave , at time ; !,! the ones tested positive; and !,! , ! the average density of recovered predicted by the model. !,! are considered to be statistically independent.

Additional results
In this section we provide additional numerical results mentioned in the main paper.

Sensitivity analysis
In this section we provide a set of numerical results where we evaluated the sensitivity of the model predictions to variations of model parameters or data estimates.