Predicting performance of naïve migratory animals, from many wrongs to self-correction

Calculation of flight-step headings and movement

Terms defining flight-step movement, precision and geophysical orientation cues are listed in Table 1. Since seasonal migration nearly ubiquitously proceeds from higher to lower latitudes, it is convenient to define headings clockwise from geographic South (counter-clockwise from geographic North for migration commencing in the Southern Hemisphere). Assuming a spherical Earth, a sequence of N migratory flight-steps with corresponding headings, α_i, i = 0,…, N−1, the latitudes, ∅_i+1, and longitudes, λ_i+1, on completion of each flight-step can be calculated using the Haversine Equation⁷⁶, which we approximated by stepwise planar movement using Eqs. (1) and (2). For improved computational accuracy and to accommodate within flight-step effects, we updated simulated headings and corresponding locations hourly. A migrant’s flight-step distance ({R}_{{{mathrm {step}}}}=3.6{V}_{{mathrm {a}}}{cdot n}_{{mathrm {H}}}/{R}_{{{mathrm {Earth}}}}) (in radians), depends on its flight speed, V_a (m/s) relative to the mean Earth radius R_Earth (km), and flight-step hours, n_H. With a geomagnetic in-flight compass, expected hourly geographic headings are modulated by changes in magnetic declination, i.e., the clockwise difference between geographic and geomagnetic South^10,32.

Formulation of compass courses

For simplicity, we consider the case of a single inherited or imprinted heading. This can be extended to include sequences of preferred headings. Expected geographic loxodrome headings remain unchanged en route, i.e.,

Relative to geographic axes, expected geomagnetic loxodrome headings remain unchanged relative to proximate geomagnetic South, i.e., are offset by geomagnetic declination on departure (updated hourly in simulations)

As described and illustrated in detail by Kiepenheuer¹³, the magnetoclinic compass was hypothesized to explain the prevalence of “curved” migratory bird routes, i.e., for which local geographic headings shift gradually but substantially en route. A migrant with a magnetoclinic compass adjusts its heading at each flight-step to maintain a constant transverse component, γ′, of the experienced inclination angle, γ, so that error-free headings are (see Fig. S5 in ref. ³⁴)

In a geomagnetic dipole field, the horizontal (B_h) and vertical (B_z) field, and therefore also inclination, each depends solely on geomagnetic latitude, ∅_m:(gamma ={{{tan }}}^{-1}left({B}_{{mathrm {z}}}/{B}_{{mathrm {h}}}right)={{{tan }}}^{-1}left(2{{sin }}{phi }_{{mathrm {m}}}/{{cos }}{phi }_{{mathrm {m}}}right)={{{tan }}}^{-1}left(2{{tan }}{phi }_{{mathrm {m}}}right).) The projected transverse component, therefore, becomes

which can be substituted into Eq. (7) to produce a closed formula for magnetoclinic headings in a dipole as a function of geomagnetic latitude

with the expected initial heading, ({bar{{{{{{rm{alpha }}}}}}}}_{0}), and initial geomagnetic latitude, ∅_m,0, being constants. Equations (7) and (8) have no solution when inclination increases en route, which could occur following substantial orientation error or in strongly non-dipolar fields. We followed previous studies in allowing magnetoclinic migrants to head towards magnetic East or West until inclination decreased sufficiently^33,34,46, but also included orientation error based on the modelled compass precision.

To assess sun-compass sensitivity algebraically, and also to improve computational efficiency, we used a closed-form equation for sunset azimuth, θ_s (derived in Supplementary Note 3 and see ref. ²³),

where δ_s is the solar declination, which varies between −23.4° and 23.4° with season and latitude²³. Sunset azimuth is the positive and sunrise azimuth is the negative solution to Eq. (9) (relative to geographic N–S).

Fixed sun-compass headings represent a uniform (clockwise) offset, ({bar{{{{{{rm{alpha }}}}}}}}_{{mathrm {s}}}) to sunrise or sunset azimuth, θ_s,i (calculated using Eq. (9))

where the preferred heading on commencement of migration, ({bar{{{{{{rm{alpha }}}}}}}}_{{mathrm {s}}}={bar{{{{{{rm{alpha }}}}}}}}_{0}-{theta }_{{mathrm {s}},0}), is presumed to be imprinted using an inherited geographic or geomagnetic heading^2,10,30.

