Study species
Seven adult Japanese greater horseshoe bats (Rhinolophus ferrumequinum nippon, three males and four females) were used for the behavioral experiment. The bats of this species emit pulses consisting of a short initial FM component, a constant frequency component (CF), and a terminal FM component. These pulses are accompanied by overtones. Among the overtones, the second overtone, which has a CF component of approximately 68–70 kHz, has the highest sound pressure [11].
Experimental setup and procedure
The experiment was conducted in a corridor [4.5 m (L) × 1.5 m (W)] bordered by chains within a flight chamber [9 m (L) × 4.5 m (W) × 2.5 m (H)]. Three acrylic panels [1 m (W) × 2 m (H)] were placed in the corridor next to the chain-walls, alternating on the left and right sides, with each panel being spaced 1 m apart from each other (Fig. 5a). The seven bats were allowed to fly one after the other from the starting position through the obstacle course to a net behind the starting position 12 times. At the time of their first flight, they were completely unfamiliar with the obstacle course. In addition, the experimenter carried each bat to the starting position of the flight while covering the bat with his hands to prevent it from collecting information about its surroundings. After a bat had flown through the course once, the experimenter immediately captured the bat using a 50 × 50-cm insect net. The bat was then allowed to drink water from a plastic pipette and was brought back to the starting point by the experimenter who was covering them with their hands while carrying. During the experiment, only infrared lights were used in the flight chamber. Two high-speed video cameras (MotionPro X3; IDT Japan, Inc., Tokyo, Japan; 125 frames per second) recorded the flight paths of the bats, and microphones arranged in an array around the flight chamber recorded the emitted pulses of the bats to calculate the direction of pulse emission. We also calculated the pulse emission timing of the bats by measuring the pulses they emitted using a telemetry microphone attached to the bat’s head. Please refer to the publication of Yamada et al. [14] for further methodological details.
Finite-difference time-domain (FDTD) method
We used the two-dimensional FDTD method to simulate the echoes that return to the positions of the right and left ears of bats during their flight. The FDTD method is one of the numerical methods developed by Yee to solve Maxwell’s equations, which are the governing equations of electromagnetic waves [41]. It has been widely used in the field of acoustics since Madariaga adapted it for elastic waves [42]. The method uses the governing equations of motion and the continuity of sound pressure as given in Eqs. (1) and (2).
$$\frac{{\boldsymbol\partial}{\boldsymbol p}} {{\boldsymbol\partial}{\boldsymbol t}}+{\boldsymbol\rho}{{\boldsymbol c}_{\mathbf0}}^{\mathbf2}{\boldsymbol\nabla}\cdot{\boldsymbol u}=\mathbf0$$
(1)
$${\displaystyle \begin{array}{c}\frac{\boldsymbol{\partial u}}{\boldsymbol{\partial t}}+\frac{\mathbf{1}}{\boldsymbol{\rho}}\mathbf{\nabla}\boldsymbol{p}=\mathbf{0}\end{array}}$$
(2)
where p is the sound pressure, u is the particle velocity vector, ρ is the medium density, and c0 is the speed of sound. The following equations are obtained by discretizing the above governing equations on a staggered grid.
$${\displaystyle \begin{array}{c}{\boldsymbol{p}}_{\boldsymbol{i},\boldsymbol{j}}^{\boldsymbol{n}+\mathbf{1}}={\boldsymbol{p}}_{\boldsymbol{i},\boldsymbol{j}}^{\boldsymbol{n}}-\frac{\boldsymbol{\rho} {{\boldsymbol{c}}_{\mathbf{0}}}^{\mathbf{2}}\Delta \boldsymbol{t}}{\Delta }\left({\boldsymbol{u}}_{{\boldsymbol{x}}_{\boldsymbol{i}+\frac{\mathbf{1}}{\mathbf{2}},\boldsymbol{j}}}^{\boldsymbol{n}+\frac{\mathbf{1}}{\mathbf{2}}}-{\boldsymbol{u}}_{{\boldsymbol{x}}_{\boldsymbol{i}-\frac{\mathbf{1}}{\mathbf{2}},\boldsymbol{j}}}^{\boldsymbol{n}+\frac{\mathbf{1}}{\mathbf{2}}}+{\boldsymbol{u}}_{{\boldsymbol{y}}_{\boldsymbol{i},\boldsymbol{j}+\frac{\mathbf{1}}{\mathbf{2}}}}^{\boldsymbol{n}+\frac{\mathbf{1}}{\mathbf{2}}}-{\boldsymbol{u}}_{{\boldsymbol{y}}_{\boldsymbol{i},\boldsymbol{j}-\frac{\mathbf{1}}{\mathbf{2}}}}^{\boldsymbol{n}+\frac{\mathbf{1}}{\mathbf{2}}}\right)\end{array}}$$
(3)
$${\displaystyle \begin{array}{c}{\boldsymbol{u}}_{{\boldsymbol{x}}_{\boldsymbol{i}+\frac{\mathbf{1}}{\mathbf{2}},\boldsymbol{j}}}^{\boldsymbol{n}+\frac{\mathbf{1}}{\mathbf{2}}}={\boldsymbol{u}}_{{\boldsymbol{x}}_{\boldsymbol{i}+\frac{\mathbf{1}}{\mathbf{2}},\boldsymbol{j}}}^{\boldsymbol{n}-\frac{\mathbf{1}}{\mathbf{2}}}-\frac{\Delta \boldsymbol{t}}{\boldsymbol{\rho} \Delta }\left({\boldsymbol{p}}_{\boldsymbol{i}+\mathbf{1},\boldsymbol{j}}^{\boldsymbol{n}}-{\boldsymbol{p}}_{\boldsymbol{i},\boldsymbol{j}}^{\boldsymbol{n}}\right)\end{array}}$$
(4)
$${\displaystyle \begin{array}{c}{\boldsymbol{u}}_{{\boldsymbol{y}}_{\boldsymbol{i}+\frac{\mathbf{1}}{\mathbf{2}},\boldsymbol{j}}}^{\boldsymbol{n}+\frac{\mathbf{1}}{\mathbf{2}}}={\boldsymbol{u}}_{{\boldsymbol{y}}_{\boldsymbol{i}+\frac{\mathbf{1}}{\mathbf{2}},\boldsymbol{j}}}^{\boldsymbol{n}-\frac{\mathbf{1}}{\mathbf{2}}}-\frac{\Delta \boldsymbol{t}}{\boldsymbol{\rho} \Delta }\left({\boldsymbol{p}}_{\boldsymbol{i}+\mathbf{1},\boldsymbol{j}}^{\boldsymbol{n}}-{\boldsymbol{p}}_{\boldsymbol{i},\boldsymbol{j}}^{\boldsymbol{n}}\right)\end{array}}$$
(5)
where ∆ is the grid interval and ∆t is the time resolution. The terms \({\boldsymbol{p}}_{\boldsymbol{i},\boldsymbol{j}}^{\boldsymbol{n}}\), \({\boldsymbol{u}}_{{\boldsymbol{x}}_{\boldsymbol{i},\boldsymbol{j}}}^{\boldsymbol{n}},\) and \({\boldsymbol{u}}_{{\boldsymbol{y}}_{\boldsymbol{i},\boldsymbol{j}}}^{\boldsymbol{n}}\) are the sound pressure and particle velocity at a position (x, y) = (i∆, j∆) and time t = n∆t, respectively. The sound propagation is calculated by alternately solving Eqs. (3), (4), and (5). In addition, because FDTD is a time-domain and not a frequency-domain acoustic simulation, diffuse attenuation is included, while absorption attenuation is not.
Visualization of echo incidence points
Simulation space setup
In the behavioral experiment, the obstacle walls (acrylic plates) were placed from the floor to the ceiling, so that the bats avoided them by changing their horizontal flight path. In other words, vertical (z-direction) flight had almost no effect on the avoidance behavior of this experimental system, and in fact, little change in altitude was observed during flight. Therefore, considering the amount of computation, we created a two-dimensional space (x-y) at the height level of array microphones (z = 1.2 m) to reproduce this setting in the simulation space. The FDTD simulation space was located inside the corridor (4.5 m × 1.5 m) as shown in Fig. 5c, and the absorption boundary condition was the Mur second-order boundary condition [43]. As presented in Additional file 4: Table S1, the Courant–Friendrichs–Lewy (CFL) number was 0.57. The spatial resolution (dx) was 0.3 mm. Moreover, the densities of air and acrylic were 1.29 and 1.18 kg/m3, with a bulk modulus of 142.0 × 103 and 8.79 × 109 Pa, respectively.
Echo simulation
To simulate the echoes that returned to the position of the bats’ ears, we used information on the bats’ pulse emission position and direction obtained from the behavioral experiment. A sinc function signal with a wide frequency band and high time resolution was used in the simulation space as the emitted pulse. It was flat up to 110 kHz and the duration was 0.073 ms (Fig. 5d). The sinc function signal was emitted from two source positions (simulating the bat’s nostrils), which were set at 1.25 mm each on the left and right side from a given pulse emission position and in the direction of emission which was both obtained from the behavioral experiment (Additional file 5: Fig. S1a; note that 1.25 mm is based on the half wavelength of the CF2 frequency of R. ferrumequinum nippon (68 kHz)) [44, 45]. The directivity of the sound source is created by two point sources that are close to each other. Receiver positions (simulating the bat’s ears) of the returning echoes in the simulation were set at 10 mm left and right from the pulse emission points (note that 10 mm is based on the distance between the ears of R. ferrumequinum nippon) (Additional file 5: Fig. S1a). In the simulation, the impulse responses (echoes) at the two receiver positions were calculated using the FDTD method. Then, the bat’s pulse was convolved with these impulse responses to obtain echoes. The bat’s pulses in the simulation were the downward FM components created up to the third harmonic, which has been reported to be used for distance discrimination [46]. The downward FM component has a duration of 2 ms and the first harmonic decreases linearly from 34 to 25 kHz, corresponding to the terminal FM component of the echolocation pulse in R. ferrumequinum nippon (Additional file 5: Fig. S1b). The first and third harmonics were set to a sound pressure level of -40 dB based on the second harmonic [47]. The sampling frequency was 2 MHz (1/dt). When pulses are emitted from the two point sources, they propagate in a forward and backward direction (Additional files 1, 2, and 3: Movies S1, S2, and S3). Therefore, in order to exclude echoes from the back side of the bat, the simulation was performed by excluding obstacles located behind the source positions.
Calculation of echo incidence points
Bats estimate the distance to an object based on the time delay between the emitted pulse and the echoes [24, 25]. In the simulation, this echo delay was calculated by cross-correlating the simulated left and right echoes with the bat’s pulses. The peak value of the autocorrelation of the pulse was normalized, and the peak time above the threshold of the cross-correlating signal was extracted for the left and right echoes (Additional file 6: Fig. S2). The threshold value was set to 0.02 to detect almost all the peaks. The left and right peak times were matched within a window of the maximum time difference between the two receiving points (2.5 mm/340 m/s = 7.4 μs), and combinations of the left and right peak times ([tr, tl]) were obtained. The echo incidence points ((xecho, yecho)) were then calculated by solving the following two elliptic equations determined by the calculated peak times using the pulse emission positions obtained in the behavioral experiments (i.e., the centers of the two source positions in the simulation) and the left and right receiver positions in the simulation space as the focal points. The pulse emission position was set to y = 0 and the direction of pulse emission was set to x = 0.
$${\displaystyle \begin{array}{c}\frac{{\left({\boldsymbol{x}}_{\boldsymbol{echo}}-{\boldsymbol{x}}_{\boldsymbol{r}\mathbf{0}}\right)}^{\mathbf{2}}}{{{\boldsymbol{a}}_{\boldsymbol{r}}}^{\mathbf{2}}}+\frac{{\left({\boldsymbol{y}}_{\boldsymbol{echo}}-{\boldsymbol{y}}_{\boldsymbol{r}\mathbf{0}}\right)}^{\mathbf{2}}}{{{\boldsymbol{b}}_{\boldsymbol{r}}}^{\mathbf{2}}}=\mathbf{1}\end{array}}$$
(6)
$$\begin{array}{l}\frac{\left({\boldsymbol x}_{\boldsymbol{echo}}-{\boldsymbol x}_{\boldsymbol l\mathbf0}\right)^{\mathbf2}}{{\boldsymbol a}_{\boldsymbol l}^{\mathbf2}}+\frac{\left({\boldsymbol y}_{\boldsymbol{echo}}-{\boldsymbol y}_{\boldsymbol l\mathbf0}\right)^{\mathbf2}}{{\boldsymbol b}_{\boldsymbol l}^{\mathbf2}}=\mathbf1\\{\boldsymbol a}_{\boldsymbol r}=\frac{{\boldsymbol t}_{\boldsymbol r}\cdot{\boldsymbol c}_{\mathbf0}}{\mathbf2},{\boldsymbol b}_{\boldsymbol r}=\sqrt{{\boldsymbol a}_{\boldsymbol r}^{\mathbf2}-\left({\boldsymbol x}_{\boldsymbol r}-{\boldsymbol x}_{\boldsymbol r\mathbf0}\right)^{\mathbf2}}\\{\boldsymbol a}_{\boldsymbol l}=\frac{{\boldsymbol t}_{\boldsymbol l}\cdot{\boldsymbol c}_{\mathbf0}}{\mathbf2},{\boldsymbol b}_{\boldsymbol l}=\sqrt{{\boldsymbol a}_{\boldsymbol l}^{\mathbf2}-\left({\boldsymbol x}_{\boldsymbol l}-{\boldsymbol x}_{\boldsymbol l\mathbf0}\right)^{\mathbf2}}\end{array}$$
(7)
where xr0 is the center position between the pulse emission position and the right receiver position, and xl0 is the center position between the pulse emission position and the left receiver position (see Additional file 7: Fig. S3). An animation of the process of echo incidence point visualization using acoustic simulation is shown in the Additional file 3: Movie S3.
Statistical analysis
Effect of spatial learning on the echo incidence point distribution
We were interested in testing whether the positions of echo incidence points on the obstacle walls changed depending on the spatial learning status of the bats (first vs. last flight). Since the bat’s task in the behavioral experiment was to avoid the inner sides of the obstacle walls, we considered those echo incidence points that were located on the inner half of the walls (n = 1058) and calculated the distance of these points to the inner edge. We modeled this data as a function of the bat’s spatial learning status (first vs. last flight) using generalized linear mixed effect models (function glmmTMB, package glmmTMB_1.0.2.1) [48] assuming a negative binomial error distribution (nbinom1) due to overdispersion. Because several echo incidence points can result from one pulse and several pulses were used per bat, we included a random effect with a pulse-ID nested within the bat-ID. The quality of the model fit was graphically examined (function in package DHARMa_0.3.3.0) [49] and its overall significance was determined by comparing it to the respective null model that contains only the random effect via a χ2-test (function anova in stats) [50]. The significance of the factor coding for the bat’s experience was derived from a type-II-Wald- χ2-test (function anova, package car_3.0-10) [51] while the factor-levels were compared based on least-square-means (function lsmeans, package emmeans_1.5.4) [52]. The statistical significance level was set at p = 0.05.
Effects of spatial learning on flight path planning
We calculated the turn rate at 1 ms intervals from the acquired flight paths. The turn rate is the time derivative of the flight path. To investigate the relationship between the turn rate and the pulse and echo direction, respectively, we shifted the turn rate data by a time lag of τ in 10 ms steps from − 100 ms (to the left) up to + 600 ms (to the right) relative to the pulse and echo direction, respectively, and calculated the corresponding correlation coefficients. The 95% confidence intervals for the correlation coefficients were determined by a Fisher transformation of the correlation coefficients to the following ranges:
$${\displaystyle \begin{array}{c}\pm {\boldsymbol{z}}_{\boldsymbol{\alpha} /\mathbf{2}}\frac{\mathbf{1}}{\sqrt{\boldsymbol{n}-\mathbf{3}}}\end{array}}$$
(8)
where n represents the number of data and zα/2 is 1.96 for 95% confidence interval.
Then, we extracted the τ-values that were associated with the highest correlation coefficients for each bat and each category (first vs. last flight; n = 28) and modeled them on the scale of seconds as a function of the degree of spatial learning (first vs. last flight) in interaction with the factor describing the type of information used by the bat (echo vs. pulse) using linear mixed effect models (function lmer, package lme4_1.1-26) [53]. We added the bat-ID as a random effect to the model due to the repeated sampling of the same individuals. The quality of the model fit, the significance of factors within the model, and the comparisons between factor-levels were conducted using the same functions as mentioned above. The overall model significance was tested against the respective null model using parametric bootstrapping (function PBmodcomp, package pbkrtest_0.5-1.0) [54].
