Classification of anomalous diffusion in animal movement data using power spectral analysis

Black winged kite. An individual black-winged kite (Elanus caeruleus), residing in the Hula Valley, Israel, was tracked using ATLAS, an innovative reverse-GPS system. ATLAS localises extremely light-weight, low-cost tags [42], where each tag transmits a distinct radio signal which is detected by a network of base-stations distributed in the study area. Tag localisation is computed using nanosecond-scale differences in signal time-of-arrival to each station, alleviating the need to retrieve tags or have power-consuming remote-download capabilities [32, 42]. The kite was tracked for 164 consecutive days in years 2019–2020, with a mostly constant tracking frequency of 0.25 Hz, see reference [32] for more details. Our analysis is limited to data collected during the activity hours (omitting the nights for the diurnal kite) and we further excluded data collected in proximity to the observed nest to focus on local search behaviour.

White stork. An adult white stork (Ciconia ciconia) was tracked between May 2012 and July 2020 using the configuration described in reference [43]. Here the GPS location and speed were recorded during the day at a frequency of 1/300 Hz when solar recharge was high (92% of the time), and at 1/1200 Hz otherwise. We omit days with lower frequency ($< 1\%$ of tracked days) and only include localisations that occur between the first recorded velocity of $> 4$ m s⁻¹ [44].

In this work we focus on empirical evidence for 1/f noise and ageing of the PSD in ecological movement data. In several cases it is possible to relate our findings to known physical models. Here, we list those models, for which analytical results for the PSD have been reported in recent years.

2.2.1. Brownian motion

For an overdamped, one-dimensional Brownian trajectory x(t), with t ∈ [0, t_m], the dynamics is given by the Langevin equation [45]

$$\begin{equation}\frac{\mathrm{d}x(t)}{\mathrm{d}t}=\xi (t),\end{equation} \tag{ 4 }$$

