In the following subsections the metrics are defined and their theoretical properties are described. A summary is proposed in Table 1. Considering two individuals named A and B, the position of A (resp. B) at time t is denoted by \(X_{t}^{A}\) (resp. \(X_{t}^{B}\)). The distance between A at time t1 and B at time t2 will be referred to as \(d_{t_{1},t_{2}}^{A,B}\). When the distance between two individuals is regarded at simultaneous time, this will be shortened to \(d_{t}^{A,B}\). Whenever possible, metrics introduced by different authors but that are actually very similar in their definition, are grouped under a unified name and a general definition.
Proximity index (Prox)
The proximity index (Prox in [5]) is defined as the proportion of simultaneous pairs of fixes within a distance below an ad hoc threshold (Fig. 1). Other metrics in the literature are actually analogous to Prox: the coefficient of association (Ca) [12] and the IAB index [4]. Denoting by T the number of pairs of fixes in the dyad, we propose a unified version of those metrics using a kernel K (formula 1):
$$ Prox_{K,\delta} = \frac{1}{T} \sum\limits_{t=1}^{T} K_{\delta}\left(X_{t}^{A}, X_{t}^{B}\right), $$
(1)
where δ is a distance threshold parameter.
Choosing \(K_{\delta }(x,y) =\mathbbm {1}_{\{ \| x-y \| < \delta \}}\) (\(\mathbbm {1}_{\{\}}\) represents the indicator function) as a kernel leads to the Prox metric in [5], denoted by \(Prox_{\mathbbm {1},\delta }\) henceforward. Instead, choosing Kδ(x,y)= exp(−∥x−y∥2/(2δ2)) gives the IAB index. Regarding Ca, for simultaneous fixes, its definition becomes exactly the same as \(Prox_{\mathbbm {1},\delta }\) (using Ca’s adaptation to wildlife telemetry data shown in [38]).
Most of the proximity-related metrics are based on symmetric kernels and depend only on the distance between A and B; therefore, the formula notation (1) can be simplified as:
$$ Prox_{K,\delta} = \frac{1}{T} \sum\limits_{t=1}^{T} K_{\delta}\left(X_{t}^{A}, X_{t}^{B}\right) = \frac{1}{T} K_{\delta}^{+}. $$
(2)
If the distance between two individuals is below the threshold δ during their whole tracks, \(Prox_{\mathbbm {1},\delta }\) will be 1 (and 0 in the opposite case). \(Prox_{\mathbbm {1},\delta }\) might be interpreted as the proportion of time the two individuals spent together. This interpretation is, of course, threshold dependent. The IAB index provides a smoother measure of the average proximity between two individuals along the trajectory. Proximity is thus dependent on the choice of a δ parameter and of a kernel function. Graphical examples illustrating the differences in \(K_{\delta }(x,y) =\mathbbm {1}_{\{ \| x-y \| < \delta \}}\) and Kδ(x,y)= exp(−∥x−y∥2/(2δ2)) are in Additional file 1.
Coefficient of Sociality (Cs)
The Coefficient of Sociality (Cs) [26] compares the mean (Euclidean) distance between simultaneous pairs of fixes (DO) against the mean distance between all permutations of all fixes (DE).
$$ Cs = \frac{D_{E} - D_{O}}{D_{E} + D_{O}} = 1 - 2 \frac{D_{O}}{D_{E} + D_{O}}, $$
(3)
where
$$D_{O} = \left(\sum\limits_{t=1}^{T} d_{t}^{A,B} \right)/T, $$
and
$$D_{E} = \left(\sum\limits_{t_{1}=1}^{T}\sum\limits_{t_{2}=1}^{T}d_{t_{1},t_{2}}^{A,B}\right)/T^{2}. $$
Kenward et al. [26] stated that Cs belongs to [−1,1], and it has been used as a symmetrical index since. Nevertheless, that is not true. Cs equals 1 if and only if DO=0 and DE≠0, which occurs only when the two individuals always share the exact same locations. However, Cs equals −1, if and only if DE=0 and DO≠0, which is impossible (it could asymptotically approach to 1 for very large series when DO approaches infinity). Cs equals 0 when DO=DE.
If all simultaneous fixes are very proximal but not in the same locations, Cs would approach 1 (how close to 1 would depend on the value of DE as illustrated in the right hand side of Eq. 3). Moreover, only if DE<DO, Cs can take a negative value. For Cs to take a largely negative value, the difference in the numerator should be very large compared to the sum in the denominator; in Additional file 2 we show how implausible that situation is and how sensitive it is to the length of the series. The latter makes Cs from dyads of different length difficult to compare, because their real range of definition would differ. This fact is neither evoked in the work that introduced the metric [26] nor in the ones that evaluated this and other metrics [38, 41], despite the fact that in those works no value lower than −0.1 was obtained.
Indeed, [26] assumed that the permutation of all fixes is a way to represent locations of independent individuals. While this is questionable, some modified versions, as the one proposed by [62], use correlated random walks as null models and simulated independent trajectories under these models to replace DE by a more realistic reference value. Thus, a generalized version of Cs would be:
$$ Cs = \frac{ D_{chance} - D_{O}}{D_{chance} + D_{O}}, $$
(4)
where Dchance is defined through a user-chosen movement model for independent trajectories.
The Half-weight Association Index (HAI)
The Half-weight Association Index (HAI) proposed by [10] measures the proportions of fixes where individuals are close to each other (within a user-defined threshold). By that definition, HAI is exactly the same as \(Prox_{\mathbbm {1},\delta }\). However, HAI was popularized by [2] in another form that did not consider all fixes for the computation of the metric, but used counts with respect to a reference area (called overlapping zone in the original paper):
$$ HAI = \frac{K_{\delta}^{+}}{K_{\delta}^{+} + \frac{1}{2}(n_{A0}+n_{0B})} $$
(5)
where nAB (resp nA0; n0B; n00) is the number of simultaneous occurrences of A and B in the reference area SAB (resp. simultaneous presence of A and absence of B; simultaneous absence of A and presence of B; simultaneous absence of A and absence of B), and where \(K_{\delta }^{+}\) is computed over the reference area.
It is worth noticing that the HAI adaptation proposed by [2] does not correctly account for spatial joint movement, as would do a \(Prox_{\mathbbm {1},\delta }\) version constraint to the reference area; i.e. the denominator should be equal to nAB + nA0 + n0B, which is the total number of simultaneous fixes where at least one individual is in the reference area.
The dependence to the definition of an overlapping zone or reference area is discussed in the following subsection dedicated to LixnT, which also relies on the definition of a static reference area.
If the individuals remain together (i.e. in the reference area and closer than δ) all the time, HAI is close to 1, and 0 in the opposite case. An example of the computation of HAI under [2]’s definition is given in Fig. 2.
Coefficient of Interaction (L
ixn and L
ixn
T)
Minta [42] proposed a Coefficient of Interaction (Lixn) that assesses how simultaneous are the use and avoidance of a reference area SAB by two individuals:
$$ L_{ixn} = \ln{\left(\frac{n_{AB}/p_{AB} + n_{00}/p_{00}}{n_{A0}/p_{A0} + n_{0B}/p_{0B} }\right)}, $$
(6)
where pAB is the probability, under some reference null model, of finding A and B simultaneously in SAB (the same interpretation as for n when a subscript is 0; see HAI subsection). Attraction between individuals would cause greater simultaneous use of SAB than its solitary use, which would give positive values of Lixn. Conversely, avoidance would translate into negative values of Lixn, since use of SAB would be mostly solitary. A logistic transformation of the metric (LixnT) produces values between 0 (avoidance) and 1 (attraction), making the interpretation easier:
$$ L_{ixn}T = logistic(L_{ixn})= \frac{1}{1+e^{-L_{ixn}}}. $$
(7)
Minta [42] proposed two different approaches for computing the associated probabilities conditionally to the fact that the reference area is known (see examples in Fig. 2 and Table in Additional file 3). In both cases, the probabilities are estimated under the assumptions of independence in movement among the individuals and of uniform utilization of the space. Indeed this latter assumption can be relaxed and pAB can be derived from any kind of utilization distribution (see for instance [20] for the estimation of utilization distribution).
HAI and LixnT (thus Lixn as well) rely heavily on a static reference area – either known or estimated – and on the probabilities of presence within this reference area. The static reference area could be defined, for instance, as the intersection of the respective home ranges of A and B. However, there are many approaches for estimating home ranges, each one relying on particular assumptions about the spatial behaviour of the studied populations [9]. Thus, SAB is not a simple tuning parameter. The way it is defined may completely modify the output. If the reference area is equal to the whole area of movement of the two individuals, then both the numerator and the denominator in the logarithm are equal to infinity and LixnT cannot be derived. That problem could arise for extremely mobile individuals, such as tuna, turtles and seabirds [8], or fishing vessels [6], and avoiding it would require the computation of multiple dynamic reference areas. Therefore, LixnT may be better used for specific cases where the definition of the reference area relies on a deep knowledge of the spatial behaviour of the populations.
Joint Potential Path Area (jPPA)
Long et al. [39] computed the relative size of the potential encounter area at each time step of two individuals’ tracks. Assuming a speed limit ϕ, the potential locations visited between two consecutive fixes define an ellipse (Additional file 4). Then, the potential encounter area corresponds to the intersection between the ellipses of the two individuals (at simultaneous time steps; see Fig. 3). The overall potential meeting area is given by the spatial union of all those potential encounter areas. This area is then normalized by the surface of the spatial union of all the computed ellipses to produce the joint Potential Path Area (jPPA) metric ranging from 0 to 1 (see formula in Table 1). jPPA values close to 0 indicate no potential spatio-temporal overlap, while values close to 1 indicate a strong spatio-temporal match.
Several issues can be discussed here. First, no movement model is assumed and therefore the method confers the same probabilities of presence to every subspace within the ellipse regions. This is clearly unrealistic as individuals are more likely to occupy the central part of the ellipse because they cannot always move at ϕ, i.e. maximal speed. Second, the computation of the ellipses relies strongly on the ϕ parameter. If ϕ is unrealistically small, it would be impossible to obtain the observed displacements and the ellipses could not be computed. By contrast, if ϕ is too large, the ellipses would occupy such a large area that the intersected areas would also be very large (hence a large jPPA value). Alternatively, [36] proposed a dynamic computation of ϕ as a function of the activity performed by the individual at each fix. Within this approach, additional information or knowledge (i.e. other data sources or models) would be required for the computation of ϕ.
Cross sampled entropy (CSE and CSEM)
Cross sampled entropy (CSE) [51] comes from the time series analysis literature and is used for comparing pairs of motion variables [3; 18, e.g.]. It evaluates the similarity between the dynamical changes registered in two series of any given movement measure. Here we present a simplification of the CSE for simultaneous fixes and position series. A segment of track A would be said to be m-similar to a segment of track B if the distance between paired fixes from A and B remain below a certain threshold during m consecutive time steps. If we define Nm as the number of m-similar segments within the series, then CSE can be defined as (the negative natural logarithm of) the ratio of Nm+1 over Nm and might be understood as (the negative natural logarithm of) the probability for an m-similar segment to also be (m+1)-similar. Formally, CSE is defined as:
$$ \begin{aligned} CSE_{\delta}(m) &= -\ln\left\{\frac{ \sum\nolimits_{t=1}^{T-m} \mathbbm{1}{\left\{ \left(\max_{k \in [0,m]} \left| X_{t+k}^{A} - X_{t+k}^{B} \right| \right) < \delta \right\} } }{\sum\nolimits_{t=1}^{T-m} \mathbbm{1}{\left\{ \left(\max_{k \in [0,m-1]} \left| X_{t+k}^{A} - X_{t+k}^{B} \right| \right) < \delta \right\} }}\right\}\\ &=-\ln\frac{N_{m+1}}{N_{m}}, \end{aligned} $$
(8)
A large value of CSE corresponds to greater asynchrony between the two series, while a small value corresponds to greater synchrony.
CSE relies on an ad hoc choice of both m and δ. In practice, it is expected that the movement series of A and B will not be constantly synchronous and that, for a large value of m, Nm could be equal to 0, in which case CSE would tend to ∞. Therefore, the largest value of m such that Nm>0, i.e. the length of the longest similar segment, could be an alternative indicator of similarity between the series (do not confuse with the longest common subsequence LCSS; see [60]). We propose to use this measure (standardized by T−1 to get a value between 0 and 1) as an alternative index of joint movement (formula 9), which we denote by CSEM. An example of a dyad and the computation of its CSEs and CSEM is shown in Fig. 4.
$$ CSEM = \frac{\max \left\{m; N_{m} >0\right\} }{T-1}, $$
(9)
with the convention that max{∅}=0.
Correlations (r
V)
Pearson and Spearman correlations between variables such as longitude, latitude, distance, velocity, acceleration and turning angles from pairs of tracks, have been used as measures of synchrony in several studies (e.g. [16],). Correlations are easy to interpret. Pearson correlation coefficients (Table 1) assess linear correlations, while Spearman correlation coefficients based on ranks statistics capture any functional correlation. The correlation in a given V variable between dyads is denoted by rV.
Dynamic Interaction (DI, D
I
d and D
I
θ)
Long and Nelson [37] argued that it is necessary to separate movement patterns into direction and displacement (i.e. distance between consecutive fixes or step length), instead of computing a correlation of locations [55] which may carry a mixed effect of both components. To measure interaction in displacement, at each time step, the displacements of simultaneous fixes are compared (formula 10).
$$ g_{t}^{\beta} = 1- \left(\frac{\left| d_{t,t+1}^{A} - d_{t,t+1}^{B} \right|}{d_{t,t+1}^{A} + d_{t,t+1}^{B}}\right)^{\beta} $$
(10)
where β is a scaling parameter meant to give more or less weight to similarity in displacement when accounting for dynamic interaction. As β increases, \(g_{t}^{\beta }\) is less sensitive to larger differences in displacement. Its default value is 1. When \(d_{t,t+1}^{A} = d_{t,t+1}^{B}, g_{t}^{\beta }=1\); and when the difference in displacement between A and B at time t is large, \(g_{t}^{\beta }\) approaches zero. For \(g_{t}^{\beta }\) to be 0, one (and only one) of the individuals in the dyad should not move; for a sum of \(g_{t}^{\beta }\) to be equal to zero, at every time t one of the two individuals should not move.
Interaction in direction is measured by
$$ f_{t} = \cos\left(\theta_{t,t+1}^{A} - \theta_{t,t+1}^{B}\right) $$
(11)
where θt,t+1 is the direction of an individual between time t and t+1. ft is equal to 1 when movement segments have the same orientation, 0 when they are perpendicular and −1 when they go in opposite directions.
Long and Nelson [37] proposed 3 indices of dynamic interaction: 1) DId, dynamic interaction in displacement (average of all \(g_{t}^{\beta }\)); 2) DIθ, dynamic interaction in direction (average of all ft); and 3) DI, overall dynamic interaction, defined as the average of \(g_{t}^{\beta } \times f_{t}\) (Table 1). DId ranges from 0 to 1, DIθ from -1 to 1, and DI from -1 (opposing movement) to 1 (cohesive movement). Figure 5 shows an example of the three indices.
Conclusions on the theoretical properties of the metrics
Practical use (C1): While each metric concerns a concrete aspect of joint-movement behaviour, some of them, such as Cs and DI, are harder to interpret. DI mixes up the coordination in displacement and direction. When DI is close to 1, it is certainly explained by high values in both components. When it is close to −1, it is an indication of overall high displacement coordination but in opposite directions. With values around zero, however, it is impossible to know if it is because of displacement or direction or both. For Cs, because obtaining values close to −1 is extremely rare, values around zero and, more particularly, slightly negative values are difficult to interpret. In addition, the maximum attainable value depends on the length of the series, which is likely to vary from dyad to dyad (Additional file 2).
Dependence on parameters (C2): Almost every metric depends on the ad hoc definition of a parameter or component, as summarized in Table 1. This is consistent with the fact that, since there is no consensus on the definition of behaviour [34], and much less on that of collective behaviour, its study depends heavily on the definition that the researcher gives to it. It should be noted that behind each choice of a parameter value, there is also an underlying assumption (e.g. that a distance below a δ value means proximity); the difference is that parameters can be tuned, and a variety of values can be easily tested. HAI and LixnT make a critical assumption of a static reference area, and its definition, which may be tricky for highly mobile individuals, is a key issue for the computation of both metrics. On the other hand, rV and DIθ are the only metrics that do not depend on parameter tuning or assumptions for its derivation; except for the assumptions of correlations being linear, or of linear movement between two successive positions when deriving directions, respectively.