Forecasting Seasonal Habitat Connectivity in a Developing Landscape

Connectivity and wildlife corridors are often key components to successful conservation and management plans. Connectivity for wildlife is typically modeled in a static environment that reflects a single snapshot in time. However, it has been shown that, when compared with dynamic connectivity models, static models can underestimate connectivity and mask important population processes. Therefore, including dynamism in connectivity models is important if the goal is to predict functional connectivity. We incorporated four levels of dynamism (individual, daily, seasonal, and interannual) into an individual-based movement model for black bears (Ursus americanus) in Massachusetts, USA. We used future development projections to model movement into the year 2050. We summarized habitat connectivity over the 32-year simulation period as the number of simulated movement paths crossing each pixel in our study area. Our results predict black bears will further colonize the expanding part of their range in the state and move beyond this range towards the greater Boston metropolitan area. This information is useful to managers for predicting and addressing human–wildlife conflict and in targeting public education campaigns on bear awareness. Including dynamism in connectivity models can produce more realistic models and, when future projections are incorporated, can ensure the identification of areas that offer long-term functional connectivity for wildlife.


Introduction
In a developing world, providing connectivity for wildlife is widely recognized as an important component of successful conservation and management plans [1,2]. Connectivity allows wildlife to access resource patches during day-to-day movements, facilitates dispersal, encourages gene flow, and allows for range shifts and range expansion [2,3]. Connectivity, the implementation of which ranges from a single road crossing structure to large landscape corridors, is being included as a key part of many conservation plans [4], but well-connected protected areas are still relatively rare and globally fall well short of the Aichi Target 11 of the Convention on Biological Diversity [5]. More connectivity plans are therefore encouraged-and one of the components of successful connectivity plan implementation is the inclusion of one or more focal wildlife species [4].
Connectivity for wildlife is typically modeled in static environments and estimated for a single snapshot in time (e.g., [6,7]). However, wildlife responses to landscapes are inherently dynamic, as Land 2020, 9,233 2 of 20 are the landscapes themselves (e.g., [8]). Dynamism can also exist among individuals. For example, individuals from the same population may have differing responses to the same landscape features [9,10], driven by differences among males and females, age classes, spatial locations, or temperament (e.g., boldness) [10][11][12]. There may also be differing behaviors and responses to landscape features throughout a diel period and across seasons [11,13]. In addition, interannual variability may exist due to disturbance-succession processes, climatic variation, human development, and other landscape changes [14]. All these types of dynamics may interact to affect functional habitat connectivity for a species.
Incorporating dynamism into wildlife connectivity models is important for several reasons. First, it allows for more realistic modeling and accurate estimates of connectivity. Dynamic connectivity models can increase our understanding of temporal changes in connectivity and population processes across space and time [15,16]. Simulation studies have found that static connectivity models can underestimate connectivity by an average of 30% compared with models that include space and time [16], and that dispersal-limited species may have a lower extinction threshold when modeled with dynamic landscapes compared to static ones [17]. Empirical studies have found that the importance of connectivity for wildlife population dynamics may be masked or misunderstood in static landscapes and therefore underestimated as a conservation need [15,18,19]. Furthermore, corridors derived from static connectivity models may differ substantially from those derived from dynamic models, resulting in a reduction of functional connectivity and a potential misallocation of conservation resources.
Dynamic connectivity can be modeled using different connectivity algorithms such as least-cost paths [20,21], resistant kernels [22], circuit theory [23,24], and, increasingly, individual-based movement models (IBMMs; e.g., [25,26]). In cost-path, kernel, and circuit theory approaches, resistance surfaces are typically derived for different temporal periods and connectivity is analyzed separately for each (e.g., [13,27]). The use of IBMMs is attractive because they provide a more mechanistic representation of functional connectivity by simulating the sequential movements of individuals in response to context-specific environmental conditions [28,29]. This allows simulated individuals to respond to a dynamic environment within a single model, resulting in a single summary output. Simulating sequential movements of individuals may be important for understanding movement response to future environmental conditions [30]-especially spatiotemporal dynamics of distributional shifts (e.g., range expansion, reintroductions) where the leading edges of species' distributions are conditional upon individual locations at a given point in time [26,31]. IBMMs can also be parameterized with empirical data allowing for realistic simulations of movement and behavior for a population of interest [32]. Furthermore, the emergence of higher order patterns (e.g., home range size) from individual-level behaviors allows simulated patterns to be compared against empirical estimates of those patterns as a form of model calibration [33,34]. Connectivity estimates from IBMMs can be obtained by calculating the number or density of individuals moving through every pixel on the landscape [25,26,35].
We incorporated individual, daily, seasonal, and interannual dynamics into an IBMM for black bears (Ursus americanus) in the Commonwealth of Massachusetts, USA, to predict statewide functional connectivity. The black bear range in the state is expanding toward the greater Boston metropolitan area (MassWildlife unpublished data: https://www.mass.gov/service-details/learn-about-black-bears) and our goal was to predict movement routes and colonization in this expanding range. We simulated present-day and future movement on an annual basis from the current time step to the year 2050 and incorporated projected changes in human development into our models. We summarized density of movement from the IBMMs and identified areas where movement is likely to be concentrated in the future. Because black bears in Massachusetts are expanding their range towards more populated areas of the state, this information will be useful to managers for predicting and responding to human-wildlife conflict, as well as for targeting education campaigns.

Data and Methods
All analyses were conducted in R software [36].

Black Bear Data
Empirical data for parameterizing the IBMMs were from GPS-collared female black bears in central and western Massachusetts. Black bears were fitted with Telonics Gen 3 or Gen 4 collars programmed to acquire a fix every 45 min. Collar data was collected from 34 individuals between 2012 and 2017, totaling 65 bear-years (collar duration, number of fixes, and other collar information are provided in Appendix A). The GPS data were assessed for positional accuracy and fixes from the Gen 3 collar were removed if the position dilution of precision (PDOP) was > 5. Fixes from the Gen 4 collars were removed if the fix was unresolved, if the fix had a PDOP > 5 and was uncertain or was a 2D fix, or if the fix had a PDOP > 20 and was a certain or 3D fix. This filtering was performed to minimize the locational error to < 100 m and resulted in a mean data loss of 3.93% [37,38]. More capture, handling, and data processing procedures are detailed in Zeller et al. [11].

Resistance Surfaces
The IBMMs were modeled across resistance surfaces representing the ease or difficulty of bear movement across each grid cell in the study area. Individual black bears have varying movement responses to landscape features during different seasons and different times of day, which may translate into different resistance surfaces for each season/diel period [11]. Bears in our study area hibernate in the winter, and during the nonhibernation period they move the least in the spring and the most in the summer [11]. Movement response to different land cover types also changes seasonally [11]. During a 24-h period, bears tend to move more during the day than at night. However, daily movement patterns in our study area may change with season and for different land cover types [11]. To incorporate these movement dynamics into the IBMMs, we used six previously developed resistance surfaces reflecting bear movement: (1) spring/day, (2) spring/night, (3) summer/day, (4) summer/night, (5) fall/day, and (6) fall/night (Appendix B).
The resistance surface for each season/diel period was derived through the following four-step procedure (full details in Zeller et al. [11]). First, a step selection function was developed for each individual bear using only data for that season/diel period. Second, the step selection function was used to project the probability of movement across the study area. Third, the movement probability surfaces were combined across bears in a spatially-weighted approach following Osipova et al. [10] so that spatial differences in bear behavior and movement were maintained across the study area. Specifically, we created a Euclidean distance surface for each bear during each season/diel period from the centroid of that bear's season/diel minimum convex polygon home range. We then took the inverse distance values and normalized them to sum to 1 across each surface and used the resulting surface to weight the individual movement probabilities at each grid cell. We summed the weighted movement probability surfaces to create a combined movement probability surface. Fourth, we created a resistance surface by calculating the linear inverse of the combined movement probability surface. These resistance surfaces have low resistances in grid cells with a high probability of movement and high resistances in grid cells with a low probability of movement. Resistance surfaces were rescaled from 0-1 to 1-10 for use in the IBMMs.

Step Length Distributions
The movement kernels for our IBMMs were based on empirical step length distributions from the collared bears. Bears in our study area move most in the summer and during daytime periods compared with other seasons and nighttime periods [11]. To reflect these differences, we calculated movement distributions for the same six seasons/diel periods as the resistance surfaces. Black bear step lengths followed a Pareto distribution and the shape and scale parameter of the Pareto distribution were estimated for each season/diel period with the gpd.fit function of the gPdtest package ( [39]; Appendix D; Figure A2).

Individual-based Movement Model
The Massachusetts Division of Fisheries and Wildlife estimates there are approximately 5000 black bears in the western and central parts of the state (unpublished data). The highest density of bears is to the west of the Connecticut River, with densities diminishing to the east. An 'expanding' range has also been identified, which has a low number of recently established bears ( Figure 1). To reflect these density changes, and the fact that we are only modeling female bears, we distributed 2000 simulated individuals to the west of the Connecticut River, 990 in the established range to the east of the Connecticut River, and 10 in the expanding range. Start points of these 3000 individuals were sampled probabilistically on the spring/day probability of movement surface ( Figure 1). We assume that, by simulating females from home range data, our results will be conservative in terms of movement propagation. However, the focus on female movement is important at the expanding edge of ranges since they determine population establishment and whether reproduction is occurring in these areas.

Individual-based Movement Model
The Massachusetts Division of Fisheries and Wildlife estimates there are approximately 5,000 black bears in the western and central parts of the state (unpublished data). The highest density of bears is to the west of the Connecticut River, with densities diminishing to the east. An 'expanding' range has also been identified, which has a low number of recently established bears ( Figure 1). To reflect these density changes, and the fact that we are only modeling female bears, we distributed 2,000 simulated individuals to the west of the Connecticut River, 990 in the established range to the east of the Connecticut River, and 10 in the expanding range. Start points of these 3,000 individuals were sampled probabilistically on the spring/day probability of movement surface ( Figure 1). We assume that, by simulating females from home range data, our results will be conservative in terms of movement propagation. However, the focus on female movement is important at the expanding edge of ranges since they determine population establishment and whether reproduction is occurring in these areas. We used a resistance kernel-based approach [22] to simulate individual bear movements that combined our empirical distributions of step lengths with our dynamic resistance surfaces to estimate the probability of an individual moving to any given pixel based on the intervening resistance values. At the beginning of each time step we first built a resistance kernel at the start point of each individual. We determined the individual's maximum movement potential (i.e., the extent of the kernel) during a single time step in resistance distance units defined as the 97.5th quantile of the Pareto distribution derived from the parameters of the appropriate season/diel period. We hereafter refer to this quantity as an individual's 'bank account. ' We then calculated the resistance distance of every grid cell in the study area from the start point using the rawspread function in the GRIDPROCESS package [40] (Figure 2a). The extent of the resistance kernel stops when the bank account reaches a value of 0. In areas of low resistance, the bank account will be spent more slowly allowing the kernel to spread farther than in areas of high resistance. The output from the rawspread function provides a We used a resistance kernel-based approach [22] to simulate individual bear movements that combined our empirical distributions of step lengths with our dynamic resistance surfaces to estimate the probability of an individual moving to any given pixel based on the intervening resistance values. At the beginning of each time step we first built a resistance kernel at the start point of each individual. We determined the individual's maximum movement potential (i.e., the extent of the kernel) during a single time step in resistance distance units defined as the 97.5th quantile of the Pareto distribution derived from the parameters of the appropriate season/diel period. We hereafter refer to this quantity as an individual's 'bank account. ' We then calculated the resistance distance of every grid cell in the study area from the start point using the rawspread function in the gridprocess package [40] (Figure 2a). The extent of the resistance kernel stops when the bank account reaches a value of 0. In areas of low resistance, the bank account will be spent more slowly allowing the kernel to spread farther than in Land 2020, 9, 233 5 of 20 areas of high resistance. The output from the rawspread function provides a surface of proximity values to the start point. Therefore, we subtracted the values from one to obtain a resistance distance kernel.
This point then became the new start point and the procedure was repeated. We simulated steps for each bear at 45-minute intervals from April 1 to November 15, the average den emergence and entry dates for our sample of bears, for a total of 7,296 steps per individual. We started simulations on the spring/day resistance surface. We then alternated between the night resistance surface and Pareto parameters and the day resistance surface and Pareto parameters every 16 time-steps. We replaced the resistance surface and Pareto parameters from the spring with those of the summer at the 2,369th time step and those of the summer with those of the fall at the 4,172th time step. These changes corresponded with ecologically identified seasons for our collared population of black bears and with the dates of the GPS data used to estimate the seasonal resistance surfaces (see [11]). Figure 2. Procedure for simulating one step for the black bear individual-based movement models (IBMMs). From a start point on the resistance surface, (a) we applied and took the inverse of the rawspread function [40] which calculates resistance distance from the start point. We then (b) applied the Pareto density function of the appropriate diel period and season to the resistance distance surface to create a movement probability kernel and transformed the kernel to sum to 1. We then (c) probabilistically sampled a destination (end) point on the movement probability kernel.
We also describe the IBMM procedure with the standardized Overview, Design Concepts, and Details (ODD) protocol recommended by Grimm et al. [42,43] in Appendix D. Figure 2. Procedure for simulating one step for the black bear individual-based movement models (IBMMs). From a start point on the resistance surface, (a) we applied and took the inverse of the rawspread function [40] which calculates resistance distance from the start point. We then (b) applied the Pareto density function of the appropriate diel period and season to the resistance distance surface to create a movement probability kernel and transformed the kernel to sum to 1. We then (c) probabilistically sampled a destination (end) point on the movement probability kernel.
We converted these resistance distance values to relative probability values using the Pareto distribution parameters of the appropriate season/diel period. This produced our final movement probability kernel (Figure 2b). The destination point (end point) for that time step was sampled from the movement probability kernel with the strata function of the sampling package (Figure 2c; [41]). This point then became the new start point and the procedure was repeated. We simulated steps for each bear at 45-min intervals from April 1 to November 15, the average den emergence and entry dates for our sample of bears, for a total of 7296 steps per individual. We started simulations on the spring/day resistance surface. We then alternated between the night resistance surface and Pareto parameters and the day resistance surface and Pareto parameters every 16 time-steps. We replaced the resistance surface and Pareto parameters from the spring with those of the summer at the 2369th time step and those of the summer with those of the fall at the 4172th time step. These changes corresponded with ecologically identified seasons for our collared population of black bears and with the dates of the GPS data used to estimate the seasonal resistance surfaces (see [11]).
We also describe the IBMM procedure with the standardized Overview, Design Concepts, and Details (ODD) protocol recommended by Grimm et al. [42,43] in Appendix D.

Model Validation and Calibration
We calculated the step lengths and minimum convex polygon home range sizes of our simulated individuals and compared them with our empirical sample of bears (described above) to determine how well our IBMMs represented movement and space use of black bears in our study area. Our first simulations resulted in comparatively short step lengths and small home range sizes (Appendix E). This is expected because the probability distribution (i.e., our Pareto distributions) used with a resistance kernel represents the maximum movement potential of an individual in a theoretical landscape with resistance = 1. In contrast, the observed step lengths used to estimate our Pareto distribution parameters occurred within real heterogeneous landscapes with varying resistances. Parameterizing a resistance kernel with an observed distribution of step lengths will therefore underestimate an individual's movement potential and the extent of the movement probability kernel when applied to real resistant landscapes. Because we were unable to estimate the distribution of black bear step lengths in a nonresistant landscape, we increased the scale parameter of the Pareto distributions by 10% to allow our resistance kernels to more closely approximate a black bear's maximum daily movement potential. This resulted in more realistic step lengths and home range sizes. Increasing the scale parameter above 10% resulted in unrealistically large home range sizes (Appendix E).

Future Projections
We used data from the Designing Sustainable Landscapes project (DSL, [44]) to incorporate future landscape change and predict bear movement on the expanding edge of their range in Massachusetts. The DSL urban growth model simulated human development across the northeastern US at 10-year time steps from 2020 to 2080 [45]. Resultant development categories were low-, medium-, and high-density development. We used the years 2020-2050 for our analysis. In the current time step, low-, medium-, and high-density development measured 5.5%, 3.7%, and 1.4% of the study area. In 2050, these development categories increased to 6.8%, 4.7%, and 2.5%, a percent increase of 24%, 27%, and 79% respectively. Over this time period, the human population in the state is projected to grow by at least 6% [46].
Most of our individual bear step selection function (SSF) models had development variables that matched the DSL categories. In those individuals, we used the DSL output directly in our future projections. However, approximately 30% of bear SSF models identified percent impervious surface as the most influential development variable. Therefore, we also transformed the DSL projections into percent impervious surfaces by reclassifying the low-, medium-, high-density development to 25%, 50%, and 100% impervious surface respectively. We combined these new percent impervious surfaces with our current (2019) percent impervious surface layer for each of the four future projection time steps. We predicted the SSF models for each bear for each season/diel period with these future projections, combined the projected probability of movement surfaces, and transformed them into resistance surfaces as described above.
We considered our first year of simulation to be 2019 and we simulated bear movement annually from 2019-2050. We used the last point from each simulated individual from the previous year as the start point for simulations in the next year. The 2020 resistance surfaces were used for the years 2020-2029, the 2030 resistance surfaces for 2030-2039 and so on. We summed the number of individual simulated paths crossing a pixel for each year and across all years to develop a movement density surface. Paths were derived by connecting all consecutive simulated points for each individual in each year.  Therefore, step lengths in our simulated bears are shorter, but home ranges are slightly larger than the empirical data. However, we did not include memory, territoriality, or interactions with conspecifics in our IBMMs and would therefore expect slightly shorter step lengths and larger home range sizes. Movement patterns and space use were visually similar among our simulated and empirical bears (e.g., Figure 3).

Our
Land 2020, 9, x FOR PEER REVIEW 7 of 21 the empirical data. However, we did not include memory, territoriality, or interactions with conspecifics in our IBMMs and would therefore expect slightly shorter step lengths and larger home range sizes. Movement patterns and space use were visually similar among our simulated and empirical bears (e.g. Figure 3).

Discussion
We used an IBMM to model dynamic habitat connectivity for female black bears in Massachusetts, USA and found that black bears have high future movement potential within and beyond their expanding range in the state. We were able to incorporate individual, daily, seasonal, and interannual dynamics into our IBMM. Individual black bear differences were captured in the resistance surfaces used in the IBMMs by using the spatially-weighted approach presented by Osipova et al. [10]. This approach allows responses to landscape features to change depending on location. For example, bears in more developed areas of the state show lower avoidance of human development than more rural bears [11] and these differences are captured in these weighted surfaces. We also used resistance surfaces and movement parameters for two times of day and three seasons to capture movement and behavioral differences during these different periods. Lastly, we incorporated interannual dynamics by modeling IBMMs annually and incorporating future human development projections at 10-year time steps. By combining these four levels of dynamism, we were able to create realistic movement models and predict future connectivity.
IBMMs are well suited for modeling dynamic connectivity. As we demonstrated, they allow individuals to respond to changeable landscape features at different time steps and can integrate daily, seasonal, and interannual changes within a single model. More fundamentally, IBMMs allow population-level responses to emerge from the decisions made by thousands of individuals in response to local environmental conditions. This provides a more mechanistic approach for estimating dynamic landscape patterns of resource use and movement and, ultimately, functional connectivity. [32,47]. IBMMs have been used in this way to model dispersal corridors for 55 generations of little owls (Athene noctua) in Germany and Switzerland [26], connectivity under different land management scenarios for bighorn sheep (Ovis canadensis) in British Columbia, Canada [48], and connectivity for coral reef larvae using daily ocean current velocity data and different dispersal capabilities [25].

Discussion
We used an IBMM to model dynamic habitat connectivity for female black bears in Massachusetts, USA and found that black bears have high future movement potential within and beyond their expanding range in the state. We were able to incorporate individual, daily, seasonal, and interannual dynamics into our IBMM. Individual black bear differences were captured in the resistance surfaces used in the IBMMs by using the spatially-weighted approach presented by Osipova et al. [10]. This approach allows responses to landscape features to change depending on location. For example, bears in more developed areas of the state show lower avoidance of human development than more rural bears [11] and these differences are captured in these weighted surfaces. We also used resistance surfaces and movement parameters for two times of day and three seasons to capture movement and behavioral differences during these different periods. Lastly, we incorporated interannual dynamics by modeling IBMMs annually and incorporating future human development projections at 10-year time steps. By combining these four levels of dynamism, we were able to create realistic movement models and predict future connectivity.
IBMMs are well suited for modeling dynamic connectivity. As we demonstrated, they allow individuals to respond to changeable landscape features at different time steps and can integrate daily, seasonal, and interannual changes within a single model. More fundamentally, IBMMs allow population-level responses to emerge from the decisions made by thousands of individuals in response to local environmental conditions. This provides a more mechanistic approach for estimating dynamic landscape patterns of resource use and movement and, ultimately, functional connectivity. [32,47]. IBMMs have been used in this way to model dispersal corridors for 55 generations of little owls (Athene noctua) in Germany and Switzerland [26], connectivity under different land management scenarios for bighorn sheep (Ovis canadensis) in British Columbia, Canada [48], and connectivity for coral reef larvae using daily ocean current velocity data and different dispersal capabilities [25]. Modeling connectivity in a static environment may mask the importance of connectivity for maintaining wildlife population dynamics and may result in identifying corridors that are only periodically functional [15,16,18]. Because of the vast resources needed to conserve some wildlife corridors, it is not only important to model these corridors in a way that reflects behavioral and landscape dynamics, but also to ensure the longevity of these corridors into the future. To this end, we included future development projections in our IBMM and summarized connectivity at our final time step. If connectivity is being modeled at each time step separately (as with cost path, resistant kernel, or circuit theory approaches; e.g., [10,13]), identifying areas of consistent connectivity regardless of development or climate change is key to prioritizing conservation for long-term persistence of wildlife movement (e.g., Jennings et al. this issue).
Our IBMM was built using empirically derived resistance surfaces and movement step lengths and simulated reasonably realistic patterns of black bear movement and space use. Nevertheless, we recognize several areas in which our model could be improved. Our empirical population was comprised of females and all observed movements were within-home-range movements. Including males, especially young dispersers, would have likely resulted in longer step lengths and more movement. Female bears do not typically disperse far from their natal ranges [49,50]. However, this philopatry may be density dependent and females may disperse further in recently established black bear range [51]. Given that we could only model female home range movement, our model is a conservative approximation of movement and connectivity for our study area. Incorporating male bears and data from dispersing individuals would increase our insights into bear movement for the entire Massachusetts black bear population. A more direct incorporation of individual differences could also have been included in the IBMM by calculating step lengths and Pareto parameters for each individual bear and then assigning these parameters to each simulated individual based on a random draw. Other improvements to our model include incorporating births, natural deaths, and hunting take as well as different behavioral states, territoriality, and motivation [52]. Incorporating these parameters may have resulted in step lengths and home range sizes that were further aligned with our empirical bears, especially if step lengths were modeled as a function of land cover type. A more accurate understanding of spatial variation in current black bear density across our study area, particularly in the expanding portion of their range, may also improve our understanding of future landscape use of black bears. However, given the computational demands of our IBMM, we were unable to incorporate these model additions and are confident that our results are a reasonable approximation of black bear movement and connectivity in Massachusetts.

Conclusions
Incorporating dynamism into wildlife connectivity models increases their complexity, but also their realism. Using GPS collar data on black bears as well as future projections of human development, we were able to include four levels of dynamism-individual, daily and seasonal, and interannual dynamics -within a single connectivity model. The increasing availability of animal tracking data as well as historical and future geospatial time series data increases our ability to model dynamic connectivity for wildlife. Ensuring realism and future longevity of wildlife corridors will be important for the success of corridor conservation initiatives in providing long-term functional connectivity for wildlife. With our IBMM, we predicted black bears will increase their colonization within the current boundary of their expanding range in Massachusetts and will move beyond this boundary into uninhabited areas of the state. The model developed here might be combined with other information to determine the location, amount, and spatial configuration of public education needed to reduce human-bear conflicts (e.g., [53]), or where ordinances may be implemented to remove food attractants. The model may also be used to identify road-crossing hotspots that may be the target of mitigation measures. As black bears are expected to continue to move towards and colonize more populated areas around the greater Boston metropolitan area, information on where this movement Appendix A Table A1. Female bear ID, year, and date range of GPS collars used in the analyses. Number of fixes outside the hibernation period are the number of GPS locations acquired after cleaning data for inaccurate fixes. Number of steps used in the step selection function analysis for each season are also provided. Steps were only used in the analysis if the start and end points of a step were at consecutive 45-min intervals, indicating that there were no missing fixes between the two locations. Steps were only used for a season if a full season's worth of data were available for that bear year.

ID
Year

Entities, state variables, and scales
The model was comprised of one lower hierarchical level-the individual-and one higher hierarchical level-the environment. Individuals (i.e., entities) were all female black bears of breeding age in the state of Massachusetts. Movement characteristics of individuals were derived from within home-range movements of our sample of GPS collared bears in the state. Because 45-min GPS data were used, each time step in the model was 45-min. For a single year (i.e., the non-hibernation period: April 1st to November 15th), there were 7296 time steps. These dates correspond to the average den exit and entry rates for the sample population. Three thousand individuals were simulated for each year of the simulation.
The individuals moved through their environment, in which each 30 m pixel was assigned a resistance to movement value from 1-10. The model used six different resistance surfaces: spring/day, spring/night, summer/day, summer/night, fall/day, fall/night. The resistance surfaces were derived from step selection functions (SSFs) that predicted the probability of movement for each pixel in the study area. GPS data from 76 bear-years were combined with the following environmental variables to create the SSFs: land cover, ruggedness, slope, roads, development, and water. Details on the SSF and probability of movement surfaces are provided in Zeller et al. [11]. The linear inverse of the probability of movement surfaces were re-scaled to a range of 1-10 to obtain the resistance surfaces.
The first year of the simulation was 2019 and was considered to be the current environment. To incorporate future development, we used data from the Designing Sustainable Landscapes Project urban growth model, which provides development projections at 10-year intervals starting in 2020 [45]. In our 2020 model, these projected development variables were incorporated into the resistance surfaces to reflect a growing human footprint. The 2020 resistance surfaces were used for 10 years of the simulation. This process was repeated again for 2030, 2040, and 2050. The extent of the modeling environment was the entire state of Massachusetts plus a 10 km buffer to prevent edge effects.

Process overview and scheduling
The movement model began in the year 2019 and 7296, 45-min time steps were run for each year. For each bear in each year, simulations started on the spring day surface. The resistance surfaces and movement parameters (determined by the Pareto distribution, see below) were changed within the modeling environment at the corresponding time steps. Every 16 time steps the night surfaces and movement parameters replaced the day ones. The summer surfaces and movement parameters replaced the spring ones at step 2369 and the fall surfaces and movement parameters replaced the summer ones at step 4172. This process was repeated annually for 31 more years to the year 2050. The years were run sequentially and the last location for an individual in a given year was used as the starting location for that same individual the following year.

Design concepts
Basic principles: The model used individual movements to simulate population-level patterns of movement and connectivity across the study area. The model incorporated the principle that individuals will move in a manner that minimizes the costs of movement (e.g., energetic costs, mortality risks) and used resistance surfaces to quantify these costs. Furthermore, the model explicitly accounted for changes movement costs among individuals, between diel periods (i.e., day and night), and across seasons. The model did not directly evaluate changes in movement costs but rather cumulative effects on population-level movement processes.
Emergence: The movement of thousands of individuals over many decades allowed for the emergence of population-level movement and connectivity across a large study area. Using movement rules that change an individual's response to the environment during different seasons and times of day while incorporating a changing environment allowed for the realistic emergence of concentrated movement pathways connecting habitat patches.
Adaptation: Individual movements were simulated as a function of landscape resistance such that individuals were more likely to move across low-resistance pixels. Resistance was inversely related to habitat suitability as estimated through step selection functions such that low-resistance pixels represented more suitable habitats (e.g., better foraging opportunities, suitable shelter sites), lower mortality risk (e.g., reduced human influence, lower risk of road mortality), or some combination thereof.
Objectives: The endpoints following each time step were probabilistically sampled using resistance surfaces so that individuals were more likely to cross pixels with low resistance values.
Learning: None Prediction: By simulating future individual movements, the model predicted how individuals may move in response to environmental features and growing human development. This allowed for the identification of concentrated areas of movement and connectivity along the expanding edge of the black bear's range in Massachusetts and areas of potential range expansion.
Sensing: Individuals could sense (i.e., had a nonzero probability of moving to) pixels within some distance of the starting point at each time step that was based on observed step lengths for a given time of day and season. This maximum distance was defined as the 97.5th quantile of a Pareto distribution fit to observed step lengths across bears for a given time of day and season. However, the probability of moving to a given pixel decreased with increasing distance from the starting point according to a cumulative cost distance function as reflected in the resistance kernel (see below).
Interaction: None Stochasticity: Each simulated step length was sampled probabilistically from a resistance kernel surface (see below).
Collectives: None Observation: The location of the endpoints of all movement steps for each individual were the basis for observations from the model. From these points, step lengths and home range sizes were calculated. Paths were also constructed for each individual in each year by connecting consecutive points. Paths were used to calculate our connectivity metric which was the number of paths crossing a pixel.

Initialization
The study area was divided into three sections: west of the Connecticut River, east of the Connecticut River, and the expanding range. Unpublished data from the Massachusetts Division of Fisheries and Wildlife indicates that bear densities are highest in the western section, lower in the eastern section, with just a few bears currently living in the expanding range. To reflect these differences, 2000 individuals were initialized in the western range, 990 in the eastern range, and 10 in the expanding range.
The location of the start points was determined by sampling probabilistically on the spring/day movement surfaces so that start points tended to be in areas with higher movement probabilities. At the first time step in the first year, movement was initiated from each of these locations on the spring/day resistance surface with movement parameters drawn from the Pareto Distribution parameterized with GPS data from spring daytime locations of the collared sample of bears.

Input
There were three data inputs to the model. One was the initial start locations of the 3000 individuals. The second were the season-and time-specific resistance surfaces. The third were the season-and time-specific parameters for estimating the Pareto Distribution. To derive these Pareto parameters, the GPS collar data from our sample of female black bears were parsed into six data sets matching the seasons and times of day used in the SSFs and the resistance surfaces. From these GPS data, step lengths were calculated, a Pareto distribution was fit to the step length distribution, and the two parameters of the Pareto, shape and scale, were estimated.

Submodels
The procedure for simulating each movement step for an individual was as follows: A resistance kernel was built around the start point. The shape and spread of this kernel was determined by the Pareto distribution for that season/time of day (after model calibration the scale parameter for the Pareto was increased by 10% to allow for more realistic simulated movements, see main text and Appendix E). The kernel spread was bounded at the 97.5th quantile of the Pareto distribution and this boundary determined the maximum movement potential for a step. The rawspread function in the gridprocess R package [36,40] was used to calculate the resistance distance from the start point to every grid cell in the study area on the resistance surface for that season/time of day. The spread stopped when the boundary value was reached. This output was subtracted from one to obtain a resistance distance kernel around the start point.
The Pareto distribution parameters were again used to convert the resistance distance values around the start point to relative probability values. The destination point was then sampled from this probability surface with the strata function of the sampling R package [41]. The start point and end point make a movement step and the procedure was then repeated with the end point becoming the start point for the next step. Table A2. Minimum convex polygon home range sizes and step lengths of our sample of empirical and simulated bears. For the simulated bears, we calibrated the individual-based movement models by increasing the scale of the Pareto distribution. We selected a 10% increase in the scale of the Pareto distribution for our final models since this resulted in the most realistic combination of home range sizes and step lengths for female black bears in Massachusetts.