Delayed Response and Biosonar Perception Explain Movement Coordination in Trawling Bats

Animal coordinated movement interactions are commonly explained by assuming unspecified social forces of attraction, repulsion and alignment with parameters drawn from observed movement data. Here we propose and test a biologically realistic and quantifiable biosonar movement interaction mechanism for echolocating bats based on spatial perceptual bias, i.e. actual sound field, a reaction delay, and observed motor constraints in speed and acceleration. We found that foraging pairs of bats flying over a water surface swapped leader-follower roles and performed chases or coordinated manoeuvres by copying the heading a nearby individual has had up to 500 ms earlier. Our proposed mechanism based on the interplay between sensory-motor constraints and delayed alignment was able to recreate the observed spatial actor-reactor patterns. Remarkably, when we varied model parameters (response delay, hearing threshold and echolocation directionality) beyond those observed in nature, the spatio-temporal interaction patterns created by the model only recreated the observed interactions, i.e. chases, and best matched the observed spatial patterns for just those response delays, hearing thresholds and echolocation directionalities found to be used by bats. This supports the validity of our sensory ecology approach of movement coordination, where interacting bats localise each other by active echolocation rather than eavesdropping.


Introduction
Group movement patterns are dependent on perceptual inputs, governed by cognitive mechanisms, and constrained by the motor abilities of the interacting agents [1][2][3][4]. The intrinsic sequential occurrence of these three stages gives rise to delayed movement responses [5]. Current collective movement interpretations acknowledge such delays [6][7][8][9], but to explain observed patterns the existence of generic, albeit plausible, interaction rules are postulated. Since individual decisions depend on the information about the actions of other conspecifics [10][11][12][13][14] and movement response are shaped by an animal's sensory-motor abilities, it is however biologically most realistic to explain interaction exclusively from perception biases and delayed responses without inferring social forces from the observed movements.
Even though various models are capable of capturing many of the emerging features of coordinated movement patterns, the individual behavioural rules are heuristics of the underlying mechanisms. There is awareness that more biologically oriented experimental analysis and modelling is necessary [15]. Recent experimental studies on golden shiners (Notemigonus crysoleucas) [16] have indicated that speed regulation along the direction of travel appears a more dominant component of the interaction compared to alignment, whereas studies in mosquitofish reported the lack of any justifiable alignment rule [17]. From three-dimensional stereographic images of starling flocks [18], evidence of a topological rather than a metric interaction rule was brought to the fore, i.e. individuals interacted with an average of six or seven neighbours rather than every conspecific within a certain distance. Other network topologies have also been proposed recently based on animal communication across a 'Voronoi shell' of nearest-neighbours [19] or the network of visible neighbours estimated from a ray-casting algorithm around a fish eye [20]. Other investigations have also questioned the assumptions about how individuals perceive and integrate information from their neighbours. Experiments on saithe (Pollachius virens) pointed to the importance of the use of lateral line compared to vision on school structure and dynamics [21].
All these recent efforts point to the importance of sensory inputs that shape an animal's movement decisions. Animal species that use an active sensory modality, such as echolocation, offer an opportunity to tap directly into the source of this sensory information. Active-sensing, compared to passive, requires the emission of a signal whose spatial propagation is governed by quantifiable physical laws [22]. Therefore, once the call intensity and echolocation directionality of an echolocating species are known, it is possible to infer with sufficient accuracy over what spatial range individuals can detect and track each other, either by active sensing or, over greater distances, by eavesdropping on the others' signals [23][24][25].
Here we introduce and test the biosonar movement interaction (BSMI) hypothesis, whereby bats interact with conspecifics, to perform coordinated flight, based on the information collected by their biosonar and not by eavesdropping. While foraging, bats send out sonar search pulses several times per second. Whenever a sound pulse hits an object an echo is generated, which then travels back to the bat who interprets the echo for target detection, recognition and localisation. Because real objects do not reflect all the impinging sound back to the echolocator (target strength < 1) and because the sound has to travel the bat-object distance twice incurring double propagation losses, the biosonar survey range is limited to several meters at best (see e.g. [24,26,27]). The spatial profile of biosonar is therefore determined by the amplitude of the biosonar pulses and by the directionality of sound emission as well as the directionality of bat hearing, where the biosonar detection range is greatest on axis and drops off towards either side. The core of the BSMI hypothesis is thus that the directional spatial profiles of the individuals' biosonar fields determine the animals' mutual positions and headings during movement interactions, with the information obtained via eavesdropping having a minor role in movement coordination. Note that by eavesdropping we mean locating conspecifics by listening to their echolocation calls [23] and not locating oneself by trailing a conspecific and extracting information from the echoes resulting from their calls [28].
Sound field emission profiles of echolocating bats have been documented in various species [29] including Daubenton's bat, Myotis daubentonii, our study animal. This widespread medium-sized species [30] has been selected for two reasons. On one hand, the animal's cognitive abilities allow for rich behavioural patterns. On the other hand, its habit to forage low over still water surfaces to glean arthropod prey [31] reduces the spatio-temporal complexity of the interactions to just a two-dimensional plane.
To test our BSMI hypothesis we collected movement trajectories of Daubenton's bat pairs foraging freely in the wild. Based on these observations, we proceeded in three steps: (i), (ii), and (iii). First, we built a behavioural classifier (i) based solely on the paired movement paths, which segments these movement trajectories into different behaviours, and determines actorreactor roles by extracting response delays. We then tested the BSMI hypothesis by comparing the relative positions of the interacting bats with the calculated biosonar perception field (ii). Finally, we constructed and tested a BSMI model of interacting bats (iii) whereby the mechanisms of interaction are based exclusively on aerodynamic constraints. i.e. actual speed and acceleration, and sensory inputs, i.e. the known spatial extent of the sonar perception field. The quality of the BSMI model is judged against its ability to reproduce the observed spatio-temporal interaction dynamics of the bats.

