Abstract

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 forjust 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.

Author Summary

What information individuals use for movement coordination is, however, very difficult to know, except for echolocating bats. These flying mammals perceive their surroundings by emitting loud and high-pitched biosonar calls and listening for the returning echoes. Because bat biosonar imaging is much sparser in information than vision, we can accurately measure the biosonar calls of interacting bats with a group of microphones and then calculate what each of the individuals perceived. When observing pairs of Daubenton’s bats foraging low over water for stranded insects, we found they have intriguing ‘traffic rules’—they chase each other, perform tandem turns and even slow down to avoid collision. When we then modelled their biosonar vieW of the surroundings during these interactions we discovered that one simple trick suffices to create all their interactive behaviours: once another individual is close enough for your biosonar to pick up its echo, copy this individual’s flight direction Within 4—5 of your own Wingbeats.

Introduction

The intrinsic sequential occurrence of these three stages gives rise to delayed movement responses [5]. Current collective movement interpretations acknowledge such delays [6—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 conspecif1cs [10—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.

There is awareness that more biologically oriented experimental analysis and modelling is necessary [15]. Recent experimental studies on golden shiners (Notemigonus cry-soleucas) [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 mosquito-fish 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 neigh-bours rather than every conspecif1c within a certain distance. Other network topologies have also been proposed recently based on animal communication across a ‘Voronoi shell’ of near-est-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].

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—25].

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 in-dividuals’ 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].

This widespread medi-um-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.

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 actor-reactor 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-tempo-ral interaction dynamics of the bats.

Results

The experimental setup 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 15th, 22nd, 29th; July 1“, 7th 8th, 9th, 15th; and August 19th and 25th. Overall we collected nearly 70,000 individual data points (fixes) including 70 interactions as displayed in Table 1.

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.

Table 1 . Summary of the entire experimental dataset. Breakdown of the original spatio-temporal data into the individual and paired flight trajectories together with mean and standard deviation values of the temporal and spatial length as well as speed of all flights.

The presence of a black dot with coordinates (t, T) in the Cid(t, T) plot shows that bat j has copied the heading of bat i with a delay 1‘ 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).

Our analysis allowed us to classify each recorded data point from all seventy paired trajectories as noninteracting (74%) and interacting behaviour (26%), with the latter further subdivided into chases (4%) and coordinated flights (22%). Noninteracting 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. 81 Fig. panel (0).

This allowed us to distinguish movement paths among noninteracting 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 81 Fig. panels (a)-(f) for a graphical representation of the 6 steps of this procedure).

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 noninteracting bats tended to fly closer than random.

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.

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 1‘, directionalities of the echolocation field—asymmetry parameter A in Equation (11) below—and hearing threshold—parameter 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 i 5 dB and A = 16 i 1 (root mean standard error RMS 2 0.0179), see Fig. 4. We used this latter value to calculate sound fields displayed in Figs. 2A and 3.

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 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.

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 be-haviour 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.

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. noninteracting, 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 the low number of observations above -|135|° relative heading (Fig. 2C) in unclassified flights, which indicates that even bats outside mutual echolocation range avoid head-on flights.

2A and 2D when lflijl > 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].

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 conspecif1cs 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 conspecif1cs to the locations and timings of foraging success of others, which they inadvertently indicate by their feeding buzzes.

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 conspecif1c provides a moving reference frame with which to determine when, where and how individuals respond to the actions of a conspecif1c and represents a powerful tool in the arsenal of computational techniques for movement ecology [41, 42].

Methods

Experimental methods

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 rangeflnder DLE5O 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 m2 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 (SS 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

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

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 twith bat j at time t + 1‘ expressed by which varies between 1, indicating perfect alignment, 0, when the heading of the bats is orthogonal, and -1 when in opposite directions.

Integrating vz-J-(t, 1‘) 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

Note that the horizontal axis is relative to bat i. To interpret time relative to bat j one ought to compute le(t, T), but only the exchange 1' —> j together with t—> t

To determine the effects of the choice of correlation threshold on our results, we studied how the extracted delays change as a function of vc. Since the appropriate 1/6 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 VC 2 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 vc. 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 g VC 3 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 1/6 = 0.95. The choice 1/6 = 0.95 was hence used to identify contour levels in the TDDC plot where interactions might have occurred (e.g. in Fig. 1B, 81 Fig. panel (b), 82 Figs. panels (c) and (d)). The time-dependent delayed separation. To identify chase flight behaviour in the dataset, we considered Rij(t, T), the spatial separation between bat i at time t and bat j at time t+ T averaged over a time-window w, namely

and t— T + 2At. To avoid classifying trajectories where bats crossed each other’s paths with arbitrary headings and delays as chases, the set of paired values (If, T) where interaction occurred were first extracted from the analysis of C,-j(t, T). Rij(t, T) was then computed at the specific (t, T) values and the interaction was classified as chase if R,-j(t, T) 3 dc. Hence chase flights were a subcategory of coordinated flights, as bats were both well aligned, Cij(t, T) 2 1/6, and had small spatial separation, Rij(t, T) 3 dc, in the delayed reference frame. The threshold used was dc = 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.

For each time tin the TDDC and TDDS plots Within the respective regions of interaction, i.e. Cz-J-(t, T) 2 vc and Rl-j(t, T) 3 dc, it is possible to associate different delay values T. A procedure was thus developed to extract a single interaction delay T, 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].

