BatTool: an R package with GUI for assessing the effect of White-nose syndrome and other take events on Myotis spp. of bats

Background Myotis species of bats such as the Indiana Bat and Little Brown Bat are facing population declines because of White-nose syndrome (WNS). These species also face threats from anthropogenic activities such as wind energy development. Population models may be used to provide insights into threats facing these species. We developed a population model, BatTool, as an R package to help decision makers and natural resource managers examine factors influencing the dynamics of these species. The R package includes two components: 1) a deterministic and stochastic model that are accessible from the command line and 2) a graphical user interface (GUI). Results BatTool is an R package allowing natural resource managers and decision makers to understand Myotis spp. population dynamics. Through the use of a GUI, the model allows users to understand how WNS and other take events may affect the population. The results are saved both graphically and as data files. Additionally, R-savvy users may access the population functions through the command line and reuse the code as part of future research. This R package could also be used as part of a population dynamics or wildlife management course. Conclusions BatTool provides access to a Myotis spp. population model. This tool can help natural resource managers and decision makers with the Endangered Species Act deliberations for these species and with issuing take permits as part of regulatory decision making. The tool is available online as part of this publication.


Background
Bats in the Myotis genus, including the Little Brown Bat (Myotis lucifugus) and Indiana Bat (M. sodalis), face population-level threats in the eastern United States and Canada. The emerging fungal disease White-nose syndrome (WNS) has caused massive decreases in population sizes and is predicted to contribute to further declines as the disease spreads farther west across North America [1]. The Little Brown Bat was one of the most common bat species in the eastern United States until the arrival of White-nose syndrome. The drastic decrease in Little Brown Bat populations has led the U.S. Fish and Wildlife Service to consider listing the species under the Lakewood, CO, USA Full list of author information is available at the end of the article Endangered Species Act [2]. Conversely, the Indiana Bat was one of the first species listed under the Endangered Species Act [3]. In addition to WNS, these two species face other threats from anthropogenic activities such as wind energy development [4,5].
Population models have emerged as one method to understand and manage wildlife populations in light of uncertainty [6]. These models may include biologically important attributes such as different life stages (e.g., juveniles and adults). Decision makers and resource managers use these models to explore different scenarios. Possible scenarios might include no management (status quo) or different management approaches. Possible stressors that might be included within the models include harvest (e.g., hunting or fishing) or other takes such as energy development or habitat loss. These models may also address variability and uncertainty through the inclusion of stochasticity. Models may include variability relating to http://www.scfbm.org/content/9/1/9 small population sizes (demographic stochasticity), variability associated with environmental conditions (e.g., droughts vs wet years; environmental stochasticity), and uncertainty in parameter estimates (e.g., 2 births and 1 death per year vs 10 births and 9 deaths per year) [7].
Thogmartin et al. [8] developed a population model for studing the effects of WNS on Myotis spp. The original model was written in Matlab (MATLAB and Statistics Toolbox Release 2012b, The MathWorks, Inc., Massachusetts, United States), but the source code was not included as part of the publication nor easily usable by decision makers at agencies such as the U.S. Fish and Wildlife Service. We developed this model into an R [9] package to assist decision makers in using the code. R was chosen because it is Open Source and freely available to interested users. The model we present within this manuscript contains two different components: 1) a command-line deterministic and stochastic model and 2) a graphical user interface (GUI). The command line option allows R-savvy users to include the model as part of their own script. The GUI was specifically developed for U.S. Fish and Wildlife Service decision makers desiring a tool specifically implementing the model presented by Thogmartin et al. [8].

