Understanding the spatiotemporal pattern of grazing cattle movement

Understanding the drivers of animal movement is significant for ecology and biology. Yet researchers have so far been unable to fully understand these drivers, largely due to low data resolution. In this study, we analyse a high-frequency movement dataset for a group of grazing cattle and investigate their spatiotemporal patterns using a simple two-state ‘stop-and-move’ mobility model. We find that the dispersal kernel in the moving state is best described by a mixture exponential distribution, indicating the hierarchical nature of the movement. On the other hand, the waiting time appears to be scale-invariant below a certain cut-off and is best described by a truncated power-law distribution, suggesting that the non-moving state is governed by time-varying dynamics. We explore possible explanations for the observed phenomena, covering factors that can play a role in the generation of mobility patterns, such as the context of grazing environment, the intrinsic decision-making mechanism or the energy status of different activities. In particular, we propose a new hypothesis that the underlying movement pattern can be attributed to the most probable observable energy status under the maximum entropy configuration. These results are not only valuable for modelling cattle movement but also provide new insights for understanding the underlying biological basis of grazing behaviour.


Introduction
Animal movement is a highly complex process driven by various random and deterministic mechanisms involving a large number of causing factors [1,2].It has been proposed that spatiotemporal patterns in movement may arise from moving strategies that evolve to optimise foraging efficiency [3,4], decision-making processes in response to external stimuli [5], environmental conditions or landscape features [6,7,8], collective dynamics and social interactions [9,10], memory and home-return behaviour [11], just to name a few.Fully unravelling the complexity of animal movement as well as sorting out the intricate relations between the observed spatiotemporal pattern and various underlying causing factors remains a difficult scientific challenge.
For over a century, our attempts to understand animal movement have been limited to a qualitative level due to the lack of high-quality data that can provide fine-grained spatiotemporal description of movement [2].Recently, new tracking technologies such as the Global Positioning System (GPS) have been deployed in animal tracking to obtain continuous time-resolved moving trajectories with high spatiotemporal resolution.The emerging high-quality movement data enables the application of quantitative analysis and mathematical characterisation on mobility patterns at different spatiotemporal scales.
One common approach to analyse animal movement is to represent the time-resolved trajectory as discrete moving steps under the framework of random walk [12,13,14].
In this context, the dispersal kernel p(r) in space, which characterises the general distribution of step length r in the trajectory, is considered to be a significant footprint of movement [15,16].The detailed functional form of p(r) is indicative of a specific type of random walk and the underlying dynamics of movement.For example, an exponential kernel function p(r) ∼ e −r/rc with r c being the characteristic length scale is the signature of the classic Brownian walk that obeys the central-limit theorem and exhibits a normal-diffusive pattern.A scaling dispersal kernel characterised by a power-law function p(r) ∼ r −γ with γ being the scaling exponent is the signature of the Lévy walk which exhibits high heterogeneity and super-diffusive pattern.Much effort has been devoted to study the dispersal kernel p(r) for different animal species using real movement data, from small insects like honey bees [17], marine life like jelly fish and whales [18], birds like albatross [19,20], to mammals like monkeys [21] and human [16,22,23].One controversial topic that attracts tremendous attention in this line is about whether the observed suggested strong evidence for the existence of Lévy walks in animals, it has been argued that this evidence may come from statistical artefacts or inappropriate manipulation of data [24,25,26].Indeed, the dispersal kernel p(r) only provides partial information on the spatial pattern and does not fully capture all important movement aspects such as the temporal spectrum that depicts the switch among different activity modes over time.
Recently it has been found that scaling phenomena in movement can also arise in the waiting time distribution p(τ w ) that characterises the time span of non-moving period, or the inter-event time distribution p(τ e ) that characterises the time between two successive moving activities [16,11,27].These findings suggest there is a need for more detailed investigation on the spatiotemporal pattern in movement beyond the dispersal kernel.
Here we use a dataset of high-frequency GPS samples to study the movement of grazing cattle.The continuous time-resolved trajectories allow us to gain comprehensive insight into the spatiotemporal pattern.We use a two-state 'stop-and-move' model to describe the mobility pattern, dividing the trajectory into alternate moving and nonmoving states (see Material and Methods).Specifically, the non-moving state indicates that the animal remains within a radius ∆r in space for at least ∆t in time, where ∆r represents the spatial resolution limit in the observation and ∆t is a tuning threshold parameter to specify the minimum time span.The non-moving segment in the trajectory can be viewed as a single point in space, which we call waiting location, with a length τ w ≥ ∆t in time, which we call waiting time correspondingly.On the other hand, the moving state indicates that the animal is in a transition from one waiting location to another, which can be described as a trip (l, τ m ) in the trajectory with l being the distance between the two waiting locations and τ m being the time elapse of the trip.The representation of mobility pattern in this approach is shown in the schematic diagram of Fig. 1.
Under this representation, we observe a very interesting spatiotemporal pattern in which the two activity states are of unique statistical characterization.In particular, we find that the dispersal kernel or trip length distribution p(l) is best described by a hybrid exponential distribution, which indicates that the trajectory has a two-level hierarchical structure in space and each level appears to follow a Brownian walk.This is in contrast to the widely-observed Lévy walk patterns in other species.Despite the absence of scaling law in the spatial dispersal, we find that the waiting distribution p(τ w ) in the time domain is best described by a truncated power-law.Possible underlying mechanisms and ecological implications accounting for this phenomena are discussed (see Discussions).
Our results provide new quantitative insights into grazing cattle movement that are lacking in most previous work.Indeed, understanding grazing/foraging animal movements is not only a critical issue in biological science but also of fundamental importance to many practical issues such as farm and livestock management [28], the maintenance of biodiversity in ecosystems [29] and developing better tracking [30] and virtual fencing technologies [31].

