Modelling and tracking the flight dynamics of flocking pigeons based on real GPS data (small flock)

Ecological Modelling 344 (2017) 62–72

Contents lists available at ScienceDirect

Ecological Modelling journal homepage: www.elsevier.com/locate/ecolmodel

Modelling and tracking the flight dynamics of flocking pigeons based on real GPS data (small flock) Ayham Zaitouny a,b,∗ , Thomas Stemler a , Michael Small a,b a b

School of Mathematical and Statistical Science, The University of Western Australia, Crawley, Western Australia, Australia Commonwealth Scientific and Industrial Research Organisation (CSIRO), Mineral Resources, Kensington, Western Australia, Australia

a r t i c l e

i n f o

Article history: Received 23 June 2016 Received in revised form 18 November 2016 Accepted 19 November 2016 Keywords: Animal behaviour Animal movement Flocking Nonlinear model Tracking GPS

a b s t r a c t Improvements in our ability to collect and process real-time data from animal behaviour have lead to increasing interest in the quantitative analysis of collective motion in various species. This motion can be characterised as movement of a group of animals that interact among themselves and with the surrounding environment to preserve the cohesive form of the group. Examples of this motion include a flock of birds, a fish school and a group of deer or sheep. The problem becomes particularly challenging as the group dynamics and the individual dynamics need to be disentangled for biological analysis. Within this framework we describe a new mathematical technique that may be applied in order to understand animals group movement and behaviour. In particular, we model and track the behaviour of a small flock of pigeons and verify using real data that the dynamics of each individual are a combination of two different mechanisms (local and global). © 2016 Elsevier B.V. All rights reserved.

1. Introduction Understanding animal movement is a challenging mathematical exercise as this movement is driven by complex requirements of the individual, such as interactions between animals or even interactions with different animal groups. The specific question of understanding animal behaviour has been approached using different methods. The research so far falls into two main categories: The first looks at the problem from the point of view of biology. The problem of this approach is that it requires a large set of biologically meaningful parameters (Ballerini et al., 2008; Hemelrijk and Hildenbrandt, 2011), and, of course, is reliant on the correctness and completeness of the underlying biological model. The second approach deals with the problem statistically, the two main models are Hierarchical state-space model (Jonsen et al., 2006; Patterson et al., 2008) and Hidden Markov model (Franke et al., 2004; Patterson et al., 2009). However, in the last few years new approaches have been proposed to unravel this problem. These new approaches use ideas of complex networks and nonlinear dynamical systems (Walker et al., 2010; Kattas et al., 2012a,c,b; Xu et al., 0261).

∗ Corresponding author at: School of Mathematical and Statistical Science, The University of Western Australia, Crawley, Western Australia, Australia. E-mail addresses: [email protected], [email protected] (A. Zaitouny). http://dx.doi.org/10.1016/j.ecolmodel.2016.11.010 0304-3800/© 2016 Elsevier B.V. All rights reserved.

In this paper we use a nonlinear dynamical model to solve the problem of understanding the collective motion of animals, specifically, the behaviour of a flock of pigeons. We verify our work using GPS real data provided by a study by Nagy et al. (2010). Numerical investigations on this data set allowed us to explore some essential features that derive the flock behaviour. Consequently, we initially used a simple model that described the flock dynamics accurately – with the exception of isolated singular points in the data. These singularities arise because of atypical behaviour of particular pigeons together with the artifacts of using real data set. Hence, we improved our model by introducing a hyperbolic tangent indicator and ultimately overcame these singularities. The purpose of this framework is to demonstrate, based on real data, that the behaviour of each individual pigeon in the flock is generated by a combination of local and global dynamics and to distinguish between periods when each mechanism dominates the total behaviour. While we have grounded our work with the specific example of pigeon flock dynamics, the methods are entirely generic and, we expect, applicable to a wide range of similar problems. 2. Flock dynamics overview The data set provides dynamic information (position, velocity and acceleration) for a small flock of pigeons. It includes two different behaviours of the flock (free-flight and homing flight). Since the data we are using is a GPS data set, it was necessary for us to verify its reliability. Therefore we analysed and filtered the provided data

A. Zaitouny et al. / Ecological Modelling 344 (2017) 62–72

63

Fig. 1. Flock dynamics overview.

using a new tracking methodology based on the idea of shadowing filter (Stemler and Judd, 2009; Zaitouny et al., 2016b,a; Judd, 2015).1 Using our tracking technique we were able to filter the position trajectory for each individual pigeon, consequently using our tracking model we estimated the corresponding velocities and accelerations based on Newtonian mechanics. The results have been introduced and extensively describe in our previous paper (Zaitouny et al., 2016c). The data analysis we have made indicate that our results are more robust, accurate and complete using the filtered data set than the original one.2 In our investigations we use a free flight file from the data set. This file includes eight pigeons (A-F-G-H-I-J-K-L). In order to have an overview on the behaviour of the flock we plot the projections of the positions on each plane. A previous study by the authors Zaitouny et al. (2016c) shows that there are two distinct periods “Flying and Nonflying”, besides the atypical behaviour of some pigeons (pigeon H), as depicted in Fig. 1a. In addition, if we look more closely at the flying period (see Fig. 1b and in particular the period between datum 8500 and 10,000) we notice clearly that some pigeons leave or join the flock randomly. Therefore, it is necessary to define an indicator to determine the status of each pigeon at each time point. As has been shown in our previous data analysis (Zaitouny et al., 2016c), the criterion we will use for flying or nonflying behaviour is the velocity magnitude. In Fig. 1c, which shows a plot of the velocity magnitude of pigeon A along the whole time series, we can see clearly the two different behaviour, flying and nonflying. Moreover we notice

