Memories of Migrations Past: Sociality and Cognition in Dynamic, Seasonal Environments

Seasonal migrations are a widespread and broadly successful strategy for animals to exploit periodic and localized resources over large spatial scales. It remains an open and largely case-specific question whether long-distance migrations are resilient to environmental disruptions. High levels of mobility suggest an ability to shift ranges that can confer resilience. On the other hand, a conservative, hard-wired commitment to a risky behavior can be costly if conditions change. Mechanisms that contribute to migration include identification and responsiveness to resources, sociality, and cognitive processes such as spatial memory and learning. Our goal was to explore the extent to which these factors interact not only to maintain a migratory behavior but also to provide resilience against environmental changes. We develop a diffusion-advection model of animal movement in which an endogenous migratory behavior is modified by recent experiences via a memory process, and animals have a social swarming-like behavior over a range of spatial scales. We found that this relatively simple framework was able to adapt to a stable, seasonal resource dynamic under a broad range of parameter values. Furthermore, the model was able to acquire an adaptive migration behavior with time. However, the resilience of the process depended on all the parameters under consideration, with many complex trade-offs. For example, the spatial scale of sociality needed to be large enough to capture changes in the resource, but not so large that the acquired collective information was overly diluted. A long-term reference memory was important for hedging against a highly stochastic process, but a higher weighting of more recent memory was needed for adapting to directional changes in resource phenology. Our model provides a general and versatile framework for exploring the interaction of memory, movement, social and resource dynamics, even as environmental conditions globally are undergoing rapid change.


INTRODUCTION
Seasonal migrations are widespread among terrestrial, aquatic, avian and invertebrate species (Dingle, 2014). For many species, migration is an extremely successful strategy, allowing a far greater number of individuals to inhabit landscapes which might not otherwise be able to support large numbers year round (Fryxell et al., 1988). The evolutionary stability of a migratory strategy essentially relies on the fitness benefits of accessing seasonal resources, whether for energetic gain, predator avoidance, or a suitable environment for reproduction, outweighing the energetic and survival related costs of migration (Avgar et al., 2014).
Proximate causes, drivers and mechanisms for migration vary widely across and even within species (Berthold, 1999;Shaw, 2016). Some migrants follow a "green wave" of spring vegetation as it flowers across altitudinal or latitudinal gradients (Bischof et al., 2012;Kölzsch et al., 2015;Merkle et al., 2016). These migrations can be considered "tactical, " as they can occur-as an extreme simplification-purely as response to local conditions. Other migrants perform long-distance migrations in anticipation that critical resources will be available at the time of arrival at the end point of migration (Abrahms et al., 2019). This second behavior involves the greatest trade-off between the costs and benefits of accessing those highly seasonal and localized resources. This approach can be considered "strategic" in the sense that it is driven not by immediate cues but by an anticipation based on prior experience (Bracis and Mueller, 2017;Merkle et al., 2019;Bauer et al., 2020).
Migration can be a very successful strategy, with migratory ecotypes of the same species often outnumbering non-migratory conspecifics. Migratory caribou and reindeer Rangifer tarandus, for example, are several orders of magnitude more abundant than non-migratory woodland, mountain and forest ecotypes (Festa-Bianchet et al., 2011;Uboni et al., 2016). However, the question of whether migratory animals are more or less resilient to environmental disruptions in the environment remains open and largely case-specific (Moore and Huntington, 2008;Hardesty-Moore et al., 2018;Xu et al., 2021). On the one hand, migratory species may be more vulnerable as disruptions in either of the seasonal ranges or along a migratory corridor can have significant negative impacts (Wilcove and Wikelski, 2008;Seebacher and Post, 2015;Kauffman et al., 2021). On the other hand, migratory species might be more resilient due to their general wide-ranging mobility (Robinson et al., 2009). The resilience of a migratory population depends on the plasticity and adaptability of the population, which can take multiple forms, reflecting variation in where, when and whether the migration occurs Xu et al., 2021).
Cognitive processes, in particular spatial memory, have been shown to be important mechanisms for the reinforcement and maintenance of migration (Merkle et al., 2019;Bauer et al., 2020). Similarly, sociality and social learning are likely essential to maintaining migration (Guttal and Couzin, 2010;Fagan et al., 2011;Berdahl et al., 2018;Jesmer et al., 2018). However, the interacting role of sociality and spatial memory for the plasticity of migration and the resilience of the behavior when faced with a changing environment are generally unknown, though it has been hypothesized that the importance of these cognitive processes depend on the predictability of these resources (Riotte-Lambert and Matthiopoulos, 2020). Because the scenarios underlying migration are manifold and complex, mathematical modeling may provide some insights and help clarify where, when, and under what conditions we might expect migration behavior to emerge, to be adaptive, to be maladaptive, or to collapse.
Here, we develop a diffusion-advection model with sociality and memory to explore the resilience of a migratory population under various dynamic, seasonal resource distributions. In formulating the model, our goal was to identify the minimum set of movement and memory parameters required to generate an adaptive, migratory behavior. This includes the ability to learn to migrate from non-migratory initial conditions, simulating the release of naive animals in a seasonal environment (Jesmer et al., 2018); to lose the propensity to migrate if the resource distribution does not require it, also a commonly observed phenomenon (Wilcove and Wikelski, 2008); and to assess the resilience or fragility of a migratory population against changing resource distribution dynamics, including both stochasticity and trends in spatial and temporal distributions, mirroring the effects of climate change (Park et al., 2020).
We anticipated that under many conditions a blending of tactical (i.e., direct response to resource availability or perception) and strategic (i.e., memory-driven and forwardthinking) behavior will help foragers navigate dynamic, seasonal environments. Over-reliance on either strategy should be maladaptive. We further anticipated that a shorter-term memory updating is needed to navigate trends in resource spatial distribution and temporal distribution (phenology), but that a longer-term reference memory is needed to navigate resource distributions that are stochastic (Lin et al., 2021). Similarly, we anticipated that a balance between very low sociality and extreme sociality would lead to the most resilient migratory process.

