- Methodology article
- Open Access
Bivariate Gaussian bridges: directional factorization of diffusion in Brownian bridge models
© Kranstauber et al.; licensee BioMed Central Ltd. 2014
- Received: 18 November 2013
- Accepted: 5 February 2014
- Published: 1 March 2014
In recent years high resolution animal tracking data has become the standard in movement ecology. The Brownian Bridge Movement Model (BBMM) is a widely adopted approach to describe animal space use from such high resolution tracks. One of the underlying assumptions of the BBMM is isotropic diffusive motion between consecutive locations, i.e. invariant with respect to the direction.
Here we propose to relax this often unrealistic assumption by separating the Brownian motion variance into two directional components, one parallel and one orthogonal to the direction of the motion.
Our new model, the Bivariate Gaussian bridge (BGB), tracks movement heterogeneity across time. Using the BGB and identifying directed and non-directed movement within a trajectory resulted in more accurate utilisation distributions compared to dynamic Brownian bridges, especially for trajectories with a non-isotropic diffusion, such as directed movement or Lévy like movements. We evaluated our model with simulated trajectories and observed tracks, demonstrating that the improvement of our model scales with the directional correlation of a correlated random walk.
We find that many of the animal trajectories do not adhere to the assumptions of the BBMM. The proposed model improves accuracy when describing the space use both in simulated correlated random walks as well as observed animal tracks. Our novel approach is implemented and available within the “move” package for R.
- Dynamic Bivariate Gaussian bridge
- Dynamic Brownian bridge movement model
- Utilisation distribution
- Animal tracking
- Home range and space use modelling
The availability of global positioning system (GPS) and affordable satellite telemetry has revolutionised the study of animal movement, allowing users to estimate the location of individuals at high spatial and temporal resolution. Consequently, the amount of highly-resolved data has increased by orders of magnitude providing potential improvements and new challenges for the analysis of animal movement. For many applications, among others habitat and home range analysis, it is important to convert these trajectories into spatially explicit probability distributions. Previously, methods for estimating animal space use, such as minimum convex polygons (MCP) or kernel home range estimates, ignored the temporal sequence of the locations and instead concerned with the statistical independence of successive locations of the tracked animals. Now however, novel approaches can take profit from richer data sets by incorporating time (i.e. temporal autocorrelation) into the modelling of space use.
Brownian bridge movement models (BBMM) estimate space use intensity by stochastically modelling the movement of animals between any two consecutive locations [1, 2]. Brownian bridges connect two consecutive locations by conditional Brownian random walks that start at a given location and end at the following location with a duration equal to the observed time lag between the two locations. Brownian random walks are one of the simplest random walk models and assume movement in any direction to be equally likely, movement steps are drawn from a bivariate Gaussian distribution . The interpolation, or bridge, incorporates both an empirically derived tendency to diffuse and an estimation of the error associated to the true locations. The BBMM has been extended to account for changes in the movement behaviour of animals (dynamic Brownian Bridge movement models: dBBMM ) by allowing the diffusion parameter of the Brownian motion to change according to changes in the behaviour of the animal along its trajectory [5, 6].
Brownian bridges have the benefit of being efficiently calculated and can be fitted straightforwardly to existing trajectories. Further, these bridges can convert a set of locations into a well-defined and spatially-explicit probability density function that models space use intensity (i.e., a normalized and rasterized spatial probability distribution). Such a probability surface can easily be associated to the corresponding landscape information (e.g. [7, 8]). As in any modelling procedure, however, there are constraining assumptions, i.e. the movement between locations is assumed to be diffusive (normal diffusion) and isotropic. In its dynamic version (i.e., the dBBMM), the normal diffusion assumption is somewhat relaxed, as that the dispersion parameter is allowed to vary over time (different diffusion coefficients can be estimated), accounting for rather complex multi-scale, composite-Brownian type of motions [9, 10]. Nevertheless, the movement between any two locations is equally likely to happen in all directions and is quantified as an isotropic diffusive process.
Directional bias has been incorporated differently in widely different movement models (e.g., correlated random walks, biased random walks, multi-scale random walks) (e.g. [11, 12]). Benhamou  introduced the idea of directional bias in bridging models. Biased Random Bridges  rely on the advection-diffusion model, which shows both a diffusive and an advective component, and can generate long-term drift. Essentially Benhamou assumes that the role of anisotropy in generating a “direction in diffusion” is minor and that most of the effect should come from the advection. It happens that the advection component disappears in bridge calculations. Because of this, classical Brownian bridge formulations can be used within Biased Random Bridges, except for diffusion anisotropy. In biased random walks the emergence of diffusion anisotropy depends on different movement parameters that shape the turning angle distributions and can contribute to directional bias . Given that in Biased Random Bridges the advection component is lost a relevant question remains: does taking diffusion anisotropy into consideration in bridge models improve space use estimations?
Bivariate Gaussian bridges
Here we present a novel approach to modelling animal movement as a generalisation of the BBMM: the Bivariate Gaussian bridges (BGB), which allows us to factorize diffusion (i.e. Brownian variance estimates ( )) in two elementary directions: the direction toward the next location and the direction orthogonal to it. In the case of the BBMM (and dBBMM) the probability density around the mean position of the animal between any two locations is assumed to be isotropic. In BGB, however, we factorize the movement variance in a parallel and an orthogonal component to the straight line (with constant speed) connecting two consecutive locations (a segment). This results in two normally distributed probability densities Pr p and Pr o .
Finally, all observation techniques have an error associated with estimating the true position of the animal. We add this location error distribution, which, for the sake of simplicity is assumed to be a bivariate normal distribution (standard deviation δ), to the probabilities Pr p and Pr o to account for this measurement error.
The index can vary between −1 and 1, where 0 means σ m,p and σ m,o are equal in size corresponding to Brownian motion. In case I d →1 there is no orthogonal diffusion an thus all movement is along the straight line, whereas I d →−1 indicates there is no parallel diffusion. Various studies have found turning angles to be informative measure for inferring behavioural state [15–17]. The correlation in heading has also been investigated in various contexts such as search efficiency  and orientation analysis  and has been measured by the sinuosity index . We therefore think that an index capturing heading fluctuations of the animal between succssessive locations could be indeed a very informative measure of behaviour, that in combination with other measures could determine behavioural modes classification. For the variance estimation we used the dBBMM framework , accounting for changes in the diffusion variance across time. Hence BGB incorporate changes in the movement characteristics along long journeys and, since I d can change continuously, provide additional information on changes in the directionality of the movement. Bivariate Gaussian bridge models should not only derive space use quantification more accurately (by estimating movement anisotropies), but should also allow to explore temporal patterns of I d , that can be used to elucidate changes in movement behaviour.
In the following, we formally introduce the BGB framework and assess its improvement (compared to the dBBMM) in predicting the locations of simulated correlated random walks with varying degrees of correlation and step lengths. Finally, we use a series of empirical trajectories of different animal species to assess how realistic the isotropy assumption is by investigation of I d and comparing the estimated space use derived from the BGB (anisotropic) with the one derived from Brownian (isotropic) bridge movement models.
Finally, since changes in the behaviour of tracked animals can lead to changes in magnitude and relative proportion of σ m,p and σ m,o over time, we allow the variance estimation in the BGB to dynamically take different values, resulting in dynamic Bivariate Gaussian bridges (dBGB). Using the Bayesian Information Criterion (BIC) in conjunction with estimation of variance within and between sections of the trajectory we identify significant changes in the composition of σ m,p and σ m,o as introduced by Kranstauber et al.  based on the work of Gurarie et al. . This allows the model to track changes in a trajectory along both, scale of movement and in the directionality. The algorithm uses a window of a set number of locations in which the best values for σ m,p and σ m,o are searched for. Within this window, both σ m,p and σ m,o are calculated over the entire length of the window as well as for any combination of subsections before and after possible break points. A likelihood value is calculated for each description of the window that contains either no break, a single break in σ m,p , a single break in σ m,o , or a break in both. These likelihoods are then compared using the BIC. For each window we choose the optimal descriptor based on the BIC. The window is then moved 1 step through the trajectory and the likelihood optimization is repeated. For each segment we thus obtain multiple estimates of σ m,p and σ m,o that are averaged using the variances. For the dBGB case, σ m,p and σ m,o are allowed to change independently. This allows in total for 4 different change scenarios: no change, a change in either σ m,p or σ m,o , or a change in both. Both, the size of the window and the margin, are parameters of the algorithm that define the granularity of behavioural change detection in the dynamic estimation of σ m,p and σ m,o in dBGB, or σ m in the dBBMM. Margins are the minimal number of locations used for variance estimation, and thus define the number of possible changes within the window. Larger windows have more power to identify changes in behaviour but come at the cost of not being able to describe frequent changes in behaviour. Since more variables (σ m,p and σ m,o vs. σ m ) are involved, it may be sensible to use larger windows and margins for a more accurate approximation of the true values. It is important to explore various parameter setting and assess if the results make sense and describe the expected behavioural changes. One possible way to do this is by using cross-validations .
Utilisation density calculation
With σ m,p and σ m,o for every segment known, the probability of utilisation of an area (e.g., a map grid cell) is defined by the orthogonal (Δ o ) and parallel distances (Δ p ) from the grid cell center to the expected location on the connecting line of the segment using a normal distribution. To calculate the integrated probability density, each segment is split into equally sized integration steps proportional to the time interval between the two consecutive locations defining the segment. For each time step the expected location can be calculated, using this expected location Δ p and Δ o to each grid center can be calculated. These functions are implemented in R  and will be available within the move package . For details we refer to the Methods Section at the end of the manuscript.
We evaluate the dBGB first on a set of simulated trajectories and then apply it to GPS tracks of a selection of different species.
To explore the behaviour of σ m,p and σ m,o , we simulated a variety of random trajectories using a correlated random walk function ( implemented by Calenge et al. ). We varied both the amount of correlation within the random walk from 0.5 to 0.999 as well as the scale of step lengths from 0.1 to 10. Each track consisted of a total of 249 (odd number required for the leave-one-out approach) locations whereby a circular Gaussian error with a standard deviation of 0.01 was added to each location.
The dBGB variances as well as the dBBMM variances were estimated on these trajectories. Since the simulated tracks did not contain any behavioural change, window and margin sizes, which define the granularity of behavioural change detection in the dynamic estimation of σ m,p and σ m,o in dBGB as well as σ m in dBBMM, were chosen to be large (margin: 31 and window size: 71). In order to evaluate the performance of the estimation of the utilisation distribution (UD), we compared the predicted probability densities of the dBGB to those of the dBBMM for these trajectories. To test performance of each method we did cross validations by excluding 35 locations, that were regularly spaced in time. The UD was calculated using each method without the selected locations. We then examined the performance by extracting the UD values at the omitted locations and calculating an index of performance. The performance index was calculated as the geometric mean of the probability densities of dBGB divided by the probability densities of dBBMM (following ) where values greater than 1 indicated better performance of dBGB.
An overview of the different tracks for which we evaluated the dBGB model
I d (>95 % )
I d (<95 % )
Eurasian eagle-owl (Bubo bubo)
Fisher (indiv #1) (Martes pennanti)
Fisher (indiv #2) (Martes pennanti)
Straw-coloured fruit bat (Eidolon helvum)
Turkey vulture (Cathartes aura)
Waved albatross (Phoebastria irrorata)
White stork (Ciconia ciconia)
The variances of the tracks were calculated using a window size of 39 and margin of 15. We also simulated 1000 Brownian motion tracks with 39 locations (same length as the window size), as a confidence interval for the directionality index. We then evaluated how many of the segments fell outside the expected 95% interval. If substantially more than 5% of the I d values fall outside the confidence interval we interpreted this as an indication that the trajectory as a whole did not correspond to Brownian motion.
To compare the performance of dBGB against dBBMM also in real animal trajectories, we again calculated the performance index by omitting 50 locations, which were used for a cross-validation. For this calculation we used the same window size and margin of 39 and 15.
Finally, we present the track of one individual straw coloured flying fox (Eidolon helvum) to visually highlight the obtained UD contours from dBGB in comparison with those obtained from the dBBMM. This African fruit bat roosts in the colony during the day and moves in a very directed manner to individual fruiting trees to forage during the night. We excluded segments during daytime where no movement was recorded.
Investigating the variances of different tracks showed that large parts of the tracks deviated from what would be expected under a Brownian motion assumption. All 7 empirical trajectories in our study had a median I d larger than 0 and a large proportion (on average 0.66) of the segments fell above the 95% percentile of the distribution of the simulated Brownian tracks (Table 1). If all tracks adhered to the assumption of Brownian motion, then we would expect that only 2.5% of the locations would fall above this interval. The index over time revealed that the tracks showed extended periods of time with a high I d values interspersed by bursts of low I d values (Figure 4). The cross validation on the trajectories showed that the dBGB resulted in either a similar or better fit compared to the UD derived from dBBMM (Table 1).
Our method demonstrates the potential accuracy improvements when computing the UD, and highlights the constraints of initial Brownian bridge models  that assume isotropic and homogeneous diffusion across time. The method presented here relaxes both assumptions by: i) locally characterizing diffusion, thus becoming variable in time, and ii) analysing the diffusion properties across two major and orthogonal axes of motion, thus becoming variable in space. The directional bias in trajectories can be a product of various processes (e.g. correlated random walks, biased random walks, landscape features). Although we only formally investigate correlated random walks we think that dBGB performs equally well if the directional bias has a different cause because the variance estimates are adjusted to the observed trajectory.
Biased Random Bridges  assume an advective and a diffusive component, which can incorporate the effect on the directional bias attributable to diffusion anisotropy. More mechanistic insight can be extracted from a bridging model by decoupling advection from diffusion, and thus, Biased Random Bridges represent a clear conceptual improvement with respect to standard BBMM . Nevertheless, Biased Random Bridges do not modify nor improve the estimated UD. Instead, BGB can approximate the idea of advection/diffusion by considering the diffusive non-isotropic process, and thereby improve the accuracy of space use and utilisation density distributions. We did not conduct a direct comparison between the dBGB and the Biased Random Bridges because for the latter no dynamic version is defined. It would be hard to assess where to attribute the differences to.
Worth noting are the jagged contours that appear when the directionality index (I d ) is high, such as the high probabilities right in front and behind the observed locations. Within strongly directed movement periods these probabilities overlap with the previous and next segment. It is likely that these probabilities are visible in contour lines when an animal starts or stops a directed movement period, where the transition to a stop causes extension of the contours behind the observed stop location. These jagged contours can for example be observed in Figure 5 in the upper right corner where shuttling between the localities occurs. Other bridge models suffer from similar artifacts but they might be less conspicuous because the resulting contours are smoother. In any case, it is important to note, that despite these artefacts, UD estimates are overall more accurate than with former models, and that the previous models also contain biases even though perhaps less conspicuous.
In general, the dBGB predicts the location of the observed trajectories better and only in some cases performs slightly worse than the dBBMM. Better dBBMM performance over dBGB may be due to the fact that the estimation of more variables (σ m,p and σ m,o vs. σ m ) increases the noise. The I d provides a measure of directedness independent of the step sizes of the movement. It could also be used as an indication of where the largest differences between dBBMM and the dBGB are to be expected: essentially where the largest deviations of I d from 0 are observed. When the time interval between observations increase, the directional persistence of the correlated random walk decreases . We would therefore expect that I d →0 if the time interval between observations increase, this would mean the difference between the dBGB and dBBMM decreases. When the location error is of the same order of magnitude as the movement variances, the division between movement variance and location error becomes more difficult, resulting in estimates where one of the variances becomes 0. Also the directionality index (I d ) becomes scale-dependent when location errors are high.
Although we do not investigate it here, it is likely that the estimated I d , σ m,p and σ m,o contain relevant intra and inter individual variation. For example, migrating versus sedentary herds of caribou have a far narrower turning angle distribution  which would result in a higher I d . Hence, studying I d , σ m,p , and σ m,o spatiotemporal patterns across individuals, ecological contexts, or species, could provide more mechanistic insights into animal home range and space use behaviour.
It is clear that many observed trajectories do not adhere to the assumptions of isotropic, homogeneous Brownian motion. Our model had the highest performance gain for correlated random walks with high directional correlation. Further work defining analytical descriptions of bridge functions for frequently used random walk models (e.g., correlated random walks, Lévy walks, or continuous-time random walks) is needed since the dBGB does not formally describe the probability density of any of these random walks.
For notational convenience we followed the notation and variable definitions of Horne et al. .
Probability density function
where z is any location in the space.
Including location error
The source code for the methods are available in the R package “move”.
The International Max Planck Research School for Organismal Biology for support of BK; F.B. acknowledges Grants RyC-2009-04133, and BFU2010-22337 from the Ministry of Science and Innovation (Spain), and the Human Frontier Science Program (Ref. RGY0084/2011); Rechenzentrum Garching of the Max Planck Society for support and maintenance of the computational cluster; Jordi Pages (CEAB-CSIC) for initial discussions about variance decomposition and Biased random walks with FB, and Riek van Noordwijk and Scott LaPoint for comments on the manuscript; Track owners for providing sample data: Fisher data, Scott LaPoint, supported by NSF (Grant 0756920 to Roland Kays) and National Geographic Society Waitt Grant Program (Grant W157-11 to SL), Eagle owl data, commissioned by Paul Voskamp (province of Limburg, the Netherlands) collected by René Janssen (Bionet natuuronderzoek), Stork Data, Martin Wikelski, HUJ-MPIO (Hebrew University Jerusalem, PI Ran Nathan) and MPIO joint Deutsch-Israelisches Project Grant and the Storchenhof Lohburg, Turkey Vulture data, Keith Bildstein and the Hawk Mountain Sanctuary, Waved Albatross data, Sebatian Cruz and the Galapagos National Park Service, Straw-coloured fruit bat, Dina Dechmann, Jakob Fahr and Martin Wikelski. We also thank three anonymous reviewers for their valuable comments.
- Bullard F: Estimating the home range of an animal: a Brownian bridge approach . Master’s thesis, The University of North Carolina 1999.Google Scholar
- Horne JS, Garton EO, Krone SM, Lewis JS: Analyzing animal movements using Brownian bridges . Ecology 2007, 88:2354–2363. 10.1890/06-0957.1View ArticleGoogle Scholar
- Codling EA, Plank MJ, Benhamou S: Random walk models in biology . J R Soc Interface 2008, 5:813–834. 10.1098/rsif.2008.0014View ArticleGoogle Scholar
- Kranstauber B, Kays R, LaPoint SD, Wikelski M, Safi K: A dynamic Brownian bridge movement model to estimate utilization distributions for heterogeneous animal movement . J Anim Ecol 2012, 81:738–746. 10.1111/j.1365-2656.2012.01955.xView ArticleGoogle Scholar
- LaPoint S, Gallery P, Wikelski M, Kays R: Animal behavior, cost-based corridor models, and real corridors . Landsc Ecol 2013, 28:1615–1630. 10.1007/s10980-013-9910-0View ArticleGoogle Scholar
- Pages JF, Bartumeus F, Hereu B, Lopez-Sanz A, Romero J, Alcoverro T: Evaluating a key herbivorous fish as a mobile link: a Brownian bridge approach . Mar Ecol Prog Ser 2013, 492:199–210.View ArticleGoogle Scholar
- Takekawa JY, Newman SH, Xiao X, Prosser DJ, Spragens KA, Palm EC, Yan B, Li T, Lei F, Zhao D, Douglas DC, Muzaffar SB, Ji W: Migration of waterfowl in the East Asian flyway and spatial relationship to HPAI H5N1 outbreaks . Avian Dis 2010, 54:466–476. 10.1637/8914-043009-Reg.1View ArticleGoogle Scholar
- Sawyer H, Kauffman MJ, Nielson RM, Horne JS: Identifying and prioritizing ungulate migration routes for landscape-level conservation . Ecol Appl 2009, 19:2016–2025. 10.1890/08-2034.1View ArticleGoogle Scholar
- Petrovskii S, Mashanova A, Jansen VAA: Variation in individual walking behavior creates the impression of a Lévy flight . Proc Natl Acad Sci USA 2011, 108:8704–8707. 10.1073/pnas.1015208108View ArticleGoogle Scholar
- Reynolds AM: Distinguishing between Lévy walks and strong alternative models . Ecology 2012, 93:1228–1233. 10.1890/11-1815.1View ArticleGoogle Scholar
- Bartumeus F, da Luz MGE, Viswanathan GM, Catalan J: Animal search strategies: a quantitative random-walk analysis . Ecology 2005, 86:3078–3087. 10.1890/04-1806View ArticleGoogle Scholar
- Gautestad AO, Mysterud I: Intrinsic scaling complexity in animal dispersion and abundance . Am Nat 2005, 165:44–55. 10.1086/426673View ArticleGoogle Scholar
- Benhamou S: Dynamic approach to space and habitat use based on biased random bridges . PLoS One 2011, 6:e14592. 10.1371/journal.pone.0014592View ArticleGoogle Scholar
- Codling EA, Bearon RN, Thorn GJ: Diffusion about the mean drift location in a biased random walk . Ecology 2010, 91:3106–3113. 10.1890/09-1729.1View ArticleGoogle Scholar
- Bailey H, Thompson P: Quantitative analysis of bottlenose dolphin movement patterns and their relationship with foraging . J Anim Ecol 2006, 75:456–465. 10.1111/j.1365-2656.2006.01066.xView ArticleGoogle Scholar
- Fryxell JM, Hazell M, Börger L, Dalziel BD, Haydon DT, Morales JM, McIntosh T, Rosatte RC: Multiple movement modes by large herbivores at multiple spatiotemporal scales . Proc Natl Acad Sci USA 2008, 105:19114–19119. 10.1073/pnas.0801737105View ArticleGoogle Scholar
- Gurarie E, Andrews RD, Laidre KL: A novel method for identifying behavioural changes in animal movement data . Ecol Lett 2009, 12:395–408. 10.1111/j.1461-0248.2009.01293.xView ArticleGoogle Scholar
- Bartumeus F, Catalan J, Viswanathan GM, Raposo EP, da Luz MGE: The influence of turning angles on the success of non-oriented animal searches . J Theor Biol 2008, 252:43–55. 10.1016/j.jtbi.2008.01.009View ArticleGoogle Scholar
- Benhamou S: Detecting an orientation component in animal paths when the preferred direction is individual-dependent . Ecology 2006, 87:518–528. 10.1890/05-0495View ArticleGoogle Scholar
- Benhamou S: How to reliably estimate the tortuosity of an animal’s path: straightness, sinuosity, or fractal dimension? J Theor Biol 2004, 229:209–220. 10.1016/j.jtbi.2004.03.016View ArticleGoogle Scholar
- R Development Core Team: R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing; 2012. ISBN 3–900051–07–0Google Scholar
- Kranstauber B, Smolla M: Move: visualizing and analyzing animal track data . 2013.Google Scholar
- Kareiva PM, Shigesada N: Analyzing insect movement as a correlated random walk . Oecologia 1983, 56:234–238. 10.1007/BF00379695View ArticleGoogle Scholar
- Calenge C, Dray S, Royer-Carenzi M: The concept of animals’ trajectories from a data analysis perspective . Ecol Inform 2009, 4:34–41. 10.1016/j.ecoinf.2008.10.002View ArticleGoogle Scholar
- Wikelski M, Kays R: Movebank: archive, analysis and sharing of animal movement data . World Wide Web electronic publication Accessed in 2014, http://www.movebank.org
- LaPoint S, Gallery P, Wikelski M, Kays R: Data from: animal behavior, cost-based corridor models, and real corridors . 2013.http://dx.doi.org/10.5441/001/1.2tp2j43g Google Scholar
- Mandel JT, Bildstein KL, Bohrer G, Winkler DW: Movement ecology of migration in turkey vultures . Proc Natl Acad Sci USA 2008, 105:19102–19107. 10.1073/pnas.0801789105View ArticleGoogle Scholar
- Dodge S, Bohrer G, Weinzierl R, Davidson SC, Kays R, Douglas D, Cruz S, Han J, Brandes D, Wikelski M: The environmental-data automated track annotation (Env-DATA) system: linking animal tracks with environmental data . Movement Ecol 2013, 1:3. 10.1186/2051-3933-1-3View ArticleGoogle Scholar
- Cruz S, Proaño CB, Anderson D, Huyvaert K, Wikelski M: Data from: The Environmental-Data Automated Track Annotation (Env-DATA) System: linking animal tracks with environmental data . 2013.http://dx.doi.org/10.5441/001/1.3hp3s250 Google Scholar
- Bergman CM, Schaefer JA, Luttich SN: Caribou movement as a correlated random walk . Oecologia 2000, 123:364–374. 10.1007/s004420051023View ArticleGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.