Disentangling the biotic and abiotic drivers of emergent migratory behavior using individual-based models

Understanding the drivers of movement, migration and distribution of individuals is important for insight into how species will respond to changing environmental conditions. Both abiotic and biotic factors are thought to influence migratory behavior, but their relative roles are difficult to disentangle. For migratory marine predators, both temperature and prey availability have been shown to be significant predictors of space use, though often researchers rely on physical proxies due to the lack of data on dynamic prey fields. We generated spatially explicit individual-based movement models to evaluate the relative roles of abiotic (sea surface temperature; SST) and biotic (prey availability) factors in driving blue whale (Balaenoptera musculus) movement decisions and migratory behavior in the eastern North Pacific. Using output from a lower trophic ecosystem model coupled with a regional ocean circulation model, we parameterized a blue whale movement model that explicitly incorporates prey fields in addition to physical proxies. A model using both SST and prey data reproduced blue whale foraging behavior including realistic timing of latitudinal migrations. SSTand prey-only population models demonstrated important independent effects of each variable. In particular, the SST-only model revealed that warm temperatures limited krill foraging opportunities but failed to drive seasonal foraging patterns, whereas the prey-only model revealed more realistic seasonal and interannual differences in foraging behavior. Our individual-based movement model helps elucidate the mechanisms underlying migration and demonstrates how fine-scale individual decision-making can lead to emergent migratory behavior at the population level. Moreover, determining the relative effects of the physical environment and prey availability on the movement decisions of threatened species is critical to understand how they may respond to changing ocean conditions.