Results
Dispersal kernel in moving state.We first turn attention to the spatial dispersal kernel or the trip length distribution p(l) over the whole population.Obtaining the functional form of p(l) directly from empirical data requires a binning process, which has been known to have statistical distortion for data with a broad distribution [26].To avoid the disadvantage of binning, we use the complementary cumulative distribution P (l) ≡ l p(l)dl for statistical analysis.We process the data using three different parameter sets with ∆r = 5m corresponding to the resolution limit of positioning device and ∆t = 1, 2, 3 mins respectively [33].To describe the dispersal kernel P (l) shown in Fig. 2 a-c, we consider four commonly-used candidate models [32]: (1) power-law; (2) truncated power-law with exponential cut-off; (3) exponential; (4) mixture exponential.Using the maximum-likelihood estimation (MLE) to fit the candidate models and the Akaikeinformation-criterion (AIC) for model selection [24] (see Supporting Information), we find that the best model to describe P (l) is the mixture exponential where l 1 and l 2 are the characteristic lengths in each mixture component, q is a parameter specifying the mixture proportion, and l min = ∆r = 5m is the lower bound in observation.
Another significant statistical feature of the moving state is the trip time distribution p(τ m ), which is also best described by the mixture exponential model, as shown in Fig. 2 d-f.This is consistent with our expectation that the trip time τ m is strongly correlated with the trip length l.
The mixture exponential here indicates that the spatial pattern of grazing cattle is governed by two different Brownian-like dynamics with different characteristic scales, suggesting a hierarchical structure of the movement.If we consider that the landscape is formed by a number of patchy areas, the first exponential distribution will represent the short-range movement that occurs within a patch, and the second exponential distribution with a larger characteristic length will represent inter-patch movements.To better reveal this hierarchy structure in mobility pattern, we perform clustering on the waiting locations using a density-based clustering algorithm DBSCAN which is efficient in discovering significant clusters with irregular shape from noisy data points, as shown in Fig. 3.
After grouping the waiting locations into clusters, the trips in the mobility pattern fall into two categories, intra-cluster trip and inter-cluster trip.We find that the trip length distribution for each of these two types of trips can be well described by a single exponential distribution, as shown in Fig. 4. It is worth noting that the mixture exponential still renders the highest AIC weight among the four candidate models for the intra-cluster trip length distribution, while the single exponential has the highest AIC weight without the mixture exponential.However the difference between the two components in the mixture model is comparably small (l 1 = 11.20,l 2 = 22.05, q = 0.44, ∆t = 2 mins), suggesting that the two components are not strongly distinguishable and a single exponential is a reasonable alternative model in this scenario.We also observe that some long-distance inter-cluster trips are of high similarity, indicating that transitions from one cluster to another are not completely random and spontaneous, but could be driven by a deterministic process such as memory or herding.
Waiting time distribution.To characterise the waiting time distribution, we compare three different models, namely exponential, power-law and truncated power-law with exponential cut-off (see Supporting Information).We observe that the waiting-time distribution p(τ w ) is best described by a truncated power-law distribution with γ being the scaling exponent and τ c w being the structural cut-off.As shown in Fig. 5, varying the threshold parameter ∆t in data processing does not affect the emergence of scaling phenomena.The scaling law in the waiting time distribution is indicative of the heterogeneous grazing dynamics of cattle, which could be related to the landscape hetero-geneity, the complex decision-making dynamics or the energy management of movement (See Discussion).We also find that the waiting time distributions in the main clusters discovered by the DBSCAN algorithm are all well described by a truncated power-law, suggesting the scaling behaviour in waiting time distribution is invariant at the clusterlevel.
Without taking into account correlation between activities, the temporal spectrum of mobility pattern can be approximated as a two-state renewal process where the time span of the alternate moving and non-moving activities are randomly drawn from the distribution functions p(τ m ) and p(τ w ) respectively.To test the validity of this approximation, we measure the pairwise Pearson correlation coefficient between the time span of consecutive activity segments in the following four situations: (1) the non-moving segment and the next moving segment (r = −0.0344,p = 0.0363); (2) the moving segment and the next non-moving segment (r = −0.0616,p = 0.000174); (3) two consecutive non-moving segment (r = 0.0787, p = 1.58 × 10 −6 ); (4) two consecutive moving segment (r = 0.0719, p = 1.24 × 10 −6 ).We find that none of these shows significant correlation.
The result suggests that short-range correlation does not exist in the temporal spectrum, i.e. the time span of the previous activity has little influence on the time span of the next activity, and the temporal dynamics can be approximately described by a two-state renewal process without considering long-range correlation.
Individual mobility pattern.The population-based statistics presented above are not necessarily representative of the individual patterns.It has been suggested that the characteristics of population statistics may differ from their individual counterpart after being aggregated over population.For example, the observed Lévy walk pattern in population may arise from individual heterogeneity [34].To test whether the individual pattern is consistent with the population-based statistics, we use the same model selection procedure to fit the individual statistics (∆r = 5m, ∆t = 2 mins).We find that the trip length and trip time distribution for each individual is best described by the hybrid exponential, with only one exception in the trip length distribution.On the other hand, the waiting distribution for each individual is best described by power-law or truncated power-law (see Supporting Information).This suggests that the composite Brownian walk in space as well as the scaling law in waiting time distribution are not an statistical artefact due to the mixture of different individual patterns, but they appear to be universal for all individuals.Although all individual spatiotemporal patterns are best described by the same distribution functions, the fitted parameters vary from individual to individual.For example, the exponents of the truncated power-law for waiting distribution estimated by the maximum-likelihood method range from γ = 1.6 to γ = 2.5.This indicates that the internal properties encapsulated by the scaling exponent γ are different among individual cows, although their activities appear to be governed by the same dynamics.