1

A MATLAB code of our tracking approach is provided in Zaitouny (2016). This is due, in a large part, to the fact that our methods requires an estimate of acceleration, and the filtering technique describe in Zaitouny et al. (2016c) was able to provide this in a more robust fashion. 2

that a pigeon is flying when its velocity magnitude is greater than 4ms−1 , therefore we can define our indicator as follows:

 I i (t) =

1,

if vi (t) > 4 ms−1

0,

otherwise

(1)

where i refers to pigeon i, t to the time point, and v to the velocity magnitude. Furthermore, based on our investigations (Zaitouny et al., 2016c) regarding the position projections on the different directions, we found that the variation in the x and y-directions are much higher than the variation in the z-direction. Table 1 confirms this observation, where one can observe that the velocity variation in the z-direction is only 4 ms−1 , while in the x and y-directions it is 30 ms−1 . Similarly, we can observe the different behaviour in the x and y accelerations compared to the z acceleration. Consequently, we infer that the pigeons are flying approximately on the xy-plane when they are flying as a cohesive flock. Of course, this is a peculiarity of this particular system. In genuinely three-dimensional systems, we could integrate all available information for all degrees of freedom.

Table 1 Free flight: ranges of x, y and z components of position, velocity and acceleration. Position (m) Range Velocity (ms−1 ) Range Acceleration (ms−2 ) Range

x 2920 → 3060 x˙ −15 → 15 x¨ −8 → 8

y 1760 → 1940 y˙ −15 → 15 y¨ −8 → 8

z 110 → 145 z˙ −2 → 2 z¨ −2 → 2

64

A. Zaitouny et al. / Ecological Modelling 344 (2017) 62–72

3. Flock dynamics: initial simple model Initially we used a simple model (only using the flying–nonflying indicator) to understand the flock behaviour, this model matches well the real data set besides of some singularities caused by some atypical behaviours of some pigeons (see below). However, in the next section we introduce an improved version of our model which enhances our results and overcome these singularities. The essential idea of our model is that the dynamics of each pigeon in the flock are derived by the combination of two different mechanisms. The first results from the global behaviour of the flock, while the other is produced from the local interactions among the individuals inside the flock. Fig. 2a illustrates this idea. In essence, one can see how the flock moves as a single compact entity and this forms the global mechanism. On the other hand the individual pigeons move and change their relative positions within the flock. This describes the local mechanism which is caused by the interactions among the individuals. Consequently, the acceleration of each pigeon at each time can be defined as: ai (t) = ag (t) + ail (t),

(2)

where a is the total acceleration, ag is the global acceleration, al is the local acceleration, i = 1, 2 . . . 8 refers to different pigeons and t is the time. Our model is based on two essential assumptions: • The flock we are modelling is small in size (only 8 pigeons) and each individual inside the flock interacts with all other pigeons in the flock (fully connected network). This assumption is made because our model is driven from real data, and the data set we have only provides information about a small numbers of pigeons.3 • Our first assumption leads us to the second one, that is, since the flock we are dealing with is small, no splitting will occur in the flock, that is, we never have two flocks with two distinct mass centres. This is a reasonable assumption as we model a flock of only few pigeons. However, this assumption does not contradict with the fact that some pigeons join or leave the flock randomly as shown above in the real data set. 3.1. Global dynamics The global dynamics may be interpreted as the alignment bahaviour with the neighbours, following the flock leader (if such an individual exists) or an external force acting on the flock caused by some pigeon requirements, for example searching for food or reaching their home in the case of homing flight. However, in our situation for a free circular flight since we have real GPS data, the global behaviour can be estimated from the dynamics of the centre of mass of the flock: ag (t) = amc (t),

(3)

