Methods for computing movement parameters in the Brownian bridge movement model
We first discuss how various movement parameters are calculated in the BBMM and similar models. We then provide details on the specific methods used in the two case studies. The BBMM assumes that an entity exhibits Brownian motion between measured locations. A Brownianbridge is the distribution of this process conditioned on the locations of both endpoints. To model uncertainty in the measured locations and to avoid a degenerate probability distribution at the time of a measurement, the locations are often assumed to be normally distributed around the measured location. All of the following calculations are performed for individual Brownian bridges and only use the directly adjacent measurements. Note that in the presence of measurement errors the sequence of observations does not satisfy the Markov property [14], and any Brownian bridge actually depends on more than just the adjacent measurements. Thus, we need to assume that the measurement error is small relative to the diffusion coefficient.
If we assume that we have two locations x
i
, x
i+1 measured at times t
i
, t
i+1 with variances \({\delta _{i}^{2}}\) and \(\delta _{i+1}^{2}\) respectively, the position X
t
at a time t∈[t
i
,t
i+1] follows a circular bivariate normal distribution with parameters
$$\begin{array}{@{}rcl@{}} \boldsymbol{\mu}(t) &=& (1-\alpha) \boldsymbol{x}_{i} + \alpha \boldsymbol{x}_{i+1}, \\ \sigma^{2}(t) &=& (t_{i+1}-t_{i}) \alpha(1-\alpha) D + (1-\alpha)^{2} {\delta_{i}^{2}} + \alpha^{2} \delta_{i+1}^{2}. \end{array} $$
Here, \(\alpha = \frac {t-t_{i}}{t_{i+1}-t_{i}}\) is a variable that linearly moves from 0 to 1 as t moves from t
i
to t
i+1 and D is the diffusion coefficient of the Brownian motion, which is often estimated by a maximum likelihood method [9]. When the trajectory contains different movement states over time, it may be appropriate to vary the diffusion over time rather than to keep it constant [10].
Given these probability distributions, derived parameters such as distance or speed (relative to a time scale) can be determined [16]. These parameters are important building blocks for the detection of many movement patterns. We summarize the results on the distributions of these parameters here, for full derivations we refer to [16] and online Additional file 4. Note that the derivation of velocity in [16] does not handle all possible dependencies and is superseded by the derivation in Appendix 1.
If the positions of two animals A and B at time t have independent circular normal distributions with means μ
A
(t) and μ
B
(t) and variances \({\sigma _{A}^{2}}(t)\) and \({\sigma _{B}^{2}}(t)\) respectively, the distance between A and B has a Rice distribution with parameters |μ
A
(t)−μ
B
(t)| and \(\sqrt {{\sigma _{A}^{2}}(t) + {\sigma _{B}^{2}}(t)}\). The average velocity over a time interval [t
1,t
2] is given by the difference between two (generally not independent) circular normal distributions, for X(t
2) and X(t
1). The velocity has a circular normal distribution with mean \(\frac {\boldsymbol {\mu }(t_{2})-\boldsymbol {\mu }(t_{1})}{t_{2}-t_{1}}\), while the expression for the variance depends on the number of location measurements that were obtained between t
1 and t
2.
Let t
s
, t
i
and t
f
be the time stamps of three consecutive observations with location variances \({\delta _{s}^{2}}\), \({\delta _{i}^{2}}\) and \({\delta _{f}^{2}}\) respectively, chosen such that t
s
≤t
1<t
i
. The observation at t
f
is only needed in the calculations if t
i
<t
2≤t
f
. The variance of the velocity is:
$$ {\sigma_{V}^{2}} (t_{1}, t_{2}) =\left\{\!\! \begin{array}{ll} \frac{{\delta_{s}^{2}} + {\delta_{i}^{2}}}{(t_{i}-t_{s})^{2}} + \left(\frac{1}{t_{i}-t_{s}} + \frac{1}{t_{2}-t_{1}}\right) & \text{if}\; t_{1} < t_{2} \leq t_{i},\\ \vspace*{9pt} \frac{\sigma^{2}(t_{1}) + \sigma^{2}(t_{2}) - 2\left(\frac{t_{1}-t_{s}}{t_{i}-t_{s}}\right)\left(\frac{t_{f}-t_{2}}{t_{f}-t_{i}}\right){\delta_{i}^{2}}}{(t_{2}-t_{1})^{2}} & \text{if \(t_{i} < t_{2} \leq t_{f}\),}\\ \vspace*{6pt} \frac{\sigma^{2}(t_{1}) + \sigma^{2}(t_{2})}{(t_{2}-t_{1})^{2}} & \text{otherwise.} \end{array} \right. $$
Let μ
V
and \({\sigma _{V}^{2}}\) be the parameters of the velocity distribution over a time interval [t
1,t
2]. Speed (the absolute value of velocity) over this interval then has a Rice distribution with parameters |μ
V
| and σ
V
. The direction of this velocity has a distribution with density
$$\begin{aligned} f(\gamma) =&\, \frac{e^{-\frac{\nu^{2}}{2}}}{2\pi} + \frac{\nu\cos\eta}{2\sqrt{2\pi}} \exp\left(\frac{\nu^{2}\left(\cos^{2}\eta-1\right)}{2}\right)\\ &\times\left(1+\text{erf}\left(\frac{\nu\cos\eta}{\sqrt{2}}\right)\right), \end{aligned} $$
where \(\nu = \frac {|\boldsymbol {\mu }_{V}|}{\sigma _{V}}\) is the noncentrality of the velocity distribution and η=atan2(μ
V
)−γ is the angle between the direction of the mean and the direction under consideration.
To obtain spatial distributions of speed, we consider the speed over a time interval [t+Δ
t
s
,t+Δ
t
f
], after fixing the position at one time t to a fixed location. If the time interval contains the time at which the position is fixed, i.e. Δ
t
s
≤0 and Δ
t
f
≥0, the position distributions at both endpoints of the interval are independent. The conditioned velocity and speed distributions are then determined from these two distributions. The spatial distribution of speed and the effect of the choice of the time scale (Δ
t
f
−Δ
t
s
) is illustrated in Fig. 5 by the example of the data used in the second case study.
We do not give the details about these position distributions here, but refer to Appendix 1. Let μ
s
, μ
f
, \({\sigma _{s}^{2}}\) and \({\sigma _{f}^{s}}\) represent the respective means and variances of the conditioned positions at both endpoints of the interval. Then by independence of the positions the velocity distribution conditioned on X
t
=x is given by
$${\fontsize{9.3}{6}\begin{aligned} \boldsymbol{V}_{\boldsymbol{x}; t} (t+\Delta t_{s}, t+\Delta t_{f}) &= \frac{\boldsymbol{X}_{t+\Delta t_{f}} - \boldsymbol{X}_{t+\Delta t_{s}}}{\Delta t_{f} - \Delta t_{s}}\\ &\sim \mathcal{N}\left(\frac{\boldsymbol{\mu}_{f} - \boldsymbol{\mu}_{s}}{\Delta t_{f} - \Delta t_{s}}\right)\left(\frac{{\sigma_{s}^{2}} + {\sigma_{f}^{2}}}{\left(\Delta t_{f} - \Delta t_{s}\right)^{2}}\right). \end{aligned}} $$
As discussed before, the speed has a Rice distribution. We determine the average speed at a particular location by computing a weighted average over time of the mean speed. The weight is given by the probability density of the animal’s position at the given time and location. That is,
$$ S(\boldsymbol{x}) = \frac{1}{\int\! f_{\boldsymbol{X}{_{t}}}({\boldsymbol{x}})\,dt} \int f\boldsymbol{X}_{t}({\boldsymbol{x}}) \mathbb{E}\left[|\boldsymbol{V}_{\boldsymbol{x};t}(t+\Delta t_{s}, t+\Delta t_{f})|\right] dt. $$
((1))
Methods for the analysis of the movement speed of vervet monkeys in relation to their environment
Vervet monkeys are group-living primates that are abundant throughout most of sub-Saharan Africa [28]. They occur in stable, mixed-sex groups of typically 25-30 animals that consist of multiple adult males and females along with their offspring. Patterns of home range selection and general space use are strongly affected by external environmental factors such as primary productivity and vegetation structure [29] as well as the distribution of food, surface water and perceived predation risk [30].
In order to investigate whether the movement speed of animals is similarly affected by external variables, the data used in this case study were collected on a wild group of vervet monkey ranging freely in their natural habitat in Kwazulu-Natal, South Africa, during December2010. A digital telemetry collar (e-obs Type 1A, 69 gper unit, equivalent to just over 2 % of the tagged animal’s body weight; All work at the Inkawu Vervet project was approved by the relevant local authorities (the ethical boards of Ezemvelo KwaZulu-Natal Wildlife and the University of Cape Town, South Africa), and complies with EU-directive 2010/63/EU on the protection of animals used for scientific purposes) was deployed on a single adult female within the group and programmed to obtain GPS-fixes at hourly intervals throughout the daily activity phase of the animals (05:00 – 19:00). Given that vervet monkey groups typically move as coherent units through the landscape, GPS-coordinates obtained from the tagged female were taken to represent the movement of the entire group. Local vegetation density was estimated from a multi-spectral, high-resolution (0.50 × 0.50 m 2 pixel size) satellite image (WorldView II, DigitalGlobe Inc.) obtained over the study-period. From this image, we calculated the Normalized Difference Vegetation Index (NDVI) [31], a well-established spectral correlate of primary productivity and vegetation structure.
In our dynamic BBMM calculations, we did not consider bridges at the beginning of the day that stayed very close (≤50m) to the starting location, as this indicated the monkeys had not commenced moving yet, and similarly at the end of the day near the final location. On the remaining bridges the method by Kranstauber et al. [10] was used to estimate the diffusion coefficient (using a margin of 3 and a window size of 7). The average speed distribution presented in the Results section, was computed at a time scale Δ
t of 5 minutes. Mean speed was computed as defined in Equation 1, over two time intervals relative to the focal point: one directly preceding it (i.e. Δ
t
s
=−Δ
t, Δ
t
f
=0), and one directly following it (i.e. Δ
t
s
=0, Δ
t
f
=Δ
t). If we had used only one of these intervals, we would not have been able to compute a speed near the beginning or end of the daily activity period, which could have resulted in missing values in the distribution. For the analysis with only one diffusion coefficient we used the method by Horne et al. [9]. The R scripts that were used in this analysis are provided as Additional file 5.
Methods for migration of European bee-eaters
This case study deals with the northward migration of the European bee-eater through the Arava Valley in southern Israel. The species is a very common passage migrant during both autumn and spring throughout the entire country [32]. In the 2005 and 2006 spring migration seasons, a total of 11 bee-eaters were trapped, marked and tagged with radio transmitters. Using portable systems, birds were followed over a total of 810 km during which their flight mode was established throug h both wing flap signals and the unique signature of circling flight in the recorded transmission (for details see [21,22]; Bee-eater trapping permits were obtained from the Israeli Nature and Parks Authority (permits 2005/22055, 2006/25555) and the experimental procedure was approved by the Animal Care and Use Committee of the Hebrew University of Jerusalem (permit NS-06-07-2)). Trajectories were annotated with simulated atmospheric conditions at appropriately short and small scales using the Regional Atmospheric Modeling System (RAMS; [33,34]). The relationship between bird flight mode (flapping, soaring-gliding and mixed flight) and atmospheric conditions are described in [22]. That study confirmed that turbulence kinetic energy (TKE, in m 2/s 2), as an indicator of convective updraught intensity in the atmosphere, facilitates soaring and gliding. In the current study, the relationship between bird flight mode and the movement path was estimated by calculating the effects of bird flight mode on the animal diffusion coefficient in the BBMM [16].
The relation between TKE and flight mode as well as between flight mode and the diffusion coefficient was determined by considering only movement stretches with flapping and pure soaring-gliding modes (hence omitting the mixed flight modes). The mixed flight mode is highly variable and biomechanically not as well defined as flapping or soaring-gliding flight.
A univariate logistic model was fitted to estimate the fraction of soaring flight (s) as a function of TKE. Model and parameter significance was tested for this model (using a 0.05 significance level), as well as the overall classification error. Subsequently, the BBMM was fitted to segments with flapping flight and soaring-gliding flight separately, resulting in estimates for the diffusion coefficients for each of these flight modes. Next, the diffusion coefficient for the mixed flight mode was estimated by weighting the two diffusion coefficients with the fraction of time spent in each flight mode:
$$D_{m} = (1-s) D_{f} + s\cdot D_{s}, $$
where D
m
, D
f
and D
s
refer to the diffusion coefficients of respectively mixed, flapping and soaring-gliding flight. The fraction s is obtained from the aforementioned logistic model. Using this parameterisation, the complete flight trajectories are estimated per bird by the BBMM.
In addition to the estimated model coefficients, the results of this analysis are presented in the form of probability maps of movement for selected individuals, showing not only the most likely movement path but also the uncertainty in this as a function of distance between observation points and flight mode (as illustrated in Fig. 4). The R scripts that were used in this analysis are provided as Additional file 6.