Discussion
In this study we have found that under a two-state 'stop-and-move' representation the spatiotemporal pattern of grazing cattle exhibits a hierarchical structure in space and an asymmetric temporal spectrum, which can be described by a composite Brownian walk interspersed with power-law distributed non-moving periods.This finding is in contrast to the patterns observed in human [27] and T-cell [35] mobility where the moving and non-moving states are both characterised by a scaling law.Since detailed statistical characterisation on free-range animal movement based on high-frequency GPS trajectories is still largely missing, this finding can provide new perspectives to our understanding for grazing animals movement and useful leads to the underlying ecological basis of grazing behaviour.
A simple deterministic scenario that can give rise to the observed scaling law in waiting time distribution is that the environment is structured according to the same heterogenous statistics.We can consider that different location r in the landscape is of different quality or resource abundance described a quality function Q( r).If the cattle simply spend their time for feeding on one location proportional to the quality Q( r) at that location, i.e. a 'greedy' strategy, Q( r) would be the observed waiting distribution.
Stochastic processes and spontaneous behaviour can also account for the observed spatiotemporal pattern.Recently, a plausible decision-based queueing process in which the animal executes activities from a stochastic priority list has been used to interpret the scaling law observed in the waiting time of marine predators [18].This model was originally proposed to explain the power-law distributed inter-event time observed in the communication pattern in human dynamics [22].Specifically, the model assumes that the animal performs the two activities waiting and moving with probability x 1 and x 2 = 1−x 1 at a regular basis, where x 1 and x 2 are the priority of the activity drawn from a random distribution p(x).If the animal moves, it changes its context and therefore its likelihood to move or stay also changes.As a result, the priority will be redrawn from the random distribution, representing the change of state due to the movement.This model can generate the power-law distribution in waiting time as well as the exponential distribution in step-size.By introducing a deterministic component to the decision probability, the model can be also tuned to generate different scaling exponent γ accounting for the various scaling phenomena in different species.The model is recast in a dynamic prey-predator environment where the moving probability x 1 can be interpreted as the likelihood of finding a prey in the vicinity.
We can also consider the movement as a two-state point process, in which the probabilities that the animal switches its state are q A (from moving to non-moving) and q B (from non-moving to moving) [38,39].It is well known that the state duration is exponential distributed when the switching probability is constant and independent of time [39].
Recently, it has been suggested that the power-law distributed duration can be attributed to the reinforcement dynamics, such that the switching probability is proportional to the time that animal has spend in its current state, i.e. the longer the animal stays in its current state the less likely it will change it [36,37,38].
Another explanation is to associate the movement pattern with the energy state of the animal using a maximum entropy approach.In this context, each moving and nonmoving activity is associated with a certain amount of energy loss E l or energy gain E g .According to the maximum entropy principle, the distribution of E l and E g over all activity segments should follow a Boltzmann distribution p(E l,g ) ∼ e −E l,g (See Material and methods).The validation of the maximum entropy approach is mainly subject to two conditions: (1) each individual activity is independent and has no influence on others; and (2) the energy intake and expenditure is maintained by two different mechanisms and can be treated as two isolated systems.The first condition is supported by our test on the correlation between consecutive activities, while the second is intuitively understandable.
Following this formulation, it is straightforward to derive that when E g ∝ log τ w and E l ∝ τ m the observed scaling law in waiting time as well as the exponential distribution in trip time can be reproduced.That is to say, the energy intake increases logarithmically as grazing time increases, while the energy expenditure due to moving increases linearly with the moving time or distance.It is interesting to note that the logarithmic energy intake function has been suggested for grazing animals before, and the linear energy expenditure or cost function has been widely observed in many single-mode movements of human transportation activities [41,40].It is well known that energy status can affect animal movement, but a quantitative understanding of their relation is still unclear.
Our proposed maximum entropy approach can potentially fill this gap by establishing a connection between the energy function and the observed mobility pattern, suggesting that the detailed energy intake or expenditure as a function of time in different activities can be inferred from statistical features of the macroscopic mobility patterns such as waiting time or step-length distributions.The conjectured relation can be tested in future experiments by measuring detailed energy intake or consumption using laboratory techniques.

