Asymmetry hidden in birds’ tracks reveals wind, heading, and orientation ability over the ocean

See allHide authors and affiliations

Science Advances  27 Sep 2017:
Vol. 3, no. 9, e1700097
DOI: 10.1126/sciadv.1700097


Numerous flying and swimming animals constantly need to control their heading (that is, their direction of orientation) in a flow to reach their distant destination. However, animal orientation in a flow has yet to be satisfactorily explained because it is difficult to directly measure animal heading and flow. We constructed a new animal movement model based on the asymmetric distribution of the GPS (Global Positioning System) track vector along its mean vector, which might be caused by wind flow. This statistical model enabled us to simultaneously estimate animal heading (navigational decision-making) and ocean wind information over the range traversed by free-ranging birds. We applied this method to the tracking data of homing seabirds. The wind flow estimated by the model was consistent with the spatiotemporally coarse wind information provided by an atmospheric simulation model. The estimated heading information revealed that homing seabirds could head in a direction different from that leading to the colony to offset wind effects and to enable them to eventually move in the direction they intended to take, even though they are over the open sea where visual cues are unavailable. Our results highlight the utility of combining large data sets of animal movements with the “inverse problem approach,” enabling unobservable causal factors to be estimated from the observed output data. This approach potentially initiates a new era of analyzing animal decision-making in the field.


Animal navigation has interested scientists for many decades but is now seen as an emerging field owing to animal tracking technology. This technology has enabled the movements of numerous species to be recorded over larger spatiotemporal ranges at high resolution using the Global Positioning System (GPS) or satellite tracking (1). However, the behavioral and cognitive mechanisms involved in long-distance animal orientation have yet to be resolved (2). Investigating the effects of extrinsic factors on navigational decisions by analyzing up-to-date, high-resolution animal trajectories requires a state-of-the art modeling approach.

Navigational decisions and environmental forces that determine movement trajectories can be estimated from additional sensors [for example, recording animal body orientation (3) and brain activity (4)] or remote sensing [for example, wind or current conditions (5, 6)]. These methods are undoubtedly useful; however, their applicability and/or resolution are often limited. For example, long-term recordings of animal body orientation continue to be technically demanding (3), and/or satellite-based wind information is spatiotemporally coarse [spatial resolutions are tens of kilometers, and temporal resolutions are 3 to 12 hours (5)]. These limitations hinder the ability to precisely quantify the navigational response of animals to their surrounding environments, because the external disturbance and the response of animals change dynamically and simultaneously at very small time scales (2).

Here, we used an “inverse problem approach” to estimate the causal factors affecting animal movement by using their tracks to quantify navigational decisions and environmental disturbance (for example, wind). The inverse problem approach is widely used in physics and engineering to estimate unobservable causal factors from observed output data (7). However, the use of this approach remains uncommon in biology (8, 9), especially in animal ecology, because it requires sufficiently large data sets. Another important hindrance is the lack of relevant underlying theory that describes how the model and its parameters that are to be inferred relate to the data sets. The biologging technique, whereby data loggers are attached temporarily to free-living animals in the wild, and the high-resolution data it produces now enable the use of the inverse problem approach in this research field. Some studies have used this framework implicitly to estimate either physiological factors [for example, the body density (10), air volume (11), and heat conductance (12) of animals] or external factors [for example, prey patch quality (13) and wind (1416)] by applying biomechanical models or optimal theory to biologging data, but these studies have not yet achieved the estimation of the internal process related to the animal’s decision-making in response to these factors.