where amc is the acceleration of the flock’s centre of mass. To find this acceleration, first we need to find the trajectory of the mass centre positions as follows:

 i I (t)P i (t)) i  Pmc (t) = i (

(

i

I (t))

(4)

3 Here we take note one simplification in the behaviour of small flocks. With larger flocks we cannot so easily consider a fully connected network between the individuals, rather that each pigeon is assumed to be influenced only by the closest neighbours.

where Pmc (t) ∈ R3 refers to the mass centre position at time t, Pi (t) ∈ R3 is the position of pigeon i at time t and I is the “flying–nonflying” indicator function defined in Eq. (1). Along this trajectory uncertainties will occur because of the random leaving or joining by some pigeons. We managed to overcome this problem by using the tracking filter methodology, which was introduced earlier by us Zaitouny et al. (2016a,b), Judd (2015) and has been implemented successfully to track a moving object (Zaitouny et al., 2016a,c). In this study we use this methodology to filter the estimated trajectory of the centre of mass of the flock in order to obtain more reliable results. Fig. 2b shows the mass centre trajectory estimated by Eq. (4) (in blue) as well as the estimate obtained by the aforementioned tracking filter technique (in red). One can see how the initial mass centre position trajectory includes discontinuous jumps. However, these jumps have been treated sufficiently well by our tracking filter (for example, between datum 600 and 800 in Fig. 2b where during this jump one pigeon has joined the flock). Additionally, our tracking technique (Zaitouny et al., 2016b,c; Judd, 2015) enables us to estimate the corresponding acceleration of the centre of mass based on its trajectory which is our main aim as the acceleration is the basic factor influencing the global dynamics. The x-coordinate of the mass centre acceleration, which is presented in Fig. 2d, will be considered as the global acceleration of each individual pigeon in the flock according to our model Eq. (3). In addition, Fig. 2c shows the filtered position trajectory of the centre of mass in three-dimensional space. It is an approximately elliptical4 path which we henceforth consider as the global trajectory of each pigeon in the flock. 3.2. Local dynamics: attraction and repulsion As mentioned above, the local dynamics are derived by the interactions among the individuals inside the flock. In our model, the essential factor that captures the intra-flock behaviour is based on the distance between individuals. That is, we propose an attraction–repulsion model to describe the local mechanism in order to avoid any collision between the individuals as well as preserve synchronously the cohesive form. In Fig. 2e we illustrate schematically our idea about the local dynamics: At any time point, each pigeon inside the flock is affected by the other pigeons in the flock by an attractive or a repulsive force based on the distance between them. In other words, at time t each pigeon is affected by a resultant force formed from the sum of attraction and repulsion forces. To estimate these local forces we use a potential function (Helbing, 2001). Since the main factor involved in the local dynamics is the pairwise distance, we propose a potential driven by this distance factor. An overview of the cumulative pairwise distance distribution between the pigeons along the whole time series (Fig. 2f) leads us to propose a potential of the form shown in Fig. 3a. We notice that this potential distinguishes between two different regions (repulsion and attraction) separated by an optimal distance. Based on numerical investigations of different potential functions, we found an appropriate potential formed by a combination of reciprocal and exponential functions, defined by: V (D) = n

1

D

−

1 m



e−D/k ,

D>0

(5)

where D is the distance variable, n, m and k are positive real parameters. The exponential term provides decaying attraction as a function of distance, while the hyperbolic coefficient (1/D − 1/m) provides an optimal separation. Studying this potential function shows us that it has one positive turning point, which occurs when its derivative equals zero. This

4

See Appendix A.

A. Zaitouny et al. / Ecological Modelling 344 (2017) 62–72

65

Fig. 2. Illustration of the model ideas and estimations of the mass centre dynamics (global dynamics).

condition is satisfied when D2 − mD − mk = 0, which means that the turning point is given by:

D=

m+



m2 + 4mk 2

(6)

Obviously one of these two turning points is positive while the other is negative. Since we only consider positive distances, the negative turning point will be ignored. In fact, this critical value of distance, D* , is the desired distance that a pigeon prefers in order to be sufficiently far from other pigeons. In other words, it is the critical distance where the force will be transformed from attractive to repulsive or vice versa. Using this potential, the force acting on pigeon i affected by pigeon j is defined as the negative of the potential function gradient: F(Dij ) = −∇ V (Dij ),

(7)

 where Dij =

(xj − xi )2 + (yj − yi )2 + (zj − zi )2 is the distance

between pigeons i and j. Obviously, this force is acting on the straight line in three-dimensional space passing through pigeons i and j. The direction will be determined depending on whether the force is attractive or repulsive. Consequently, the resultant local force acting on pigeon i at time t is defined as: F i (t) = I i (t) ∗



− I j (t) ∗ V  (Dij (t))

j= / i

∗

(xj (t) − xi (t), yj (t) − yi (t), zj (t) − zi (t)) Dij (t)

(8)

where I is the flying–nonflying indicator defined in Eq. (1). In the resultant force formula, this indicator is used to obtain zero force for nonflying pigeons and to exclude the effect of nonflying pigeons on others. Dij (t) is the distance between pigeons i and j at time t. xi (t), yi (t), zi (t) are the position components of pigeon i at time t.

66

A. Zaitouny et al. / Ecological Modelling 344 (2017) 62–72

Since “Force=mass × acceleration” and since we use a scalar parameter n in our potential, this local force produces the local acceleration. This is reasonable as the mass factor in the forceacceleration relationship is implicitly represented by the scaling factor n of our potential function. Therefore: ail (t) = F i (t).

(9)

3.3. Model optimisation and results As we have seen above the local dynamics model includes three parameters n, m and k, therefore these parameters need to be optimized in order to obtain results as close as possible to reality. Since we are interested in a particular real application, we will optimise and verify our model with that data. Therefore, we examine the difference between our model accelerations and the real data accelerations, and try to minimize this difference as much as possible. In practice, we used the filtered acceleration instead of the provided GPS data acceleration as our analysis indicates that the data we found through the filtering are more reliable than the original one (this has been explained in detail in our paper (Zaitouny et al., 2016c)). Accordingly, the minimization problem may be formulated as:



min (

i

aidata (t) − (aig (t) + ail (t))2 ).

(10)

t

In practice, we used the optimisation tool lsqcurvefit from MATLAB library to perform this minimisation. In general one might want to use a bounded norm instead of the Euclidean when dealing with real data. For our particular application the Euclidean norm was sufficient. When applying this minimisation algorithm, we are faced with some expected problems which we managed to overcome (these problems and the minimisation approach are explained in detail in Appendix B). We implemented the minimisation algorithm described above for 1500 points (points between 8500 to 10,000 from the data file), which roughly corresponds to 5 min of continuous flying. The outcomes of the successful minimisation approaches provided an approximation of the optimal values of our three parameters (a ≈ 0.23, m ≈ 0.0013, k ≈ 999), giving a sharp and deep potential function as shown in Fig. 3a. Accordingly, using Eq. (6) we can conclude that the optimal distance where a pigeon tries to stay far from others while they are flying is: D ≈ 1.14 m, which is a reasonable result. That is, in a flock when two pigeons are less than 1.14 m apart from each other, the local interaction between them will be a repulsive force to avoid any collision. While if the distance exceeds 1.14 m, then they attract each other trying to preserve the cohesive construction of the flock. Fig. 3b shows a comparison between the local acceleration derived from the observational data (in blue) and the local acceleration estimated using our model Eq. (8) (in red); this comparison has been done for pigeon A. The local acceleration has been calculated from the observational data acceleration after subtracting the estimated global acceleration. These model estimations match well with the real data acceleration along the whole time series, except at the points marked with green segments. This is not surprising as these points refer to the time when some pigeon leaves or joins the flock randomly, therefore singularities are expected. In Fig. 3c we present a comparison of the total acceleration (global + local) of pigeon A. The acceleration denoted in light blue is the original acceleration provided from the data set. The acceleration in dark blue is the filtered acceleration we estimated using our tracking technique in order to avoid unrealistic sharp changes. One can clearly see how this method improves and smooths out the acceleration peaks and

makes them more realistic (Zaitouny et al., 2016b,c). The red depicts our model estimations (global + local acceleration), which are very close to the filtered data estimates. Fig. 3d and e shows similar comparisons for pigeon G. This pigeon joined the flock after about 400 points (i.e. more than 1 min). Again, we notice that our model estimations match quite well with the observational data, except for the singularities we mentioned above caused by leaving or joining pigeons (Fig. 3d). Moreover, Fig. 3e shows another mismatch at the beginning of the period when pigeon G started flying (about 100 points). This is also expected, since pigeon G needed about 20 s to join the flock. In fact, during this time period, pigeon G was trying to catch up with the rest of the flock, therefore it had different behaviour than the flock’s global behaviour. In summary, we can say that our model provides satisfactory results, however, two problems have been highlighted from these calculations: • Big jumps appear after calculating the mass centre’s trajectory because of the random leaving or joining pigeons, leading to singularities in estimating the global acceleration (the centre of mass acceleration). Consequently, corresponding singularities appear in the local acceleration (highlighted by green segments and ellipse in Fig. 3b and d). This issue occurs because of the flying–nonflying condition (tough threshold). • When a pigeon attempts to join the flock, singularities of a different nature occur during the first few seconds. Just after the start of flight time, and before catching up with the flock, during this period the pigeon is attracted to the collective but is not following the global behaviour yet (highlighted by a purple rectangle in Fig. 3d). To overcome these problems, we propose a new smoother indicator to improve our model results in the following section. 4. Improving the model As mentioned in the previous section, the flying–nonflying criterion was too sharp, therefore two significant problems have arisen due to the random and sudden joining or leaving of the flock. To further improve our model, we propose a much smoother condition using the hyperbolic tangent (tanh) function. The aim of using this condition is to consider a gradual joining or leaving the flock instead of the strident flying–nonflying indicator. The tanh indicator is defined after some appropriate transformations as follows: Q i (t) =

1 2



1 − tanh

1

K1



(Di (t) − K2 )

,

(11)

where Di (t) is the distance between pigeon i and the centre of mass at time t. K1 , K2 ∈ R are two parameters resulting from the transformations applied to the tanh function. K1 describes the dilation in the x-direction which affects the slope of the indicator while K2 determines the translation in the same direction. These two parameters can be roughly determined by looking at the distances of the pigeons from the centre of mass along the flying interval (Fig. 4a). This figure shows the distances from the centre of mass when most of the pigeons have already joined the flock. It can be shown that the distances vary between 0 and 12 m with the majority distributed between 0 and 5 m. Consequently, based on some numerical investigations we can choose the transformations parameters as K1 = 5 and K2 = 15 in order to obtain a consistent indicator with these observations. Hence, we choose the tanh indicator shown in Fig. 4b. This indicator function can be interpreted as a proportional percentage of how much a pigeon is involved in the flock, in other words, how much a pigeon is influenced by the global behaviour

A. Zaitouny et al. / Ecological Modelling 344 (2017) 62–72

67

Fig. 3. Optimal potential and model outcomes on pigeons A and G.

of the collective. From this figure we can see that pigeons with distances about 5 m or less from the mass centre are involved in the flock with percentage more than Q = 95%, however, for distances between 5 m to 10 m the Q indicator will gradually decrease from approximately 95% to 85%, then for larger distances this indicator decreases gradually to reach almost 0 when the distance from the centre of mass is greater than 28 m (i.e. the pigeon has no impact on the flock anymore). This proposed smoother indicator will be used, as shown in the following subsections, to estimate the global and local accelerations, aiming to improve our model outcomes and overcome the previous problems. 4.1. Global dynamics Global dynamics are driven by the behaviour of the centre of mass. However, here we will estimate the mass centre trajectory using the smoother Qi indicator so as to avoid the singularities we

faced previously due to the joining or leaving the flock by some pigeon. The idea of how to use this indicator is to consider the region around the centre of mass as an effectiveness diminishing region, that is, the farther the pigeons move away from the mass centre, the less impact they have. In other words, the impact of some pigeon on how the mass centre’s position is modified in time is proportional to the difference between its position vector and the mass centre’s position vector by the indicator factor Qi . Fig. 4c illustrates schematically this idea. Thus, we estimate the mass centre trajectory using the Qi indicator as follows:

 Pmc (t) =

(

i

Q i (t)(P i (t) − Pmc (t − 1))



(

i

Q i (t))

+ Pmc (t − 1).

(12)

To better understand this equation, we show a simple representative case of the impact of two pigeons on the mass centre dynamics during a single step (Fig. 4d). In this figure, R is a reference state. In the equation, the term Pi (t) − Pmc (t − 1) is the

68

A. Zaitouny et al. / Ecological Modelling 344 (2017) 62–72

Fig. 4. Illustrating the ideas of how we improve our model.

vector of the impact direction of pigeon i on the centre of mass at time t (black arrows in Fig. 4d). We can observe how pigeon Pj has a smaller impact on the mass centre than pigeon Pi (blue arrows) as pigeon Pj is farther away, thus, it has a smaller Q indicator.

Using this new indicator, we estimate the mass centre trajectory, and obtain satisfactory results. In fact, with this smoother threshold we managed to overcome the jumps in the mass centre position trajectory caused by some pigeon randomly joining or leaving the flock. Fig. 4e shows a comparison of the mass centre

A. Zaitouny et al. / Ecological Modelling 344 (2017) 62–72

69

position trajectory estimated by Eqs. (4) and (12). It is evident how using the tough flying–nonflying condition results in large jumps when some pigeon joins or leaves the flock. However, replacing this tough condition with the smoother Qi threshold provides us with better results and helps avoid these kind of unwanted jumps. Consequently, as expected, estimating the global acceleration based on this alternative mass centre position trajectory formulation (using Qi ) renders improved results. Using the Qi indicator, the global acceleration for some pigeon i at time t is defined as follows:

the flock on other pigeons far from them. In fact, in such a situation this pigeon will be strongly attracted by the pigeons inside the flock to join the cohesive collective. Furthermore, we give a lower boundary to the Qi indicator in order to avoid dividing by very small values which results in irregularly large estimations of the local acceleration.

aig (t) = Q i (t) ∗ amc (t).

In order to obtain the optimal potential, the same minimisation procedure described in Section 3.3 has been run for our advanced model. The outcomes are consistent with our previous ones and the optimal parameters are almost the same. The local and global dynamics have been estimated to all pigeons and the outcomes are consistent with the observed data set. Fig. 5a shows our improved model estimations of the local and total (local+global) accelerations of some pigeons. In comparison to the observed data, we can see that our estimations are close to the actual data. Our calculations demonstrate a marked improvement, particularly at the initial points of the interval and around the time points when other pigeons join or leave the flock (compare with Fig. 3). Here it is evident how we manage to overcome the singularities of the first type. Although an irregular spike has appeared in the estimations of the local acceleration around the point 300, this can be explained as a strong repulsive interaction with pigeon I due to a very short distance6 . In fact, in order to avoid any collision our model shows stronger repulsive forces against such very short distances. However, it seems that the pigeons, and areodynamic forces, are smarter than the model and they can react with smaller forces. Furthermore, in order to show how our ultimate model effectively overcomes the singularities of the second type (those appear at the periods when some pigeon is joining the flock), we show the estimations of the local and total accelerations of pigeons G and J in Fig. 5b and c, respectively. One can clearly see that our estimations of both the local and total accelerations match quite well with the data set along the entire time interval. Comparing with our previous results in Fig. 3d and e, we can see that our advanced model enhanced significantly our results particularly during the time period when the pigeons were far away from the flock and trying to join the rest of the collective. In other words, our model now can capture the attractive dynamics of a pigeon when it is trying to join the flock. Moreover, in Fig. 5a, c just after the time point 800, it seems that pigeons A and J behave atypically and spikes have appeared in

(13)

In other words, the global acceleration is now not exactly the mass centre acceleration as before, but rather it is proportional to the mass centre acceleration by the factor Qi . This is more sensible physically and biologically as when some pigeon is trying to join the flock, it will not follow the flock’s behaviour instantly. It will first be attracted to the flock, then gradually reach the flock and will finally follow the global behaviour once close enough (Fig. 4g).5 Fortunately, using this modified global acceleration to estimate the local acceleration, enables us to remove the singularities that have appeared previously in the local acceleration. Fig. 4f shows a comparison of the local acceleration of pigeon A estimated by using the two different methods of extracting the global acceleration. It is easy to notice how the results of using the Qi indicator (in red) succeed in avoiding these singularities due to random joining/leaving. On the other hand, as shown before when the Ii indicator is used (in blue), the singularities are large at these points. However, as can be seen in Fig. 4f, apart from these points, the two estimations using the different indicators are almost identical. Which means we have successfully overcome the first problem. 4.2. Local dynamics In order to further improve our model’s performance, we try to overcome the second aforementioned problem which appears in the periods when some pigeon is trying to join the flock. We aim to capture the behaviour of a pigeon which is far from the flock and trying to join the collective. The idea here is when a pigeon is away from the flock, the local dynamics are more significant than the global one. While the pigeon moves towards the flock, the global dynamics will gradually take over the local one, Fig. 4g shows a simple representative case of this idea. Accordingly, while the global dynamics of a pigeon are proportional to the Qi indicator, the local force is inversely proportion to it. Hence, we modify our local acceleration using the hyperbolic tangent indicator Qi (Eq. (11)), as follows:

F i (t) =

4.3. Results

⎧ i (xj (t) − xi (t), yj (t) − yi (t), zj (t) − zi (t)) I (t)  ⎪ ⎪ − I j (t) ∗ Q¯ ji (t) ∗ V  (Dij (t)) ∗ ; Q i (t) ≤ 0.15 ∗ ⎪ ⎨ 0.15 Dij (t) j= / i

(xj (t) − xi (t), yj (t) − yi (t), zj (t) − zi (t)) I i (t)  ⎪ ⎪ − I j (t) ∗ Q¯ ji (t) ∗ V  (Dij (t)) ∗ ; otherwise. ⎪ ⎩ Q i (t) ∗ Dij (t)

(14)

j= / i

Same variables and parameters as before, where Q¯ ji is the same tanh indicator defined in Eq. (11) but with a variable such as Dij , the distance between pigeons Pi and Pj , and parameters as K1 = 5, K2 = 15 + meanallj (Dij (t)). The purpose of using this new indicator in estimating the local forces is to construct another impact diminishing region but now around the pigeon under consideration. Which relatively means that the closer a pigeon is to another pigeon, the more impact it has on its local dynamics. The reason why we defined K2 as above is to avoid excluding the impact of the pigeons inside

5

We will discuss this further in the local dynamics subsection.

the data of the local acceleration, however, these spikes have been captured by our model estimations. To explain what is happening over this period we looked at the distances of these pigeons from the centre of mass, Fig. 6a shows the distances of pigeons A,G and J from the centre of mass as well as the corresponding Qi indicator. One can see that the distances of pigeons A and J from the mass centre around this time are about 20 m which imply smaller Qi indicators (as shown), which means these pigeons split away from the flock at this time and were attracted to join the flock again. This

6

See Fig. (C.12) in the Appendix C

70

A. Zaitouny et al. / Ecological Modelling 344 (2017) 62–72

Fig. 6. Pigeons A, G and J: distances from the centre of mass and the corresponding Qi indicator with the intervals where global and local dynamics dominate the total behaviour.

Fig. 5. Pigeons A, G and J: our model results of the local acceleration and the total acceleration compared with the filtered observed data.

explains many other such peaks along the time series of the local acceleration. In summary, it is evident how our advanced model has overcome the two problems we faced in our initial model. Our advanced model captures the dynamics and the behaviour of the pigeons along the whole time series even though they fly continuously within the flock or just decide to join or leave the group. In fact, the results support our assumption that the dynamics of a pigeon are composed of two distinct mechanisms (local and global). Moreover, our model reveals an interesting conclusion, that is, when a pigeon is flying one of these two mechanisms dominates the total behaviour. In fact, when a pigeon flies within the flock then the global dynamics dominate the total behaviour. On the other hand, when the pigeon is far away from the flock then the local dynamics are the dominant behaviour. In order to highlight this important result, we compare the magnitude of the global acceleration with

the magnitude of the local acceleration along the entire time series of pigeons A, G and J. In practice, to illustrate the dominant dynamics of the total behaviour, we estimated log(ag (t) / al (t) ) along the whole interval, Fig. 6b shows these calculations for pigeons A,G and J, when this amount is positive then the global dynamics dominate the pigeon’s behaviour, while when it is negative then the local dynamics control the behaviour. For pigeon A the global dynamics dominate its behaviour most of the time except at the beginning of the time interval and just after the point 800 where the local dynamics take over the global. Pigeons G and J, when they started flying and tried to join the flock, their local dynamics were controlling their behaviours, while when they were within the flock the global dynamics dominated their behaviour most of the time. Pigeon J, similar to Pigeon A, shows around the point 800 a period where the local dynamics took over the global dynamics. This is an important result, that is, previous studies (Zhang et al., 2014; Small and Xu, 2012) have modelled the flock behaviour using alignment condition (can be considered as the global dynamics in our case) combined with attraction–repulsion condition. According to our result it is important to take into account the dominance of either behaviour at a given time. 5. Conclusion In this study, using real GPS data, we managed to model and track the behaviour of a small flock of pigeons. Our model proves that the dynamics of a pigeon inside a flock are composed of two essential mechanisms: The first is the global mechanism (collective level) and this mechanism can be understood as the behaviour of the alignment with neighbours inside the flock or following a leader of the flock if it exists. The second mechanism is the local one

A. Zaitouny et al. / Ecological Modelling 344 (2017) 62–72

(individual level) which is caused by the intra-interactions among the individuals within the flock. Moreover, Our model led us to an important result, that is, according to the position of the pigeon within the flock one of these two components will dominate the total behaviour. In fact, this study provides a strong understanding of how a small flock preserves its cohesive construction, and how the individuals within the collective interact among each other. Our model could be a robust base for many real-life applications, such as building or designing a group of flying objects (automated drones, for example) where we can only control one of them (the leader) and the rest could follow the leader and preserve the cohesive flock without any collision.

71

results. Therefore, we implemented our minimisation on the single pigeons separately along the whole time interval under consideration. In practice, the minimisation procedure worked significantly well with those pigeons flying most of the time interval. However, the minimisation algorithm failed with the other pigeons which join or leave the flock randomly, in these cases it seems that the procedure stuck in a local minimum instead. When the minimisation algorithm worked, the results we gained were roughly similar and successfully worked even when we implemented the same estimated potential function on the other pigeons that leave or join the flock. Appendix C. Results for all pigeons

Acknowledgements The authors wish to thank Kevin Judd for valuable discussions and comments. They also would like to thank, M. Nagy, Zs. Akos, D. Biro and T. Vicsek, the authors of Nagy et al. (2010) for providing the data. Finally the authors thank UWA for funding Mr. Zaitouny to do his Ph.D.

Here we show our improved model results for all pigeons from the data set. We show comparisons for both local and total

Appendix A. Elliptical mass centre trajectory Initially one might say that the flock is flying in circular paths. Hence, we tried to estimate the radius of the mass centre trajectory (Fig. 2c). We performed this calculation by first estimating in the xand y-directions the differences between each pair of consecutive extrema and then by obtaining the means of these differences. We found that the radius of the flock trajectory in the x-direction is about 65 m while the radius in the y-direction is approximately 87 m. Based on this result, we can say that the flock is not flying along a circular trajectory, but rather along an elliptical path.

Fig. C.7. Pigeon A: our model results of the Local acceleration and the total acceleration compared with the filtered observed data.

Appendix B. Optimisation approach As mentioned in the manuscript, we are faced with some expected problems when applying the minimisation algorithm. The first issue appears because our problem is three-dimensional. Therefore, when we run the minimization approach including the z-component, undesired results were obtained. In fact, this is consistent with the conclusion we reached from the data analysis we performed in a previous study (Zaitouny et al., 2016c). There we showed that when the pigeons are flying as a flock, they fly approximately on a two-dimensional plane. This implies that the dynamics of the flock on the xy-plane are totally different from the dynamics in the z-direction. Indeed, physically the dynamics in the z-direction are explained by aerodynamics laws (gravity and the force of the wings). As a result, we exclude the z-component completely from our model algorithm, in other words reducing our problem from a three-dimensional to a two-dimensional one. Unfortunately, while this solves the previous issue, another appears. This second issue is raised due to the exclusion of the z-component from the local force formula given by Eq. (8). This caused uncertainties throughout our calculations and messed up the minimization. These uncertainties occurred when two pigeons are exactly or almost above each other, which appears as a collision in the two-dimensional projection, leading to a very large repulsion force. This is obviously an unrealistic situation in three-dimensional space. To overcome this problem, we estimate the local force in Eq. (8) as a three-dimensional value, but when we perform the minimisation in Eq. (10), we only do so for the x and y coordinates. This improved the results but did not resolve the problem. The reason for this is that when we did the minimisation over all pigeons in the flock, some uncertainties occurred because of the random leaving and joining of some pigeons. These uncertainties decreased the quality procedure of the

Fig. C.8. Pigeon G: our model results of the Local acceleration and the total acceleration compared with the filtered observed data.

Fig. C.9. Pigeon J: our model results of the Local acceleration and the total acceleration compared with the filtered observed data.

72

A. Zaitouny et al. / Ecological Modelling 344 (2017) 62–72

Fig. C.10. Pigeon F: our model results of the Local acceleration and the total acceleration compared with the filtered observed data.

Fig. C.14. Pigeon L: our model results of the Local acceleration and the total acceleration compared with the filtered observed data.

dynamics for all pigeons. We can see how our model shows high satisfactory results for all pigeons in the flock, including those pigeons which behave atypically by joining or leaving the flock randomly. See Figs. C.7–C.14. References

Fig. C.11. Pigeon H: our model results of the Local acceleration and the total acceleration compared with the filtered observed data.

Fig. C.12. Pigeon I: our model results of the Local acceleration and the total acceleration compared with the filtered observed data.

Fig. C.13. Pigeon K: our model results of the Local acceleration and the total acceleration compared with the filtered observed data.

Ballerini, M., Cabibbo, N., Candelier, R., Cavagna, A., Cisbani, E., Giardina, I., Orlandi, A., Parisi, G., Procaccini, A., Viale, M., Zdravkovic, V., 2008. Empirical investigation of starling flocks: a benchmark study in collective animal behaviour. Anim. Behav. 76, 201–215. Franke, A., Caelli, T., Hudson, R.J., 2004. Analysis of movements and behavior of Caribou – Rangifer tarandus – using hidden Markov models. Ecol. Model. 173, 259–270. Helbing, D., 2001. Traffic and related self-driven many-particle systems. Rev. Mod. Phys. 73, 1067–1141. Hemelrijk, C.K., Hildenbrandt, H., 2011. Some causes of the variable shape of flocks of birds. PLoS ONE, 6. Jonsen, I.D., Myers, R.A., James, M.C., 2006. Robust hierarchical state-space models reveal Diel variation in travel rates of migrating leatherback turtles. J. Anim. Ecol. 75, 1046–1157. Judd, K., 2015. Tracking an Object with Unknown Accelerations Using a Shadowing Filter. arXiv:1502.07743 [cs.SY]. Kattas, G.D., Prez-Barbera, F.J., Small, M., Xu, X.-K., Walker, D.M., 2012a. Unraveling complexity in interspecies interaction through nonlinear dynamical models. Acta Ethol., 1–10, http://dx.doi.org/10.1007/s10211-012-0134-0. Kattas, G.D., Xu, X.-K., Small, M., 2012b. Dynamical modeling of collective behavior from pigeon flight data: flock cohesion and dispersion. PLoS Comput. Biol., 8. Kattas, G.D., Xu, X.-K., Small, M., 2012c. Generating self-organizing collective behavior using separation dynamics from experimental data. Chaos, 22. Nagy, M., Akos, Z., Biro, D., Vicsek, T., 2010. Hierarchical group dynamics in pigeon flocks. Nature 464, 890–893. Patterson, T.A., Basson, M., Bravington, M.W., Gunn, J.S., 2009. Classifying movements behavior in relation to environmental conditions using hidden Markov models. J. Anim. Ecol. 78, 1113–1123. Patterson, T.A., Thomas, L., Wilcox, C., Ovaskainen, O., Matthiopoulos, J., 2008. State space models of individual animal movement. Trends Ecol. Evol. 23 (2), 87–94, http://www.sciencedirect.com/science/article/pii/S0169534707003588. Small, M., Xu, X., 2012. The Instability of Followers and Emergent Vorticity in Flocking Behaviour for an Experimental Interaction Rule, arXiv preprint arXiv:1207.5290. Stemler, T., Judd, K., 2009. A guide to shadowing filters for forecasting and state estimation. Physica D 238, 1260–1273. Walker, D.M., Carmeli, C., Perez-Barberia, F.J., Small, M., Perez-Fernandez, E., 2010. Inferring networks from multivariate symbolic time series to unravel behavioural interactions among animals. Anim. Behav. 79, 351–359. Xu, X.-K., Kattas, G.D., Small, M., 2012. Reciprocal relationships in collective flights of homing pigeons. Phys. Rev. E 85 (February), 026120, http://dx.doi.org/10. 1103/PhysRevE.85.026120. Zaitouny, A., 2016. One-Dimensional Tracking Methodology Using Shadowing Filter, Mathworks/Matlab Central/File exchange. Zaitouny, A., Stemler, T., Algar, S., 2016a. Optimal Shadowing Filter Tracking Methodology with Application to Chaotic Dynamics of Moving Objects (submitted for publication). Zaitouny, A., Stemler, T., Judd, K., 2016b. Tracking Rigid Body with Limited Information: Using a Shadowing Filter Based on Newtonian Dynamics (submitted for publication). Zaitouny, A., Stemler, T., Small, M, 2016c. Tracking a Single Pigeon Using a Shadowing Filter Algorithm (submitted for publication). Zhang, H.-T., Chen, Z., Vicsek, T., Feng, G., Sun, L., Su, R., Zhou, T., 2014. Route-dependent switch between hierarchical and egalitarian strategies in pigeon flocks. Sci. Rep., 4.