Predatory Behavior is Primary Predictor of Movement of Wildland-Urban Cougars

While many species have suffered from the detrimental impacts of increasing human population growth, some species, such as cougars (Puma concolor), have been observed using human-modified landscapes. However, human-modified habitat can be a source of both increased risk and increased food availability, particularly for large carnivores. Assessing preferential use of the landscape is important for managing wildlife and can be particularly useful in transitional habitats, such as at the wildland-urban interface. Preferential use is often evaluated using resource selection functions (RSFs), but RSFs do not adequately account for the habitat available to an individual at a given time and may mask conflict or avoidance behavior. Contemporary approaches to incorporate landscape availability into the assessment of habitat preference include spatio-temporal point process models, step-selection functions, and continuous-time Markov chain (CTMC) models; in contrast with the other methods, the CTMC model allows for explicit inference on animal movement. We used the CTMC framework to model speed and directionality of movement by a population of cougars inhabiting the Front Range of Colorado, U.S.A., an area exhibiting rapid population growth and increased recreational use, as a function of individual variation and time-varying responses to landscape covariates. The time-varying framework allowed us to detect temporal variability that would be masked in a generalized linear model. We found evidence for individual- and daily temporal-variability in cougar response to landscape characteristics. Distance to nearest kill site emerged as the most important driver of movement at a population-level. We also detected seasonal differences in average response to elevation, heat loading, and distance to roads. Motility was also a function of amount of development, with cougars moving faster in developed areas than in undeveloped areas.