Introduction
Animal migration is at once a critically important and a globally threatened ecological process (Bauer and Hoye, 2014;Wilcove and Wikelski, 2008). Migration increases survival by enabling animals to exploit seasonally available resources, transports nutrients, propagules and pathogens between ecosystems, and facilitates trophic interactions by moving large numbers of predators or prey over wide ranges (Bauer and Hoye, 2014). Migratory species are also among those most at risk from human-induced rapid environmental change (Horns and Şekercioğlu, 2018;Wilcove and Wikelski, 2008), necessitating an understanding of the ecological and environmental underpinnings of this behavioral phenomenon. Though migration for many species is often understood to be a means for reaching distinct breeding and foraging ranges to enhance survival and reproduction (Dingle and Drake, 2007), a range of drivers can influence finer-scale foraging patterns within a larger-scale seasonal migratory behavior. On a daily scale, migratory animals have to decide whether to remain in a foraging hotspot or move to a new foraging ground based on present and past conditions (Abrahms et al., 2019a;Palacios et al., 2019). Mechanisms influencing migratory behavior and space use range from endogenous factors such as memory (Abrahms et al., 2019a;Bracis and Mueller, 2017;Scott et al., 2014) and navigational capacity (Lohmann et al., 2008), to exogenous factors such as social information (Guttal and Couzin, 2010;Jesmer et al., 2018;Mueller et al., 2013) or predation pressure (Fryxell and Sinclair, 1988).
Bottom-up processes including responses to environmental cues and tracking of resource availability have also been shown to profoundly affect an animal's movement decisions during migration (Merkle et al., 2016;Thorup et al., 2017). These environmental drivers can be either biotic, such as forage availability, or abiotic, such as temperature, snowfall, or water availability, and, importantly, can differ in their importance by the spatial scale of inquiry. For instance, forage quality (biotic) and water availability (abiotic) were shown to drive broad-scale spatial distribution patterns in wildebeest migrations, though biotic factors alone drove distributions at finer spatial scales (Holdo et al., 2009). Similarly, primary production and water temperature have been shown to be significant predictors of broad-scale distribution patterns for many migratory marine predators (Block et al., 2011;Hazen et al., 2013), but prey availability is hypothesized to dominate space use over abiotic factors at local scales (Benoit-Bird et al., 2013;Hazen et al., 2009;Kenney et al., 2001).
Given the challenges in monitoring migratory species and their environments over large spatial extents and timeframes, the relative roles of biotic and abiotic factors in driving the fine-scale movements and timing of migratory megafauna are difficult to disentangle in many systems. This is especially true for marine migrants, for which habitats can shift at daily to weekly timescales (Steele and Henderson, 1994) and physical proxies are often used in absence of dynamic prey fields (Abrahms et al., 2018;Croll et al., 2005;Owen et al., 2018;Scales et al., 2014). Nevertheless, identifying the dominant bottom-up processes and thresholds for decisions driving foraging migration is key to adequately capture seasonal movement patterns, interannual variability in migration timing, interaction with anthropogenic threats, and anticipating the responses of migratory species to environmental change.
We developed a spatially-explicit, individual-based model (IBM) to tease apart the biotic and abiotic drivers of fine-scale movements leading to emergent migratory behavior for a marine predator, the blue whale (Balaenoptera musculus). Blue whales are the largest animal to have lived on the planet, and their large energetic demands and requirement for dense patches of krill have spurred the evolution of their physiology and foraging behavior (Goldbogen et al., 2011;Goldbogen et al., 2015;Slater et al., 2017). In order to maximize their energetic gain and reproductive success, blue whales like many baleen whales migrate from tropical breeding grounds to temperate or arctic feeding grounds annually, timing their migration to reliable peaks in their prey, krill (Abrahms et al., 2019a;Fossette et al., 2017). In the eastern North Pacific, blue whales perform latitudinal migrations between winter breeding grounds off the coast of Central America and in the Gulf of California and summer feeding grounds that extend from Baja California to Washington state in the California Current System (Bailey et al., 2009;Irvine et al., 2014).
Previous models aimed at understanding blue whale movements have focused on species-habitat relationships, using environmental proxies for prey to explore where and when the whales occupy the California Current. These have included static models based on shipboard sightings (Becker et al., 2016;Redfern et al., 2019), photo-ID work , global models based on sea surface height alone (Pardo et al., 2015), foraging models  and dynamic models based on changing oceanography (Abrahms et al., 2019b;Hazen et al., 2017). Many of these models have found that sea surface temperature, mixed-layer depth, and chlorophyll-a were most correlated with blue whale habitat, often tracking the seasonality of upwelling dynamics. Chlorophyll-a represents a direct proxy for productivity and krill aggregations (Abrahms et al., 2019a;Bailey et al., 2012;Suryan et al., 2012), whereas temperature may indicate aggregative features that can concentrate prey Scales et al., 2014). Mixed-layer depth is often an ocean-model derived variable that is used instead of chlorophyll-a with similar predictive capacity (Becker et al., 2016). In addition, stochastic dynamic programming techniques have been used to develop models that integrate habitat patch quality and energy demands of lactating females to identify optimal movement decisions (Pirotta et al., 2019;Pirotta et al., 2018). However, these studies lack explicit information on the distribution of prey aggregations, which are known to significantly shape blue whale space use and foraging decisions .
Individual based models (IBMs) provide an alternative framework that treat individuals as autonomous agents. Each individual is prescribed a set of deterministic or stochastic rules that inform and update the behavior, movement, or interaction of agents. Ensembles are created by the simulation of many agents, allowing for inference of population statistics, geographic locations, and emergent behaviors, while still retaining the discreteness and variability of individuals and small populations. A strength of IBMs is the ability to test individual decision processes and mechanisms underlying emergent behaviors and self-organization. Examples include flocking of starlings (Hildenbrandt et al., 2010), stripe formation in zebra fish (Volkening and Sandstede, 2015), and aggregation of locusts into destructive hopper bands (Bernoff et al., 2020). From a population modeling standpoint, IBMs have successfully been used to explore the distributions and behaviors of marine and terrestrial animals, including sea lions (Fiechter et al., 2016), juvenile chinook salmon , sardine and anchovy Rose et al., 2015), Icelandic capelin (Barbaro et al., 2009), and caribou (Latombe et al., 2014). In the context of blue whale populations, emergent processes are a result of individual decisions and include both behavioral and spatiotemporal aspects such foraging rates and the timing of migrations.
The objective of this paper is two-fold: first we introduce an IBM framework that reproduces observed behaviors of blue whales in the California Current System, and then explore model ensembles to examine contributions from abiotic (temperature) and biotic (krill availability) factors on migratory behavior. Including explicit realistic movement of individuals adds information about the temporal accessibility of habitats by considering transit speeds and foraging opportunities; details ignored by species distribution models, other correlative techniques, and the stochastic dynamic programming models. Specifically, we compare the independent and combined effects of sea surface temperature and krill availability on (1) foraging behaviors and latitudinal distributions of the blue whales during their northward migration, and (2) transition to southward migrations. This comparative approach allows us to evaluate the relative importance of biotic and abiotic fields on population-level behaviors. Moreover, this is one of the first studies to use direct prey fields in place of environmental proxies on a broad scale to examine the drivers of blue whale migration.