Memory Movement Model
In designing our study, our goal was to develop a minimal heuristic in which the following processes were explicitly modeled: (1) Random or exploratory movement, (2) attraction to resources, (3) sociality in the movements, (4) a long-term (or reference) memory of large-scale movement behavior, and (5) a short-term (or working) memory that updates movement behavior based on recent experience.
A diffusion-advection equation provided a computationally efficient and versatile framework for examining just such a system. We consider a population moving in one dimension in a constrained domain D and distributing itself according to the following equation: (1) where u represents the population distributed in time and space. The first term is the diffusion term, capturing the fast time-scale exploration and "random" movements of individuals, with ε is the diffusion rate.
The second term represents the attraction to a dynamic resource h, with the proportionality of the advection to the gradient of the resource given by the parameter α (note, the population and resource distributions are functions of both space and time u(x, t) and h(x, t) -we omit the dependent variables in the notation for brevity). This is the well-studied standard chemotaxic resource-following behavior. We borrow the general notation from earlier related work expanding diffusion-advection models to incorporate non-local information (Fagan et al., 2017) and behavioral switching .
The third term captures the collective or social advection term of the population via a non-local, density dependent function v s (u, x). If this function takes the form of a convolution around a non-local kernel k, i.e., v s (u) = k(x) * u(x), and if that kernel is odd, an attractive or "swarming" behavior can be generated (Mogilner and Edelstein-Keshet, 1999). We use the kernel analyzed by Mogilner and Edelstein-Keshet (1999): The convolution of u with this kernel has the property of pushing the population in a positive direction when x < u , and in a negative direction when x > u , where u is the mean location of the population. The parameter λ is a length scale of sociality, roughly one-half the size of the swarm, and β is a parameter that quantifies the overall strength of sociality. Finally, the last term captures the direct advection that emerges from a memory-driven migratory behavior. This term evolves with a set of parameters θ y that slowly change each year y ∈ {0, 1, 2, ...}, i.e., the count of periods τ : y = ⌊t/τ ⌋. The migration is specified by six parameters θ : the timing of the start and duration of two anticipated seasons (e.g., summer and winter) t 1 , t 1 , t 2 , t 2 , and the spatial coordinates of the population centroid for each season x 1 and x 2 . The remembered migratory speed term is a simple step function given by: t > t 1 and t ≤ t 1 + t 1 s 12 ; t > t 1 + t 1 and t ≤ t 2 0; t > t 2 and t ≤ t 2 + t 2 s 21 ; t > t 2 + t 2 or t ≤ t 1 where the migration speeds s 12 and s 21 from the respective ranges are set such that they arrive at x 1 at t 1 , depart at t = t 1 + t 1 , arrive at x 2 at t = t 2 , and depart at t 2 + t 2 . Thus, s 12 = (x 2 − x 1 )/(t 2 −(t 1 + t 1 )) and s 21 = (x 1 −x 2 )/(t 1 −(t 2 −τ + t 2 )). This step-like migration function is a one-dimensional version of the migration parameters estimated for individuals  and populations  in empirical studies.
We consider these six parameters to be the known or remembered determinants of the migratory behavior, with an initial set θ 0 determining the reference migration behavior. This reference migration is updated each year by the experience of the population. To perform this updating, we estimate a new set of parameters θ y after each year, and combine these new parameters with the reference parameters according to the following weighted mean: where each of the six parameters is updated according to Equation 3 identically. The estimates θ y are obtained via a least-squares minimization of the migration track (m(t, θ ) = t 0 v m (t ′ , θ y ) dt ′ ) against the spatial mean of the population process in year y (i.e., u(t) = X u y (t, x)dx). The parameter κ ∈ (0, 1) captures the reliance on that long-term memory. When κ = 0, all of the actionable memory is from the preceding year. When κ = 1, the actionable memory is entirely the reference memory.
The model is confined to a one-dimensional bounded domain [−χ, χ], with no flux outside of the boundaries. Formally, this no-flux condition means the following conditions must be met In practice, the design of our resource space (see below) and other parameterization lead to 0 or near 0 values of both h(x) and u(x), and the simpler ∂u(−χ, t)/∂t = ∂u(χ, t)/∂t = 0 boundary condition provides a good approximation.
As there are no birth or death processes, the total population remains fixed and constant, for convenience integrating to 1. Furthermore, the parameters remain constant throughout time, with no adaptation or mutation-selection process. Our interest is in the ability of a fixed set of movement and memory parameters to navigate an intra-and interannually dynamic, seasonal environment.