To understand this, consider that bat 1' acts at time t and bat j reacts at time (t + 1‘) then a subsequent action by bat i at time t + At can be reacted to by bat j at time (t + At + T + AT), Where AT represents a change in delay. If AT < —At then (t + At + T + AT) < (t + T), and bat j performed the second reaction before the first reaction at time (t + 1‘). Such situations are deemed nonsensical and hence changes in delay should be restricted to AT/At 2 —1.

The actual response delay across the trajectory is obtained by computing the cumulative value of Cl-J-(t, T) from left to right, i.e. from the first to the last observation time, according to three requirements: Cl-J-(t, T) < vc, the change of delay along a path must fulfil the time-ordered requirement, and the actual Cz-J-(t, 1‘) value along a possible path gets reduced if delay changes over time. 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 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 in-teraction—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 con-specific’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 AT/At 2 —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.

Any region of the TDDC map exceeding the threshold vc and with the corresponding TDDS region less than or equal to dc pertains to chase-flight, therefore the interaction path is constructed by minimis-ing 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.

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.

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 Gil(t, T) 2 1/6 for all (t, 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.

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 the. These segments subsequently transition to parallel motion for both individuals beginning at time tp. There is obviously a large range of possible delays for the parallel segment, [Ii IT], and a small range, [rim 1%] , for the high-curvature segment, and a highly dynamic range 1” p for the transition region the g t g tp. However due to our condition, AT/At 2 —1, if the change larly 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, [rim, 1;]. 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 t1 and t2, 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.

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. 81 Fig. panel (d) and 82 Fig. panel (0). Hence the uncertainty in the chase flight delays is quantified simply by the (small) range of possible delays.

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 88 Fig. panel (i). There, around 1‘ = 1.3 s the maximum of Cz-J-(t, 1‘) (green line) suddenly jumps from T 2 0.2 s to T = 1 s. This jump occurs because bat is 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

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 conspecif1c 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—12) and S9 Fig. for further details). The simulations were run in discrete time, for 10 s and At = 20 ms corresponding to the resolution of the observations, and continuous space by generating initial random positions and headings inside a region of 30 x 30 m2.

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 yi(t). The actual heading at time t— At and time 1‘, taken by bat i once the aerodynamic constraints have also been accounted for, are represented, respectively, by gbl-(t — At) and gbz-(t). The angular rate of change fii(t) is limited by aerodynamic constraints corresponding to a maximum lateral acceleration ac of 4 g (4 times the gravitational acceleration of 9.81 ms‘z), providing a maximum given by fii(t) = aC/vl-(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 yi(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 conspecif1c 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 conspe-cif1c. 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.

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(f, K) represents a von Mises distribution with mean direction 5 and concentration parameter K. The concentration parameter K = 557, was obtained by a fit (R2 = 0.997) to the turning angles of unclassified flying bats.

Interaction 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 Tin the past. In other words, at time tthe desired heading of bat i is drawn from a von Mises distribution centred around the average between bat is heading at time t — At and the heading that bat j had at time Tin 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 (R2 = 0.999) to the observed chases, whose concentration parameter was K = 2473 (S9 Fig. panel (b)). The desired heading of bat i is thus

Sound field modelling

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 a, eXpressed in dB, at angle C from the emission direction and returning from a bat at distance p in meters, is given by

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 i90°, respectively. For simplicity we also fitted a cosine function giving an asymmetry parameter equal to 2 dB (R2 = 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.

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 9,7-(1‘) and 91-,(1‘) and their relative distance d(t) = dij(t) = dji(t), and computing a,j(d(t), 90(0) for the sound emitted by bat i and aji(d(t), 917(0) for the sound emitted by bat j. If 05,]- or afi 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 1‘. Combining Equation (9) and (10), the desired heading of bat 1' becomes and similarly for bat j.

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).

Supporting Information

Sample spatio-temporal paired trajectory. Example of a paired flight trajectory, 81 Fig. Summary of the data analysis procedure. Graphical representation of the sequential steps involved in extracting behavioural states to identify actors and reactors (af) and their rel-82 Fig. Representative sample of coordinated flight and chase. The top panels show two different recorded trajectories from the dataset: (a) is an example of coordinated flight behaviour and (b) is an example of chase flight behaviour. Panels (c) and (d) show contour level plots of the corresponding TDDC function Cl-J-(t, 1‘). Similarly panel (e) and (i) show the corresponding TDDS function Rij(t, T). In panels (a) and (b) the similarity in the TDDC plots is evident. It is possible to distinguish clearly the bats’ behaviour only When comparing their respective

Symbols indicate positions of the reactor relative to the location (centred) and heading (upwards) of the actor. For behaviour deemed unclassified, the individual at the centre is picked at random for each pair. Parameters of emission directionality, source level and sound absorption are the same as in Fig. 2A. Red lines are isocontours of the emitted sound field. (TIFF)

Foraging flight trajectory and cumulative emission directionality of Myotis daubentonii over approximately a 30 m diameter pool. (a) Top vieW and (b) lateral vieW of flight trajectory. Open circles are 126 positions of individual sound emissions as determined by acoustic tomography. Filled symbols represent the positions of recording microphones. (c) Measured emission directionality With 15° resolution based on calibrated microphone recordings, and taking transmission losses of sound on the way from the bat to the microphone for the temperature and humidity [22] at the time of recording into account. Flight direction is to the top. Emission amplitude expressed in dB is normalized to the highest observed value. Thin solid lines: maximum and minimum relative amplitude. Bold solid line: mean relative amplitude. Bold dotted line: cosine function fitted to relative emission amplitudes. A is the parameter of the fitted cosine function, namely A[cos(C) — 1], resulting in

Two polar histograms showing all the recorded headings of the individually flying bats (a) and the paired flying bats (b) captured during the experiment. 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°, (i33.1°), for the individual flying bats and 348.0° and 168.0°, (i36.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)

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 1/6. The effect of changing vc 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 3 VC < 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 vc. Trajectory 17, 20, 40 and 64 were depicted, respectively, in 82 Fig. panels (a) and (b), and

A schematic correlation map demonstrating the effect of trajectory curvature on delay extraction. Within the plot the correlations are distinguished by those below a selected threshold (whitespace) and those above the threshold (shaded area), separated by the correlation threshold itself (black-line). Two distinct flight behaviours are drawn: high correlation flight where the bats are well aligned but turn frequently, and parallel flight where the bats fly in straight lines. Times and delays have been marked to demonstrate the effect of the time ordering condition (AT/At 2 — 1) on the uncertainty of delay extraction. (TIFF)

In panels (a) and (b) two sample trajectories are drawn with their corresponding TDDC functions in panels (d) and (f), respectively. The maximum of the TDDC function averaged over all times is shown in panel (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-or-dered procedure between 1 s < t < 2 s are extracted by considering the TDDS plot (displayed in 81 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.

Panels (a) and (b) compare the distributions of turning angles for the observed unclassified and chase flights, respectively. A fit of the histograms with a von Mises distribution is also shown as well as the maximum lateral acceleration imposed on the model animals. Panels (c) and (d) show two examples of trajectories that emerge from the model in which the bats act independently and interactively. At the end of each trajectory we show the bats’ echolocation range obtained by selecting hearing threshold of 10 dB and an emission directionality with an asymmetric parameter A = 16 in Equation (11). Panels (e) and (i) show the bats’ headings corresponding to the coloured trajectories in (c) and (d), respectively. because the animals’ separation distance was at the limit of their hearing range. (TIFF)

Acknowledgments

We thank Dr Dieter Vanderelst for his comments on earlier draft of the ms and for his help on hearing directionality, Ieremy Cusack for data collection and him and Dr Holger Gorlitz for the Video analysis and trajectory generation. We also thank Bristol Water plc for support and site access.

Author Contributions

Performed the experiments: LG MH. Analyzed the data: TIM LG MH. Contributed reagents/materials/ analysis tools: MH LG. Wrote the paper: LG MH.