Environmental and prey data
Environmental and prey data are simulated using an implementation of the Regional Ocean Modeling System (ROMS) coupled with the biogeochemical NEMURO model (Fiechter et al., 2016;Fiechter et al., 2015;Kishi et al., 2007;Rose et al., 2015) for years 2000-2010. NE-MURO describes nutrient, phytoplankton, and zooplankton dynamics in each grid cell, which are advected and diffused in space with ROMS. The ROMS-NEMURO output is generated prior to the IBM calculations and provides environmental conditions including daily sea surface temperature (SST) and relative krill density ρ at 3 km horizontal resolution over a domain ranging from 116 to 128 ∘ W and 32-44 ∘ N (Fiechter et al., 2018). Krill relative density is a near-surface unitless quantity and is a real number between 0 (low) and 1 (high). The ROMS-NEMURO domain serves our primary focus and study area as blue whales are known to exhibit extensive foraging in this region from May-October and the model provides physical and prey fields at sufficiently fine spatial and temporal resolution to inform behavior in the IBM. Fig. 1 shows quarterly averages for SST and krill density over the full computational domain for the year 2008; a year with ENSO neutral conditions in the California Current region and relatively high krill abundance in the coastal upwelling zone. In April-September, the upwelling along the coast is evidenced by cooler nearshore waters and locally enhanced krill concentrations.

Model description
Spatially-explicit individual-based models (IBMs) used in population dynamics seek to capture the motion of individuals in a heterogeneous landscape through a series of discrete updates. Here, the IBM is formulated as a semi-Markov state switching model with N-states. Each state N is correlated with a behavior and characteristic state movements are described by step-length and turning angle distributions. Locations and states for all individuals are updated at each time step. States are selected probabilistically with elements p j,k of an N × N state transition matrix τ define the probability of transitioning from state j to k. Ensembles are created by the simulation of many agents, allowing for inference of population statistics and behaviors, while still retaining the discreteness of individuals.
IBMs allow us to program behavior and see which aspects most closely reproduce observed migration patterns and the probabilistic updates naturally incorporate variations in each decision. To explore how environmental and prey influences impact whale movements throughout one seasonal migration period (May-December), we define and analyze ensembles from the two IBMs described below.

Transit-Forage model
We used published movement parameters in our model to test the hypotheses of what drives blue whale foraging migrations. Previous studies of blue whale tagging data (Bailey et al., 2009) have identified two distinct behavioral states, corresponding to transit and foraging. In transit, whales tend to move in a straight line, whereas slower speeds and higher turning angles are observed in foraging. As such, we employed a two-state model, with states 1 (transit) and 2 (forage) used to evaluate how SST and prey availability drive whale movements, behaviors, and distributions.

North-South model
Previous correlative models have shown higher success rates of seasonal versus annual models when predicting blue whale migration patterns (Abrahms et al., 2019b;Hazen et al., 2017). Blue whales likely rely on different biological and physical processes during their northward migration to feeding grounds versus southward migrations to breeding grounds, with the transition likely triggered by bioenergetic or breeding demands. Additionally, a higher percentage of transit is reported during the southward migration with individuals only stopping to forage if the prey conditions are sufficiently dense (Bailey et al., 2009). Northward and southward states are therefore split for the migration study to account for the stricter foraging requirements while returning to southern breeding grounds. Mechanisms for transition to southward migrations are therefore studied with a four-state model; states 1 and 2 describing transit and forage during northward migration and states 3 and 4 for southward migration.
We use a migration strategy based on the individual foraging rate averaged over the previous 10 days, which is implemented in the transition probabilities from 1 and 2 (north) to 3 (south). Low foraging rates increase the probability of transitioning to southward migration, and the time-average provided sufficient opportunity for agents to travel between and seek out new krill patches. Again, we test the effect of SST and prey density on this transition. Additional migration strategies including a krill satiation threshold in which southward transitions are initiated when whales consume sufficient levels of krill, and natal breeding homing were tested for transitioning to southward migration and are described in the Supplement.