Seasonal Resource
We ran this model on a spatial domain x ∈ [−100, 100], and a periodicity τ = 100 (i.e., 100 day years). We were interested in an approximately periodic resource dynamic, i.e., one in which h(x, t) ≈ h(x, t − τ )). We generated two types of resource distributions. A "non-surfable" resource (island resource), and weakly surfable resource (drifting resource). Both are characterized by a peak in time and space centered at m x at m t , and −m x at τ − m t (for example, locations 30 and −30 at times 25 and 75, respectively). These pulses have a shared time scale of duration s t and a spatial scale of extent s x , the standard deviation in the time and space dimension, respectively. The island resource is simply two uncorrelated bivariate normal distributions where is the bivariate Gaussian distribution function, and the normalizing constant K is selected such that the average total amount of resource throughout the year is 1, i.e., 1 τ T X h(x, t)dx dt = 1. The drifting resource differs from the island resource in that the total amount of resource at any given time X h(x, t)dx = 1. This property is attained by distributing the resource as a rescaled beta distribution, where the shape and scale parameters vary sinusoidally in such a way as to make the standard deviations and means match the desired values of m x , m t , s x , s t (see Supplementary Materials for details). Both types of resources are illustrated in Figure 1.
Within a given year, the resource is entirely symmetric: h y (x, t) = h y (−x, τ − t). However, in scenarios exploring climate change we allow the peaks to vary with directional trend and stochasticity according to: m x (y) ∼ N(µ x + γ x y, σ x ) and m t (y) ∼ N(µ t + γ t y, σ t ), where the µ, γ and σ terms are the mean, slope and variance, respectively, for the location and time duration of the pulse. Thus, if γ = 0 and σ = 0, the conditions are constant across years and if γ x > 0 there is a shift of the resource toward the extremes of the domain. While we did not explore phenological shifts in timing, those can readily be modeled as well. These trends model the pole-ward shift of peak resources and the earlier spring phenology occurring with a warming global climate (Renner and Zohner, 2018). The spatial and temporal scales of the resource peak (s x and s t ) remain constant in all of our simulations.

Metrics
The main metrics we were interested in are migration mismatch, foraging efficiency and adaptation to directional trends.
Migration mismatch captures the combined difference between the migration phenology and the resource phenology in time and space. Spatial mismatch MM x is the absolute difference between the migration targets and the resource peaks: Temporal mismatch is the difference between the arrival time and the peak of the resource if arrival is post-peak, the difference between the departure time and the peak of the resource if departure is pre-peak, and 0 if the seasonal duration spans the peak, i.e., MM t = max{t 1 − m 1 , m 1 − (t 1 + t 1 ), 0} + max{t 2 − m 2 , m 2 − (t 2 + t 2 ), 0}. Thus, the total mismatch is the sum of these: TM = MM x + MM t . A mismatch of less than 1 is essentially perfect, a mismatch of 1-5 we consider excellent, and beyond 50 the system can be said to have failed to keep track of the resource dynamics.
To quantify the foraging efficiency, i.e., the organisms' ability to track the distribution of the resources over space and time, we use a continuous form of the Bhattacharyya coefficient (Bhattacharyya, 1943) which quantifies the similarity between two distributions. We compute this coefficient at every time point in a given year, and take the mean across the equilibrium year to determine foraging efficiency (FE). Thus, the foraging efficiency index is: where the spatial integral is taken over the domain. This metric is constrained to be between 0 and 1. For simulations with a constant resource, we ran the model until a quasi-equilibrium (stationary) state was achieved, i.e., where the Bhattacharya index of the population distribution across subsequent years reached a value of 0.99999. Once stationarity was attained, we computed the migration mismatch and foraging efficiency metrics, as well as the number of years required to reach stationarity.
For numerical runs with climate change, we first run a simulation with a given parameter set until stationarity, as above, and then begin shifting the location of the resource poleward with a slow, moderate or rapid trend (γ x = 0.25, 0.5, and 1, FIGURE 1 | Examples of various seasonal resource distribution functions, contrasting short duration, but wide pulses (σ t = 3, σ x = 12; left panels), long duration but spatially concentrated pulses (σ t = 12, σ x = 3; right panels), and isolated resource pulses (upper panels) from the weakly drifting resource (lower panels). The total amount of resource is identical across all scenarios. In the weakly drifting resources, the total amount is constant at all times, and uniform in the middle of the phase (time = 0, 50, and 100).
Frontiers in Ecology and Evolution | www.frontiersin.org respectively), and/or by adding stochasticity (spatial standard deviation 3, 6, 9, or 12). For stochasticity analyses, we compare foraging efficiency across a range of the reference memory parameter κ. For analyses that included directional trends, with or without stochasticity, we quantified the ability of the system to keep track of climate change with a spatial adaptation (SA) index. This index is the ratio of the slope of the memorybased migration location over time, i.e., SA = γ /γ x where the adaptation slope estimate is the regression coefficient of the spatial coordinate of the migration against time (i.e., m x,i = γ x i + m x,0 , where i is the year), and γ x is the rate of drift of the resource peak (Table 1). An SA equal to 1 suggests that the process is keeping up with climate change, an SA of 0 indicates that the process is not responding at all to climate change. Values greater than 1 (super-adaptation) are possible, as are values less than 1, which correspond to a loss of migration behavior. All movement model parameters, resource parameters, and metrics are summarized in Table 1.