With a TCSC, preferred headings relative to sun azimuth are adjusted according to the time of day. In the context of sun-compass use during migration, Alerstam and Pettersson²² related the hourly “clock-shift” induced by crossing bands of longitude (∆h = 12 ∆λ/π), to a migrant’s time-compensated adjustment given the rate of change (i.e., angular speed) of sun azimuth close to sunset

resulting in a “time-compensated” offset in heading on departure ((varDelta bar{{{{{{rm{alpha }}}}}}}cong varDelta {{{{{rm{lambda }}}}}},sin phi), which Eq.(4)). Equation (4) results in near-great-circle trajectories for small ranges in latitude, ∅, until inner clocks are reset. The feasibility of TCSC courses over longer distances (latitude ranges) relies on two critical but little-explored assumptions: (1) time-compensated orientation adjustments are presumed to follow the angular speed of sun azimuth (Eq. (11)) retained from the most recent clock-reset site, and (2) to negotiate unpredictable migratory schedules, migrants are presumed to retain their preferred geographic heading on arrival at extended stopovers²².

Regarding the first assumption, time-compensated adjustments could also be influenced by proximate speeds of sun azimuth even when inner clocks are not fully reset. We, therefore, use distinct indices to keep track of “reference” flight-steps for clock-resets (c_ref,i) and time-compensated adjustments (s_ref,i). TCSC flight-step headings can then be written as

where θ_s,i represents the sunset azimuth on departures, c_ref,i specifies the most recent clock-reset site (during which geographic headings are also retained, i.e., ({bar{{{{{{rm{alpha }}}}}}}}_{i}={{{{{{rm{alpha }}}}}}}_{i-1})), and s_ref,i specifies the site defining the migrant’s temporal (hourly) rate of “time-compensated” adjustments (Eq. (11)). For TCSC courses as conceived by Alerstam and Pettersson²², reference rates of adjustment to sun azimuth are reset in tandem during stopovers, i.e., ({s}_{{{mathrm {ref}}},i}={c}_{{{mathrm {ref}}},i}), but we also considered a proximately gauged TCSC, where migrants gauge their adjustments to currently experienced speed of sun azimuth, i.e., ({s}_{{{mathrm {ref}}},i}=i).

Regarding the second assumption, retaining geographic headings on arrival at stopovers is not consistent with ignoring geographic headings between consecutive nightly flight-steps, and may be difficult to achieve while landing. We, therefore, examined a more parsimonious alternative (Fig. 7d, Supplementary Fig. 3) where migrants retain their (usual) TCSC heading from the first night of stopovers, i.e., as if they would have departed on the first night. This alternative also simplifies Eq. (12) to

where the index t_i−1 here represents the departure date from the previous flight.

Sensitivity of compass-course headings

Sensitivity was assessed by the marginal change in expected heading from previous (imprecise) headings, (partial {bar{alpha }}_{i}/partial {alpha }_{i-1}). When this is positive, small errors in headings will perpetuate, and therefore expected errors in migratory trajectories will grow iteratively. Conversely, negative sensitivity implies self-correction between successive flight-steps. Geographic and geomagnetic loxodromes are per definition constant relative to their respective axes so have “zero” sensitivity, as long as cue-detection errors are stochastically independent.

For magnetoclinic compass courses in a dipole field, sensitivity can be calculated by differentiating Eq. (8) with respect to previous headings:

All three terms in the denominator indicate, as illustrated in Fig. 3b, that magnetoclinic courses become unstably sensitive at both high and low latitudes, and any heading with a significantly East–West component.

Sensitivity of fixed sun compass headings is non-zero due to sun azimuth dependence on location (Eq. (9)):

The sine factor on the right-hand side in Eq. (14) causes the sign of (partial {bar{alpha }}_{i}/partial {alpha }_{i-1}) to be opposite for East to West or West to East headings, and tan θ_s also change sign at the fall equinox (due to solar declination changing sign). The azimuth term in the denominator indicates heightened sensitivity closer to the summer or winter equinox and at high latitudes, and, conversely, heightened robustness to errors closer to the spring or autumnal equinox (since ({{tan }}{theta }_{{mathrm {s}},0}to pm infty)). This seasonal and directional asymmetry is illustrated in Fig. 3c, e.

TCSC courses (Eq. (12)) involve up to three sensitivity terms, due to dependencies on sun azimuth, longitude and latitude:

The first square-bracketed terms in Eqs. (15a, b) are identical to the fixed sun compass (Eq. (14)), reflecting seasonal and latitudinal dependence in sun-azimuth. For headings with a Southward component (α₀ < 90°), the second bracketed terms are always negative, i.e., sensitivity-reducing, resulting in a broad range in latitude and headings with self-correcting headings (Fig. 3c–f). The third bracketed term in Eq. (15b) (for a proximately-gauged TCSC) is also negative, and in fact increasingly so until clocks are reset, but remains small in magnitude compared to the second term.

Spatiotemporal migration model

We wrote a model in MATLAB to simulate and assess the feasibility and robustness of each compass course to spatiotemporal effects on a global scale, based on our compass course formulations (Eqs. (5)–(12)). For the species simulations, we also incorporated spatiotemporally dynamic geomagnetic data (MATLAB 2020b package igrf)⁵¹, assuming a default season, fall 2000. For the generic migrant simulations, we assumed a geomagnetic dipole Earth, i.e., ignored variation in magnetic declination. Sunset azimuth was computed using Eq. (9) (this was, e.g., two orders of magnitude faster than the routine further requiring time of day and longitude used in ref. ³⁴). In all cases 10,000 individuals were simulated, until they either arrived in goal areas, passed 1000 km South of the goal latitude or exceeded the maximum number of steps, N_max. To avoid migrants overshooting narrow goal areas within a single flight-step, we assumed they could identify goal areas in flight (checked once per decile of flight-step durations; Table 2).

For the species simulations, optimal inherited headings were determined using the MATLAB nonlinear solver fminbnd, among candidate initial headings clockwise from East (−90° clockwise from S) to West (90°). For sun compass courses, which can potentially begin with Northward headings^34,64, we tested initial headings between NE (−145°) and NW (145°). Error scenarios were assessed for both 5°–60° directional precision among flight-steps, and biologically relevant scenarios with 5°–40° compass precision in 5° intervals (assuming equivalent precision in cue detection, transfers and maintenance), and both in the absence of and including 15° within-flight drift. To result in 15° expected drift per flight step, hourly concentration in drift was adjusted as an autocorrelated process with lag 1. We further assumed a default between-individual variability of 2.5° (κ = 525, circular length = 0.999). For the marsh warbler uncertainty analysis (Fig. 6), we varied this between 0° and 10° (κ = 33, circular length = 0.985) in 1° interval, and also tested 20° within-flight drift.

For the generic migrants, we simulated migration between 65°N–0°N and between 45°N–25°N to a goal with a radius of 500 km, for biologically relevant scenarios with 0°–60° compass precision in 1° interval, and both in the absence of and including 15° within-flight drift. To obtain all possible routes between these latitudes (Fig. 5a), we varied initial (inherited) headings in 0.5° intervals.

All relevant model parameters are listed in Table 2. These were chosen to match known studies and migration patterns. In some cases (e.g., marsh warbler and monarch butterfly), goal areas reflect plausible destinations from which migrants presumably use other cues to home or pin-point to even narrower known winter or passage areas^2,3. When flight-step distances and stopover durations were less known or certain, these were chosen to ensure modelled migration was consistent with known departure and arrival dates.

Generic migrants departed on September 15 ± 5 (mean and standard deviation), in 8-h flight-steps at speeds of 12.5 m s⁻¹. Sequences of 5 consecutive flight-steps were interspersed by stopovers of 5 ± 2 days.

Formulating migratory performance

Performance (arrival probability) of independent stepwise movement on a plane to a (circular) goal area of radius R_goal will approximate a cumulative normal distribution (erf function), modulated by the expected number of steps and the angular breadth to the goal area, which for long-distance routes is (beta ={{{tan }}}^{-1}left({R}_{{mathrm {{goal}}}}/{R}_{{{mathrm {mig}}}}right)cong {R}_{{{mathrm {goal}}}}/{R}_{{{mathrm {mig}}}}). Assuming uniform population headings and a sufficiently large number of flight-steps with directional precision ({sigma }_{{{mathrm {step}}}}), a first approximation of performance is

where the expected number of steps, (hat{N}), scales with the minimum (error-free) number of steps,

N₀, multiplied by a ratio of Bessel functions (Supplementary Note 2). From Eq. (16) we see that within the planar and high-precision limits, performance will increase with ({beta }_{{{mathrm {adj}}}}=beta sqrt{{N}_{0}}), which is the length-adjusted goal breadth (Eq. (3)).