State transition functions and movement distributions
Behavioral states are selected probabilistically and are a function of current environmental covariates X (SST), krill density ρ, and model parameters θ. The probabilities encode how environmental and prey conditions influence the state transitions. The 2 × 2 state transition matrix for the transit-forage model is shown in Eq. (1) as an example, with the 4 × 4 matrix for the north-south model given in the Supplement Here, is defined as the probability of foraging and s t is the behavioral state at time t. The probability of foraging is a weighted sum of environmental and krill transition functions and is given by with the weights w 1 and w 2 controlling the relative importance of SST and krill and = + z w w 1 2 to normalize probabilities to 1. The environmental transition probability is the simple linear response function that selects for colder temperatures. Krill transition probabilities are set by the logistic selection function where probability of foraging increases with krill density. Model parameters are provided in Table 1. We assume future states depend only on current conditions and are independent of the past and current state, so that τ is constant on the columns and rows sum to 1. In the model, behavioral states are updated every 6 h. Fig. 3 shows the probability of foraging computed for the full domain using monthly SST and krill inputs; figures of daily probabilities are provided in the supplement.
Across both the transit-forage and north-south models, the forage and transit states have the same movement distributions, regardless of north-or southward migration.
Step lengths are positive real numbers in units of meters and follow gamma distributions with probability density function (Bailey et al., 2009;Scales et al., 2016) where Γ(a) is the gamma function. Turning angles are sampled from a von Mises distribution with probability density function given by where I 0 (κ) is the Bessel function of order 0. A 0 ∘ turning angle corresponds to straight movement in all states except southward transit , 3 where 0 ∘ is set to south. Parameters used in step length and turning angle distributions for transit and forage states are given in Table 1 and movement distributions are shown in Fig. 2. The gamma distribution shape and scale parameters yield average step lengths of 22.2 ± 12.9 km (transit) and 6.3 ± 5.8 km (forage) for a 6 h time period, consistent with (Bailey et al., 2009;Scales et al., 2016).

Model implementation
For its implementation in the IBM, the semi-Markov model formalism is combined with an area restricted search (ARS). The ARS adds intelligence by allowing individuals to locally sense areas of high foraging potential. The algorithm is as follows, and is further outlined in Fig. 2. Let Y w, t be the geographic location of whale w on time step t. Individual locations are initialized on a given day uniformly at random within a pre-defined region of the domain. The probability of foraging P w, t is calculated for each location Y w, t and used to select the next behavioral state + s w t , 1 . An area-restricted search is conducted on the perimeter of a box of radius R to find the location Y w with highest probability of foraging. If such a location exists, the turning angle is centered around the direction from Y w, t to Y w .
Step lengths and turning angles are sampled from the state distributions, and used to update the locations to + Y w t , 1 . Positions and states are updated at 6 h intervals using daily ROMS output.
Ensembles of 2000 whales for years 2000 -2010 were generated using the transit-forage and north-south models. To understand the independent roles of SST and krill in driving the observed movement patterns of the northward migration in the transit forage-model, we additionally generated model ensembles using only SST or krill as inputs. These one-input ensembles can be thought of from a modeling perspective as probing the independent effects of each variable, or from a biological perspective as populations with distinct strategies that base their fine-scale decisions on only one of the two variables. We analyzed the effects of these strategies on broad-scale migration patterns, seasonal and interannual differences in prey consumption, and transitions to southward migration. The same environmental and krill transition functions were used in all cases, and turning a driver off or on is equivalent to adjusting the weights in the state transition matrix. Results are additionally measured against a correlated random walk (CRW) with no environmental or prey preferences. The CRW has the same behavioral states and movement distributions, with the probability of foraging fixed at 0.4, in alignment with published foraging rates for May -November (Bailey et al., 2009). Krill densities and SST are recorded in the CRW model at whale locations for comparison purposes.
To measure and compare the foraging success between population ensembles, we introduce a foraging potential metric Ω. The value of Ω for each individual w provides a sum of the krill density encountered while in a foraging state and is given by Here, ρ(Y w, t ) is the krill density at the location of whale w, is an indicator function that takes the value one if whale w is foraging (in state 2) and zero otherwise, and the summation is over all time steps t.