Simulation Studies
We explored this model using numerical differencing of a system of ordinary differential equations (ODE's) approximating the PDE in Equation (1) with the Runge-Kutte algorithm using the deSolve (Soetaert et al., 2010) and ReacTran (Soetaert and Meysman, 2012) packages in R. We additionally used the nlsLM function in package minpack.LM (Elzhov et al., 2016) for robust and fast annual estimation of the migration parameters. The complete code is available as an R package (memorymigration) available on GitHub at https://github.com/EliGurarie/memorymigration and as an interactive Shiny application at https://spot3512.shinyapps.io/ memorymigrationshinyapp/.
We assessed a wide range of parameter values and resource geometries and dynamics with the goal of answering the four main questions: (1) Can this model adapt to a discrete shift in peak resource location and timing? What is the relative role of memory and sociality for adaptation? (2) Can this model acquire a migratory behavior from a non-migratory initial condition? (3) What is the role of a reference memory for dealing with stochastic resource dynamics? (4) Can this model adapt when the resource peaks shifts in space? Details of parameter combinations and reported metrics are provided in respective results sections.

Adaptation to Resource Phenology
The ability of this system to attain a stable, migratory state that matches the dynamics of the resource is illustrated in Figure 2.
In the illustrated scenario, it takes nearly 40 years to attain an equilibrium, and the eventual steady state is one where the centroid of the migration lines up exactly with the centroid of the resource, and the arrival timing coincides with the peak of resource availability. Notably, the path to this equilibrium is somewhat indirect, with the later winter range taking more time to stabilize than the earlier summer range. The eventual steady state is one where the foraging efficiency is relatively high, near 0.6 compared to an initial value of 0.3. However, the increase in the foraging efficiency was not entirely monotonic, as the system moved through some slightly sub-optimal stages in adjusting its migration behavior. We ran this process for 8,100 parameter combinations crossing different values of the movement process (α, β and λ) and resource characteristics (σ x and σ t ), and present the total mismatch (TM) against all those combinations in Figure 3. In all of these simulations, memory was entirely recent (κ = 0), since there can be no benefit to relying on a sub-optimal reference memory. We compared a set of diffusion rates ε between 1 and 8, but only illustrate results for ε = 4.
A well-matched migration phenology (TM < 5) occurred under very many combinations of parameter values, but all parameters play interacting roles. Among the more intuitive results are that greater values of α (resource following) lead to an improved ability to match the migration. Resource peaks with larger spatial extent (higher σ x ) are generally better for migration matching.
Less intuitive was the high importance of the sociality parameters, in particular the spatial scale of the swarming. Higher levels of social attraction (β) led to improved migration matching except in those cases where the sociality scale λ was high. Thus, for example, at λ = 20, no simulations at β ≥ 200 managed to acquire or maintain a matched migration. However, at λ = 50 or 100, the migration was slightly better matched at high values of β (Figure 3). The spatial extent of the swarm was a remarkably significant variable. Smaller swarms were able to match migration only at low values of social attraction (β = 200), and relatively high values of resource attraction (α ≥ 600). Random forest analyses, whether on the log of total mismatch or on the classification of a perfect match, uniformly show that the most important variables (Breiman, 2001) were α and λ (4.14 and 4.02 proportional increase in MSE), and the least important was σ t , with a 0.5 proportional increase.
Overall, foraging efficiency was strongly correlated with migration matching, as expected. At high mismatch (> 50), foraging efficiency was low (mean 0.29, s.d. 0.16) compared to the near-perfect matching migrations (mean 0.58, s.d. 0.14). However, somewhat higher mismatch (1 to 5) showed an even higher overall foraging efficiency (mean 0.62, s.d. 0.18-see also Figure 4). Figure 5 illustrates the ability of the model animals to learn to migrate in a weakly drifting resource environment with a narrow pulse of resource peaking at 30 and -30 (at days 25 and 75), but a uniform distribution of resource at times 0 and 50. In order to learn to migrate, the system needed to have a higher exploratory impulse (higher diffusion constant ε), a stronger resource advection (higher α) and somewhat weaker sociality (lower β). The qualitative behavior of this process was to start drifting toward the summer resource, while slowly developing a weak pulse toward the winter resource as well. After first locking in on the summer resource, the winter migration, driven both by . The color scheme reflects the total mismatch, i.e., the sum of the absolute differences between the migration timing and locations from the resource peak. The white squares represents parameter combinations where the PDE could not be solved for artifactual numerical reasons, that all correspond to a failed adaptation (high mismatch). high diffusion and high resource following, slowly extended itself until both narrow peaks of resource were consistently reached.

Learning to Migrate
The model had, in general, a difficult time learning migration from a non-migratory initial condition. Out of 4,047 successful runs, only four attained mismatch below 1, and 130 below 5. Conditions that were more conducive to learning migration were pulses of longer duration (high σ t ), but smaller in scope (low σ x ), suggesting that the feedback that encourages migration needs to be compact in space but long enough in duration to lock in to the memory.

Directional Climate Change
To assess the ability of the system to adapt to a trending climate, we generated scenarios with slow, moderate and fast outward directional shifts in the resource peak (0.25, 0.5, and 1 units/year, respectively). We then assessed 40 parameter combinations for FIGURE 5 | Example of model learning to migrate. The resource is a "weakly drifting" resource and the initial (year 0) condition is non-migratory. The simulation was run for 100 years, and a sampling of those years (labeled) are presented in (A): all years from 0 to 10, followed by 20, 40 and 100. Otherwise, panels are as in Figure 2. Additional parameter values were ε = 5, α = 500, β = 50 and λ = 40. each of those scenarios, high and low values of resource following (α = 400 and 100), high and low values of sociality (β = 400 and 100) and 10 values of the spatial scale of sociality (λ = 20 to 200). The spatial and temporal scale of the resource pulses were fixed to σ x = 12 and σ t = 6, a combination which analyses in section 3.1 showed were generally "easy" to adapt to. We computed the adaptation index and foraging efficiency for each of the 120 runs (Figure 6). We were interested in the dynamics against λ due to the consistently high importance of this parameter for matching migration in steady states. Our main index of interest was the spatial adaptation (SA) to trends.
As Figure 6 shows, higher values of resource following (α = 400; orange circles) are nearly universally better for keeping up with climate change (SA values near 1). Furthermore, when combined with high sociality (β = 400; right panels), nearly all parameter combinations do a good job keeping up with climate change (SA values ranging between 0.53 and 0.85 for a swarm size greater than 50). However, that maximum value is still less than 1, suggesting that truly matching a steadily drifting trend is very difficult. Smaller social spatial scales (λ < 50) have a very hard time adapting when the social attraction is high, but do fairly well when social attraction is low. Larger sized swarms do progressively worse across more parameterizations, e.g., in the most rapid climate change scenario, the SA drops from 0.83 to -0.13 as the swarm increases in size from 40 to 200 (encompassing, essentially, the entire spatial domain). FIGURE 6 | Adaptation to a steadily drifting resource. In three scenarios, the spatial coordinates of the resource drift by 0.25, 0.5, and 1 unit per year (top to bottom, respectively). The y-axis is the spatial adaptation index (SA), i.e., the trend of the memory-driven migration divided by the resource drift trend. Values near 1 indicate a behavior that keeps up with climate change, values near 0 indicate no change in migration behavior, and negative values indicate a trend that is opposite to the climate trend. We compare across spatial scales of sociality (λ-x-axis), for low and high values resource following (α = 400 and 100-orange and blue dots) and low and high values of sociality (β = 100 and 400, left and right panels). The size of the circles is proportional to the foraging efficiency of the resulting parameter combinations. The bottom-right boxplots indicate the final year foraging efficiency against SA; purple and blue boxes indicate the highest values, orange and gray lower values.
A rather more dramatic pattern is visible for the lower foraging attraction scenario (α = 100; blue circles). Notably, no parameter combination at this value comes close to keeping up with the rapid climate change (SA range -0.64 to 0.13). For slower climate change, however, there is a window of values for the swarm size between 40 and 80, where the SA exceeds 1, but then crashes quite rapidly to negative values of SA as that swarm size increases. These "super-adaptive" processes indicate a unique sweet spot where a swarm is large enough to capture and adapt to the drifting resource, but not so large that the information gathered in a given year is too weak to adjust the migratory behavior in a following year.
As anticipated, better adaptation to the drifting resource correlated strongly with higher foraging efficiency (inset boxplots).

Reference Memory and Stochasticity
While recent memory can be helpful for adapting to a single novelty or a smoothly changing conditions, we hypothesized that a more conservative approach that relies on a reference memory may be beneficial when conditions change stochastically. We tested this hypothesis by solving a set of models across a range of κ values fom 0 (all recent memory) to 1 (all reference memory). In these scenarios, we ran the system for as many years as needed with no stochasticity to acquire a stationary state (i.e., similarity index greater than 1-1e-6). We then used the stationary state as the reference memory, and then ran the process for an additional 50 years with a stochasticity (i.e., standard deviation in peak location of the resource) ranging from 0 to 12, and present the resulting average foraging efficiency (Figure 7).
Overall, as expected, the greater the stochasticity, the lower the foraging efficiency. Further, as we predicted, highest level of κ can significantly help foraging efficiency, with some variation across the spatial scale of sociality, especially in more highly stochastic scenarios. When that scale of sociality is high enough (λ = 120, blue colors) there is greater probability of overlap with a stochastic resource, and a conservative, stable migratory regime is much more beneficial in the long run.

Stochasticity and Trends
We added 30 years of directional trends to the variously stochastic process described above, and assessed the adaptation index against the reference memory parameter κ (Figure 8). Over-reliance on reference memory (κ = 1) by definition does not allow the system to keep up with climate change, leading to an adaptation index of 0. However, in many cases a balancing of recent and reference memory (κ value between 0.6 and 0.8) in many cases was slightly but significantly better than relying entirely on recent memory. The smaller spatial scale (in the selected parameter space) does a generally better job than the larger spatial scale at lower stochasticity. At higher level of stochasticity, however, the larger spatial scale outperforms the smaller spatial scale, which completely loses track of the climate change.

DISCUSSION
Animals navigate complex, dynamic and patchy environments. When there is a strongly localized and seasonal component to the resource dynamics, movement strategies limited to straightforward resource-following taxis necessarily fail to efficiently exploit available resources. It is in these cases, quite common in the natural world, that seasonal migration FIGURE 7 | Foraging efficiency (FE) across various values of reference memory κ (x-axis) for increasing amounts of interannual stochasticity (ψ x , panels left to right), and two values of sociality spatial scale λ = 80 and 120. For the processes with non-zero stochasticity (ψ x > 0), the process was run 90 times for values of κ. Points represent the average of the FE's across all 50 years and 90 replicates. In these scenarios, the resource following parameter α = 100, the social attraction β = 400 and diffusion ε = 4. Note that the y-axis scale is constrained over a relatively narrow range (0.6-0.7).
FIGURE 8 | Role of reference memory in adapting to climate change for increasingly stochastic resource dynamics. We ran the model with a moderate rate of climate change (mean shift: 0.5/year) at five increasing levels of stochasticity (inter-annual standard deviation of resource peak 0, 3, 6, and 12, left to right panels). For non-zero stochasticity, we ran the process 30 times and present the mean and standard error of the spatial adaptation index across various values of the reference memory parameter κ: where κ = 0, the system modifies its migration based entirely on recent experiences; at κ = 1, the memory never changes from the reference memory. Other parameter values are resource following α = 100, social attraction β = 400, and diffusion ε = 4. becomes a viable, even necessary, strategy. However, when resources start shifting in space and time-as is occurring at an accelerated pace with recent global climate change-the migration phenology itself must exhibit some plasticity. It is our conjecture that this plasticity is facilitated by a memory-driven process in which recent experiences inform strategic behaviors in subsequent years.
By allowing a population to adjust its migratory behavior based on recent experiences with the resource location, the model we presented here emulated (a) the successful navigation of an environment with temporally and spatially isolated seasonal resource patches, (b) the emergence of a migratory behavior from an essentially resident or naive initial condition, and (c) some intrinsic robustness to changes in those environmental resources, whether steadily shifting trends or inter-annual stochasticity. The relatively simple, social and memory-driven mechanism was able to adapt to long-term changes in resource dynamics, even with inter-annual stochasticity, and may thereby provide a framework with which the interaction of memory, movement, social and resource dynamics can be further explored. Importantly, our model was in no ways evolutionary, as it contained no birth-death processes or selection pressures. Thus, we used foraging efficiency as a convenient metric of the utility of migration, though this was not a measure explicitly maximized by the model. Other metrics, such as foraging efficiency in a given season, or probability of survival or reproduction relative to resource availability (Bauer et al., 2020) may respond differently across model parameters and could be useful in understanding the relative success of alternative migratory strategies in different contexts. However, the overall annually averaged foraging efficiency metric provided the broadest linkage between resource dynamics and animals' locations and was consistent with the minimal biological assumptions and generality of our framework.

Adaptation and Resiliency
Our goal was to understand the combinations of factors that lead to a resilient migration behavior. The model we describe was a final iteration of a sequence of models which failed to develop or maintain social migration behavior. For example, in earlier versions memory was modeled as an attractive advection mathematically identical to the resource attraction, but with the attractor being the location of the population in previous years. These models proved to be inefficient at generating a consistent social migratory behavior, i.e., only under very specific parameter combinations and "easy" conditions was a migratory equilibrium attained, and that equilibrium was highly unstable to perturbations. Only a clear, directed advective process with an explicit seasonal signal (i.e., the remembered migration timing, rates, and targets which were remembered in our model) could generate the patterns we aimed to capture. This suggests, somewhat indirectly, that migration behavior is unique as a fundamental, long-term, and risky strategy, profoundly different from the kind of tactical resource response which governs shorter-scaled animal redistributions.
Similarly, iterations of the model that did not have some amount of social cohesion tended to diffuse away without establishing a consistent, migratory stationary state. In fact, sociality parameters-in particular, the spatial scale λ-were, unexpectedly among the most important parameters for determining the resiliency of the process. Populations with small spatial scales tended to have a more difficult time locking in to an adaptive migratory pattern, and only when social attraction was relatively weak. On the other hand, overly large spatial scales compromised the ability of the process to track climate change, due to a dilution of the population's ability to concentrate over available resource patches and remember the corresponding benefits.
The ability to adapt a migration also depended strongly on properties of the resource dynamics. In particular, the reinforcement of memory and foraging is strongest when patches are concentrated in time, but relatively large in space. Interestingly, in most stable patterns, the eventual targeted migration arrival time coincided with the peak, rather than the beginning, of the resource dynamic. This indicates that the longdistance social migration behavior may be particularly reinforced when the targeted resource is very sudden. This is the case for the rapid green-up that occurs in high latitudes as snow recedes in tandem with extended day lengths leading to an intense greenup period (Park et al., 2020) or, for example, when resources are linked to the short-duration early blooming phenology of very particular plants (Post and Forchhammer, 2007;Renner and Zohner, 2018).
Even with no strong intrinsic propensity to migrate and a weak phenological resource pulse to follow, our model captured the ability to acquire a strong and adaptive migration behavior (Figure 5). Learning migration, however, requires a very strong resource attraction, higher levels of exploratory behavior (e.g., diffusion, and larger spatial scale of sociliaty), and-oftenmany more years, findings that echo empirical observations (Jesmer et al., 2018).
Despite the ability of the process to adapt under many stable conditions, our migration model (and, perhaps, migration behaviors in general) can also be considered somewhat fragile. Under many shifting conditions, e.g., increasing stochasticity, rapidly shifting resources, a shift in some of the system parameters, or even a shift in the spatial and temporal extent of resources, migration can collapse and turn into a non-migratory, residential behavior (Figure 3). This sensitivity may explain why partially migratory populations are so common and, apparently, evolutionary stable (Berthold, 1999;Chapman et al., 2011), as well as the wide range of migration plasticity shown in wild populations, even within a species (Xu et al., 2021).

Biological Interpretation of Parameters
Diffusion-advection models of animal movement and redistribution are grounded in the general idea that animal movements, somewhat like movements of physical particles, combine random (diffusive) components with directed (advective) components (Skellam, 1951;Turchin, 1998;Okubo and Levin, 2001). While direct relationships between diffusion models and movement data are somewhat tenuous (Gurarie and Ovaskainen, 2011;Potts and Schlägel, 2020), as a theoretical tool for exploring processes they are invaluable for their versatility and the relative ease of numeric computation of the partial differential equations (PDEs) that describe them mathematically.
Despite its evident abstraction, our goal was to develop a model where all parameters have well-defined biological interpretations. The diffusion (ε) captures short time-scaled randomness of movement, reflecting exploratory and short-term dispersive behavior. The foraging advection strength (α) captures the attraction of the population to better quality resources at a relatively large scale. These two parameters, the basic ingredients in diffusion-advection models of animal movement, have direct parallels to empirically estimated properties of animal behavior: diffusion is closely related to families of random walk models (Gurarie and Ovaskainen, 2011) while the advective taxis is related to the step and resource selection functions that are routinely estimated from movement data (Potts and Schlägel, 2020). The spatial scale of the social group (λ) captures the spatial extent of the population, i.e., a population-level home range (Noonan et al., 2019). Diffusion-advection models can also be interpreted as a probabilistic description of a single individual's movement. In this case, λ would correspond to an individual home-range and β would be an individual's tendency to be drawn to the center of that home range, akin to an individual migratory Ornstein-Uhlenbeck process .
The sociality parameter (β) quantifies the strength of an individual's desire to approach the center of the social group. While this parameter is not typically measured, it may in principle be possible to estimate in a manner analogous to a step-selection function by replacing environmental variables with presence of conspecifics as a covariate. The ratio between α and β can be interpreted as the relative importance of foraging to social cohesion, which appears to be important in predicting the resilience of migration.
Migration timing, rate, and seasonal range location parameters can be straightforwardly estimated from movement data (Cagnacci et al., 2015;Gurarie et al., 2019) and synchrony of migration timing and site fidelity are well-documented for many migratory species (Joly et al., 2021). Thus, for example, Gurarie et al. (2019) explicitly estimated the ranging area, timing, and seasonal range locations for migratory caribou, identifying the kind of inter-annual variation that is reflected in the stochastic scenarios explored here, as well as trends in timing.
The reference memory parameter κ is, of course, impossible to observe directly. Our model does, however, allow us to explore in an heuristic way the conditions under which a strong cultural tendency to migrate with certain fixed patterns can help a population hedge against stochasticity (Abrahms et al., 2019;Fagan, 2019). An extremely conservative behavior is the best way to hedge against stochasticity with no directional changes (high κ values in Figure 7), as there is no benefit to change behavior based on recent experiences if they provide no information about future conditions. However, this extreme conservatism is, by definition, incapable of adapting when there is a consistent shift in resource distribution (Figure 8). In cases where both processes are occurring, we did see a slight improvement in adaptability when long-term reference memory was balanced against a strong response to recent experience (see peaks in Figure 8).
Clearly, our exploration of the model was not exhaustive. We did not explore, for example, the resilience of the migration process to changes in resource timing, which would correspond to the widely observed earlier onset of spring as measured by green-up and flowering phenology (Cleland et al., 2007). We hope that making the model available, including via the interactive interface, will facilitate further independent exploration of these processes.

Social Learning and Collective Knowledge
Models have shown that collective knowledge is important, if not essential, to the evolution and process of migration (Guttal and Couzin, 2010;Shaw and Couzin, 2013;Berdahl et al., 2018). Many migratory organisms are social, and social learning is an acknowledged, non-genetic method for transmitting information (Kashetsky et al., 2021). Furthermore, the general role of social learning for improving a population's ability to track resources has been studied not just in animal systems, but in synthetic systems inspired by social behavior of animals such as optimization heuristics algorithms and the study of swarm robotics (Şahin, 2005;Brambilla et al., 2013). Because our model is not individual-based, we can not identify any specific mechanism (e.g., leader-follower) of social information transfer. But, in a generic way, our model assumes that migration is driven by a collective decision for the timing and locations of seasonal ranges, consistent with the known social and exogenous (e.g., daylength related) triggers for migration. Further, the underlying assumption of the migration "urge" is consistent with the strong endogenous programs to migrate, e.g., the seasonal restlessness known as Zugunruhe exhibited by many birds (Berthold, 1999;Helm, 2006). However, in its generic diffusionbased approach to randomness, our model indirectly captures individual-level variation in migration parameters, an inevitable property of any population-level process .
In contrast to the many individual-based models of the evolution of migration (e.g., Guttal and Couzin, 2010;Anderson et al., 2013;Shaw and Couzin, 2013), our model did not include any selection, inheritance or birth or death processes. For example, Anderson et al. (2013) explored the resilience of a population under selective pressure under persistent trends and increased stochasticity of a drifting optimal resource window, showing that a certain amount of heritable phenotypic plasticity is necessary to adapt successfully to climate change even at the cost of efficiency. Our model underscores the fact that some level of resilience and adaptability can be attained with a purely cognitive process that balances sociality with long and short term collective memory. Importantly, this knowledge can be transmitted through social and cultural, rather than genetic, pathways. The high level of sociality among migratory animals, as well as multi-annual parent offspring bonds, are an evident pathway for that kind of transmission. As with those evolutionary models, however, it is clear that when changes are too rapid, no amount of cognition can help entirely mitigate against adverse outcomes. Furthermore, if behaviors are not sufficiently plastic (i.e., if κ is too close to 1), then adaptation is very difficult.
Given the slow scale of fitness selection and the constant change in environmental conditions, it is possible that certain inherent properties of populations, for example the "conservatism" captured by the κ parameter, are themselves selected for to maximize resilience over a long time scale in stochastic environments. The structure of the reference memory in our model was a rather simplistic approach to introduce conservatism or lag to the shifting migration parameters. In our model that reference memory is eventually entirely forgotten, whereas a more sophisticated approach would separate a slowly varying cultural memory, perhaps that is transmitted genetically or culturally, i.e., on the scale of generations, against shorter-scaled responses. In an evolutionary model, we might hypothesize that the overall rates of long-and short-term memory shifts would be related both to the scales of short and long-term fluctuation of the resource, i.e., the auto-correlation scale, strength of trends, and stochasticity of the resource dynamics.

Summary
Rapid environmental change, both global warming and increased anthropogenic development, is causing severe and dramatic impacts to the widespread and generally successful strategy of seasonal migration for many taxa, and the fate of many animal migrations is a topic of increasing concern (Wilcove and Wikelski, 2008;Kauffman et al., 2021). The ability of animals to respond to these changes depends deeply on their behavioral plasticity and cognitive abilities. The importance of those abilities is in direct proportion to the difficulty in studying them directly. By quantitatively exploring the properties of a heuristic model that distill many of the main properties of wild populations in dynamic and seasonal environments, we hope to have identified some broad patterns that might guide further empirical exploration of the cognitive underpinnings of adaptability and resilience.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
EG and WFF provided the original idea. EG and SP developed and ran models and analysis. All authors contributed to the article and approved the submitted version.

FUNDING
EG, WFF, GCC, and RSC were supported in part by NSF award DMS 1853478. EG and WFF were further partially supported by NSF grant IIBR 1915347.