In the equation for flight-step longitude (Eq. (2)), the secant factor (cosine of latitude in the denominator) reflects the poleward convergence of longitudinal meridians. This means that for any compass course, orientation errors at higher latitudes will exert a greater influence on overall longitudinal error:

Aggregated across entire migration routes, the effective longitudinal error will scale approximately as in a Mercator projection⁵⁸:

where ∅₀ and ∅_A are the initial (natal) and arrival latitude, respectively. To include latitudinal contribution to error, we modulated the multiplicative factor L according to route-mean orientation, (bar{alpha }),

For TCSC courses (bar{alpha }) was estimated as the average between the initial and final great circle bearings between the natal and goal locations. The spherical-geometry factor, G, is largest for purely Eastward or Westward headings (G = L > 1) and nonexistent for North–South headings (G = 1, reflecting no longitude bands being crossed). We expected this factor to affect compass courses differentially according to their error-accumulating or self-correcting nature.

We further modified the effective goal-area breadth to account for a (geographically) circular goal area on the sphere, i.e., effectively modulating the longitudinal component of the goal-area breadth at the arrival latitude, ∅_A:

To account for differential sensitivity among compass-courses, we generalized the normal many-wrongs relation between performance and number of steps, (1/{hat{N}}^{eta }), from η = 0.5 (Eqs. (3) and (16)) to

where b < 0 reflects iterative augmentation of errors and b > 0 self-correction, and s represents a modulating exponential damping factor, consistent with the limiting circular-uniform case (as κ → 0, i.e., ({sigma }_{{{mathrm {step}}}}to infty)), where no (timely) convergence of heading is expected with an increasing number of steps.

In assessing performance, we also accounted for seasonal migration constraints via a population-specific maximum number of steps, N_max (Table 2; this became significant for the longest-distance simulations with large expected errors, i.e., small ({{{{{{rm{kappa }}}}}}}_{{{mathrm {step}}}}=1/{sigma }_{{{mathrm {step}}}}^{2})). The probability of having arrived at the goal latitude can be estimated using the Central Limit Theorem:

where I_j is the modified Bessel function of the first kind and order j⁵³, and σ_C (the standard deviation in the latitudinal component of flight-step distance) can be calculated using Bessel functions together with known properties of sums of cosines^53,77 (Supplementary Note 2).

Regression-estimated performance

We fit the parameters in the spherical-geometry factor (Eq. (18)) and many-wrongs effect (Eq. (20)) according to expected performance, estimated as the product of sufficiently timely migration (Eq. (21)) and sufficiently precise migration, now generalized from Eq. (16), i.e.

This resulted in up to four fitted parameters for each compass course

1. i.

   an exponent, g, to the spherical-geometry factor (Eq. (19)), i.e., G^g, reflecting how growth or self-correction in errors between steps further augments or reduces this factor,

2. ii.

   a baseline offset, b₀, to the “normal” exponent η = 0.5, which mediates the relation between the number of steps and performance (Eq. (20)),

3. iii.

   an exponent s reflecting how decreasing precision among flight-steps dampens the many-wrongs convergence (Eq. (20)),

4. iv.

   for TCSC courses, a modulation, ρ, to the offset, b₀, quantifying the extent to which self-correction increases with increased flight-step distance R_step, i.e., ({{b={b}_{0}R}_{{{mathrm {step}}}}^{{prime} }}^{rho }) in Eq. (20), where ({R}_{{{mathrm {step}}}}^{{prime} })is the flight-step distance scaled by its median value among species.

Parameters were fit using MATLAB routine fitnlm based on compass course performance among species and seven error scenarios (5°, 10°, 20°, 30°, 40°, 50°, and 60° directional precision among flight-steps), for all combinations (including or excluding the four parameters). The most parsimonious combination of parameters was selected using MATLAB routine aicbic, based on the AICc, the Akaike information criterion corrected for small sample size⁵⁷. Null values for the spherical-geometry parameter were set to g = 1, and for the parameters governing convergence of route-mean headings b₀ = 0, s = 0, and, for TCSC courses, ρ = 0 (for loxodrome courses, ρ = 0 by default, i.e., was not fitted).

Statistics and reproducibility

Our simulation results, regression fitting and AICc-model selection are reproducible using the MATLAB scripts (see the section “Code availability”).

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Source: Ecology - nature.com