Modeling assumptions
We summarize the model assumptions that were made in the development of these individual based models to study how daily to monthly SST and prey variability contribute to seasonal whale movements and foraging migrations.
• Model whales based their fine-scale decision only on SST and/or krill availability. While other environmental factors may contribute to the process, SST is often highly correlated with whale presence (Gill et al., 2011;Hazen et al., 2017;Palacios et al., 2019) and is one of the most accessible variables for use in models. The direct importance of prey in constraining distributions and migrations of marine megafauna has not yet been evaluated.
• Future behavioral states only depend on current prey and environmental conditions and are independent of current and past states. Model runs are initiated at the start of the migration season and whales are not allowed to learn from prior migrations.  • Upon entering the domain, all agents are in transit state 1 with the turning angle 0 ∘ aligned with due north.
• In the north-south model, the transition probability from 1 and 2 to 4 is 0. That is, whales only transition from north transit-forage states to the south transit state 3 . Transition to state 3 was restricted until day of year 150 to allow the agents to explore the domain and find foraging grounds.

Parameter selection and sensitivity analysis
When possible, model parameters were determined from previously published resources to reduce the number of free variables. The environmental transition probability function is parameterized to select for colder temperatures associated with upwelling. The step length and turning angle distributions for the transit and forage states are derived from the speeds and turning angles of the transit and area restricted search behavioral states in (Bailey et al., 2009) and scaled to a 6 h time step. The 6-hour time steps were selected to best correspond with the 3 km spatial resolution of the ROMS data. The average step lengths associated with a 6-hour interval allow the whales to transit and forage within the local domain without stepping over and missing potential foraging opportunities. Area restricted search is designed as a local search algorithm, thus the search radius was set to match this average step length. Thus, the free parameters and model variables are reduced to β and ρ * used in the krill transition probability, the day of year individuals enter the domain (start date), and starting location.
A sensitivity analysis indicates that ensemble spatiotemporal distributions and population level transit-forage behaviors are robust to reasonable changes in the free parameters. Results of the sensitivity analysis and additional justification for parameter selection are provided in the Supplement. The start date and starting locations of individuals were identified as the most influential parameters on the population's northward migration and foraging success, whereas changes to the parameters β and ρ * have minimal impact. The model sensitivity to start date aligns with recent studies of blue whale call patterns which reveal variability in the timing and success of the northward migration (Szesciorka et al., 2020). To account for the increased sensitivity these variables add, the start date is uniformly distributed from May 1, -June 1, and the starting location is selected uniformly at random from a wide box in the southern part of the domain.

Results
We aggregated population ensembles for all years to analyze broadscale patterns. Further, we considered specific years to highlight transient and interannual trends.

Northward migration
Using only fine-scale decisions, the transit-forage model produced realistic broad scale northward migration patterns, timings, and distributions. Fig. 4 shows the aggregated results for all years 2000-2010. Specifically, the ensembles reach a northernmost mean latitude by early October (Fig. 4a) and peak abundance of blue whales in the Monterey Bay region (defined by latitudes 36.6-370N) between July and October (Fig. 4b). These distributions align with previous reports on migration characteristics, estimates based on satellite tagging data ( Fig. 4a; Abrahms et al., 2019a;Irvine et al., 2014), and surveys of shipboard sightings (Croll et al., 2005;Fossette et al., 2017). In addition, the ensembles show a higher level of foraging during May -October, consistent with empirical observations (Bailey et al., 2009;Irvine et al., 2014), with a sharp decrease in foraging rates (Fig. 4c) as krill abundance begins to decrease October in the domain, in concordance with the end of upwelling.
The first panel in Fig. 3 displays July 2008 aggregates for the probability of foraging and foraging locations from the transit-forage model, and emphasizes that these two variables are not directly correlated. In the IBM foraging locations depend on the probabilities, but also explicitly on the current location of each agent and the travel times to productive habitats. In July, the majority of the population remains in southern regions of high foraging potential and has yet to reach the rich northern foraging grounds, despite the high probability of foraging above latitude 40 ∘ N.