Material and methods
Dataset description.The dataset consists of continuous 0.5-Hz GPS samples for 34 individuals covering an observation period of over 50 hours.We select the data of 31 individuals in which there is no discontinuity in GPS samples and we choose a continuous 30-hour observation window during which the animals were grazing in a confined 600m × 400m rectangular area.The trajectory for each individual cow can be denoted by a sequence L = {p i }, where p i = (x i , y i , t i ) represents a GPS sample with (x i , y i ) being the position coordinates and t i being the timestamp.We use moving average filtering to reduce the noise and smooth the trajectory with a 10 sec moving window, such that Classification of mobility pattern.We define the non-moving segment of a trajectory as a set of consecutive points L w = {p k , p k+1 , . . ., p k+m−1 }, which satisfies the following three conditions: (1) the distance d k,j from the starting point p k to any other point p j of the segment must be smaller than a certain threshold ∆d, i.e. max k<j<k+m d k,j ≤ ∆r; (2) the distance from the starting point p k to the point following the ending point of the segment p k+m must be larger than ∆r, i.e. d k,k+m > ∆r; (3) the time span of the segment must be longer than a certain threshold ∆t, i.e. t k+m − t k > ∆t.In this definition, the first two constraints are made to identify consecutive points that are likely to represent an identical position within in a certain proximity.The third constraint imposes a minimum time span of the non-moving segment that can be tuned to exclude some very-short random activities such as a pause when encountering an obstacle, as well as making the extracted non-moving segments more representative of meaningful activities such as grazing or resting.After extracting the non-moving segments, we simply define the points between two non-moving segments as the moving segments.The approach here is in analogy to the definition of staying points for continuous GPS samples in most spatiotemporal analysis of human mobility [27,42].The value of ∆t is suggested to be 2 − 3 mins [33].
Maximum entropy principle.The maximum entropy principle originates from statistical mechanics, which assumes that the configuration of microscopic states of a complex system (e.g. the energy of each particle) leading to the macroscopic observation is the one that maximise the entropy of the system.Suppose the system consists of N noninteracting particles and has a total energy U , such that N = n i and U = n i E i where n i denotes the number of particles at a specific energy state E i .Then the ensemble that represents all possible configurations of the system is called the canonical ensemble, and the probability p(E) that a particle has a specific energy state E is denoted by P (E) ∝ e −E/ Ē .Here we assume that the alternate moving and non-moving activities in the mobility pattern operate in two independent systems, while the individ-

Figure 1 :
Figure 1: A diagram of the spatiotemporal pattern under the two-state 'stopand-move' representation.(a) The temporal spectrum of activities illustrated in a spike train.Colour segments on the time-axis represent alternating waiting (white) and moving (red) activities over time.(b) The raw trajectory of an individual cow before processing.(c) The spatial pattern extracted from the raw trajectory in panel (b) can be projected as a transition graph, where the waiting locations for non-moving segments are represented by red dots and the trips are represented by blue solid lines.

Table 4 :
The AIC weights for the individual trip length distribution Cow ID Powerlaw Truncated powerlaw Exponential Hybrid exponential

Table 6 :
The AIC weights for the individual waiting time distribution