Results
We recorded flight interactions of foraging Daubenton's bats flying in the South-West corner of Barrow Tanks reservoir number three near Barrow Gurney in England. The experimental set-up allowed for a 20 ms temporal resolution of flight tracking (see Methods section). The recordings were taken over 10 evenings in the summer of 2009: June 15 th , 22 nd , 29 th ; July 1 st , 7 th , 8 th , 9 th , 15 th ; and August 19 th and 25 th . Overall we collected nearly 70,000 individual data points (fixes) including 70 interactions as displayed in Table 1.
To identify interacting flight behaviour we used two mathematical tools: the time-dependent delayed correlation (TDDC) and the time-dependent delayed separation (TDDS). For a given point in time, the TDDC compares the heading of one of the bats with the heading of the other at all other times, thus allowing us to establish whether one bat was copying the direction of the other (coordinated flight). Similarly the TDDS compares the position of one of the bats at one time with the position of the other bat at all other times, thus allowing us to specify if one bat was copying the position of the other (chase flight). The remaining data, when no indication of interaction emerged, were deemed unclassified.
In Fig. 1A we show a sample trajectory (shown also in S1 Video), whose reaction delays have been extracted from the time-dependent delay directional correlation (TDDC) plot displayed in Fig. 1B-the corresponding time-dependent delay separation (TDDS) plot is not shown in this example as it reveals absence of chase behaviour. The presence of a black dot with coordinates (t, τ) in the C i,j (t, τ) plot shows that bat j has copied the heading of bat i with a delay τ at time t. A straight black line segment indicates that the bats have maintained a constant reaction delay over a certain period of time, whereas an inclined line, e.g. in the intervals 3.14 s < t < 3.20 s and 3.50 s < t < 3.88 s, illustrates an interaction with changing delay. If the extracted delay turns negative, bat i now copies the heading of bat j and actor-reactor roles have swapped. In this sample trajectory, this swap has occurred in two stages. At t = 1.4 s bat i first terminates the nearly parallel flight it was maintaining with bat j, and subsequently, at t = 1.8 s, tries to align its heading with that of bat j (see caption of Fig. 1 for more details).
We analysed those flights which had more than one bat in the field of view of the camera at one time. Our analysis allowed us to classify each recorded data point from all seventy paired trajectories as non-interacting (74%) and interacting behaviour (26%), with the latter further subdivided into chases (4%) and coordinated flights (22%). Non-interacting is the definition given to flight segments where there is low correlation between the headings of the individuals. Coordinated flight is defined as a flight segment with high correlation between the headings of the individuals, and chase flight is a flight segment where one individual occupies the previous positions of another. Note that chase flights are a subset of coordinated flights. As shown in Fig. 1C chase flights and most coordinated flights were characterized by similar response delay values of up to 500 ms, which corresponds to 5-7 biosonar calls, i.e. information updates, at the known average call interval of 70-100 ms [32]. One marked distinction between chases and coordinated manoeuvres was in their leader-follower dynamics, synonym here of actor-reactor. During chase flights, the individual ahead was clearly the leader [33], as the follower was retracing the movement path of the individual in front. This leader-follower hierarchy was not equally prominent when bats were performing coordinated flights, for example the reactor was not necessarily behind the actor (see e.g. S1 Fig. panel (f)).
In summary to realise step (i) in the validation of our BSMI hypothesis, we reconstructed movement trajectories from videogrammetry data, analysed the corresponding TDDC and TDDS plots, and extracted the associated delay values. This allowed us to distinguish movement paths among non-interacting flights, chases and coordinated flights. For the last two categories, the sign of the extracted delay values provided the sought actor-reactor roles as a function of time for each interacting paired trajectory (see S1 Fig. panels   To quantify alignment dynamics of the interacting bats, in Fig. 2 we plot bat separation distance, relative heading and relative position [19]. Chases were characterised by low separation distances (Fig. 2B), high alignment ( Fig. 2C) and narrow frontal exposure angles (Fig. 2D), while coordinated flights occurred over larger separations, wider range of relative headings and only a moderate preference for frontal exposure angles (see S2 Fig. for the difference in the TDDS and TDDC plots between a coordinated flight and a chase). Unclassified data showed no preference in relative headings and exposure angles, yet comparison with randomised pairing of trajectories (dashed lines) revealed that unclassified data points were overrepresented at separations below 3.5 m and underrepresented above approximately 7 m. This indicated that even non-interacting bats tended to fly closer than random.
With the identified actor-reactor roles, it became possible to superimpose the bat's echolocation field and test whether a plausible explanation for the observed interacting flights is the existence of an acoustic response threshold-step (ii) in our BSMI hypothesis. A match between the actual spatial patterns of interactions with these superimposed sound fields would indicate a potential functional link. In Fig. 2A, where we plotted the actor's position relative to the heading and location of the reactor, we added isocontours of the calculated biosonar field, i.e. the amplitude of the echo reflected from the actor's body returning to the reactor. We noticed that 99% of observations were above 0 dB and 97% above 10 dB, i.e. almost no reaction  occurred when the actor's echo was below the likely range for the reactor's hearing threshold of 0-10 dB [34]. In short, there was no reaction when the actor was outside the reactor's echolocation field.
Step (iii) in our BSMI hypothesis consists of investigating whether alignment dynamics with delay in combination with the perceptual bias created by the shape of the echolocation field suffice to produce the observed behavioural interactions. For that purpose we constructed a mathematical movement model (BSMI model) of echolocating bat pairs capable of generating interacting trajectories that mimic the patterns of chases and coordinated flights observed in the field data (see Methods). From the spatio-temporal patterns emerging in the BSMI model, we generated relative position plots in Fig. 3 for different delay response τ, directionalities of the echolocation field-asymmetry parameter A in Equation (11) below-and hearing thresholdparameter B in Equation (12) below. By varying these parameters in the model we produced actor-reactor and chaser-chasee spatial relationships that can be compared to the observations in Fig. 2A. We found that the best fit was achieved for response delays between 100 ms and 500 ms, and with B = 10 ± 5 dB and A = 16 ± 1 (root mean standard error RMS = 0.0179), see Fig. 4. We used this latter value to calculate sound fields displayed in Figs

Discussion
Movement analysis has a long tradition in animal ecology [35,36], however, the segmentation and classification of trajectories of co-moving individuals is an area of research still in its and heading (upwards) of the other: blue for chases, green for coordinated flights, and grey for unclassified behaviour. In the latter case, the individual at the centre is picked at random for each pair. Reactor is in the centre with upward heading with symbols showing actor's positions and with red lines representing isocontours of the reactor's calculated echolocation field. Parameters: cosine directionality with front-rear difference of 36 dB accounting for both emission and hearing directionality (see Equation (11) below). Source level is 110 dB [29] and absorption is 1.28 dB m −1 [22]. Inset indicates separation distances d ij , relative headings ϕ ij , i.e. the angular difference between the two velocity vectors, and exposure angles θ ij , i.e. the angular position of the actor with respect to the heading of the reactor. Panel infancy [37]. Here we have shown that a use of the velocity and separation delayed correlation maps, respectively the TDDC and TDDS functions, allowed us to classify interaction patterns in the movement trajectories of coflying bat pairs. Considering only those reactions that are strictly ordered in time, it is possible to reconstruct the delay with which one bat copies the heading of the other bat as a function of time.
With our automatic classifier we have found that bats foraging in the field move in coordination five times more often than chasing the bat in front. Although a chaser has a spatial advantage when either the leader fails to capture a prey item (2nd attacker advantage) or when a patch of resources is found (shared exploitation), a chase might also become an aggressive behaviour if the chaser directs unpleasant or even harmful loud signals at the chasee. On the other hand, coordinated flight could improve search efficiency because both individuals hear when and where the other has searched for and found prey, which can allow each member of the pair to survey larger areas of the water surface per unit of time through eavesdropping.
Empirical evidence strongly suggests that active echolocation was used by the reacting individual for alignment during interactions: 99% of all interactions were within the echolocation field and the entire echolocation field down to echo amplitudes of 0-10 dB was used ( Fig. 2A). The possibility that eavesdropping per se might have produced the observed patterns was tenuous. We reached that conclusion by drawing the reactor's position relative to the actor (see S3 Fig.). In that case no match of spatial interaction patterns with the received sound appeared obvious as interactions did not follow any sound exposure isocontour and most interactions happened in a direction opposite to the main lobe of emission. However, the fact that unclassified, i.e. non-interacting, bats flew at closer than random distances (Fig. 2B), often outside the maximum calculated echolocation range but within eavesdropping range (S3 Fig.), suggests the idea that it is beneficial for a bat to remain in the eavesdropping region of others, thereby possibly profiting from overhearing the other's search success. Eavesdropping might also account for In contrast to other species, where alignment occurs by copying the heading of the individual in front, e.g. in Surf scoters [38], we observed several coordinated flights where the reactor was in front of the actor ( Fig. 2A and 2D when jθ ij j > 90°). Although our dataset is limited in size, these large exposure angles indicate that bats may be able to respond to an individual behind them because echolocation can provide information from rear angles e.g. in Fig. 2A an interaction, represented by a row of green points, tracing just outside the -40 dB contour from a position of 3 o'clock (lateral) to 6 o'clock (directly behind). Also note that in contrast to other movement coordination studies foraging Daubenton's bats search for ephemeral and depletable prey and therefore their direction of motion is not targeted at one specific location of a shareable resource. As a result their respective position and orientation could change considerably facilitating actor-reactor swaps within an interaction.
Our BSMI model was able to create interactive behaviours (chases and coordinated flights) given that the model is based on biosonar perception with no input from movement interaction data other than observed flight speeds and turning angle distributions. Most importantly, the model output can reproduce the observed spatial interaction patterns only when the emission directionality A, hearing threshold B and delay were at their natural values. The best fit directionality value A = 16 is in the middle of the known range of emission directionality for this species (A = 7.3 for an independent measurement shown in S4 Fig., and A = 25.7 [29]), and the hearing threshold of B = 10 dB is also within the expected range [34].
In this study we have provided an answer to how individuals react to each other while they are foraging together, but it was not possible to determine why. As at least some data were collected at dates after juvenile bats had taken flight, it is not entirely impossible that some flight interactions were training flights between mothers and their offspring. Increased foraging success might be a driver for coordinating flight manoeuvres. This might be possible because bats that coordinate their movement at short distance have each other survey zones overlap and may thus enter into coordinated exploitations in presence of locally rich prey patches. Alternatively, with small prey density a bat may enter agonistic interactions forcing conspecifics away to monopolise a diminishing resource patch. The chases we observed might well contain such an agonistic component [39]. With the help of collective movement models of delayed interaction [40] future research could address some of these issues by testing responses of conspecifics to the locations and timings of foraging success of others, which they inadvertently indicate by their feeding buzzes.
In conclusion, we corroborate our BSMI hypothesis of foraging bats by (i) classifying bat interactions in the wild from highly resolved trajectory data, (ii) explaining the individuals' relative positions based on the calculated biosonar perception field, and (iii) fitting the experimental observations with a bio-inspired movement model of interacting bats. In the process we have constructed a procedure to identify actor and reactor based not on their relative distance but on the delay dynamics of their relative orientation (defined through relative headings and exposure angles), which is independent of the animal species studied and can be applied to all types of coordinated movement data. The emerging actor-reactor relations and the associated delays with which an animal copies the heading of a conspecific provides a moving reference frame with which to determine when, where and how individuals respond to the actions of a conspecific and represents a powerful tool in the arsenal of computational techniques for movement ecology [41,42].

Experimental methods
Animal movement trajectories were obtained by videogrammetry using two 525EX CCIR (Watec Inc., Newburgh, NY) with Cinegon 1.8/4.8 IR lens (Jos. Schneider Optische Werke GmbH, Bad Kreuznach, Germany) on tripods at heights between 4.13 and 5.66 m above the water surface as measured with Laser rangefinder DLE50 professional (Robert Bosch GmbH, Gerlingen, Germany). Pitch towards the water surface and lateral roll angles of the cameras were measured using Clinometer GeoMaxiclin (Geo Supplies Ltd, Sheffield, UK) to the nearest half degree. The recording area of approximately 900 m 2 was illuminated with 850 nm IR light by a Raymax 200 (Raytec, Ashington, UK) positioned between the cameras at the shoreline. The video outputs were recorded synchronously on one XM-DVR PRO S153 (Datatoys, Mequon, MI) solid state dual-channel video recorder. Videos with bat flights were selected manually. Bat positions in the videos were selected by hand in every frame and the corresponding spatial positions calculated in custom written scripts (M. Holderied and H. Goerlitz) in MATLAB (The MathWorks Inc., Natick,MA). Each camera was calibrated and pixel based ray bundles for each camera were obtained using VMS software (Mark Shortis and Stuart Robson). Fields of view of both cameras did overlap only slightly and bat positions were not triangulated by stereo-imaging but rather by intersection of the pixel bundles originating at position, height and orientation of the single respective camera with a 2D plane 20 cm above the water surface, which was taken as the average flight height of this species [43]. The slight left-right asymmetry in the dataset indicates a preference for clockwise turns and a slight tendency to move parallel to the shoreline (S5 Fig.). This resulted in fewer trajectories moving away from the cameras, and thus making the (projected) reconstruction of the flight paths with similar level of errors. Videos were recorded in VGA format at 25 Hz frame rate. Each VGA frame is based on two half frames (odd and even lines taken in alternation). We analysed each half frame independently and interpolated the respectively missing odd or even lines such that the effective rate was 50 half frames per second, allowing for a temporal resolution of flight tracking equal to Δt = 20 ms. This is about a factor four higher than the average call interval, i.e. the information update rate of Daubenton's bats [32].

Data segmentation
The dataset of paired bat trajectories was analysed to determine the fraction of data where interactions had occurred. For those data deemed to represent interaction we identified leader and follower dynamics in pairs of co-flying bats based on the sign of the automatically extracted delay values as a function of time with the help of the time-dependent directional correlation (TDDC) function, shown in Equation (2) below, and the time-dependent distance separation (TDDS) function, displayed below in Equation (4). Analysis of the TDDC and TDDS plots allowed us to segment paired movement trajectories into different behaviours.
The time-dependent delayed directional correlation. A good indicator of the alignment of two individuals is the normalised scalar product of the velocity of bat i at time t with bat j at time t + τ expressed by n ij ðt; tÞ ¼ṽ i ðtÞ jṽ i ðtÞ j Áṽ j ðt þ tÞ jṽ j ðt þ tÞ j ; which varies between 1, indicating perfect alignment, 0, when the heading of the bats is orthogonal, and -1 when in opposite directions. In our recorded trajectories it was clear that bats changed behaviour dynamically as time t progressed. Integrating ν ij (t, τ) over all N values of the observation times [6] would wash out many of the transient features. An average over a smaller time-window empirically adjusted to the temporal scale of the observed dynamics was more appropriate [7,44,45] and represented a time-averaged correlation between the direction of the two bats. This averaging was also necessary to compensate for high-frequency noise due to pixelation and measurement errors. For each paired trajectory we thus computed the TDDC function where v ij (t, τ) has been averaged over a symmetric time-window 2wΔt (t k = kΔt) for times t and delays τ constrained by These constraints generate the sharp boundaries (see e.g. in Fig. 1B, S1 Fig. panel (b) and S2 Fig. panels (c) and (e)) and correspond to the line of points given by τ + t = 0 s and τ + t − NΔt with Δt = 0.02 s. To obtain sufficiently smooth trajectories the window width w = 2 was chosen, corresponding to a moving average over five data points for each coordinate, two on each side of the spatial coordinate pair. Note that the horizontal axis is relative to bat i. To interpret time relative to bat j one ought to compute C ji (t, τ), but only the exchange i ! j together with t ! t − τ in C ij (t, τ) would produce the same TDDC plot. To define when two animals performed an interacting flight, it was necessary to set a threshold value v c for C ij (t, τ) at and above which the heading of the two bats were deemed sufficiently correlated. To determine the effects of the choice of correlation threshold on our results, we studied how the extracted delays change as a function of v c . Since the appropriate v c is clearly species dependent, requiring a judgement of an animal's ability to align, our choice was a compromise guided by the need to capture bats' coordination in the observed manoeuvres, e.g. the turns displayed in Fig. 1A. On one hand, we used a sufficiently high threshold to avoid misinterpreting noise features with alignment and to capture when a reacting animal attempted to copy the actions of the actor. On the other hand, we reduced the threshold value from v c = 1, corresponding to persistent and perfect alignment, to increase the number and length of the interacting trajectories across the entire dataset. In S6 Fig. we display how the length of these interacting segments changes in each trajectory for different choice of v c . Overall there is a limited increase on the number of interacting data points for threshold values smaller than 0.95. In addition Gaussian fit of the resulting delay histograms for values 0.85 v c 0.95 did not show significant differences either in the standard deviation or the 95% confidence interval from a Gaussian fit of Fig. 1C with the choice v c = 0.95. The choice v c = 0.95 was hence used to identify contour levels in the TDDC plot where interactions might have occurred (e.g. in Fig. 1B, S1 Fig. panel (b), S2 Figs. panels (c) and (d)).
The time-dependent delayed separation. To identify chase flight behaviour in the dataset, we considered R ij (t, τ), the spatial separation between bat i at time t and bat j at time t + τ averaged over a time-window w, namely In Equation (4) (x l (t), y l (t)) with l = i or j is the coordinate pair of bat i or j at time t. When this separation equals zero for positive τ, then bat j occupied all previous positions that bat i had between time t − τ − 2Δt and t − τ + 2Δt. To avoid classifying trajectories where bats crossed each other's paths with arbitrary headings and delays as chases, the set of paired values (t, τ) where interaction occurred were first extracted from the analysis of C ij (t, τ). R ij (t, τ) was then computed at the specific (t, τ) values and the interaction was classified as chase if R ij (t, τ) d c . Hence chase flights were a sub-category of coordinated flights, as bats were both well aligned, C ij (t, τ) ! v c , and had small spatial separation, R ij (t, τ) d c , in the delayed reference frame. The threshold used was d c = 0.15 m, which is smaller than 0.25 m, the average wingspan of a Daubenton's bat. It turned out that for all behaviours classified as chases the average separation in delayed time were at most 0.09 m.
Extracting delay values. For each time t in the TDDC and TDDS plots within the respective regions of interaction, i.e. C ij (t, τ) ! v c and R ij (t, τ) d c , it is possible to associate different delay values τ. A procedure was thus developed to extract a single interaction delay τ, and its associated uncertainty, for each time t. These values can then be used to identify leader and follower and subsequently to classify the trajectory in terms of the behavioural mode of interaction (chase and coordinated flight).
There are several potential methods to extract such 'interaction path' (delays as a function of time) from the TDDC plot. Two straightforward methods consist of taking the delay value which maximises the correlation: over the extent of the entire recording of a bat pair [6], or over a time-window, either forward from [44,45] or symmetric around each time-step [7].
Whilst easy to perform, these procedures are poor at reconstructing biologically realistic interactions most importantly because they violate the simple assumption that each bat's reactions are strictly ordered in time. To understand this, consider that bat i acts at time t and bat j reacts at time (t + τ) then a subsequent action by bat i at time t + Δt can be reacted to by bat j at time (t + Δt + τ + Δτ), where Δτ represents a change in delay. If Δτ < −Δt then (t + Δt + τ + Δτ) < (t + τ), and bat j performed the second reaction before the first reaction at time (t + τ). Such situations are deemed nonsensical and hence changes in delay should be restricted to Δτ/Δt ! −1.
In the TDDC plot all possible 'interaction paths' represent delay values that change over time. The actual response delay across the trajectory is obtained by computing the cumulative value of C ij (t, τ) from left to right, i.e. from the first to the last observation time, according to three requirements: C ij (t, τ) < v c , the change of delay along a path must fulfil the time-ordered requirement, and the actual C ij (t, τ) value along a possible path gets reduced if delay changes over time.
More specifically, each possible path starts from the coordinates (t 0 , τ Ã ), with t 0 = wΔt and (−N + w)Δt < τ Ã < (N − w)Δt, and finish at (t max , τ max ), with t max = (N − w)Δt. Any path is thus represented by the coupled elements (t, τ) with and the corresponding delay values, with M limited by the various constraints.
The 'optimal' path has the highest overall correlation for each pairwise trajectory given by p l ðt l Þ ¼ X NÀw n¼w C ij ðt n ; t l n Þ; if C ij ðt n ; t l n Þ ! v c and t l n ¼ t l nÀ1 ; C ij ðt n ; t l n Þ À q; if C ij ðt n ; t l n Þ ! v c and t l n 6 ¼ t l nÀ1 ; 0 otherwise; where q is a penalising cost associated with changing delay. The inclusion of this penalising cost in the algorithm for extracting delay values was imposed to account for the added cognitive requirements that an animal would incur if it were to interact by constantly changing the timeframe of reference. As the animals' relative locations fluctuate considerably during an interaction-change of speed and heading is a necessity during foraging activities-a constant delay may in fact represent a steady source of information through which to control the conspecific's actions and take decisions accordingly. The appropriate q value to use depends on the profile of the correlation map which in turn is affected by the spatio-temporal trajectories. For our bats, the value of q had to be small enough to allow for delay changes when approaching the contours delimiting the interaction region in the correlation map, and large enough to restrict the change of delays when away from those limiting contours. These goals were met with the choice q = 0.02. By setting to 0 the contributions of all paths with profiles Δτ/Δt ! −1 and favouring those with zero slope, the extracted delays represent a compromise between having the longest segments within the disconnected regions of interaction and the shallowest slope. In this way we identify interaction events as those behaviours with the longest duration and with limited variations in delay response. In other words we associate delayed interactions with actions and reactions that are persistent rather than brief or instantaneous.
The procedure for extracting delays for chase flight behaviour differs slightly. Any region of the TDDC map exceeding the threshold v c and with the corresponding TDDS region less than or equal to d c pertains to chase-flight, therefore the interaction path is constructed by minimising the delayed separation, (rather than maximising the directional correlation), again satisfying the constraint that interactions be strictly ordered in time and with a penalising cost for changing delay. If a trajectory contains both chase flight and coordinated flight the chase flight interaction path is found first and then the coordinated flight interaction path is chosen as the maximum correlation path that connects to the chase paths.
Uncertainty in delay extraction. Current methodologies to determine reaction delays from movement data are unable to specify ranges of uncertainty of the extracted values. The use of the TDDC and TDDS, however, provides the means to quantify this delay uncertainty. This uncertainty, which can be asymmetric around the extracted delay, is affected by the classification (chase or coordinated flight) as well as the curvature of the trajectories with lower uncertainty being obtained when bats make interactive turns as compared to when they move straight in coordination.
Choosing between multiple delays at a given time in the region of high correlation in the TDDC plot is a source of uncertainty. For coordinated flight the number of possible delays at each time is dependent on the trajectories themselves. Consider two extreme examples, one where the paired trajectories are parallel lines and the other where both trajectories contain many turns and are therefore curved and winding. In the first case the bats are perfectly aligned at all times and for all delays. In such a case C ij (t, τ) ! v c for all (t, τ), and identifying uniquely the delay with which the bats performed this flight is not possible. In the second case, as a result of all the turns, the delay with which the bats flew is more apparent and the number of delays that exceed the correlation threshold will converge to the correct one as the threshold is increased.
For a given threshold, one could determine the range of values around the extracted delay for each value of time in the TDDC plot. However using the relationship between the delay range and the curvature of the trajectory allows to obtain a smaller bound of the delay uncertainty. To understand this, consider S7 Fig. where we plot an hypothetical TDDC function whereby the associated paired trajectories contain flight segments with high-curvature ending at time t hc . These segments subsequently transition to parallel motion for both individuals beginning at time t p . There is obviously a large range of possible delays for the parallel segment, ½t # p ; t " p , and a small range, ½t # hc ; t " hc , for the high-curvature segment, and a highly dynamic range for the transition region t hc t t p . However due to our condition, Δτ/Δt ! −1, if the change in the minimum delay is greater than the elapsed time, i.e. ðt # p À t # hc Þ=ðt p À t hc Þ < À1, (or similarly for the maximum delays looking backwards in time), some of the parallel-flight delays cannot contribute to the interaction path and the delay choice is made from the reduced range, ½t # pÃ ; t " p . In this way knowledge from the high-curvature flight has reduced the uncertainty in choosing delays from the parallel-flight. Although we have only compared two times to illustrate the idea, the procedure consists of drawing the −1 gradient lines forward in time at each time point from the minimum delay value, and backward in time from the maximum delay value. Anywhere these gradient lines are within a threshold contour, the delay is reduced to the value along that line. For example for all times between the marked values of t 1 and t 2 , the −1 gradient line lies inside the contour and so reduces the delay-range for those times. The extracted delay uncertainty is therefore quantified as the delay range after the time ordered condition has been applied to each point along the contour.
For coordinated flights we thus have that the straighter a trajectory, the larger is the uncertainty in the delay. This however does not occur in chase flights for which the delay uncertainty is independent of the curvature of the flight because only a small delay window allows one bat to occupy the other bat's previous position (see e.g. S1 Fig. panel (d) and S2 Fig. panel (f)). Hence the uncertainty in the chase flight delays is quantified simply by the (small) range of possible delays.
A comparison between the use of the time-independent directional correlation function [6], as well as of the 'maximal' path in C ij (t, τ) [7], with our procedure to extract delay values is displayed in the 8-panel S8 Fig. In panels (d) and (f) one can see that tracking the maximum of the TDDC function gives a delay path that displays sudden jumps, at times between positive and negative values. Furthermore, it is apparent from panels (g) and (h) that this maximum delay path violates our time-ordering procedure being outside the allowed uncertainty. A clear example of the inability to identify the correct response delays by tracking the maximum of the directional correlation is shown in S8 Fig. panel (f). There, around t = 1.3 s the maximum of C ij (t, τ) (green line) suddenly jumps from τ ' 0.2 s to τ = 1 s. This jump occurs because bat i's heading after the right turn at t = 1.4 s is approximately parallel to bat j's heading after its own right turn at t = 2.3 s (see the corresponding trajectory plot in Fig. 1A).

The BSMI model
An agent based model was created to demonstrate the effects of asymmetric perception and reaction delays on the spatio-temporal trajectories of simulated bat pairs. The model contains the minimal ingredients necessary to account for the salient sensory biology features of the bats and their interaction. The individual movement statistics of a bat follows that of a correlated random walker [46,47] with turning angles and speed values obtained from the observed ones. Each animal is limited in its manoeuvrability by aerodynamic constraints and tries to align to a conspecific present within its interaction field. Each bat echolocates every 100 ms with a certain hearing and emission directionality and detects the presence of the other animal only when above a certain threshold. When another bat is detected, an individual responds by copying the other's heading with some delay (see Equations (8)(9)(10)(11)(12) and S9 Fig. for further details). The simulations were run in discrete time, for 10 s and Δt = 20 ms corresponding to the resolution of the observations, and continuous space by generating initial random positions and headings inside a region of 30 × 30 m 2 .
Flight dynamics. Bat i's heading ϕ i (t), speed magnitude v i (t) and location (x i (t), y i (t)) are updated at each time step as follows: v i ðtÞ $ Rð4:81 m s À1 ; 2:18 m s À1 Þ; if j g i ðtÞ À i ðt À DtÞ j b i ðtÞDt; x i ðtÞ ¼ x i ðt À DtÞ þ v i ðtÞ Dt cos ½ i ðtÞ; In Equation (8) the symbol * means that v i (t) is drawn from the Rician distribution [48] R(μ, σ) with non-centrality parameter μ and scale parameter σ, which was fitted (R 2 = 0.97) to the observed movement patterns of paired flying bats. The desired heading of bat i, that is the heading that a bat would choose depending only on its behavioural state (see below) at time t, is γ i (t). The actual heading at time t − Δt and time t, taken by bat i once the aerodynamic constraints have also been accounted for, are represented, respectively, by ϕ i (t − Δt) and ϕ i (t).
The angular rate of change β i (t) is limited by aerodynamic constraints corresponding to a maximum lateral acceleration a c of 4 g (4 times the gravitational acceleration of 9.81 ms −2 ), providing a maximum given by β i (t) = a c /v i (t). The 4 g limit was measured to contain 99.9% of the turning angles from the paired flight data. Examples of an interacting and independent path generated by the BSMI model are displayed in S9 Fig. The bat may be in either of two behavioural states: independent or interacting. The behavioral state determines the shape of the distribution from which an individual would draw its desired turning angle γ i (t), and it is updated in two ways. (i) Every five time steps an individual would send an echolocation call to determine the heading of the other bat. When a conspecific is detected, a bat enters an interacting behavioural state and attempts to align by selecting a turning angle that reduces the difference in relative heading. If the other bat is not detected, the individual enters an independent state and it maintains that behaviour until it detects a conspecific. Once in the independent state the bat selects a turning angle from a broader distribution. (ii) The other way a bat may change its behavioural state is by moving from the interacting to the independent state at each time step. Even though a bat is in the process of aligning to another individual detected in the last echolocation call, it may find convenient to evade alignment, e.g. for the purpose of being the first to detect a prey or to ensure exclusive access to certain areas of the water's surface. To reproduce this feature in the model an individual in the interaction state may choose, at each 20 ms time step, to enter an independent state with a small probability of 0.05.
Independent flight. When an individual is in the independent state, the animal's turning angles are drawn from a von Mises distribution [49], displayed in S9 Fig. panel (a), centred around its heading at the previous time step. The desired heading of bat i at time t is thus where M(ξ, κ) represents a von Mises distribution with mean direction ξ and concentration parameter κ. The concentration parameter κ = 557, was obtained by a fit (R 2 = 0.997) to the turning angles of unclassified flying bats.
Interacting flight. lnteraction in the form of alignment was constructed by setting a bat's desired heading as the mean of its own with the one of its conspecific [2,50] but with the characteristic feature that the heading of the conspecific is the one at time τ in the past. In other words, at time t the desired heading of bat i is drawn from a von Mises distribution centred around the average between bat i's heading at time t − Δt and the heading that bat j had at time τ in the past. The concentration parameter of the von Mises distribution was chosen to be intermediate between the turning angle distribution observed in unclassified flights and a distribution with no noise. To avoid the need to impose additional assumptions that specify under what conditions individuals would choose to enter a coordinated or chase behaviour, we used the turning angle distribution fitted (R 2 = 0.999) to the observed chases, whose concentration parameter was κ = 2473 (S9 Fig. panel (b)). The desired heading of bat i is thus Sound field modelling Echolocation directionality. A bat may enter into an interacting state if the received echo amplitude is above its hearing threshold. This echo amplitude is the result of the echolocation directionality of the animal as well as the target strength of the echo-giving object and the transmission losses during sound propagation [51]. By considering a source level SL = 110 dB [29] and a target strength TS = -10 dB, the echo amplitude α, expressed in dB, at angle z from the emission direction and returning from a bat at distance ρ in meters, is given by The third term, in curly brackets, in Equation (11) corresponds to the azimuthal echolocation directionality of the bat consisting of (i) an emission directionality with a (positive) asymmetry parameter A and (ii) a hearing directionality. The representation of the emission directionality with a cosine law is supported by the experimental results displayed in S4 Fig.  Hearing directionality was calculated from a head-related transfer function (HRTF) obtained from a 3D acoustic modelling of a Myotis daubentonii head for the frontal hemisphere and for frequencies between 35 and 45 kHz [52]. Hearing directionality was very broad with attenuation being -2 and -3 dB at ±90°, respectively. For simplicity we also fitted a cosine function giving an asymmetry parameter equal to 2 dB (R 2 = 0.44), which we have used in the third term of Equation (11). The fourth term, in square parenthesis, is (i) the sound attenuation of a spherical point source relative to a reference distance of 0.1 m, and (ii) the sound absorption losses c = −1.28 dB m −1 for 10°C, 86% relative humidity and 45 kHz following [22]. The factor of 2 in front of the square parenthesis accounts for the distance to the target and back.
Flight dynamics simulation with echolocation. In a simulation run of the model, every five time steps bats i and j emit an echolocation call with the axis of emission corresponding to the bats' headings. The attenuation of the emitted sound at the location of the other bat is then estimated by knowing the exposure angles θ ij (t) and θ ji (t) and their relative distance d(t) = d ij (t) = d ji (t), and computing α ij (d(t), θ ij (t)) for the sound emitted by bat i and α ji (d(t), θ ji (t)) for the sound emitted by bat j. If α ij or α ji are greater than the attenuation threshold B, the bats carry on independently of the conspecific. However if either attenuation is less than or equal to B the corresponding bat will enter the interacting state after a reaction delay τ. Combining Equation (9) and (10) and similarly for bat j.
To verify that the delays extracted from the movement trajectories are biologically meaningful, we ran the model using no delay, with delays within the observed range 100-500 ms, as well as larger delay values between 500 ms and 2 s. We generated spatio-temporal trajectories of co-flying pairs as a function of the emission directionality A and threshold B parameters with values -60 dB B 60 dB and 0 A 72. The model-generated trajectories were then analysed and classified with the analysis of the TDDC and TDDS plot, and the locations of the actor relative to the reactor compared to the observed spatial patterns of interaction as displayed in the relative position plot (Fig. 3). The line connecting the two recording cameras corresponds approximately to 180°and 0°. In both panels the measured headings, rather than being uniformly distributed, are biased towards mean angles. Computing the mean and standard deviations we obtain 349.4°and 169.4°, (±33.1°), for the individual flying bats and 348.0°and 168.0°, (±36.7°), for the paired flying bats. This bias might be caused by the bats trying to avoid the perimeter demarcating their foraging environment: one part along the shoreline where the cameras were mounted and the other part perpendicular to the latter and located just outside the recording region. (TIFF) S6 Fig. Sensitivity to the choice of directional correlation threshold. For each of the 70 trajectories that comprise the paired flight dataset, the percentage of that trajectory deemed interacting flight behaviour is shown for different correlation thresholds v c . The effect of changing v c is trajectory-dependent with almost no change in some cases, e.g. 14, and large variation in others, e.g. 56. An analysis of each of the 70 paired flights with thresholds smaller than the one used (0.85 v c < 0.95) reveals that only four trajectories (4,38,53,59) brings about interacting behaviour that is otherwise missed. In all the other trajectories, a lower threshold simply increases the length of the interacting paths by converting unclassified segments to classified ones at the beginning or the end of a segment already deemed interacting. The right panel shows the total fraction of interacting flight behaviour for all the data for different thresholds v c . Trajectory 17, 20, 40 and 64 were depicted, respectively, in S2 Fig. panels (a) and (b), and S8 Fig. panels (a) (and also S1 Fig. panel(a)) and (b), whereas trajectory 7, 29 and 56, plotted respectively in S10 Fig. panels (a) (c) and (e), respectively. In panel (d) and (f) we also show the delay path extracted from our time-ordered procedure (black line) and from selecting the so-called 'maximal path' [7], that is the maximum at each time of Equation (2). For panel (d) delay values (and uncertainties) for our time-ordered procedure between 1 s < t < 2 s are extracted by considering the TDDS plot (displayed in S1 Fig. panel (d)) when the interaction becomes a chase. A comparison of the delay values extracted with the different procedures from panel (a) and (b) is displayed, respectively, in panels (g) and (h) with shaded light red areas representing the uncertainty in the extracted delay from our combined analysis of the TDDC and TDDS function.