We simultaneously estimated not only an external factor (ocean wind) but also the animal’s navigational decision-making (heading; that is, the direction in which an animal is oriented) from GPS tracking data (Fig. 1). Note that the movement of flying or swimming animals is affected by the motion of the fluids (flow; that is, wind or currents) they traverse, and the direction of their movement relative to the ground, recorded via GPS, can differ from the heading (17). Accordingly, animals constantly need to control their heading in a wind or current (flow) to reach the distant goals (18). If a bird were to align its heading vector with the preferred direction (northward), then the direction of the track vector would be northeast, deviating from the preferred direction (Fig. 1A). That is, the bird would be drifting due to the flow (known as full drift). The bird would be able to reduce the drift by directing its heading vector point somewhat into the flow (Fig. 1, B and C). This is generally referred to as flow compensation. The case in which the birds’ track direction coincides with the preferred direction is known as complete compensation (Fig. 1C); other cases are considered partial compensation (Fig. 1B). Compensation for flow drift increases the reliability of the bird reaching its destination and, when the flow is constant from the start to the end of the journey, reduces the time spent traveling and thus minimizes its total energy cost (17, 19). Birds flying over land are expected to use landmarks for referencing their position relative to the ground to prevent drifting due to lateral wind. Contrary to this, birds moving over the ocean such as seabirds, where landmarks are absent, do not have any referencing beacons (20). By investigating their orientation in response to wind (for example, drifted and compensation), we can elucidate both their navigational strategy and capacity.

Fig. 1 Conceptual framework of the inverse problem approach with biologging data.

The black and red lines indicate the process by which the movement path (the line at the bottom right) of an animal is generated and how the inverse problem approach is applied, respectively. In general, the animal’s movement is affected by both external factors (for example, wind or currents) and internal factors (for example, the internal state, animal navigation capacity, and motion capacity). Here, we simultaneously estimated information relating to animal heading (internal factors) and ocean wind (external factors) by analyzing time series data recorded by bird-borne GPS loggers. (A to C) Patterns of animal response to flow. (A) Full drift. (B) Partial compensation. (C) Complete compensation. The red, blue, and black arrows indicate the heading, flow, and track vectors, respectively. The green arrow indicates the preferred direction of the animal (the direction of the destination).

Using the inverse problem approach, we propose a new method to estimate the heading and flow vectors to reveal animal decision-making in response to the flow, using tracking data only. For this, we need a model that quantitatively relates the observation data (that is, track data) to the parameters to be estimated (that is, heading and wind vectors). We constructed such a model by using a (biased) random walk model, one of the basic theories for animal movement analysis that describe the direction and length of the animal track vector using probability densities (2124). This formulation has three advantages: (Property 1) The model can deal with random fluctuation observed in animal movement paths (25). (Property 2) Variation of the probability density parameter allows the models to describe various types of movement modes, from exploratory searching to directed movement toward a particular direction or location (23, 24). (Property 3) It is easy to fit the model to the real data and find the best-fit parameter values by formulating the maximum likelihood function and/or state-space modeling framework (2628). We propose a new model that is similar to previous models with respect to the above three properties, but which differs in terms of the implementation of the random walk, that is, “inversely” estimating the environmental factor (flow vector) and the internal decision-making (heading vector) of animals from the “resultant” animal tracks. We applied the method to the homing tracks of streaked shearwaters (Calonectris leucomelas) recorded at a resolution of one location per minute and investigated their responses to wind.


Model fitting to the track data

The intuitive concept of our model is shown in Fig. 2. We consider a situation in which an animal is moving toward a particular destination through a moving fluid. Our model describes the distributions of the heading and track vectors. It predicts that the heading vector is distributed “symmetrically” along its mean vector (Fig. 2A), but the track vector should be distributed “asymmetrically” along its mean vector via the effect of the flow (Fig. 2B). We can estimate the heading and wind vectors by fitting our model to the track vectors calculated from the tracking data. The detail of the model and the conditions required for model fitting are explained in Materials and Methods and the Supplementary Notes.

Fig. 2 Heading and track vector distributions of model prediction.

(A) Example of probability distributions of heading vector pH(uH, vH), distributed symmetrically along the mean heading vector (red arrow). The radius of the red dashed circle is equal to the mean air speed. (B) Example of a probability distribution of a track vector pT(uT, vT), gained by moving pH(uH, vH) to the wind vector (blue arrow). It is asymmetric along the mean track vector (black arrow). The radius of the white dashed circle is equal to the mean speed of the track vector.