Underlying population model
Thogmartin et al. [8] previously published the population model forming the backbone of BatTool. We include a flow chart of the model (Figure 1), the equations (Equations 1,2,3,4,5,6,7,8 and 9), and variables (Table 1) within this article as well as an overview of the biology underlying the model. Additional analysis of the model was published with the original article [8].
Indiana Bats and Little Brown Bats migrate between summer maternity roost sites and winter hibernacula. Pups are born at roost sites and then migrate to hibernacula during the fall. At this point, the pups become first-year breeders (colloquially referred to as juveniles within our model). The juveniles overwinter at the hibernacula. The juveniles then migrate to summer roost sites during the spring. Our model does not directly consider spring migration mortality. A proportion of the juveniles breed. The breeding and non-breeding juveniles may have different survival rates within the model for the summer and fall seasons. The juveniles migrate back to hibernacula during the fall and become adults. The adults then overwinter and migrate in the spring to the summer roost sites. Like the juveniles, there are both breeding and nonbreeding adults. After summer, the adults migrate to the hibernacula during fall. This cycle continues until the bats die [4,5].
Our model reports the bat population size during winter that would be found at a specific hibernacula. This was done because most bat surveys are conducted at the hibernacula and the winter populations are best understood and monitored for both the Little Brown Bat and Indiana Bat. Our model only follows females within the model. The input population is divided by two and the results are multiplied by two under the assumption of an even gender distribution. This is a common assumption in population ecology because males do not limit population size in many non-monogamous species including bats [4,5,7]. Our model is a matrix model (a series of discretetime difference equations) that follows the population P through time. P(t) is a two-entry vector with the top entry being the number of juveniles and the bottom entry being the number of adults at time t, where t is the time in year. The projection matrix, A, moves the population forward one year (Table 1) The population at the next year is We decomposed the projection matrix A (Equation 1) into the seasonal projection matrices in order to facilitate seasonal "take" and include WNS mortality during the winter. Although take is formally defined under the Endangered Species Act of 1973 to include "harass, harm, pursue, hunt, shoot, wound, kill trap, capture, or collect, or to attempt to engage in any such conduct", our model considers all take as mortality-causing events.
where ⊗ is the outer product (element-wise matrix multiplication function). A is decomposed into 5 matrices (Equations 4,5,6 and 7. The spring and fall projection matrix for non-reproducing individuals becomes The summer projection matrix for reproducing individuals becomes The spring projection matrix becomes The winter projection matrix becomes This allows the seasonal take parameters (winter τ wi ; spring τ sp ; summer τ su ; and fall τ fa ) to be inserted into the projection matrix: A simple, ceiling carrying capacity, K , is also used within the model. Once K is reached, A becomes the identity matrix. The value for K can either be specified by the user or come from population survey data. BatTool also includes optional stochasticity. Environmental stochasticity is included by modifying the input parameter with a uniform distribution, parameter ± Uniform(-envs, envs), where "envs" is a user-specified value. A safeguard is also included to ensure the parameter stays within (0, 1). Demographic stochasticity may also be included within the model. When demographic stochasticity is included, a binomial distribution replaces the simple matrix calculations. As an example, the number of juveniles surviving winter would become The births are also replaced by a binomial distribution. This is appropriate because each female Myotis bat may only produce a maximum of 1 offspring per year. Another distribution would be needed if an individual could produce more than one offspring (e.g., Poisson).

Data inputs
BatTool includes several different data inputs (

Lambda table
The ratio of the population at year t + 1 compared to the year t is commonly called lambda in population ecology http://www.scfbm.org/content/9/1/9   [10]. This is because the growth rate for a linear model (such as our matrix projection model) is also the eigen value, which is commonly represented with the Greek letter lambda (λ) [11].

LambdaEstimates table
The LambdaEstimates table contains estimated lambda values for each mentioned hibernaculum. This table is populated with hibernaculum-specific population rates of change [8].

Hibernacula table
The Hibernacula table lists hibernacula names, counties, take values, and observed population counts. The hibernacula counts are plotted as part of the output. The carrying capacity, K, defaults to be 1.5 × the maximum observed population at a hibernacula. Also, the starting population within the model is the last year of the observed population counts, but this value can be changed in the GUI by the user. The take description includes the start, duration, and amount occurred in each season. We included an example table that the user can amend in their own studies.

WNS Infection tables
The WNSInfection Probability

Package installation
This package may be installed by downloading it from the journal's additional materials. We have included both the raw package ending in tar.gz (Additional file 1) and a file compiled for Windows ending in .zip (nested within Additional file 2). Additionally, File 2 is a zip file that also contains data necessary for the GUI to Run. To install the package, use the package installer included as part of R (see ?install.packages for help). Additional install directions are included as part of the of the readme.txt file located in the Additional file 2. The gWidgetstcltk package and required dependencies are needed for the GUI to work. After installing the package, use the library(BatTool) to load the tool.

Command lines tools
The two major functions within the package are the deterministic model (main_pop) and the stochastic model (pop_stochastic). To see an example of the deterministic model, use the following lines of code:

> library(BatTool) > example(main_pop)
This will produce Figure 2. In this example, the population grows until it reaches its carrying capacity. The example also shows the juvenile and adult populations.
The stochastic model runs multiple simulations and includes several different options worth noting. Running the example for the function will show 50 example population trajectories with the mean and 95% credibility interval overlayed on the plot (Figure 3). This function requires that the number of simulations (or replicates) be specified by the user. Three levels of stochasticity may be run with this model (Figure 4). The model includes parameter uncertainty for any lambda value or range of lambda values. Environmental stochasticity may be specified with a value of zero indicating no environmental stochasticity. Demographic stochasticity may be turned on. Both types of stochasticity may be included. The different levels of stochasticity are also shown with the following example for this function. > library(BatTool) > example(pop_stochastic)

Graphical user interface
The GUI is housed within a demo in the BatTool package. Models from the GUI start with the last year of observed data being year 0 (e.g., if there are observations through 2012 for a hibernacula, year 1 of the output would be 2013).

> library(BatTool) > demo(GUIcode)
Running the demo will launch the GUI ( Figure 5). Changing the hibernacula number will load data for a new hibernacula after the return key is pressed. Clicking on the "Hibernacula number:" button will launch a table that shows hibernacula information including user-contributed names corresponding to hibernaculumspecific identification numbers. The default starting population is the last population from the last observed year and the default Hibernaculum limit is 1.5× the largest observed population at the hibernacula. Two different scenarios may be run and different options may be set for each scenario. These options are listed under different tabs ( Table 3).
The default WNS Year of Infection is based upon the lookup table if the data are available. If the data are not available, the probability of infection for the specified species is used and a random year of infection is used for each simulation. Alternatively, the year of infection may be entered manually; similarly, the probability of infection occurring within a hibernaculum can be adjusted manually.
The default female WNS take parameters for each county are part of the Hibernacula table. These parameters may be changed in either the GUI or the csv file. Example hibernacula 998 contains non-trivial take parameters as an example case. The female take parameters used in the GUI only affect adults. Conversely, the simple model allows either the adult population or juvenile population to suffer take events; similarly, the probability of infection occurring with a hibernaculum can be adjusted manually.
Results from the GUI are stored in a new folder, "ResultsSingleHib\temp". The user may change the temporary folder name within the GUI prior to each simulation; otherwise, past runs will be overwritten. The user can also modify the output figure under the "Graphing and results options" tab. The figure resulting from the GUI (Figure 6) includes the means and credible intervals for two scenarios, any previously observed population data, as well as 4 horizontal lines. The horizontal line at zero represents extinction. The horizontal line at 10,000 bats represents a priority benchmark size for the winter population according to the U.S. Fish and Wildlife Service recovery plan, whereas the horizontal lines at 500 and 2,000 bats represents lower priority hibernacula sizes.
Clicking "RUN" causes the simulations to start and clicking "Restore defaults" reverts settings to their default values. User settings are reported in the results folder. http://www.scfbm.org/content/9/1/9 The default values for some parameters are loaded from parameter tables ( Table 2). The Critical Parameter tab does not contain any user editable data and is provided as a summary. The take parameters (take, begin year, and duration) may be changed for each season.

Import custom data into the GUI
Custom data may be incorporated into the GUI two different ways. First, values may be directly entered. Second, input tables may be changed. The WNS scenarios may be changed by either changing the default scenario tables or editing the Scenario 1 file (WNS_other_1.csv) or Scenario 2 file (WNS_other_2.csv) file in the working directory.

Case study Background
Population viability analysis (PVA) is a quantitative framework for understanding the effects of stressors on populations [6]. This approach allows conservation biologists, decision makers, and risk assessors to compare different management actions (or lack of action). The U.S. Fish and Wildlife Service uses an analytical framework for assessing stressors, which includes PVA as one component. Assessing the effects of wind energy development on the Indiana Bat consists of three steps: 1. Evaluating individual Indiana Bat exposure to action-related stressors and response to that exposure (i.e., likelihood of exposure to wind turbines and the likelihood of death or injury upon exposure); 2. Integrating those individual effects to discern the consequences to the population(s) to which those individuals belong (i.e., what are the effects to the reproductive potential and survival of maternity colonies and hibernacula); and 3. Determining the consequences of any population-level effects to the species at the Recovery Unit and species levels (i.e., will this action affect the likelihood of recovery at these two scales?) For our case study, we focus on Step 2. Our location is based upon an actual project, but the location has been anonymized for this case study to maintain data confidentiality.

Model settings
We conducted two different assessments. The first was for a maternity colony. The second was for a hibernaculum. All parameters were the same across the two assessments other than the initial population size and hibernaculum limit. A stationary condition (λ ∈ [0.99 − 1.01]), but slightly declining population due to model stochasticity, was used. The scenarios used for this assessment did not include White-nose syndrome. Each simulation was run for 50 years and 1,000 simulations were run. The maternity colony assessment had a starting population of 80 and hibernaculum limit of 200. Two bats were taken for 30 years each spring and fall for an annual take of 4 bats per year. This level of take would represent a small, but reasonable loss associated with a wind farm. For the hibernaculum assessment, two different take scenarios were examined. The first scenario only included the loss of 2 bats each spring and fall. This scenario results in the same pattern of take as the maternity colony take scenario. The second hibernaculum scenario included the loss of 300 bats each spring and fall for 30 years for an annual take of 600 bats per year. This level of take would represent take from multiple facilities affecting a hibernaculum. These values are take authorizations requested by wind energy generation concerns. Note that our model does not include spatial structure and this limits the use of our http://www.scfbm.org/content/9/1/9 model for studying wind energy take at the species level or other large spatial scales. This limitation occurs because the model was developed to initially assess White-nose syndrome at a hibernaculum.

Results and conclusions
The take of 4 females per year (2 during spring, 2 during fall) caused a greater population decline for the maternity colony, but not the hibernaculum (Figure 7, the left panel versus the center panel). The take of 600 females per year was sufficient to increase the rate of decline as well (Figure 7, right panel). Simply evaluating the loss of individuals at hibernaculum or larger scales failed to account for the spatial dynamics of the species. For example, take of only 4 females per year did not yield a detectable effect at the hibernaculum level, but the loss of 4 individuals could lead to the loss of an entire maternity colony if immigration is insufficient to overcome the long-term loss of breeding individuals to take from wind energy development. This impact was not detectable by simply evaluating the loss of 4 individuals from the hibernaculum population because the magnitude of loss relative to the population size was minuscule relative to the stochasticity experienced by the population. These findings indicate Figure 7 Case study figures. Figures from the case study from three take scenarios. Scenario 2 had take for each set of simulations. The solid lines are the mean outputs and the dashed lines are 95% credible intervals. See the text for differences between scenarios. http://www.scfbm.org/content/9/1/9 that efforts to minimize bat mortality (e.g., altering turbine speeds [12]) may be needed at the development site if real losses are equivalent to those tested in these simulations.

Conclusions
BatTool is an R package designed to help natural resource managers and decision makers. The package contains a population model accessible through both a GUI and command-line interface. The main functions of the command line are the main_pop model function and pop_stochastic function. These functions may be used to simulate the population-level effects of WNS and wind energy development. There is also a GUI included as part of this package allowing users who are less comfortable with a command-line interface to use and change the model inputs. Because of the ease of use of the GUI, this package can also be used as part of population ecology or natural resource management courses.