In contrast to many resource selection studies, one of the primary goals of continuous-  The flexibility of the CTMC framework can account for time-varying responses to landscape 115 drivers by allowing coefficients to vary temporally , and it can also be 116 implemented in a Bayesian hierarchical framework, allowing for inference on individual-and 117 population-level drivers. 118 Given the increasing potential for human-wildlife conflict as development permeates 119 rural and wildland areas along the Front Range and elsewhere in the West, we sought to extend 120 previous work by explicitly modeling cougar movement to identify key drivers of their behavior, 121 and in doing so, better understand their use of the wildland-urban landscape in both space and  Figure 1). We were interested in examining seasonal differences in 142 cougar movement due to temporal variability in the landscape-level covariates. For example, we 143 expected a stronger response to prey-based covariates in summer months, because mule deer 144 fawns, a primary prey source for cougars, are born in June (Pojar and Bowden 2004) and fawns 145 are at a disproportionately high risk for predation (Hornocker 1970 (Vectronics GmbH, Berlin, Germany) programmed to obtain fixes every three hours.

150
Our study area comprised a 2,700 km 2 region in the Colorado Front Range to the north-151 west of Denver ( Figure 1). The study area consisted of a matrix of private (43%) and public 152 (57%) land (Blecha 2015 We used a Bayesian hierarchical CTMC model to evaluate drivers of cougar movement; this 162 model is an extension of the model proposed by Hanks et al. (2015) and allows for inference on 163 movement rates and directional bias, as opposed to use, in continuous time and discrete space. 164 The initial step in the CTMC framework is to estimate a continuous movement path from the 165 observed data points. We used the functional movement model proposed by Buderman et al. 166 (2016) to account for observation error and predict locations every ten minutes for the selected 167 two weeks of each month. A random subset of paths from the posterior predictive distribution of 168 the movement model were spatially discretized to a latent variable formulation with a cell size of 169 100 m 2 , which was the largest cell size among the available covariates. cell, the probability that they move to a particular neighboring grid cell (directionality) is the 175 ratio between the movement rate into that cell and the movement rate out of the preceding cell.

176
Therefore, higher proportional rates indicate directional bias in movement. Thus, movement rate 177 parameters, which are a function of covariates (i.e., landscape variables that correspond to the 178 position of the cell on the landscape), control both motility and directionality. Hanks et al. (2015) 179 showed that the likelihood of the CTMC model (i.e., the product of the motility and directional 180 components) for movement can be expressed as a Poisson GLM using a latent variable 181 formulation.

182
In the latent variable formulation, each transition corresponds to four data points (the four    period. This is similar to the idea proposed by Johnson (1980), where preference was determined 209 by comparing some measure of usage and availability of a landscape feature on an individual 210 basis.

211
On average, we expected cougars to respond similarly to landscape covariates. Therefore,  we need to select a priori by scaling the 2 term to each parameter, and ~ ( , 2 ). Both 239 2 and 2 serve as regularization terms, where 2 is selected a priori and 2 = 0.1.

240
Finally, to assess whether males and females exhibited different amounts of temporal 241 variation in their response to potential movement drivers, we fit the GLM and GAM models to (increasing distance, orient away from a feature). All rasters were aggregated to a 100 m 2 266 resolution, which is within the distance that a cougar might typically move over a ten-minute 267 interval (Dickson et al. 2005). 268 We hypothesized that a number of landscape covariates may contribute to transition rates 269 and directional bias of cougars: mule deer (Odocoileus hemionus) utilization (as a proxy for 270 availability), distance to nearest potential kill site, distance to nearest structure, distance to 271 nearest road, elevation, heat insolation load index, and topographic wetness. We also used an We calculated distance to nearest structure (m) as the Euclidean distance to the nearest 298 man-made roofed structure (Blecha 2015). Distance to road was calculated using major roads 299 data (i.e., a major highway primarily for through traffic usually on a continuous route and streets 300 whose primary purpose is to serve the internal traffic movement within an area) obtained from 301 Colorado Department of Transportation. Due to increased human activity around structures and 302 roads, we expected cougars to move faster when closer to roofed structure and distance to nearest

344
There was no detectable effect of many of the landscape covariates on average motility or 345 directionality at a population-level ( Figure 2). However, distance to potential kill site emerged as 346 the primary driver of both motility and directionality in the GLM framework (Figures 2, 3). As despite a non-significant population-level response (Figures 3, 4).

366
Our results suggest that distance to nearest potential kill site was also the predominant 367 motility and directionality driver in the diel time-varying framework (H-GAM; Figures 5 and 6). performance across all months, whereas the GLM performed the worst (Figure 9). The models 400 that were a mixture of a GAM and GLM, varying by sex, were generally equivalent (Figure 9).

401
The largest difference between the two sex-varying models was observed in October, when the 402 better of the two models included time variation for males and no time variation for females 403 (Figure 9).

405
The strong observed response to distance to nearest potential kill site is likely due to returning 406 visits to the carcass and unmeasured fine-scale covariates related to landscape features that 407 increase the likelihood of a successful hunting attempt. Blecha (2015) found that kill sites, 408 compared to preceding locations, occurred more frequently in areas with higher housing 409 densities and lower topographic positions, such as drainage areas, despite drainage areas having 410 lower prey availability. We did not measure hunting success, but we did find that cougars moved 411 toward lower elevations at dusk, when cougars are likely to hunt or return to a carcass. We found

759
Additional supporting information may be found in the online version of this article at the 760 publisher's website.

Appendix A
We spatially discretize a posterior predictive continuous path from the movement model to the resolution of the rasters of interest and decomposed into two elements: c, a state sequence consisting of the sequential grid cells (of N possible grid-cells) visited by the individual, and τ , a vector of residence times that describe how long the individual spent in each grid cell. We describe the cell sequence in terms of the transition rates α where α i j is a parameter controlling movement from cell i to cell j that can be a function of spatial covariates: If we designate t as the t th observation in the state-sequence (t ∈ T ), then the residence time τ t is exponentially distributed with a rate equal to the sum of all α ij (the total transition rate): α ij e −τ k N j=1 α ij . (2) In the above notation, [τ t |β] represents the probability distribution of the random variable τ t given the parameters β; this notation will appear again. We assume that it is impossible to move directly to non-neighboring cells, and therefore α ij = 0 for all j except for the cells adjacent to cell i.
When an individual transitions to a neighboring cell, the probability of transitioning to cell c t+1 = l is Assuming independence, the joint likelihood is the product of the transition probabilities and the residence times in the state sequence c is: Using a latent variable representation, where and then the product of [z ctk , τ t |β] over all N is proportional to the likelihood of the observed transition: The above process is parameterized with a single realization from the movement model, however, we have failed to account for the uncertainty in the animal's path. To avoid computational storage limitations, we use multiple imputation to account for the uncertainty in the path and make approximate posterior predictive inference on transition rates.