Our model was applied to the homing portions of the GPS tracks of streaked shearwaters (N = 33; represented by the colored part in Fig. 3) (see Materials and Methods). Because the animal movement pattern and flow vector are expected to change spatiotemporally, we conducted the model fitting by selecting 51 min of positional data (Xt−25 to Xt+25), and the mean bird heading and wind (flow) vectors were estimated every consecutive minute at time point t. Thirty of the trips partially fitted the model [that is, fulfilled conditions 1 (the heading probability density is anisotropic) and 2 (the angle between the mean heading and the mean track vector is less than 90°)] (see Materials and Methods). Examples of tracks and estimated vectors are shown in Fig. 4 (A to C). Note that the asymmetry (or symmetry during light wind) of the track vector distribution predicted by the model (Fig. 2B) was also observed in the real data (Fig. 4, D to F). The other three trips did not fit the model at any point and were excluded from the analysis.

Fig. 3 Tracks of streaked shearwaters.

Tracks of streaked shearwaters (N = 33). In the tracks, the parts used for analysis, that is, those in which the bird is in the homing phase and not on land (the boundary is indicated by blue dashed lines), are indicated by the colored lines. The orange dot is the nesting colony on the study site (destination).

Fig. 4 Track vector distribution of real track.

(A to C) The black lines, red arrows, blue arrows, and orange dots indicate the homing tracks, the estimated heading vectors, the estimated wind vectors, and the study site, respectively. (D to F) The track vectors calculated from the position data in the 51-min time window [the part of the track colored green in (A) to (C) is plotted (green dots)]. The radius of the black dashed circles are equal to the mean speed of the track vector. The black, red, and blue arrows represent the mean track vector, the model-estimated mean heading vector, and the wind vector, respectively. The light blue dashed arrow represents the reanalysis wind data. (F) The upper right panel is a magnification of the estimated and reanalysis wind vectors in the main panel. The track vectors are distributed asymmetrically along its mean vector when a crosswind exists and almost symmetrically when the wind is weak, as predicted by the model in Fig. 2B.

Comparison between the wind vectors estimated from bird tracks and simulated by the atmospheric model

There was a significant correlation between the mean estimated wind vector and the mean reanalyzed data wind vector [generalized vector correlation coefficient of 1.29 that takes a value between 0 (no correlation) and 2 (complete correlation); N = 32; P < 0.0001] (Fig. 5A). The wind speed (P = 0.0019) (Fig. 5B) and direction (Fisher-Lee correlation coefficient, 0.20; P = 0.0013) (Fig. 5C) of the two velocities were also significantly correlated. Although the bird-based wind speed was strongly correlated with that of the reanalyzed wind data, the wind speed was 37% lower than that of the reanalyzed data (Fig. 5B). The difference between the bird-based wind direction and reanalyzed data was larger in light wind (fig. S1).

Fig. 5 Comparison between the wind vectors estimated by the model and by the reanalysis data set.

(A) Comparison of the wind vector (N = 32). The first and second rows show the mean wind vector of reanalysis data (wind that is calculated by the atmospheric simulation model) and the mean estimated wind vector (wind that is estimated using bird tracks), respectively. The length of the line indicates the wind speed, and the direction indicates the direction in which the wind blows. (B) Comparison of the wind speed. The red solid line is the regression line (y = 0.37x + 1.7), and the red dashed line represents y = x. (C) Comparison of the wind direction. The red dashed line is where y = x.

Wind compensation over the sea

Using the track vector direction and the estimated bird heading direction, we investigated the orientation strategy of streaked shearwaters over the sea. The regression method of Green and Alerstam (29) provided the degree of compensation and preferred direction (the direction in which the birds intended to fly under calm wind conditions) in three sections [northern coast (NC), offshore (OS), and western coast (WC)] (Fig. 6). When the track direction is linearly regressed against α (Fig. 6A), the slope of the line indicates the degree of compensation and the intercept indicates the preferred direction. Full drift, that is, when a bird does not compensate at all for wind drift, is characterized by a slope of 1. In all three sections, the slope was significantly lower than 1 (Fig. 6, B to D, and Table 1), indicating that the birds were compensating for wind drift irrespective of the distance from the coastline.

Fig. 6 Orientation strategies of streaked shearwaters.

