Management efficacy in a metapopulation model of white‐nose syndrome

The fungal pathogen Pseudogymnoascus destructans (Pd) causes white‐nose syndrome (WNS), an emerging disease that affects North American bat populations during hibernation. Pd has rapidly spread throughout much of the continent, leading to mass mortality and threatening extinction in several bat species. While previous studies have proposed treatment methods, little is known about the impact of metapopulation dynamics on these interventions. We investigate how the movement of bats between populations could affect the success of five WNS control strategies by posing and analyzing a two‐population disease model. Our results demonstrate that vaccination will benefit from greater bat dispersal, but the effectiveness of treatments targeting fungal growth or disease progression can be expected to diminish. We confirm that successful control depends on the relative contributions of bat‐to‐bat and environment‐to‐bat contact to Pd transmission, and additionally find that the route of transmission can influence whether interpopulation exchange increases or decreases control efficacy.

mortality. A variety of potential control methods targeting the determinants of WNS are currently under development; however, no treatments have been tested on a broad scale and the disease continues to persist unabated. While conducting large-scale experiments on affected species may be costly and impractical, mathematical modeling provides a powerful way to simulate disease dynamics. Using mathematical models, proposed treatment strategies can be tested before laboratories and government agencies commit resources to trials in the wild.
Multiple modeling studies have explored various aspects of WNS disease dynamics via both continuous-time and discrete-time models (Erickson et al., 2014(Erickson et al., , 2016Kramer et al., 2019;Langwig et al., 2017;Maslo et al., 2017;O'Regan et al., 2015;Reynolds et al., 2015;Russell et al., 2015;Thogmartin et al., 2013). Additionally, a few other modeling studies have considered the efficacy of several proposed control methods (Cornwell et al., 2019;Hallam & McCracken, 2011;Meyer et al., 2016). To our knowledge, however, no previous research has investigated how the exchange of individuals between bat populations can alter the efficacy of WNS control strategies relative to a single-population setting. Yet, interpopulation exchange is an important management consideration because hibernating North American bat species are known to travel long distances (Davis & Hitchcock, 1965;Fenton, 1969;Humphrey & Cope, 1976;Kaleigh et al., 2013). Moreover, recent work has demonstrated that bat populations exhibit cryptic connections which bridge groups and species to pathogen dynamics, indicating that interconnectivity within and between populations is an even more significant factor in WNS dynamics than previously thought (Hoyt et al., 2018).
In this paper, we investigate the impact of bat metapopulation dynamics on the long-term effectiveness of recently proposed control strategies by posing and numerically analyzing a twopatch mathematical model. We test five control strategies: fungicide intervention, microclimate intervention, soil bacteria intervention, ultraviolet light intervention, and vaccination. We adapt our model from the one proposed in Meyer et al. (2016) by extending it to incorporate movement of bats between populations, updating it to reflect recent WNS research, and considering new intervention strategies. In particular, we focus on the dynamics arising from the existence of dispersal between two autonomous populations of bats and from the possibility for different controls in each of the two locations. Our approach constitutes the first comprehensive modeling study to analyze the metapopulation dynamics of WNS together with control strategies targeting multiple routes of Pd transmission.

| Model overview
In this study, we focus on Myotis lucifugus (the little brown bat). Among North American bat species affected by WNS, little brown bats face a particularly significant threat given that most of their colonies have declined by more than 90% since the introduction of Pd (Frick et al., 2015;Langwig et al., 2012). Once one of the most widely distributed bat species in the United States and Canada, the little brown bat is now predicted to become regionally extinct within the next decade (Frick et al., 2010;Langwig et al., 2012;Reynolds et al., 2015). For reasons that are yet unknown, however, some populations of little brown bats have stabilized and are even growing in size following the invasion of Pd (Langwig et al., 2012(Langwig et al., , 2017. Hence, the little brown bat is frequently used as a model organism for studying hibernating North American bat species and has been the subject of several past WNS modeling studies (Cornwell et al., 2019;Langwig et al., 2017;Meyer et al., 2016;Reynolds et al., 2015).
Little brown bats, along with other hibernating North American bat species, follow four distinct annual phases: swarming, hibernation, roosting, and roosting with birth ( Figure 1). Since Pd transmission dynamics remain qualitatively unchanged within each phase but differ drastically between phases, our model employs phase-specific equations adapted from the hybrid continuous/discrete model structure introduced in Meyer et al. (2016). The equations in each phase track four classes of individuals: the numbers of susceptible bats (S), exposed bats (E), infected bats (I ), and free-living Pd colony-forming units (CFUs) (P). We use the ode45 function in MATLAB (The MathWorks, Inc., version R2019a) to numerically integrate the within-phase equations. All simulations begin at the start of the swarming phase, and the initial conditions for every subsequent phase are taken to be the final population values of the preceding phase.
We define susceptible bats to be those that are susceptible to Pd conidia colonization but have not yet acquired the fungus. Some data suggest the existence of a period between Pd exposure and the manifestation of clinical WNS symptoms (Langwig et al., 2017;Lorch et al., 2011). Thus, we also include the exposed class for bats that carry Pd fungus but are still asymptomatic. We assume exposed bats are not infectious and hence do not shed Pd conidia into the environment or transmit Pd via bat-to-bat contact. Infected bats, on the other hand, are both symptomatic and infectious; they die at an increased rate, shed Pd spores into the environment, and transmit Pd to other bats via physical contact. Therefore, the total size of bat population j is N S Additionally, an infectious environmental reservoir of Pd can persist in cave and mine sediments regardless of the presence of bats (Lorch et al., , 2013Reynolds et al., 2015). Hence, we include a P j class to account for the free-living pathogen within the hibernaculum of population j.
F I G U R E 1 Schematic diagram illustrating the yearly temporal arrangement of all model phases and discrete reclassification events To determine how bat metapopulation dynamics affect persistence, we compare model dynamics arising in a two-population setting with baseline dynamics from a single-population setting. That is, we first evaluate each control strategy on its own within a single-population model without metapopulation dynamics and then extend our base model to incorporate movement between two populations, which we denote A and B. Accordingly, we use the term "metapopulation" to refer to the combined grouping of population A and population B. These two populations are simulated autonomously: each population occupies its own hibernaculum and thus has its own S, E, I , and P classes, and the two populations are not necessarily subjected to the same controls. To account for the exchange of bats between populations A and B, our model moves a proportion of the bats in each population to their respective classes in the other population immediately following the end of swarming.

| Swarming phase
The swarming phase lasts approximately 2 months from mid-August to mid-October (simulation days 1 through 61). During this time, each population congregates near a hibernaculum to mate and accumulate fat stores in preparation for hibernation (Fenton & Barclay, 1980). Studies suggest that the swarming phase does not play a significant role in the transmission of Pd; despite high bat-to-bat contact, normal bat immune system function and high body temperatures preclude infection at this time (Langwig et al., 2014). Hence, we assume that during swarming exposed bats remain in the exposed class and do not become infected. Since no infectious bats are present, susceptible bats can become exposed to Pd from the environment but not from bat-to-bat contact. Following previous modeling studies (Cornwell et al., 2019;Meyer et al., 2016;Reynolds et al., 2015), we assume that Pd continues to grow logistically on cave sediments throughout the year, even in the absence of spore shedding from infected bats. We likewise adhere to precedent by modeling both environment-driven and host-driven transmission of Pd as density-dependent (Cornwell et al., 2019;Meyer et al., 2016).
The swarming phase equations are given by: During the swarming phase the susceptible class (S) decreases at a combination of two rates: a density-dependent rate ϕ of pathogen transmission from environment-to-bat contact and the natural rate of bat mortality μ. The exposed class (E) increases at the densitydependent rate ϕ of pathogen transmission from environment-to-bat contact and decreases at the natural rate of bat mortality μ. The environmental Pd reservoir (P) grows logistically at the natural free-living Pd growth rate η with carrying capacity K Pd . We assume that no infected bats are present at the beginning of the simulation and that all infected bats either die or move to the exposed class before the beginning of the swarming phase in subsequent years. Therefore, the infected class (I) is absent during swarming.

| Dispersal reclassification
Little brown bats can travel long distances throughout the roosting and swarming phases (Davis & Hitchcock, 1965;Fenton, 1969;Humphrey & Cope, 1976;Kaleigh et al., 2013), but we assume that bat populations do not mix at all during the hibernation phase. Since in our model bat-to-bat transmission of Pd only occurs during hibernation, metapopulation dynamics during other phases would have no direct bearing on the force of infection. Thus, while bats from both hibernacula might intermingle during the entirety of the roosting and swarming phases, exchange between populations only becomes epidemiologically significant when bats do not return to their hibernacula of origin at the beginning of the next hibernation phase. Therefore, we account for all possible between-population movements throughout the active phases via a single annual dispersal event. After the end of swarming but before the beginning of hibernation, our model instantaneously moves a proportion σ of each bat population class (S E , , and I ) to the respective class in the other population.
The dispersal reclassification equation is given by: where S E I , , A

| Hibernation phase
The hibernation phase lasts approximately 7 months from mid-October to mid-May (simulation days 62 through 273). During hibernation, bats lower their basal metabolic rates and remain torpid, causing their body temperatures to drop to a mere 2-8°C (Bouma et al., 2010;Boyles & Willis, 2010). Due to their lowered immune system function and body temperatures, bats are most susceptible to Pd infection at this time (Verant et al., 2012). Pd infects bat epidermal tissue and can create distinctive lesions on the muzzles and wings of infected bats, leading to a cascade of physiological issues (Verant et al., 2014;Warnecke et al., 2012). These physiological issues, in turn, cause abnormally frequent and long arousals from torpor through which infected bats can deplete their fat stores and consequently die (Reeder et al., 2012;Verant et al., 2014). In our model, after an average of 83 days (Lorch et al., 2013) exposed bats move into the infected class and become symptomatic and infectious. We assume that during hibernation susceptible bats can become exposed to Pd from contact with infected bats as well as with the environment. Pd continues to grow logistically in the hibernaculum environment, but during hibernation the rate is increased by the shedding of spores from infected bats.
The hibernation phase equations are given by: All of the equation terms present in the swarming phase also exist here, with an addition of the following: The rate at which the susceptible class (S) decreases is intensified by a densitydependent rate β of transmissive bat-to-bat contact, and consequently the rate at which the exposed class (E) increases is intensified by the same amount. Furthermore, the exposed class decreases at an additional rate τ indicative of disease progression. Here the infected class (I) appears and is made up of those individuals whose disease progressed from the exposed stage to the symptomatic-infectious stage. The infected class decreases at the WNS-induced mortality rate δ, as well as at the natural bat mortality rate μ. Finally, the free-living Pd class (P) increases at an intensified rate comprised of the sum of its natural growth rate η and the rate ω of Pd spore shedding from infected bats.

| Roosting reclassification
A bat's ability to leave the hibernaculum and relocate to a summer roosting site is likely dependent on its viability for flight and overall level of health. Therefore, some bats will die from WNS-associated conditions shortly after the end of the hibernation phase. Rather than track these individuals, for simplicity we incorporate an instantaneous reclassification event immediately before the start of the roosting phase.
The roosting reclassification equations are given by: where S E I , , f f f , and P f denote the final values of the susceptible, exposed, infected, and free-living Pd classes, respectively, at the end of the hibernation phase. Likewise, S E I , , i i i , and P i denote the resulting initial conditions for the roosting phase.
After the end of hibernation but before the beginning of roosting, a proportion sδ 1 + 1 of infected bats (I ) are assumed to be non-moribund and move to the exposed class (E), while the remaining infected bats die. Note that this proportion depends on the rate of disease-induced mortality δ and a scaling constant s estimated in Meyer et al. (2016). The remaining population classes are not affected by the reclassification.

| Roosting phase
The roosting phase lasts approximately 3 months from mid-May to mid-August (simulation days 274 through 309 and again days 331 through 365 for the nonreproductive periods). After bats leave hibernacula in mid-May, they disperse in small groups to summer roosting sites formed in trees, abandoned buildings, and under piles of wood or rock (Fenton & Barclay, 1980). Bats remain active throughout the roosting phase, thus having full immune system function and high body temperatures of over 20°C. Consequently, this phase is an important time for bats to heal and recoup body weight. Several studies have shown that no Pd infection occurs during the summer months, and that bats can gradually clear the fungus and heal damaged wings (Fuller et al., 2011;Meteyer et al., 2011). We model the recovery process by continuously moving bats from the exposed class to the susceptible class over the course of the roosting phase.
The roosting phase equations are given by: The susceptible class (S) increases at a rate a as some exposed bats (E) slowly clear their infection, and it decreases at the natural bat mortality rate μ. The exposed class decreases at the same rate a, as well as at the natural bat mortality rate μ. Since the hibernaculum environment in which free-living Pd grows remains the same year-round , the equation for the free-living Pd class (P) is the same as during the swarming phase.

| Roosting with birth
Little brown bats give birth during a 3-week period in the middle of the roosting phase (simulation days 310 through 330). We assume that apart from reproduction, bat behavior and disease dynamics during the birth subphase are equivalent to those of the nonreproductive roosting phase.
The birth subphase equations are given by: The birth subphase equations are identical to those of the nonreproductive roosting phase, except for the addition of a logistic growth term in the equation for the susceptible class (S).
The birth term causes the number of susceptible bats to increase at an additional rate γ with carrying capacity K Ml . Note that here N S E = + since the infected class (I) is absent.

| Control strategies
We investigate five promising WNS treatment methods that target different aspects of WNS dynamics (Table 1). We assume that each population can either have one control strategy actualized or none actualized at all. While strategies can differ between populations, both populations are always treated at the same intervention intensity ∈ α [0, 1], where α is a proportion of the model parameter of interest. α = 0 indicates that no control is implemented while α = 1 denotes maximum implementation intensity.

| Fungicide intervention
We consider fungicide intervention in the form of a chemical treatment that persists in the hibernaculum throughout the year, limiting the area within the hibernaculum suitable for fungal growth (Meyer et al., 2016;Reynolds et al., 2015). Alternatively, this intervention could be potentially achieved by physically removing sediment ideal for Pd growth from the hibernaculum environment (Meyer et al., 2016). We model fungicide intervention by decreasing the Pd carrying capacity K Pd by a proportion α.

| Microclimate intervention
Pd growth is known to increase with hibernaculum temperature and humidity, so bats hibernating in warmer and more humid conditions have higher fungal loads and suffer greater WNS impacts (Langwig et al., 2016;Marroquin et al., 2017). Hence, modifying hibernacula entrances to create cooler sites or restricting bats' access to the warmer portions of hibernacula could reduce the Pd load on bats and potentially slow or stop the progression of WNS symptoms (Langwig et al., 2016;Meyer et al., 2016). We note that while it is possible for microclimate intervention to reduce both the rate of WNS-induced mortality δ as well as the rate of disease progression τ, we want to compare our results directly with those of Meyer et al. (2016). Therefore, we model microclimate intervention by reducing the rate of WNS-induced mortality δ by a proportion α. Since the proportion of nonmoribund bats at the end of hibernation is determined by the expression sδ 1 + 1 , decreasing disease-induced mortality also increases the number of viable infected bats after hibernation.

| Soil bacteria intervention
Volatile organic compounds produced by several species of soil-associated bacteria, such as Rhodococcus rhodochrous, are known to exhibit anti-Pd effects . As such, soil bacteria intervention is a contact-independent control that could inhibit conidial growth both on bat tissue and within the environmental reservoir. Therefore, we implement this intervention strategy by scaling two parameters. We decrease both the rate of bat progression from the exposed class to the infected class τ and the natural growth rate of environmental Pdη by a proportion α.

| Ultraviolet (UV) light intervention
Pd is extremely sensitive to UV light since it lacks a critical enzyme in the alternate excision repair pathway, which is known to contribute to the repair of DNA lesions induced by UV light (Palmer et al., 2018). We hypothesize that artificial UV light emitted from bulbs installed in the hibernaculum would kill some Pd on bats and in the environmental reservoir, as well as inhibit the growth of surviving fungus on bats. Hence, we model UV light intervention with three model modifications. We introduce a Pd mortality term αP − into each dP d t equation to account for the loss of fungus on hibernaculum surfaces due to UV exposure. Since a lower Pd load on infected bats will potentially slow or stop the progression of WNS symptoms, we decrease the WNS-induced mortality rate δ by a proportion α. Finally, we reduce the rate of bat progression from the exposed class to the infected class τ by the same proportion α to account for nonlethal damage to fungus living on infected bats.

| Vaccination
Raccoon poxviruses expressing Pd calnexin and serine protease have been shown to induce antifungal immunity in little brown bats, which could potentially protect bat populations against WNS (Rocke et al., 2019). Practicably, these oral WNS vaccine candidates could be administered as a topical gel or paste that bats would then ingest when grooming (Rocke et al., 2019;Stading et al., 2017). While it would be easiest for natural resource managers to distribute the vaccine during the hibernation phase, when bats are congregated in the same location and are in torpor, this may not be the ideal time to vaccinate. During hibernation the susceptible pool would be at its lowest value, since many bats would already be infected. Moreover, disturbing bats during torpor might lead to more arousals and further exacerbate WNS symptoms. Thus, we choose to model vaccine administration immediately after the reproductive phase, when the susceptible class is largest and bats are more easily accessible than during swarming. We incorporate vaccination with the addition of a vaccinated class (V ) to each population. After the end of the birth subphase but before resuming the nonreproductive roosting phase, a proportion α of susceptible bats are vaccinated and moved to the vaccinated class, where they are no longer susceptible to Pd infection. In Cornwell et al. (2019) it was shown that annual administration of a vaccine that offers lifetime immunity constitutes the most effective vaccination strategy, so we assume that vaccination occurs every year and that vaccinated bats do not leave the V class except due to natural mortality. Moreover, V -class bats move between populations in the same manner as all other bats and are included in N for reproduction purposes. Appendix A displays the modified model equations with vaccination. Table 2 summarizes all of the parameters used in the model. Without the presence of disease (β ϕ = = 0) a stable disease-free equilibrium exists, but due to natural bat mortality between the annual birth subphases, the population carrying capacity K Ml must actually be higher than the initial population size. For K Ml , we chose the smallest positive integer that would cause the population to reach 15,000 after every annual birth subphase (after transiently exceeding 15,000 for the first few years of the simulation) in the absence of disease. Typically, 95% of females each give birth to one pup over a 21-day period (Fenton & Barclay, 1980), resulting in a fecundity rate of ∕ γ = 0.5(0.95 21). Banding studies have found that approximately 4% of little brown bats relocate between hibernacula and 12% relocate between summer roosting sites, for an overall 6% relocation rate any given year (Kaleigh et al., 2013). Therefore, we consider dispersal percentages ranging from 0 to 10 as biologically reasonable. In Meyer et al. (2016) it was estimated that the post-hibernation recovery probability of a viable infected bat is 0.75, so in our model it corresponds to a rate of ∕ a = 0.75 92. With the exception of the relative contributions of environment-to-bat and bat-to-bat contact to the transmission of Pd, the remaining parameters are known biological constants or were well-approximated in Meyer et al. (2016).

| Parameter estimation
Little is known about the relative contributions of environment-to-bat (ϕ) and bat-to-bat (β) contact to the transmission of Pd, even though previous studies have identified this ratio as an important determinant of disease dynamics and control efficacy (Meyer et al., 2016). Hence, we conduct our model analyses for three pathogen transmission cases: primarily environment-tobat contact, equal contributions of bat-to-bat and environment-to-bat contact, and primarily bat-to-bat contact. To select corresponding ϕ and β values for each transmission case, we used the estimation that 25% of the population survives 2 years after the initial introduction of Pd (Blehert et al., 2009). Figure 2 shows the contour of ϕ β , values corresponding to 25% survival after 2 years in a single population without control. The range of ϕ satisfying this survival criterion is approximately ⋅ (0, 6.92 10 ) −13 . We defined primarily environment-to-bat transmission to be located at 90% of the ϕ range along the contour, equal contributions of batto-bat and environment-to-bat transmission to be at 50% of the ϕ range along the contour, and primarily bat-to-bat transmission to be at 10% of the ϕ range along the contour. The three transmission cases we tested are listed in Table 3. We assume that no free-living Pd is initially present in any hibernacula, and that both populations begin only with susceptible and exposed bats. We found that the initial ratio of exposed to susceptible bats can impact the metapopulation size during the first few years, but with seasonal movement between populations all simulations across a reasonable range of ratio Not applicable values converge before 10 years (not shown). In this study we focus on long-term evaluation of management efficacy, so we disregard the initial transient behavior observed after short simulations. Henceforth, we run all simulations for 10 years and use the initial condition S E I P ( , , , ) = (14999, 1, 0, 0) 0 0 0 0 for both populations. Consequently, implementing some intervention x in population A and some intervention y in population B will produce exactly the same result as implementing intervention y in population A and intervention x in population B. For this reason, we only display the lower triangular matrix of results in Figures 4 and 6. 3 | RESULTS

| Single population
In nearly every situation, implementing a control strategy at any intervention intensity will improve bat population survival relative to no intervention ( Figure 3). The only exception is microclimate intervention when disease is primarily transmitted by bat-to-bat contact. In this case, reducing the rate of WNS-induced mortality only increases the number of secondary infections each I -class individual imparts to other bats, thereby decreasing population survival. The efficacy of each control strategy in all other situations increases with higher values of intervention intensity α. The strategies differ, however, in the slope of population survival with respect to intervention intensity. Microclimate intervention is only effective at very high intervention intensities, while vaccination can achieve 50% population survival at a low ≈ α 0.18. The other strategies fall between these extremes. Interestingly, the plots of percent survival against intervention intensity are convex for all control strategies except vaccination, indicating that improvements in intervention intensity are generally more consequential at higher values.

F I G U R E 2 Heat map displaying percent population survival for all values of ϕ β
, in a single-population setting without control. The black contour indicates 25% survival, with the three tested transmission cases marked. The population begins with one exposed bat and 14,999 susceptible bats and the simulation is run for 2 years Since we model control strategies that target multiple pathogen transmission pathways, the relative contributions of environment-to-bat and bat-to-bat transmission affect each strategy in a distinct way (Figure 3). When implementing microclimate intervention, soil bacteria intervention, or vaccination, higher population survival is always attained with higher contributions of bat-to-bat contact to Pd transmission. Fungicide intervention and UV light intervention, on the other hand, directly reduce the size of the environmental Pd reservoir. As a result, for these controls the lines denoting primarily environment-to-bat transmission and primarily bat-to-bat transmission cross twice at some intervention intensities α 1 and α 2 (where α α < 1 2 ). At intervention intensities below α 1 , a higher contribution of bat-to-bat contact to the transmission of Pd yields better population survival. The model exhibits nonlinear behavior as α increases between α 1 and α 2 , and for intervention intensities of α α > 2 a higher contribution of environment-to-bat contact now results in better survival. The values of α α , 1 2 differ between the two controls, and additional simulations (not shown) reveal that these critical points depend on the total number of years simulated.
For an example of how the primary route of pathogen transmission can affect the relative efficacies of control strategies, consider a mid-range intervention intensity of α = 0.67 (Figure 3). Under primarily environment-to-bat transmission the most effective strategy is UV light intervention, followed in order by vaccination, fungicide intervention, soil bacteria intervention, and microclimate intervention. The ranking of control strategies remains the same given equal contributions of environment-to-bat and bat-to-bat contact to Pd transmission. But if we instead consider primarily bat-to-bat transmissive contact, the ranking of most effective controls changes to the following: vaccination, UV light, soil bacteria, fungicide, and microclimate.

| Two populations
Our results demonstrate that metapopulation dynamics may impact the survival of bat populations affected by WNS and modify the rankings of control strategies for managers seeking to maximize the combined size of two populations. Whether movement between populations increases or decreases metapopulation survival, as well as the severity of its effects, depends on a variety of factors such as the chosen combination of control strategies, the intervention intensity, and the primary route of pathogen transmission. Figure 4 compares final metapopulation survival percentages after 10 years given no dispersal with final metapopulation survival percentages given 7% annual dispersal, under the assumption of equal bat-to-bat and environment-to-bat transmission. Although the differences are subtle, some pairings of control strategies are more effective in the presence of interpopulation movement while others are more effective when the two populations remain isolated from each other. Consequently, the rankings of preferred strategies differ between these T A B L E 3 Summary of transmission cases

Transmission case ϕ β
Primarily environment-to-bat transmission ⋅ 6.24 10 −13 ⋅ 6.79 10 −7 Equal contributions to transmission ⋅ 3.44 10 −13 ⋅ 3.89 10 −6 Primarily bat-to-bat transmission ⋅ 6.80 10 −14 ⋅ 9.00 10 −6 two scenarios. For example, given the absence of dispersal, implementing UV light intervention in population A and fungicide intervention in population B is preferable to the combination of vaccination in A and fungicide intervention in B. Yet, in the presence of 7% dispersal, the ranking is reversed such that vaccination in A and fungicide intervention in B is a better combination than UV light intervention in A and fungicide intervention in B. Ranking reversals can also occur when one of the populations is left untreated. For instance, if a manager faces the choice of vaccinating population A while leaving population B untreated or applying soil bacteria to A and fungicide to B, our model points to the latter option in the absence of dispersal. Incorporating 7% dispersal, however, shows that vaccinating population A and leaving B untreated is the better choice.
Other differences in pairwise rankings between these two dispersal cases are evident from Figure 4, but in every example at least one of the control strategies in at least one of the pairs is UV light or vaccination. The pairwise rankings of control combinations not including either UV light or vaccination in either pair remain the same, even though many of those control combinations do still increase or decrease in efficacy with the addition of dispersal. This principle also holds for lower intervention intensities than shown. When intervention intensities are raised above ≈ α 0.97, however, all of the control strategies become increasingly similar in effectiveness and we observe ranking reversals even among control combinations without UV light or vaccination (not shown). Figure 5 summarizes how the exchange of bats between populations affects the efficacy of each control strategy individually. As dispersal increases, the effectiveness of vaccination increases while the effectiveness of every other strategy generally decreases. UV light intervention F I G U R E 3 Each plot displays percent population survival against intervention intensity for all three transmission cases in a single-population setting. The plots differ only in the control strategy implemented. The population begins with one exposed bat and 14,999 susceptible bats and simulations are run for 10 years at mid-range intervention intensities of α 0.4 < < 0.6 is the only exception, since after initially decreasing in efficacy at low values of dispersal it begins to slowly become more effective as dispersal rises above 1.5%. When each control strategy is viewed individually, vaccination and UV light intervention are the only ones whose rankings switch due to dispersal. But as demonstrated above, pairing one strategy in population A with another strategy in population B yields more ranking reversals.
Note that some control strategies are much more effective than others ( Figure 5). Implementing a very effective control in one population at a high intervention intensity and an ineffective control in the other population at a low intervention intensity will often result in better metapopulation survival than if both controls are implemented at midrange intervention intensities (not shown). Consequently, it is usually better to manage the metapopulation by focusing all efforts on the most effective control strategy in one population instead of distributing resources to achieve mediocre intervention intensities in both populations. On the other hand, control strategies are not additive and more extreme implementations of some strategies can sometimes be as effective as the best strategies.
Whether control efficacy is greater or lower in the presence of metapopulation dynamics than in the single-population case depends on the magnitude of dispersal, the controls implemented, and the primary route of transmission. Figure 6 illustrates the potential outcomes by displaying percent changes in survival due to dispersal. Notice that the absolute value percent change in survival due to dispersal is often, although not always, greater when 8% of bats move between the populations than when 4% move. In other words, a higher proportion of bats that exchange populations can lead to a greater influence of interpopulation movement on metapopulation survival and control efficacy.
F I G U R E 4 Each cell represents a separate simulation, with the color denoting percent survival of the metapopulation. The letters denote control strategies as follows: N is no control, F is fungicide, M is microclimate, B is soil bacteria, UV is UV light, and V is vaccination. The left plot has no dispersal, while the right plot has 7% annual dispersal. Both populations begin with one exposed bat and 14,999 susceptible bats, simulations are run for 10 years, the intervention intensity is α = 0.9, and equal contributions of environment-to-bat and bat-to-bat contact are assumed Metapopulation dynamics do not affect systems with the same control strategy in both populations, since in that case the two populations are identical throughout the course of the simulation (Figure 6). When pairing vaccination with another control, the strong positive impact of dispersal on vaccination nearly always masks the weaker effects of dispersal on other controls. For example, consider vaccination in population A and microclimate intervention in population B (Figure 6). Since dispersal strongly increases the efficacy of vaccination but only slightly decreases the efficacy of microclimate intervention, the combination of vaccination in one population and microclimate intervention in the other will still be positively affected by increased dispersal. Combinations of vaccination with UV light intervention, however, can be either positively (Figure 6c) or negatively (Figures 6a,b) impacted by increased dispersal because dispersal drastically changes the efficacies of both of these controls individually.
The impact of dispersal on control efficacy is more varied among combinations not involving vaccination. Sometimes dispersal impacts the combination more negatively than it impacts each control individually. For example, notice from Figure 6c that implementing soil bacteria in population A with no control in population B results in a slightly negative percent change in survival due to dispersal, and implementing fungicide in population A with no control in population B results in a slightly positive percent change in survival due to dispersal.
F I G U R E 5 Each plot displays percent survival of the metapopulation against annual dispersal percentage. The control strategies are only implemented in population A, while population B remains untreated. The plots differ only in intervention intensity. Both populations begin with one exposed bat and 14,999 susceptible bats, simulations are run for 10 years, and equal contributions of environment-to-bat and bat-to-bat contact are assumed. Here absence of control results in 0.9% metapopulation survival, so all interventions except microclimate at α = 0.3 improve survival Nonetheless, pairing soil bacteria in A with fungicide in B yields a rather substantial negative percent change in survival due to dispersal. In other situations dispersal impacts the combination less negatively than it impacts each control individually. For an example, similarly consider the combination of UV light intervention in population A and microclimate intervention in population B (Figure 6a). The lack of a single dynamical trend for pairings of different controls suggests that interactions between bats and Pd exhibit underlying complexity due to multiple routes of pathogen transmission.
Surprisingly, the negative effects of interpopulation exchange are strongest when Pd transmission is driven primarily by environment-to-bat contact. The greatest negative percent change in survival due to dispersal occurs when soil bacteria intervention is implemented in population A and population B is left untreated, with WNS primarily transmitted by environment-to-bat contact: in this scenario dispersal of 8% causes a −26% change in survival relative to no dispersal (Figure 6a). When disease is primarily transmitted by bat-to-bat contact, however, the greatest negative change in survival due to 8% annual dispersal is approximately −15%, which results from applying UV light intervention to population A and fungicide intervention to population B (Figure 6c). The positive impact of dispersal on vaccination does not change as much across the three transmission cases.

| DISCUSSION
WNS is one of the most pressing wildlife diseases of this century (Blehert et al., 2009;Lorch et al., 2016), necessitating prompt development of control strategies to mitigate its severity and progression. Conducting field trials of proposed controls might not be prudent or cost-effective, so mathematical models can instead be used to easily test said strategies without endangering atrisk fragile ecosystems. Although bat metapopulation dynamics are known to be a key part of the chain of pathogen transmission in other infectious diseases (Blackwood et al., 2013), little is known about how movement between populations impacts proposed WNS interventions. We posed and analyzed a mechanistic model of WNS in little brown bats that incorporates five promising control strategies targeting several disease processes. Most importantly, we simultaneously considered the dynamics arising from an interplay between Pd transmission, choice of control strategies, and bat dispersal. Our results demonstrate that the most effective management decisions take exchange between populations into account, and highlight the importance of further quantifying bat metapopulation dynamics in conjunction with the transmission of Pd.
Our single-population baseline analysis confirms the conclusions of Meyer et al. (2016). In accordance with their findings, our results likewise call attention to the importance of understanding the relative contributions of environment-to-bat and bat-to-bat contact to the transmission of Pd. In our model, fungicide intervention is mathematically analogous to their Control IV (decreasing the environmental reservoir size) and microclimate intervention is mathematically analogous to Control I (thermal refugia), while the remainder of our controls are new. Since Meyer et al. (2016) only considered each control at intervention intensities sufficiently high to yield 33% or better 5-year population survival, the study did not detect that at low intervention intensities fungicide intervention yields better population survival with higher contributions of bat-to-bat contact to pathogen transmission. At higher intervention intensities, however, our results concord with theirs: fungicide and microclimate interventions are plausible for situations in which transmission is driven primarily by environment-to-bat contact, although the latter requires drastic intervention intensities.
F I G U R E 6 Each cell represents a separate simulation, with the color denoting percent change in metapopulation survival between no dispersal and annual dispersal values of 4% (left plot) or 8% (right plot). The letters denote control strategies as follows: N is no control, F is fungicide, M is microclimate, B is soil bacteria, UV is UV light, and V is vaccination. The subfigures differ only in the primary route of Pd transmission. Both populations begin with one exposed bat and 14,999 susceptible bats, simulations are run for 10 years, and the intervention intensity is α = 0.8. (a) Primarily environment-to-bat transmissive contact; (b) equal contributions of environment-to-bat and bat-to-bat contact to pathogen transmission; (c) primarily bat-to-bat transmissive contact Managers should pay especially close attention to the potential impact of bat metapopulation dynamics when choosing between control strategies that are similarly potent within a singlepopulation setting, as was the case for vaccination and UV light intervention in our model. While the difference in efficacy between some strategies might be negligible in the absence of movement between populations, single-population models cannot accurately predict control efficacies in a multi-population setting. According to our model, failing to account for 10% annual dispersal can lead to a suboptimal management decision that results in 23% lower metapopulation survival-a huge number for bat colonies near their tipping point (see Figure 5, α = 0.9).
Our simulations also reveal that the effect of metapopulation dynamics on control efficacy depends on the primary route of pathogen transmission. Previous studies established that successful implementation of controls is highly dependent on the relative role of each transmission route (Meyer et al., 2016). In our metapopulation model those results are even more pronounced. Differences in Pd transmission not only change the effectiveness of individual controls, but also determine how dispersal affects each control. For example, while fungicide and microclimate interventions are much less effective in the presence of dispersal if disease is primarily transmitted by environment-to-bat contact (Figure 6a), their efficacy remains virtually the same regardless of the dispersal amount if disease is instead transmitted primarily by bat-to-bat contact (Figure 6c). This conclusion suggests that empirical research investigating the metapopulation dynamics of bats in conjunction with the transmission of Pd is necessary to construct an accurate picture of WNS and, in turn, to develop more effective control strategies.
In general, our results indicate that WNS control strategies targeting fungal growth or disease progression should be expected to decline in efficacy with higher levels of dispersal. Of the controls we tested, vaccination was the only one to consistently benefit from the presence of bat dispersal. This accords with our simplifying assumptions that vaccination grants lifetime immunity and that vaccinated bats move between populations in the same way as nonvaccinated bats, thereby lessening the force of infection in every population they move to. In this manner vaccination can even improve the survival of hibernacula whose locations are unknown, since vaccinated bats would disperse throughout the entire regional metapopulation. Unless such a control can be developed and efficiently administered, managers should expect the exchange of bats between populations to diminish the effectiveness of most control strategies. We structured our model to investigate the metapopulation dynamics of WNS given the current understanding of the disease. Necessarily, this entailed making modeling assumptions that we believe are representative of actual disease dynamics and consistent with previous studies. Due to the rapid rate of research on WNS and North American bat species, our methodology, results, and conclusions might need to be revisited in future studies once new data become available and additional control methods are proposed. While our simulations point to some controls as more effective than others, we did not consider potential differences between difficulty of implementation for the various control strategies. The present lack of extensive laboratory and field trials for the controls we tested implies that our results should be interpreted as provisional; although we suggested a ranking of control strategies, pragmatic concerns over safety and feasibility could limit the options available to managers.
This study is the first to investigate WNS management in a metapopulation model, so we kept our model as simple as possible. In particular, we only considered two populations and assumed that dispersal does not change their sizes. To our knowledge this assumption is reasonable, since there is no evidence for the existence of hibernacula that absorb many individuals while contributing few to other hibernacula. If some species of bats do tend to favor larger (or smaller) hibernacula when moving between populations, however, then managers should consider such directionality of dispersal in their WNS treatment plans. Moreover, further research is needed to elucidate whether this disease system operates under density-dependent or frequency-dependent transmission and how metapopulation dynamics might affect WNS control strategies in a larger network of hibernacula. We suspect that treatment of well-connected hubs would be the best management strategy, and urge for both field-based and modeling studies to investigate WNS interventions within more complex metapopulation structures.
In summary, we have provided a modeling approach that can be used to examine disease dynamics arising from multiple transmission routes of Pd in a metapopulation setting. We assessed the long-term success or failure of several control strategies and, more significantly, we demonstrated that those outcomes depend on the presence and magnitude of bat movement between populations. Fundamentally, our results accentuate the importance of considering bat metapopulation dynamics together with different routes of pathogen transmission, since the emergent dynamics influence how dispersal between populations changes the effectiveness of each control strategy. A.2 | Dispersal reclassification

ACKNOWLEDGMENTS
The dispersal reclassification equation with vaccination is given by: