Drivers of polar bear behavior and the possible effects of prey availability on foraging strategy

Background Change in behavior is one of the earliest responses to variation in habitat suitability. It is therefore important to understand the conditions that promote different behaviors, particularly in areas undergoing environmental change. Animal movement is tightly linked to behavior and remote tracking can be used to study ethology when direct observation is not possible. Methods We used movement data from 14 polar bears (Ursus maritimus) in Hudson Bay, Canada, during the foraging season (January–June), when bears inhabit the sea ice. We developed an error-tolerant method to correct for sea ice drift in tracking data. Next, we used hidden Markov models with movement and orientation relative to wind to study three behaviors (stationary, area-restricted search, and olfactory search) and examine effects of 11 covariates on behavior. Results Polar bears spent approximately 47% of their time in the stationary drift state, 29% in olfactory search, and 24% in area-restricted search. High energy behaviors occurred later in the day (around 20:00) compared to other populations. Second, olfactory search increased as the season progressed, which may reflect a shift in foraging strategy from still-hunting to active search linked to a shift in seal availability (i.e., increase in haul-outs from winter to the spring pupping and molting seasons). Last, we found spatial patterns of distribution linked to season, ice concentration, and bear age that may be tied to habitat quality and competitive exclusion. Conclusions Our observations were generally consistent with predictions of the marginal value theorem, and differences between our findings and other populations could be explained by regional or temporal variation in resource availability. Our novel movement analyses and finding can help identify periods, regions, and conditions of critical habitat. Supplementary Information The online version contains supplementary material available at 10.1186/s40462-022-00351-4.

Local hour -Behavior exhibits circadian patterns associated with resource availability (i.e., seal distribution and behavior) -Access to seals is greatest during mid-day when haul-out behavior is most frequent -Low ambient light increases the reliance on olfaction The drift state will peak 0:00 h, olfactory search will be most common in the early morning and late evening, and ARS will be most common during the day Solar radiation " " Sun altitude " " Wind velocity -Low wind speed do not facilitate scent dispersion and hinder olfaction -High wind speed increases turbulence, fragments odor plumes, and hinders olfaction -High wind speed decrease seal haul-out behavior and decrease access to prey The drift state will be most common during low and high wind speeds and olfactory search and ARS will peak at moderate speeds Snow depth -Increasing snow depth increases the energetic cost of locomotion and decreases the effectiveness of olfaction Olfactory search will be negatively correlated with snow depth, and drift and ARS states will be positively correlated with snow depth Total precipitation -Precipitation hinders olfactory and visual search Olfactory search will be negatively correlated with snow depth, and drift state will be positively correlated with snow depth ARS may be negatively correlated or peak at intermediate levels of precipitation Sea ice concentration -Low ice concentration increases cost of locomotion.
-Low ice concentration decreases accessibility of prey -Still-hunting is most favorable at high ice concentration when seal access to air is more limited Olfactory search and drift will be positively correlated with ice concentration, and ARS will be negatively correlated with ice concentration.
Sea ice drift is driven by wind forcing, ocean currents, the Coriolis force, and internal ice stress. Due to the compressive strength of sea ice, local drift is typically affected by average force applied to it at a large scale. To allow for this, we summarized wind velocities at four buffers (0 km, 50 km, 100 km, 150 km). The buffer represents the radius around each GPS location from which wind velocities were extracted. The wind velocities of all raster cells that fell within the buffer were averaged to obtain the an estimate for the GPS point. For the buffer of 0 km, wind speed was estimated at the GPS location using bilinear interpolation.
The buffer with the lowest AIC was 100 km and was used to predict drift. Wind velocities were obtained from the ERA5 meteorological reanalysis project, which provides hourly global analysis fields at a 31 km resolution [Hersbach et al., 2020] and were extracted using the raster package [Hijmans and van Etten, 2016].
The motion of the GPS tracks was described using two data streams: step length l t ∈ (0, ∞) (the distance between consecutive locations) and turning angle ϕ t ∈ (−π, π] (change in bearing between consecutive steps), where t ∈ 1, ..., T represents the time of the first GPS location in a step [Langrock et al., 2012, McClintock andMichelot, 2018].
Step length l t and turning angle ϕ t were modelled as slope-intercept functions of wind velocities. We fit Weibull and gamma distributions to step length and von Mises and wrapped Cauchy distributions to turning angle and selected the combination of distributions that had the lowest AIC [McClintock et al., 2020]. Our final model for drift assumed that step length follows a Weibull distribution: where κ (l) t > 0 is the shape parameter and λ  et al., 2020]. The shape and scale parameters were each defined as slope-intercept functions of wind speed r t ≥ 0 as follows: where is β {1,2,3,4} ∈ R are the intercept and slope coefficients to be estimated.
According to AIC, our best model defined turning angle ϕ t using the von Mises distribution. Due to the Coriolis effect, sea ice drift is approximately -20 • relative to the wind direction in the northern hemisphere [Tschudi et al., 2010]. However, the precise angle depends on a number of factors and was estimated closer to -15 • in Hudson Bay [Togunov et al., 2021]. To allow for drift toward an unknown angle of bias relative to wind, we modelled the mean turning angle as a trade-off between bias parallel to the direction of the wind and bias perpendicular to wind. Specifically, mean turning angle µ (ϕ) t was assumed to follow a menotactic BRW circular-circular von Mises regression model [Rivest andDuchesne, 2016, Togunov et al., 2021]: where ψ t is the direction of wind at t relative to the movement bearing at t − 1. We hypothesized that the concentration of turning angle around µ (ϕ) t would increase as wind speeds increased. At higher wind speeds, the force is greater and the direction can be predicted more accurately. To capture this, we modeled turning angle concentration as a slope-intercept function of wind speed following: After fitting the BCRW drift model to the dropped collar data, we predicted the drift speed and direction given the estimated ECMWF wind vectors. Shape and scale parameters of step length were estimated using Eq. 2 and Eq. 3, which were used to estimate the latent drift speedl t as the mode of Weibull step length We used the estimated mode drift speed for the correction as it minimized the median drift-corrected step length l (c) t compared to the mean or median of the Weibull distribution.

Mean directionμ
(ϕ) t of drift was predicted using Eq. 4. Next, we estimated the drift-corrected paths by subtracting the estimated drift vectors from the observed drift vectors following: whereû t ∈ R andv t ∈ R are the estimated east-to-west and south-to-north components of the drift-corrected track. The estimated drift-corrected step lengthl (c) t was calculated following: Next, we estimated the drift-corrected trajectory of the track using the cumulative sum of the estimated u t and v t components following:x Otherwise.
are the estimated x and y coordinates of the drift-corrected track, respectively.
To test that we did not over-fit the data, we evaluated the drift correction BCRW performance in reducing the apparent speed of dropped collars by fitting it to a random subset of 10% of the dropped collar tracks and correcting for it in the entire dropped collar data set. When correcting for drift in the bear tracks, we fit the drift correction BCRW to the entire dropped collar data-set, and used the estimated coefficients and ECMWF wind to predict and correct for drift along the bear tracks. Step length of dropped collars before and after drift correction; (b) cumulative path of dropped collars before and after drift correction; (c) residual plot of step length for drift corrected dropped collars sorted by ordinal date. Grey represents original data, blue represents drift-corrected data based on BCRW fit to the full data set, and green represents drift-corrected data based on BCRW fit to a 10% random sample of dropped collar tracks.

B.2 Drift correction results
We expect that applying the drift correction would mostly reduce the movement of the drop collar close to zero. When fit to the full dropped collar data set, the drift correction reduced the mean step length by 29% from 1.76 km to 1.25 km ( Fig B1). When fit to a 10% random sample of dropped collar tracks, the drift correction reduced the mean step length by 28% to 1.26 km. (Fig. B1). The drift corrected tracks still appear to have a remaining directional bias, however it is less than the original track (Fig. B1). There appears to be some residuals with respect to day of the year ( Fig. B1C), however, including day of the year in the drift correction BCRW yielded only marginal model improvement with respect to reduction of step length (data not shown).

C Ice motion-based drift correction C.1 Drift correction methods
We tested using sea ice motion vectors produced by the National Snow and Ice Data Center [NSIDC;Tschudi et al., 2019Tschudi et al., , 2020 to correct for ice motion in the dropped collars as in Appendix B.
All the same procedures as described in Appendix B.1 were followed with the exception of the definition of step length l t ∈ (0, ∞) and turning angle ϕ t ∈ (−π, π] being defined in terms of sea ice drift rather than wind velocities. Specifically, the shape and scale parameters of step length were defined as slope-intercept functions of NSIDC ice motion speed r (N SIDC) t ≥ 0 as follows: where is β {1,2,3,4} ∈ R are the intercept and slope coefficients to be estimated.
Mean turning angle µ (ϕ) t was assumed to be a BRW circular-circular von Mises regression model [Rivest and Duchesne, 2016]: where ψ (N SIDC) t is the direction of NSIDC ice motion at t relative to the movement bearing at t − 1. In addition, we modeled turning angle concentration as a slope-intercept function of NSIDC ice motion speed following: After fitting the BCRW drift model to the dropped collar data, we predicted the drift speed and direction given the estimated NSIDC ice motion vectors. Shape and scale parameters of step length were estimated using Eq. 10 and Eq. 11 and were used to estimate the latent drift speedl t using eq. 6. Mean drift direction µ (ϕ) t of drift was predicted using Eq. 12. Finally, we estimated the drift-corrected paths, step lengths, and trajectory as described in Appendix B.1 using eq. 7, eq. 8, and eq. 9, respectively Step length of dropped collars before and after drift correction using NSIDC sea ice motion vectors; (b) cumulative path of dropped collars before and after drift correction; (c) residual plot of step length for drift corrected dropped collars sorted by ordinal date. Grey represents original data, blue represents drift-corrected data based on BCRW fit to the full data set, and green represents driftcorrected data based on BCRW fit to a 10% random sample of dropped collar tracks.

C.2 Drift correction results
When fit to the full dropped collar data set, the drift correction reduced the mean step length by 14% from 1.76 km to 1.51 km; Fig. C1). When fit to a 10% random sample of dropped collar tracks, the drift correction reduced the mean step length by 14% from 1.76 km to 1.51 km (Fig. C1).

D Drift correction HMM output comparison
This subsection compares broad differences in HMM output for four different methods of accounting for ice drift in telemetry data: 1) no drift correction/integration, 2) wind-based drift and tidal integration, 3) satellite-based drift vector subtraction and tidal integration, and 4) wind-based BCRW vector subtraction and tidal integration. The first method assumes there is no

D.1 No ice motion correction or integration
The first method assumes ice motion has no effect on bear movement.
Step length is modeled as interceptonly for each state (i.e., main text, equation 1). Turning angle is modeled following: where ψ (wind) t represents the directions of wind at time t relative to the track bearing at time t − 1, α 1 represents the magnitude of attraction toward ψ (wind) t for the olfactory search states, and α 2 represents the magnitude of attraction toward ψ (wind) t + π/2 for the olfactory search states [Togunov et al., 2021]. In this formulation, it is assumed that the direction of D and ARS is centered at 0 and that the direction of olfactory search states is exclusively detriment by the direction of wind.

D.2 Wind-based drift and tidal integration
In the second method, step length µ (l) S,t is modeled as follows: where β 1,S is the state-specific intercept coefficient for step length, β 2 is the slope coefficient for the effect of wind speed on drift speed, r Turning angle is modeled following: is assumed that the direction of drift only affects the drifting state, the ARS state turning angle is centered at 0, and that the direction of olfactory search states is exclusively detriment by the direction of wind.

D.3 Satellite-based drift vector subtraction and tidal integration
In the third method, NSIDC sea ice motion vectors (horizontal u and vertical v components) are first interpolated to the track location using tri-linear interpolation (first, bilinearly in space, then linearly in time). Next, the u (N SIDC) and v (N SIDC) components of drift are subtracted from the horizontal and vertical components of track displacement to obtain voluntary displacement fol owing: where u

D.4 Wind-based BCRW vector subtraction and tidal integration
The fourth method is the primary method described in the main text. First, wind-driven ice drift is estimated and subtracted from the polar bear tracks as described in Appendix A. Then remnant movement is modeled as described in main text equations 1-4.

D.5 Results
Based on the Viterbi-decoded states, the state frequencies were relatively comparable between the integration, satellite-correction, and BCRW-connection models: 17-19% in ARS, 46-51% in stationary drift, and 32-36% in olfactory search (Table D1). The HMM based on the raw data with no ice drift integration or correction exhibited a markedly lower frequency of ARS (12%) and stationary drift (22%), but nearly twice the frequency of olfactory search (66%). This is likely because the added motion of sea ice makes it much more difficult to identify stationary or low-velocity behaviors, while inflating the frequency of behaviors with directed travel.
When comparing specific decoded states, the concurrency of each model combination varied between 45-85% (Table D2). The largest discrepancies between models was between the raw HMM and the BCRW-based correction, followed by NSIDC-based correction. Among the models with drift correction and integration, concurrency ranged between 78-85% and the largest state discrepancy occurred in the classification of olfactory search. Specifically, when comparing the ice drift integration to the satellite-based correction (69% discovery rate) and BCRW correction (68% discovery rate; Table D2). ARS had the next lowest discovery rate between 80-87% among HMMs with drift correction/integration (Table D2). In every confusion matrix, both ARS and drift are most frequently classified as olfactory search, and olfactory search is classified as drift.
The difference in drift correction between the satellite-based method versus the BCRW-based method is particularly apparent in Fig. D1 (c) versus (d), wherein the segments classified as drift are more contracted in the BCRW-based method. The difference between the estimated step length and turning angle distributions are relatively minor between the different models ( Fig. D1 e-l).

D.6 Conclusions
The most confidently identifiable state is stationary drift, as it is characterized by smooth motion with constant speed and in some parts of the Bay is exhibits clockwise spirals due to tidal circulation. The most effective drift-correction methods should lead to drift segments becoming more contracted and increase the distinctness of spiraling due to tidal circulation. From the movement tracks (Fig. D1), it appears the BCRW-correction method used in this paper (described in Appendix B) is more effective than drift correction using satellite-based NSIDC ice motion estimates. This is also corroborated when comparing the results in Appendix B versus Appendix C, which show the BCRW-correction reduces the step length of dropped tags Step length (  by 29%, while the NSIDC-based correction reduced step length by 14%. Therefore, a BCRW correction is more effective than a NSIDC-based correction. As we do not have a ground-truth of behavioral state, it is not possible to know objectively from our data which of the presented methods is most effective at state identification. We did identify a notable difference between the state identification across the methods, particularly when comparing concurrency at the level of individual steps. A key challenge when modeling movement in a moving environment is that the observed step length and orientation both depend on the bear's voluntary step length and speed as illustrated in Figure   D2.  Michelot, 2018, Johnson et al., 2021]) may lead to bias both in the step length and turn angle, leading to compounding error and exaggerated misclassification. The lowest concurrency in state identification was of the olfactory search state, followed by ARS, then drift. This is in line with the prediction the integration methods are most effective for identifying stationary behaviors, followed by unbiased behaviors, and least effective for biased behaviors. Therefore, we believe the integration method's inability to estimate the voluntary orientation, leads to misclassification of ARS, and in particular, olfactory search.
As there is no model that can simultaneously estimate the effect of drift on movement and movement parameters, there is likely no perfect solution to correct or integrate motion of environmental data, assuming it cannot be accurately and precisely measured. However, the differences among the three correction/integration methods (i.e., BCRW-based, NSIDC-based, and integration) was much smaller than between any of them and the model assuming no drift (i.e., the HMM fit to the raw locations). Therefore, any method to integrate or correct for environmental motion will likely significantly improve the models accuracy and precision -failing to consider ice drift appears to have a large effect on state prediction and may lead to significant bias.

E Pseudo-design matrices
Pseudo-design matrix for step length in momentuHMM: The column correspond to the coefficients estimated using maximum likelihood by the model, and the rows correspond to the movement parameters: step length mean (µ Pseudo-design matrix used for turning angle: angleFormula is a special function that can be included in design matrices to model circular parameters (e.g., mean turning angle µ