(A) Definition of α. (B to D) Mean direction of the track vector plotted against α in three sections: NC (B), OS (C), and WC (D). The slope of the regression line corresponds to the degree of compensation, and the intercept (the value of the track direction at α = 0) corresponds to the preferred direction. The pink dashed line (slope, 1) corresponds to birds that do not compensate (full drift), and the gray dashed line (slope, 0) indicates a bird that completely compensated for wind drift. (E) (Left) Gray lines are homing tracks, and orange, light blue, and purple lines show the NC, OS, and WC sections, respectively. The orange, blue, and purple points are median locations (the location with median latitude and longitude of tracks in each section) of the NC, OS, and WC sections, respectively. Green arrows are the preferred direction, and the red point is the nesting colony (final goal). (Right) Preferred direction (green arrows) and nest direction from the median location (black dots) in each section. The green lines are 95% CI of the preferred direction. The patterns of wind response of birds in each section are shown. The red, blue, and black arrows indicate the heading, wind, and track vectors, respectively.

Table 1 Slope and intercept of the linear regression.

The directions of the colony from each of the median locations (the median locations of the tracks in each of the sections) (see also Fig. 6E) are also shown. PD, preferred direction.

View this table:

The degree of compensation and preferred direction differed among the sections. When the birds were flying relatively near the coast (NC and WC sections), the slope was lower than for offshore flight (OS section) and not significantly different from 0, indicating that the birds completely compensated for wind. On the other hand, when the birds were flying over the open sea (OS section), the slope was significantly higher than 0, indicating that the birds partially drifted, although they still compensated for wind to an extent (partial compensation). When the birds were flying far from the colony (NC section), the direction of the colony was within 95% confidence interval (CI) of the preferred direction. However, when the birds flew offshore or approached the coast (OS and WC sections), the preferred direction deviated from the direction of the colony in that it lay to the north of the colony (Fig. 6E and Table 1).


The biased random walk (BRW) model, the model that is commonly used for animal directed movement, assumed an explicitly or implicitly symmetrical distribution of the track vector along its mean vector (23). However, the track vector should be distributed asymmetrically along its mean vector, because flows distort the distribution. Rather, it is reasonable to assume that the heading vector is distributed symmetrically. We obtained the wind and birds’ heading vectors over the range across which the free-ranging birds moved, from GPS data only, by using the asymmetric distribution of the GPS track vector along its mean vector.

Two points should be considered with respect to our method. First, the estimated wind speed was lower than that of the reanalysis wind data. The reason for this was previously addressed [for example, the difference in the estimated height of ocean winds between the bird-based (less than 5 m) and satellite-based method (10 m) (15)]. Second, our method assumes constant heading speed (condition 3) and does not take into account animals adjusting their heading speed in response to the wind direction and speed for energetically optimal flight (3032) because the purpose of this study is to propose a simple model and to test whether the asymmetrical distribution of the track vector contains the information of the wind and heading vectors. Although including the airspeed adjustment in the model would complicate the model, it may be more realistic and valuable to include such parameters in future research. However, our assumption that the mean airspeed of a bird is constant irrespective of wind did not affect the validity of our findings because, first, the estimated wind vector correlated with the reanalysis wind data sets and, second, our method estimated the heading direction with good precision irrespective of condition 3, as shown in the numerical simulation (see note S2), and, third, because we used the regression method of Green and Alerstam (29) to estimate the preferred direction and the degree of compensation, we only required the track vector direction and the heading direction, which we successfully estimated (see above). Note that we do not need the airspeed, which we assume to be constant, indicating that the assumption does not affect our conclusions of shearwaters’ orientation strategy.

We showed that homing shearwaters could compensate for wind by adjusting their direction of flight such that they head into the wind depending on the degree of lateral wind. This homing behavior should be evaluated in the context of flight control and cognitive mechanism. When an animal compensates for crosswinds, the animal should head in a direction away from the preferred direction. This might be easy when landmarks are available; however, in contrast, it requires skillful flight control ability when the animal is moving over the ocean surface. In particular, shearwaters and albatrosses adopt a unique flying style—dynamic soaring, that is, the direction in which they are heading continuously changes within the scale of a few seconds (in this research, we estimated the mean values of the heading directions during a 50-min period)—which makes wind compensation more demanding than that in the case of flapping birds, whose heading directions remain within a small range. Thus, the shearwaters were likely to correctly evaluate the wind conditions they experience and control their flight direction during dynamic soaring, which results in optimal navigation toward their goal.

Our results also suggest the high cognitive ability of seabirds to solve an orientation task over the ocean. We could exclude the possibility, known as vector orientation, that the shearwaters continued heading in one particular direction, because turtles and the young birds of some species do (6, 33, 34). For vector orientation, only a sense of direction, that is, the ability to detect direction using a magnetic, sun, or stellar compass, is used. However, more complex orientation ability is required for the flow compensation the shearwaters use to return to their colony after foraging. This ability can often be termed a map sense: the ability to know where one is on the earth and the distance and direction to one’s destination without any landmarks (35). Because there are very few landmarks over the ocean, our finding that shearwaters compensate for the wind indicates that they have a map sense. Their map sense over the ocean might be based on a magnetic map, olfactory map, or wave patterns (20, 36).

Our analysis also revealed that the birds changed their orientation strategy (degree of compensation and preferred direction) during homing. Intriguingly, the degree of compensation was higher when the birds were flying over the sea near the coast (NC and WC sections) than over open sea (OS section), indicating that, in the former regions, they might be able to see the coast and obtain the corresponding navigational cues that were not available when flying over the open sea. In addition, the preferred direction changed with respect to the regions. The intended goal direction of birds was accurately toward the colony (final goal) at the beginning of homing flights (NC section). However, as they approached the colony (OS and WC sections), their preferred direction deviated from the colony, indicating that they tried to reach the coastline, which can be used as an additional navigational cue. Thus, although birds have a map sense and can compensate for wind over the open sea, they may try to use landmarks depending on their availability.

The inverse problem approach appears to be rarely used in behavioral ecology. However, as the volume and variety of biologging data has now entered the realm of big data (1), we expect the utilization of biologging data based on this approach to potentially reveal new causal influences on animal behavior. Our method is unique because it requires only tracking data, which is advantageous because animal position is the most widely used measurement parameter for biologging. We have shown the effectiveness of our model, for example, for use in testing a flexible orientation strategy of animals with respect to conditions en route; to date, such tests have been limited to a coarse time scale by the satellite data (37, 38) and/or to a short time scale by coastal radar studies (20, 39, 40). Although our method was applied to flying birds in this study, it should also be applicable to swimming animals subjected to currents. Accordingly, applications and improvements of our method can be used to reveal the orientation strategy in the flow of different species, and comparing these strategies among various species can also be expected to contribute to an improved understanding of the adaptation and evolution of their movement strategies in moving fluid.


The key idea of the model

First, we explain the key point of our model, which is based on the BRW (41). The BRW is a simple model for describing animal directed movement toward a particular goal direction or location. It describes the length (speed) and direction of a track vector with mutually independent probability densities, and the density of the direction has a peak in the goal direction. This independency between the direction and length of the track vector assumes a symmetric distribution of the track vector along its mean. What this assumption asserts is that, when you are walking in a particular direction, there should be some degree of fluctuation in each of the steps, but the fluctuation should not be different between the left side and right side relative to the mean track direction. This assumption seems intuitively natural. However, if you are walking on moving ground, such as a glacier, your track vector (movement relative to the “earth”) does not necessarily fluctuate symmetrically along its mean, but the velocity relative to the glacier should fluctuate symmetrically and thus obey the BRW. This is also the case with flying or swimming animals. “When flying or swimming animals are traveling in a particular direction, their velocity relative to the moving fluid (that is, heading vector) should obey BRW, but not the track vector.” This is the key point of our method and, as far as the authors know, is overlooked in random walk models in movement ecology. This key point eventually enables us to estimate the heading and wind vector only from track data, as explained in the following sections. The model is described in the next section. Then, we describe how to apply the model to real track data (see “The model fitting” section).

The model outline

Here, we describe the outline of our model (see note S1 for a more detailed explanation). When an animal moves while subjected to flow, its resultant track reflects the combined effects of the animal’s movement relative to the fluid and drift caused by the flow. We consider the condition in which the animal’s position in the x and y directions (Xk) is recorded at discrete time points k = 0, 1, …, N with the fixed time interval Δt. We define the animal’s track vector at time point k vT,k = (uT,k, vT,k) as the subtraction of consecutive positions divided by the time elapsed.vT,kXkXk1Δt(1)The track vector is the sum vector of the animal’s heading vector at time point k, vH,k = (uH,k, vH,k), and the flow vector at time point k, wk = (wk,x, wk,y)vT,k=(vH,k+wk)(2)

The aim of our method is to estimate the mean heading vector and flow vectors from the time series of the track vector {vT,1, vT,2, …, vT,N}, which can be easily calculated from the tracking data. For this purpose, we formulate vT,k as a stochastic variable generated from a particular probability density function characterized by the heading vector and flow vector.

We make the following assumptions about the heading vector and flow vector. First, heading vector vH,k is a stochastic variable whose probability density of the length [the animal speed relative to the fluid, airspeed for flying animals, and swim speed for swimming animals (hereafter referred to as the heading speed; sH,k = |vH,k|)] is not affected by that of the direction (heading; θH,k = arg vH,k). Second, the animal movement relative to flow is directional. That is, the probability density of the heading has a single peak at the particular mean heading φ. As probability densities that satisfy these two assumptions, we used the Weibull distribution f(sH,k|γ, ρ) for the heading speed (eq. S1), where γ is the mean heading speed and ρ is the so-called shape parameter, and the von Mises distribution gH,k|φ, κ) for the heading (eq. S2), where φ is the mean heading and κ is the concentration parameter that characterizes the variance of the distribution. From these two distributions, the probability density of vH,k = (uH,k, vH,k), represented by pH(uH,k, vH,k|γ, ρ, φ, κ), can be derived (eq. S3). Note that it is symmetrically distributed along its mean heading vector (Fig. 2A). This model is the BRW, which has been widely used for modeling the “track vector” (23). However, as we explained in the previous section, we found that the “heading vector” should obey a BRW for the animal movement in fluid. The final assumption is that flow vector wk = (wx,k, wy,k) can be approximated by a spatiotemporally constant vector w = (wx, wy) during the unit of model fitting. Under the assumption of this constant flow vector, we added the flow-induced drift effect. Using the relation vT,k = (vH,k + w), the probability density of the track vector pT(uT,k, vT,k|γ, ρ, φ, κ, wx, wy) can be derived from that of the heading vector pH(uH,k, vH,k|γ, ρ, φ, κ) by transforming its variables from (uH,k, vH,k) to (uT,k, vT,k) (eq. S4). pT(uT,k, vT,k|γ, ρ, φ, κ, wx, wy) corresponded to pH(uH,k, vH,k|γ, ρ, φ, κ) translated to w (Fig. 2, A and B). The important prediction of our model is that the track vectors should be distributed asymmetrically along the mean track vector (Fig. 2B). Note that the flexibility of random walk models (Property 2 mentioned in Introduction) enables the formulation of the animal movement in moving fluid, and the random fluctuation of animal track described by random walk models (Property 1) is crucial in this model.

The model fitting

Parameter estimation. When the track vector data {vT,1, vT,2, …, vT,N} are given, the likelihood of the model can be calculated as followsL(γ,ρ,φ,κ,wx,wy)=k=1NpT(uT,k,vT,k|γ,ρ,φ,κ,wx,wy)(3)

Using the maximum likelihood estimation (MLE), the parameters γ, ρ, φ, κ, wx and wy can be estimated (Property 3). Accordingly, the mean heading vector and flow vector can be estimated by fitting our model to the track vector data. By conducting this estimation for the simulated movement data, we tested the effect of a small sample size, the fluctuation of wind, and the measurement error on the accuracy of parameter estimation (we conducted MLE using “optim()” of R version 3.2.0). We found that, in some respects, these factors could be problematic for estimations; thus, we specified five conditions to exclude these problems (see notes S2 and S3 for detailed information).

Condition 1: We only accepted a result for which the estimated probability density of the heading vector is anisotropic and more widely distributed along the perpendicular direction relative to the axis of the mean heading vector than along the axis of this vector (specific representation is provided in eq. S7).

Condition 2: We only accepted results in which the angle between the mean heading and the mean track vector is less than 90°.