where ξ(t) is a Gaussian, delta-correlated white noise term with zero mean and $\left\langle \xi (t)\xi ({t}^{\prime })\right\rangle =2D\delta (t-{t}^{\prime })$, and D is the diffusion constant. For a particle governed by equation (4), the EA-TA- (equation (3) and EA-MSD are respectively given by $\left\langle \bar{{\delta }^{2}(\tau ,t)}\right\rangle =2D\tau$ and $\left\langle {x}^{2}(t)\right\rangle =2Dt$ [6, 37]. As both show identical scaling with time, the process is ergodic in the (weaker) Boltzmann–Khinchin sense of equality between (long) time and ensemble averages [35]. For this process, the mean of the PSD follows [41]

$$\begin{equation}\left\langle S(f)\right\rangle \simeq \frac{4D}{{f}^{2}},\end{equation} \tag{ 5 }$$

entailing β = 2 and that the amplitude of the PSD is linear in D without explicit dependence on t_m. Here, the non-integrability of the noise is generally cured by considering an unbounded process however it does demonstrate the significant difference between 'white' noise and processes where the noise is Brownian (sometimes called 'red' or 'brown' noise [8, 25]).

For any finite t_m the single-trajectory PSD will fluctuate from one realisation to the next. The distribution of the single-trajectory PSD (equation (1)) around its ensemble average is then quantified by the amplitude [8]

$$\begin{equation}\frac{S(f,{t}_{\text{m}})}{\langle S(f,{t}_{\text{m}})\rangle }\equiv A(f,{t}_{\text{m}}),\end{equation} \tag{ 6 }$$

where A is a random number that depends on the choice of the stochastic model (see examples below), and can also depend on the dimension d, the frequency f, and the measurement time t_m. The distribution of A is highly informative since it is model specific as shown analytically for a number of processes [8, 9, 46, 47]. In particular, it is important in studies where only few paths are available, as it relates the single PSD to its average. For BM, the distribution of the single-trajectory PSD around the mean (5) is given by equation (6), with [8]

$$\begin{equation}P(A)=\frac{2\sqrt{\pi }{A}^{(d-1)/2}}{\sqrt{3}{\Gamma}(d/2)}\;{\text{e}}^{-\frac{4}{3}A}{I}_{(d-1)/2}\left(\frac{2}{3}A\right),\end{equation} \tag{ 7 }$$

where I_ξ (⋅) is the modified Bessel function of the first kind and Γ(⋅) is the gamma function. Notably, the shape of the distribution depends on the dimensionality, where dimensions lower than the embedding space of the process can be achieved by projection or component-wise measurement [8, 9].

2.2.2. Fractional Brownian motion

As in the case of BM, FBM can also be defined in terms of the Langevin equation (4), replacing the delta-correlated noise term ξ(t) with zero-mean, power-law correlated fractional Gaussian noise ξ_fGn(t) defined by [35]

$$\begin{equation}\left\langle {\xi }_{\text{fGn}}(t){\xi }_{\text{fGn}}({t}^{\prime })\right\rangle \sim 2H(2H-1)\vert t-{t}^{\prime }{\vert }^{2(H-1)}\end{equation} \tag{ 8 }$$

for t ≠ t'. Here, H is termed the Hurst exponent and, similarly to the case of BM, the process is ergodic with $\left\langle \bar{{\delta }^{2}(\tau ,t)}\right\rangle =2\tilde{D}{\tau }^{2H}$ and $\left\langle {x}^{2}(t)\right\rangle =2\tilde{D}{t}^{2H}$, where $\tilde{D}$ is an effective diffusion constant [35, 48, 49]. In accordance with the definition of ξ_fGn(t), for H > 0.5 the noise is persistent leading to a positively correlated process, while for H < 0.5 the noise is antipersistent leading to a negatively correlated process. For FBM and other stationary (in increments) processes, the PSD is written in terms of the EA autocorrelation function ${C}_{\text{EA}}(\tau )=\left\langle x(t)x(t+\tau )\right\rangle$ following the Wiener–Khinchin theorem

$$\begin{equation}\left\langle S(f,\infty )\right\rangle ={\int }_{-\infty }^{\infty }{\text{e}}^{\text{i}f\tau }{C}_{\text{EA}}(\tau )\mathrm{d}\tau .\end{equation} \tag{ 9 }$$

Here, C_EA and the time averaged (TA) autocorrelation function, defined by

$$\begin{equation}{C}_{\text{TA}}({t}_{\text{m}},\tau )=\frac{1}{{t}_{\text{m}}-\tau }{\int }_{0}^{{t}_{\text{m}}-\tau }x(t)x(t+\tau )\mathrm{d}t,\end{equation} \tag{ 10 }$$

are identical at long times. In the limit of long measurement times it was shown that for FBM [9]

$$\begin{equation}\left\langle S(f,{t}_{\text{m}})\right\rangle \sim \begin{cases}{t}_{\text{m}}^{2H-1}{f}^{-2}\quad \hfill & H > 1/2,\hfill \\ {f}^{-1-2H}\quad \hfill & H< 1/2,\hfill \end{cases}\end{equation} \tag{ 11 }$$

i.e., for persistent FBM the PSD is ageing, while it is not for subdiffusive FBM. Here, it is also possible to obtain the distribution of the single-trajectory PSD around the mean; for example, for superdiffusive FBM (H > 1/2) in d-dimensions it was found that [9]

$$\begin{equation}P(A)=\frac{{A}^{d/2-1}}{{2}^{d/2}{\Gamma}(d/2)}\;{\text{e}}^{-A/2},\end{equation} \tag{ 12 }$$

see reference [9] for more details and results. Here, the distribution P(A) is in general different from that of BM, compare equation (7).

2.2.3. Subordinated FBM

Here we consider a process in discrete time steps, expressed by the number of jumps n = 1, 2, 3, ..., such that the autocorrelation function is given by

$$\begin{equation}\langle {x}_{n}{x}_{n+{\Delta}n}\rangle ={\Delta}{x}^{2}[{n}^{2\mathcal{H}}+{(n+{\Delta}n)}^{2\mathcal{H}}-{n}^{2\mathcal{H}}],\end{equation} \tag{ 13 }$$

where Δx is a scaling parameter. Furthermore, $\mathcal{H}$ is analogous to the Hurst exponent of the FBM, although here it is not directly related to the scaling of the EA-MSD on time (see below). In this process, the discrete-time FBM (13) is subordinated to the operational time of the CTRW, such that after each FBM step, the particle is immobilised for a random waiting time τ drawn from a fat-tailed distribution [50]

$$\begin{equation}\psi (\tau )\sim {\tau }^{-(1+\alpha )},\end{equation} \tag{ 14 }$$

with 0 < α < 1, such that ⟨τ⟩ diverges. By combining the correlation structure of the jump sizes of FBM with the unbounded waiting times of CTRW, this process effectively distinguishes the 'operative clock', defined by the sequence of jump events themselves, from the 'real' clock defined by t, thus allowing for both temporal correlations as well as ageing effects (note that pure FBM only allows for the former, while CTRW allows for the latter ⁸ ). Waiting time distributions similar to equation (14) have been reported in multiple empirical systems, e.g. [32, 51–54].

As subordinated FBM is a non-stationary process, C_EA and C_TA are not identical, and the Wiener–Khinchin theorem no longer applies. Similar observations in other processes have led to the development of the ageing Wiener–Khinchin theorem in recent years [16, 17], which enabled analytical calculations of the PSD of both stationary and non-stationary processes, relating the exponents β and z to underlying physical processes, see, e.g. [9, 17, 24, 46, 47, 55]. In particular, for subordinated FBM the PSD has been rigorously computed in reference [24]. The leading term of the spectrum was found to again differ between positively and negatively correlated increments. It satisfies [24]

$$\begin{equation}\left\langle S(f,{t}_{\text{m}})\right\rangle \sim \begin{cases}{t}_{\text{m}}^{2\alpha \mathcal{H}-1}{f}^{-2}\quad \hfill & \mathcal{H} > 1/2,\hfill \\ {t}_{\text{m}}^{\alpha -1}{f}^{-2+\alpha (1-2\mathcal{H})}\quad \hfill & \mathcal{H}< 1/2.\hfill \end{cases}\end{equation} \tag{ 15 }$$

In equation (15), for positively correlated increments $\mathcal{H} > 1/2$, the process does not exhibit 1/f noise (i.e., β = 2, as is the case for BM, see the discussion in [9]); however, the ageing properties depend on α: for $\alpha \mathcal{H} > 1/2$ the spectrum increases with the measurement time t_m as in superdiffusive FBM [9], while for $\alpha \mathcal{H}< 1/2$ the spectrum amplitude decreases with t_m. In contrast, for anticorrelated increments in terms of the number of jumps n $(\mathcal{H}< 1/2)$, for any α < 1 the spectrum displays 1/f noise (β < 2) and ageing, effecting a decrease with the measurement time (z < 0). Below we compare these theoretical predictions to the recorded movement paths of a kite and a stork, for which we can obtain the values of α and $\mathcal{H}$ and infer the non-stationary properties of different movement modes.

Importantly, the values of α and $\mathcal{H}$ in this model can also be obtained from independent analysis of the EA-TA-MSD, equation (3). For the subordinated process, when 0 < α < 1 the EA-TA-MSD is given by [56]

$$\begin{equation}\left\langle \bar{{\delta }^{2}(\tau ,{t}_{\text{m}})}\right\rangle \sim \frac{{\tau }^{1-\alpha +2\alpha \mathcal{H}}}{{t}_{\text{m}}^{1-\alpha }}.\end{equation} \tag{ 16 }$$

The EA-TA-MSD was found for our empirical data (kite and stork) in a recent study [44]; however, the MSD (and in particular the EA-MSD) is typically a noisy observable [24], and in some cases it was not possible to perform adequate power-law fits. Thus, the analysis of the PSD can be highly valuable when comparing the data to theoretical processes, as shown below. Finally, as the theoretical distribution of the single-trajectory PSD is yet unknown for subordinated FBM, we compare our results to stochastic simulations of the process. Here, the increments of the FBM, in discrete time n, are obtained using the python function fgn from the package fbm (here fGn refers to fractional Gaussian noise—the increment process of the FBM). The increments, δr_n , are projected onto two dimensions by generating a random number θ ∈ [−π, π] and setting δx_n = δr_n  cos θ and δy_n = δr_n  sin θ. FBM in two dimensions is given by the cumulative sums ${x}_{n}={\sum }_{m=0}^{n}\delta {x}_{m}$ and ${y}_{n}={\sum }_{m=0}^{n}\delta {y}_{m}$. The times between the (discrete-time) steps are then drawn from a power-law distribution (14) using the python package powerlaw [57]. Notably, there are also other ways to simulate FBM for d = 2 and in particular, to sample from a power-law distribution ⁹ ; yet, comparing between different methods is beyond the scope of this paper.

To compare data and theory, we fitted the averaged PSD of empirical and simulated trajectories to power laws. The PSD for each trajectory was computed in a straightforward manner from definition (1), and the PSDs were averaged over the relevant ensemble. To obtain the averaged PSD as a function of the measurement time, the averaged PSD was computed for different values of t_m, defined as the time passed from the first point of the trajectory. The resulting curves, both as a function of f and t_m, were then fitted to a straight line on a log-log curve, using SciPy library's curve-fit (nonlinear least squares method) in python 3.8. In all cases the fit was performed for small frequencies f ≪ 1, as to ensure the validity of the theory. Errors were computed as the standard deviation when fitting the averaged PSD for different times (when fitted versus f) or different frequencies (when fitted versus t_m). Direct errors of each fit (obtained using curve-fit) were, in all cases, negligible compared to the standard deviations described above.

Similarly to reference [32], the kite's tracks were segmented into two behavioural modes: local searches (area restricted search) and commutes (directed flights between local searches). Localisations were segmented by detecting switching points in the data—distinct points in which the bird switches between the two behaviours [58]. Switching points were detected using spatiotemporal criteria segmentation, such that localisations that are in proximity to one another both in space and time, were segmented together [27]. In accordance with the conclusions of reference [32] we independently analysed the time series ensemble that represents instances of searches and the time series ensemble that represents commutes. Below, we compare these analyses to an unsegmented ensemble of daily tracks, showing that ecological knowledge of the underlying processes is crucial to correctly identify the noise properties of the data.

For searches, the PSD (figure 1(b)) displays 1/f^β noise with $\left\langle S\right\rangle \sim {f}^{-1.50}{t}_{\text{m}}^{-0.57}$. These results are consistent with subdiffusive subordinated FBM with α = 0.43 ± 0.04 and $\mathcal{H}\simeq 0$, see equation (15). Notably, the value of $\mathcal{H}$ is close to zero and the error is relatively large, thus an exact estimate of its value is hard to achieve. An analysis of the MSD for the same ensemble of searches produces $\left\langle \bar{{\delta }^{2}(\tau ,t)}\right\rangle \sim {\tau }^{0.55}{t}^{-0.55}$, which is consistent with α = 0.45 ± 0.03 and $\mathcal{H}\simeq 0$, see equation (16). Here, both the analysis of the MSD and of the PSD suggest that subordinated FBM is a plausible model for the movement within search grounds of this kite. In contrast, for the commutes (figure 1(c)) we find that $\left\langle S\right\rangle \sim {f}^{-2}{t}_{\text{m}}^{0.60}$, consistent with (ergodic) superdiffusive FBM with H = 0.80 ± 0.03, see equation (11). Here, analysis of the EA-TA-MSD gives a scaling of $\left\langle \bar{{\delta }^{2}(\tau ,t)}\right\rangle \sim {\tau }^{1.66}$, suggesting H = 0.83 ± 0.03, in agreement with the analysis of the PSD.

Zoom In Zoom Out Reset image size

We compare these results to the PSD of unsegmented daily trajectories (figure 1(a)). Here, the trajectories include both the commutes and searches, each occurring at different temporal and spatial scales. We find that β = 2.10, however there is strong ageing in the process as $\left\langle S\right\rangle \sim {t}_{\text{m}}^{-0.66}$. Notably, the 1/f noise found during searches is completely skewed due to the lack of segmentation, and the process is most likely a mixture of (at least) two different processes, making these results hard to interpret. Additionally, the dependence of the PSD on t_m (figure 1(a), inset) does not provide a clear scaling form when normalised by f^−β for different f as is the case for searches and commutes (figures 1(b) and (c), insets). These findings highlight the importance of correctly accounting for different movement phases, allowing us to better interpret the 1/f noise in the data.

Some discussion on stationarity is in order here. For searches, the longer we observe the system, the smaller the amplitude A of the PSD becomes (figure 1(b), inset). This occurs because the longer we track the searching bird, the more likely we are to find it 'trapped' at a specific location for long periods; thus, the rate at which the animal moves is significantly reduced and the noise levels decrease accordingly [12]. In a recent work we have shown that this effect originates from long waiting times during searches [32]. The scaling exponent which is proportional to α (see equation (15)) allows for yet another way to quantify this effect. In contrast, for commutes the PSD increases with the measurement time (figure 1(c), inset). This process can be modelled with the use of FBM (i.e., α > 1) and is Gaussian and ergodic in nature. As these flights are highly correlated, persistent movement at increasingly long measurement times increases the amplitude of the PSD—hence the growth of the PSD.

Finally, in figure 2 we plot the distribution of A given by equation (6) for searches and commutes in d = 1 and d = 2 dimensions. For searches (figure 2(a)), we compare these distributions to simulation results for subordinated FBM with α = 0.45 and $\mathcal{H}=0.01$, which are consistent with the results obtained above. For both dimensions we find excellent agreement between the simulated and observed distributions. Interestingly, we find that for searches, the empirical distribution does not strongly differ between d = 1 and d = 2. Here, for subordinated FBM a theoretical distribution for A is yet unknown. In comparison, for the commutes (figure 2(b)) the distributions are in good agreement with the theory given by equation (12) for superdiffusive FBM.

Zoom In Zoom Out Reset image size

The relatively low tracking frequency of the stork (one recorded point every five minutes) does not allow to identify and independently analyse potential searches, as these will include a very limited number of points (e.g., a search of 1 h will include $\sim 12$ points). Thus, daily tracks are not segmented but rather considered as a single process, limiting the scope of the analysis. Nevertheless, to properly account for the PSD for different behaviours, daily paths are clustered into subsets that represent distinct behaviours in a bird's life cycle [43]. This is done using a Gaussian mixture model (GMM) [59] based on the logarithm of the maximum displacement (ΔX) that the stork performs during each day (figure 3(a)). We have found three significant movement modes which are highly correlated to the time of year: the mode with $\left\langle {\Delta}X\right\rangle =3$ km occurs primarily $( > 80\%)$ between April and July, the mode with $\left\langle {\Delta}X\right\rangle =33$ km occurs primarily $( > 80\%)$ between September and February, and the mode with $\langle {\Delta}X\rangle =211$ km occurs primarily $( > 87\%)$ during well-known migratory periods [43]. We thus refer to these three modes as breeding, wintering and migrating respectively, and perform the analysis below for each subset of daily paths separately. Note that the GMM identified another mode with $\left\langle {\Delta}X\right\rangle =0.4$ km. However, as this mode includes a small number of daily paths (60 days out of $\sim 2500$), we discard it from our analysis.

Zoom In Zoom Out Reset image size

For the breeding ensemble (966 days, figure 3(b)) we found $\left\langle S\right\rangle \sim {t}_{\text{m}}^{-0.45}{f}^{-1.9}$, which is consistent with subordinated FBM with α = 0.55 ± 0.09 and $\mathcal{H}=0.40\pm 0.05$, i.e., a non-stationary process with near 1/f² noise. For the wintering mode (1068 days, figure 3(c)) $\left\langle S\right\rangle \sim {t}_{\text{m}}^{0.07}{f}^{-2.10}$, i.e., the noise is approximately 1/f² noise and we find no significant ageing patterns. Importantly, the fact that we measure exponents β = 2 and z = 0 in equation (2) does not entail that the process is Brownian. On the contrary, subordinated FBM with any α and $\mathcal{H}$ that fulfil $2\mathcal{H}\alpha =1.07$ is a plausible model, as it is consistent with $\left\langle S\right\rangle \sim {t}_{\text{m}}^{0.07}{f}^{-2.10}$ and equation (15). Here, the knowledge of the analytical scaling of the PSD (equation (15)) is not sufficient to determine the values of both α and $\mathcal{H}$ independently, and complementary analysis is crucial. For instance, based on the conclusions of reference [44] one can determine $\mathcal{H}\simeq 0.7$ and α ≃ 0.77, which is consistent with the above relation. Below, we show that the distribution of the single-trajectory PSD further supports this. For the migrating mode (480 days, figure 3(d)) we get $\left\langle S\right\rangle \sim {t}_{\text{m}}^{1.27}{f}^{-2}$, which is consistent with the ballistic motion of a migrating stork. Here, the stork migrates approximately in a straight line towards a stationary target (breeding and wintering grounds for spring and fall migrations, respectively), accounting for the rapid increase in amplitude (figure 3(d), inset). Notably, the fact that the scaling with the measurement time is steeper than linear implies that, in addition to the highly correlated nature of this process, the animal is also accelerating, which is again consistent with the complementary analysis in reference [44]. ¹⁰ Indeed, storks roost in stopover sites during the night and tend to depart in the late morning, when soaring conditions improve, facilitating faster flights at lower energy costs [60].

In figure 4 we plot the distribution of the single-trajetory PSD around its average for breeding, wintering and migration in one and two dimensions. For breeding we compare these distributions to simulations of subordinated FBM with α = 0.55 and $\mathcal{H}=0.40$ (see above). Similarly, for wintering we compare to subordinated FBM with α = 0.77 and $\mathcal{H}=0.70$. Here, the predicted distributions (7) for BM for d = 1, 2, can be shown to be far from the empirical distributions, entailing that this process is not Brownian. For both the breeding and wintering modes the empirical distributions show good agreement with simulated subordinated FBM. Note that although the ageing properties differ between these modes (see insets of figures 3(b) and (c)), the distributions of the single-trajectory PSD (figures 4(a) and (b)) are similar, suggesting some kind of unifying features of the underlying dynamics of the stork during breeding and wintering, see reference [44] for further discussion. Finally, for migration (figure 4(c)) we find that the distribution is in excellent agreement with the theoretical distribution for superdiffusive FBM given by equation (12).

Zoom In Zoom Out Reset image size

In many applications, the initial measurement time (the time from which the process is first observed) is different from the initiation time of the process, where the latter is often unknown. For nonergodic processes, such a discrepancy can make an anomalous process appear closer to normal diffusion, as the observed statistics of the initial unmeasured period can be distinct from the statistics of the measured period [18, 53, 61–63]. We denote the time between the initiation of the process and the beginning of the measurement by t_a, and refer to it as ageing time [35, 63]. Here, we measure the exponents β and z in equation (2) for the kite during searches, as a function of t_a, see figure 5. We vary t_a between 1 and 512 by discarding all points occurring at t < t_a for each t_a and obtain the PSD from the remaining points in the range [t_a, T] (such that t_m = T − t_a), T being the time since the unknown initiation of the process. Notably, searches typically last 40 min (=2400 s), which is significantly longer than any choice of t_a.

Zoom In Zoom Out Reset image size

For the PSD of the searching kite, we find that variations in β are relatively small, with $< 10\%$ difference between its value at t_a = 0 and its value at t_a = 512. This entails that ageing the trajectory (starting the measurement at increasingly long lag time t_a instead of the initiation time of the process) does not affect the power-law shape and 1/f noise. In contrast, the exponent z significantly changes with t_a $( > 70\%)$, meaning that ageing the trajectory changes the non-stationary properties of the system. For instance, if the process is modelled with subordinated FBM (section 3.1), analysis of the PSD can give different values of α, depending on t_a. This effect has important implications for ecologists, as it is often the case that the initiation time of a motion and the initial measurement time are not the same. Importantly, similar ageing effects can appear in the analysis of the EA-TA-MSD, as can be seen by comparing equations (15) and (16). In figure 5 we show a remarkable agreement between the two methods (PSD and MSD) with respect to the ageing analysis of the empirical tracks.