Contributions of environment and prey on the northward migration
To understand the independent contributions of SST and krill to the observed northward migrations we generate and analyze model ensembles that use only SST or krill as inputs. Strategies based only on krill or SST lead to populations with distinct domain use patterns (Figs. 3-5). First, the July 2008 foraging location aggregates (Fig. 3) demonstrate that the influences of SST or krill drive the populations to feed in spatially distinct regions. In our model, the temperature cues lead to dispersed habitats with small foraging hotspots only localized at the coast. The influence of krill promoted increased offshore habitats and foraging hotspots.
Furthermore, when only SST is used as a driver, the proportion of the whale population in the foraging state remains relatively constant between 30-40% throughout the year, with the remaining individuals in the transiting state (Fig. 4c). Longer step-lengths from the transit distribution allow the SST-only population to explore the domain quickly, resulting in a fast migration reaching the northernmost mean latitudes earlier in the summer and overall a higher northernmost mean latitude. In contrast, the krill-only strategy mimicked the full model in the latitudinal distributions (Fig. 4a) and arrival of whales in Monterey Bay (Fig. 4d), suggesting that the influence of krill is responsible for the timing of the northward migration in the combined model.
Daily and seasonal foraging patterns are also remarkably more distinct in prey-driven populations. Both models with krill have at least a 3 times higher variance in daily foraging rates when compared to the SST model. The full and krill-only models additionally show a transition to lower foraging rates between day of year 275 -300 (beginning-late October) in concordance with the end of the upwelling season ).
Higher foraging rates from krill inputs naturally correspond to a higher foraging potential. Normalized distributions of Ω for each of the three strategies are shown in Fig. 4d. The krill-only ensemble shows higher average foraging potential than the SST-only and combined populations. Thus, the influence of krill drives the foraging success of the combined population. However, the SST-only model still had a higher mean foraging than the CRW model, indicating that the influence of favorable SST drives individuals to forage in locations with higher krill densities than they would find by random chance.

Interannual differences in prey consumption and foraging
The years 2003 and 2008 have lower and higher than average modeled krill densities, respectively, and higher (2003) and lower (2008) than average SST (Fig. 5a, Supplementary Table 2). Results from the full model show significant interannual differences in agreement with the increased levels of krill from the ROMS data. 2008 has a 4% average higher foraging rate averaged over all simulation days than 2003 (Fig. 5c), with up to 20% differences depending on the time of year. Prior to day of year 260, the 2008 population shows an 8.8% higher average foraging than 2003 with daily differences reaching as high as 20%. After this date, 2003 temporarily has a higher foraging rate before both populations reach equal foraging.
A comparison of the population ensembles for 2003 and 2008 provided understanding of how SST and prey factors contributed to the observed interannual differences. Both SST and krill populations have a 3-4% higher average foraging rate for 2008 than in 2003. Thus, the lower average SST in 2008 alters the population behavior and increases the overall foraging, but the addition of krill again provides more detailed temporal information. The SST population shows no significant differences in the foraging rates before and after day of year 260, whereas in the spring and summer months the 2008 krill population has a 10.7% higher rate than 2003. Furthermore, the foraging rates for 2008 are lower than 2003 just after day of year 260 despite the higher krill density in 2008.
Average foraging potential Ω is higher in 2008 across all strategies, with krill driven populations showing a secondary peak at higher levels. Independent of the yearly variances in prey density, the krill driven populations exhibited the sharp transition to lower foraging rates in October. Thus, our results suggest the influence of prey provides the fine scale daily, seasonal, and interannual details observed in the full model.