Condition 3: Although the heading was estimated accurately, the mean heading speed widely fluctuated and deviated from the value we obtained as an answer with only 50 data points. Accordingly, we assigned the mean heading speed γ a priori. Here, we used a value of 9.63 m s−1 (34.7 km hour−1), which is the mean speed relative to the ground for streaked shearwaters reported in a previous work (42).

Condition 4: The SD of the distance between the animal’s true fixed position and the observed fixed position is less than 5% of the distance between successive observed fixed positions.

Condition 5: The flow fluctuation is smaller than that of the animal heading vector.

Thus, we assumed that conditions 4 and 5 are satisfied and conducted the MLE under condition 3, fixed the mean heading speed value (γ) to 9.63 m s−1, and, if the estimated parameters satisfied conditions 1 and 2, checked the goodness-of-model fit as follows.

Goodness-of-fit tests. The goodness of fit of the estimated distribution of the data was verified using statistical tests. First, we checked whether the data related to the heading and heading speed were distributed according to the estimated distributions by using a Kolmogorov-Smirnov test. Second, we tested whether the heading and heading speed were uncorrelated by using Pearson’s correlation test. If the estimated distribution passed these tests, we accepted that the model fits the data.

Field experiment

Streaked shearwaters (C. leucomelas) breeding at Funakoshi-Ohshima Island, Japan (39.240°N, 141.590°E) (Fig. 3) were tracked from August to September in 2013 and 2014, which corresponds with the chick-rearing period. A GPS data logger was attached to each bird’s back with waterproof tape (Tesa) and glue. Two models of GPS loggers, GiPSy-2 (2013, N = 11 birds; 2014, N = 14 birds) and GiPSy-4 (2013, N = 9 birds; 2014, N = 13 birds), were used (Technosmart). GiPSy-2 was configured to record the position of the bird every minute, and GiPSy-4 every 30 s (but resampled to 1 min in the analysis to homogenize the sampling interval between GiPSy-2 and GiPSy-4). We began to recapture the birds at the colony to recover the loggers approximately 1 week after deployment. The overall mass of the logger was approximately 25 g for both models, which corresponded to less than 5% of the body mass of each bird.

Fitting the model to real track data

The foraging trips recorded by the GPS loggers were analyzed as follows. First, we identified their homing start position using a method developed in a previous study (42). The method defines the last phase of the track during which the distance from the colony continues to decrease as the homing phase, and in that phase, the first point where the approach speed to the colony exceeds 15 km hour−1 is defined as the homing start point. Then, to eliminate the possibility that the birds used land features to navigate, we included only bird tracks in which the homing start position was located more than 200 km away from the colony. In addition, only birds that started homing further east than the longitude of the colony were analyzed (Fig. 3).

Our method was applied to the homing portions of the tracks. Because the animal movement pattern and flow vector are expected to change spatiotemporally, we conducted the model fitting as follows. Fifty-one minutes of position data for each bird were selected every consecutive minute in a sliding time window. Using the position data in each time window, the mean bird heading and wind (flow) vector were estimated every consecutive minute at time point t, averaged between t − 25 min and t + 25 min. The bird position data were not analyzed when the birds were near land, that is, when the location of the bird at the center of the time window was closer to land than the blue dashed line, as shown in Fig. 3. For each time window, the track vector was derived by calculating the distance between successive GPS positions on the x and y axes and then dividing it by the elapsed time. Because there were some missing GPS data points, the track vector was excluded from the analysis when the elapsed time between points at which the GPS coordinates were fixed exceeded 1 min. When the track vector was less than 15 km hour−1, the bird was assumed to be resting on the sea surface, that is, not flying (42). We ensured that the bird heading vector and wind vector were accurately estimated by fitting the model only when at least 45 overtrack vector data points were contained in each window. Then, parameter estimation and goodness-of-fit tests for the probability density function of the bird track vector were calculated.

Comparing the wind vector estimated from bird tracks with the wind vector simulated by the atmospheric model

We evaluated the validity of the model estimation by correlating the estimated wind vectors with wind data from the Japan Meteorological Agency mesoscale model reanalysis data sets (at a height of 10 m above sea level), which was calculated every 3 hours. When the reanalysis data at time t were available, the mean vector of the reanalysis data was calculated and compared with the estimated wind vector at approximately time t. First, we selected the estimated wind vector with a time difference from t of less than 25 min. The mean vector of the selected wind vector was calculated. This vector was defined as the mean estimated wind vector at time t. Next, using the wind vectors from the reanalysis data at time t, we selected the wind vector data within 5 km of the recorded positions of the bird from t − 50 min to t + 50 min, except for the bird positions not contained in any time window that fit the model. This vector was defined as the mean reanalysis data wind vector at time t. The correlation between the mean estimated wind vector and the mean reanalysis data wind vector was examined using the generalized vector correlation coefficient (43, 44) that takes into account both wind speed and direction and takes a value between 0 (no correlation) and 2 (complete correlation). The correlation of the wind speed was examined using Pearson’s correlation coefficient, and the correlation of the wind direction was examined using the Jammalamadaka-Sarma correlation coefficient [calculated using the package “circular” of R (45)].

Testing the bird response to the wind over the sea

Using the track vector direction and the estimated bird heading direction, we calculated the degree of compensation and preferred direction by the method proposed in a previous study (29). First, we separated the tracks into three sections: NC, OS, and WC (Fig. 6). For each section and each track, we averaged the track vector direction and the angular difference between the track and heading vectors (α) (Fig. 6A). Then, we linearly regressed the average track direction on average α. The degree of compensation and preferred direction correspond to the slope and intercept, respectively, of the regression line (29). A slope of 1 signifies full drift (Fig. 1A), a slope of 0 means complete compensation (Fig. 1C), and the case 0 < slope < 1 means partial compensation (Fig. 1B).

In addition, to investigate the extent to which the preferred direction differs from that of the bird’s nesting colony (that is, the final goal), we defined the locations with median longitude and latitude of the tracks in each section as the median locations and calculated the direction of the colony from there.


Supplementary material for this article is available at

note S1. Detailed information regarding the model.

note S2. Numerical simulation test for checking the effect of realistic sample size on model fitting.

note S3. Numerical simulation test to assess the effect of variation of flow and track location quality on model fitting.

fig. S1. Relation between wind speed and difference between estimated and reanalyzed wind direction.

fig. S2. Probability distributions of heading vector.

fig. S3. Combinations of mean heading vector, mean track vector, and flow vector used in the simulation.

fig. S4. Histogram of difference between the value of the estimated and true heading.

fig. S5. Histograms of results that comply with conditions 1 and 2.

fig. S6. Dependence of model fitted ratio and estimation accuracy on flow fluctuation and location observation error.

Sample_Track_Simulation.doc: R code that outputs simulated track data (Sample_Track_Simulation.R)

Heading_Wind_Estimation.doc: R code for estimating the heading direction and wind from the track data (Heading_Wind_Estimation.R)

Real_Track_Data.csv: Homing track data of streaked shearwaters (N = 33)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank I. K. Shimatani for valuable suggestions and discussion that led to the conception of the idea on which the model is based. We also thank Y. Yonehara, T. Shiozaki, T. Abe, and T. Abe for logistic support in the field; P. Miller, A. Wada, K. Shiomi, K. Fukaya, A. Takahashi, Y. Y. Watanabe, and Y. Naito for discussion and for providing valuable advice; and F. Daunt, M. Müller, and A. Robson for their assistance in improving our manuscript. We furthermore thank four reviewers and the editors for their valuable comments. Funding: This study was financially supported by the Tohoku Ecosystem-Associated Marine Science (TEAMS), Bio-Logging Science of the University of Tokyo (UTBLS), Grants-in-Aid for Scientific Research from the Japan Society for the Promotion of Science (15J10905, 24241001, 24681006, and 16H06541), National Geographic (Asia 45-16), and Japan Science and Technology Agency Core Research for Evolutional Science and Technology (grant number JPMJCR1685). Author contributions: Y.G., K.Y., and K.S. designed the study and collected and contributed data. Y.G. analyzed the data. All authors participated in the writing of the manuscript (Y.G. wrote the first draft). Competing interests: The authors declare that they have no competing interests. Data and materials availability: Program codes for simulating artificial track data and for estimating the heading and the wind vector as well as homing tracks of shearwaters recorded by the GPS are present in the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article