Implications on southward migration
The addition of the southward migration states results in a southward shift in the latitudinal distributions for all strategies (Fig. 6a). All cases show similar timing with the latitudes peaking around day 300, however, the addition of krill influences transitions to southward migrations in realistic ways. Specifically, all models show a large early migration at days 150-200 as whales unsuccessful at finding and remaining within rich foraging grounds migrate south. However, models that include krill continue to transition throughout October and November with a secondary peak around day 300 corresponding to the seasonal reductions in population foraging rates that are observed in the transit-forage model (Fig. 6b). Using only SST, the proportion of whales foraging remains relatively constant and low through the spring and fall (Figs. 4b -5 c), which leads to time independent and early transitions. Furthermore, the cumulative percentage of whales transitioning into a southward migration plateaus around day 300 in the SSTonly model (Figure 6 c), whereas the krill and SST + krill models display continuing transitions to southward migration later in the season. With all strategies, not all whales transition to southward migrating states as individuals may leave the study area prior to the start of migration or remain in areas with high foraging potential.

Discussion
Understanding the bottom-up drivers of spatiotemporal distributions of highly migratory species is important for predicting changes in future habitat use and for implementing effective dynamic conservation policies. Here, we introduced an IBM which reproduced realistic foraging and migratory behaviors at both individual and population levels, suggesting that fine-scale individual decisions lead to emergent properties of the population (Mueller and Fagan, 2008). Furthermore, we used model populations driven by only SST or krill to help uncover the direct impact of each variable on whale behavior and migrations.
Our findings indicate that abiotic (SST) and biotic (krill abundance) factors both play key roles in individual movement decisions during round-trip blue whale migrations.
Abiotic factors such as SST are commonly used as proxies for prey in blue whale models (Abrahms et al., 2019a;Becker et al., 2019;Pirotta et al., 2019;Pirotta et al., 2018). While we found that SST indeed leads whales to forage in areas with higher prey availability at broad spatial scales compared to a simple CRW, important spatial and temporal patterns are missed by using SST without including krill contributions explicitly . Notably, the proportion of whales foraging remained nearly constant throughout May -December, thus missing important seasonal fluctuations in foraging behavior (Bailey et al., 2009). Additionally, driving populations only with SST led to more northerly average latitudes (Fig. 4a), inaccurate arrival times to Monterey Bay (Fig. 4d), and could not realistically reproduce the seasonal southward return migration (Fig. 6).
In contrast, populations influenced only by krill retained realistic seasonal behaviors, such as increased foraging during the upwelling season and denser krill years, as well as an appropriate timing of southward migrations (Abrahms et al., 2019a;Bailey et al., 2009). Furthermore, the influence of krill resulted in increased offshore habitats and foraging hotspots. However, these ensembles expressed unrealistically high levels of seasonal foraging  when compared to published estimates of monthly foraging behavior (Bailey et al., 2009).
When combined, krill acted to inform biologically feasible seasonal and yearly variations in foraging behavior, and SST reduced excessively high foraging rates, plausibly by accounting for the impact of unfavorable environmental conditions on movement and foraging decisions. While the results presented here represent statistics from large population ensembles run with a fixed set of parameters and using simulated SST and prey fields, we found general population distributions and migratory behaviors were robust to reasonable changes in model parameters (Supplement).
The growing use of IBMs has opened new avenues to test mechanistic questions of animal movement and make predictions for future scenarios (Bauer and Klaassen, 2013;DeAngelis and Mooij, 2005;Stillman et al., 2014;Tang and Bennett, 2010). Examples range from terrestrial to marine communities and include the role of memory in zebra (Bracis and Mueller, 2017) and elk (Bennett and Tang, 2006) migrations, environmental and bioenergetics factors influencing decision processes during bird migrations (Duriez et al., 2009), consequences of conservation and harvesting practices on fisheries (Ayllón et al., 2018;Jørgensen et al., 2008;O'Callaghan and Gordon, 2008), and how generic habitat characteristics can lead to migratory behavior (Shaw and Couzin, 2013). Our results here contribute to the wealth of spatiotemporal and mechanisitic information that can be learned from IBM ensembles.
As Fig. 3 shows, a high probability of foraging does not directly correspond to foraging whales. Correlative species distribution models provide information on available habitat yet have little consideration of behavioral factors that could limit use. In other words, a species' realized niche is not necessarily its fundamental niche. These models can result in high predictions in areas that are unavailable to migratory predators and falsely link animal presence with preferred habitats ( Fig. 3; Garshelis, 2000;Railsback et al., 2003). For example, Hazen et al. (2017), Abrahms et al. (2019b) predict high habitat likelihood in the northern part of the domain, but these ranges are less likely to be accessed by migrating whales due to movement limitations. The dynamic state model approaches in (Pirotta et al., 2019;Pirotta et al., 2018) seek to answer mechanistic factors driving the migration strategies for lactating female blue whales. Locations and energetic needs of the individuals are included, but movement is limited to stepping along a one-dimensional line of 100km boxes adjacent to the coast. The presented IBM overcomes these limitations by explicitly linking fully two-dimensional discrete movement processes and behavioral states to the available environmental and prey conditions, which revealed broad-scale migration patterns and mechanistic information about habitat use and foraging behaviors. Our results reveal the relative roles of the physical environment and prey availability in the movement decisions of a migratory marine predator and shed light on environmental processes influencing foraging and migratory behavior on daily to seasonal scales. Such information can help relate future changes in biotic and abiotic environmental conditions to a species' anticipated behavioral responses.
Migratory species continue to face an expanding number of threats from changing environmental conditions and interactions with humans (Horns and Şekercioğlu, 2018;Maxwell et al., 2013;Wilcove and Wikelski, 2008). Our findings also have important ecological and conservation implications, as ocean temperatures are warming globally and climate change is expected to redistribute species richness in lower trophic levels (Woodworth-Jefcoats et al., 2016). Previous work has predicted how ocean conditions may influence the distributions of top predators over the next century . However, changes in habitat quality may have disproportionate effects depending on the distance from normal migration routes. For example, increased habitat quality far from shore may not be discovered or may not be energetically feasible to exploit. Here, we show how year-to-year differences in conditions are clearly captured in the model ensembles (Fig. 5).
By using climate change projections as input, the model could generate more realistic predictions for the responses of blue whales to future scenarios. Increased water stratification from ocean warming is likely to affect the vertical distribution of krill, and krill populations may shift poleward or decline in response to predicted changes in coastal upwelling systems (Di Lorenzo et al., 2005;Palacios et al., 2004;Roemmich and McGowan, 1995;Rykaczewski et al., 2015;Sydeman et al., 2014). Given we found SST and krill were both important in determining multiple aspects of blue whale space use and foraging behavior, concurrent climate change impacts to both water temperatures and krill may influence blue whale distributions and prey consumption in synergistic ways. Our work highlights the value of including multiple forage-based drivers of animal movements and distributions to not only elucidate the mechanisms underlying migratory behavior, but also to help quantify the effects of anthropogenic risks on a threatened and highly migratory species.
Finally, the models presented here provide insight into the dynamic behavior of a highly migratory marine species, thereby laying the foundation to explore other aspects of the ecology and distributions of marine megafauna. Additionally, an individual-based model framework allows for testing the roles of other variables or processes on emergent migratory behavior. For example, animal culture and social interactions have been shown to influence the migratory capacity of other species, including ungulates (Jesmer et al., 2018) and birds (Mueller et al., 2013). Blue whales have been shown to communicate over large distances (roughly 500 -1,000 km) (Watkins et al., 2000), but the impacts of such communication on migrations have not been explored. The role of sociality in whale migrations could readily be tested in our modeling framework by allowing agents in sufficiently dense krill patches to signal their location and attract others. Optimal group sizes could also be evaluated by modeling trade-offs associated with density-dependent foraging opportunities (Markham et al., 2015). Finally, memory effects are believed to play a role in migratory behaviors (Abrahms et al., 2019a), and their role could be tested through the use of climatological inputs. Thus, the flexibility of a spatially-explicit individual-based movement model such as that presented here provides the opportunity to test a wide range of hypotheses in animal behavioral ecology.

Authors' Contributions
SJB and ELH and SD conceived the original idea and all authors contributed to study design. SD developed the model and performed analyses. SD and BA and ELH wrote the manuscript. SJB and JF contributed to editing the text.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.