The dynamical complexity of seasonal soundscapes is governed by fish chorusing
Data collectionThe acoustic recordings were collected during 2017 off the Changhua coast (24° 4.283 N/120° 19.102 E) (Fig. 5) by deploying a passive acoustic monitoring (PAM) device from Wildlife Acoustics, which was an SM3M recorder moored at a depth of 18–20 m. The hydrophone recorded continuously with a sampling frequency of 48 kHz and a sensitivity of −164.2 dB re:1 v/µPa. The acoustic files were recorded in the.WAV format with a duration of 60 minutes. The hydrophone setup prior to deployment is shown in Fig. 6. Table 2 contains the details for the monitoring period with the corresponding season and the number of hours of recordings each season used in this study. Studies have shown that the presence of seasonal chorusing at this monitoring site in the frequency range of 500–2500 Hz is caused by two types of chorusing15,38, with chorusing starting in early spring, reaching a peak in summer, and starting to decline late autumn, and silencing in winter6. Previous studies6,15,38 at this monitoring site have derived the details of two types of fish calls responsible for chorusing (Type 1 and Type 2); Supplementary Fig. 1 shows the spectrogram, waveform, and power spectrum density of the individual calls. Supplementary Table 1 tabulated are the acoustic features of the two call types. The monitoring region, Changhua, lies in the Eastern Taiwan Strait (ETS). The ETS is ~350 km in length and ~180 km wide64. The ETS experiences diverse oceanographic and climatic variations influenced by monsoons in summer and winter65 and extreme events caused by tropical storms, typhoons in summer, and wind/cold bursts occurring in winter66,67,68.Fig. 5: Study area located off the Taiwan Strait.Map of the Changhua coast located in Taiwan Strait, Taiwan depicting the deployed passive acoustic monitoring recorder at site A1. The map was produced in Matlab 9.11 (The Mathworks, Natick, MA; http://www.mathworks.com/) using mapping toolbox function geobasemap(). Full global basemap composed of high-resolution satellite imagery hosted by Esri (https://www.esri.com/).Full size imageFig. 6: Setup of the SM3M submersible recorder.SM3M recorder fastened to the steel frame (length and breadth = 1.22 m, height = 0.52 m) with plastic cable zip ties prior to deployment.Full size imageTable 2 Passive acoustic monitoring device specifications and monitoring duration during different seasons.Full size tableAcoustic data analysisThe acoustic data were analyzed using the PAMGuide toolbox in Matlab60. The seasonal spectrograms were computed with an FFT size of 1024 points and a 1 s time segment averaged to a 60 s resolution. The sound pressure levels (SPL) were computed in the frequency band of 500–3500 Hz and programmed to provide a single value every hour, thus resulting in 984, 1344, and 1440 data points in spring, summer, and winter, respectively (Supplementary Data 1).Determining the regularity and complexity with the complexity-entropy planeThe complexity-entropy plane was utilized in this study to quantify the structural statistical complexity and the regularity in the hourly acoustical and seasonal SPL time series data. The C-H plane is a 2D plane representation of the permutation entropy on the horizontal axis that quantifies the regularity, and the vertical axis is represented by the statistical complexity quantifying the correlation structure in the temporal series.For a given time series ({{x(t)}}_{t=1}^{N}), the N’ ≡ N − (m − 1) the values of the vectors for the length m > 1 are ranked as$${X}_{s}=left({x}_{s-(m-1)},{x}_{s-(m-2)},ldots ,{x}_{s}right),s=1,ldots ,,{N}^{{prime} }$$
(1)
Within each vector, the values are reordered in the ascending order of their amplitude, yielding the set of ordering symbols (patterns) ({r}_{0},{r}_{1},ldots ,{r}_{m-1}) such that$${x}_{s-{r}_{0}}le {x}_{s-{r}_{1}}le ..,..le {x}_{s-{r}_{(m-1)}}$$
(2)
This symbolization scheme was introduced by Bandt and Pompe69. The scheme performs the local ordering of a time series to construct a probability mass function (PMF) of the ordinal patterns of the vector. The corresponding vectors (pi ={r}_{0},{r}_{1},ldots ,{r}_{(m-1)}) may presume any of the m! possible permutations of the set ({{{{{mathrm{0,1}}}}},ldots ,m-1}) and symbolically represent the original vector. For instance, for a given time series {9, 4, 5, 6, 1,…} with length m = 3, provides 3! different order patterns with six mutually exclusive permutation symbols are considered. The first three-dimensional vector is (9, 4, 5), following the Eq. (1), this vector corresponds to ((,{x}_{s-2},{x}_{s-1},{x}_{s})). According to Eq. (2), it yields ({x}_{s-1}le {x}_{s}le {x}_{s-2}). Then, the ordinal pattern satisfying the Eq. (2) will be (1, 0, 2). The second 3-dimensional vector is (4, 5, 6), and (2, 1, 0) will be its associated permutation, and so on.The permutation entropy (H) of order m ≥ 2 is defined as the Shannon entropy of the Brandt-Pompe probability distribution p(π)69$$Hleft(mright)=,-{mathop{sum}limits _{{pi }}}pleft(pi right){{{log }}}_{2}p(pi )$$
(3)
where ({pi }) represents the summation over all possible m! permutations of order m, (p(pi )) is the relative frequency of each permutation (pi), and the binary logarithm (base of 2) is evaluated to quantify the entropy in bits. H(m) attains the maximum ({{log }}m!) for (p(pi )=1/m!). Then the normalized Shannon entropy is given by$$0le H(m)/{{{{{rm{ln}}}}}},m!le 1$$
(4)
where the lower bound H = 0 corresponds to more predictable signals with fewer fluctuations, an either strictly increasing or decreasing series (with a single permutation), and the upper bound H = 1 corresponds to an unpredictable random series for which all the m! possible permutations are equiprobable. Thus, H quantifies the degree of disorder inherent in the time series. The choice of the pattern length m is essential for calculating the appropriate probability distribution, particularly for m, which determines the number of accessible states given by m!70,71. As a rule of thumb, the length T of the time series must satisfy the condition T (gg) m! in order to obtain reliable statistics, and for practical purposes, Bandt and Pompe suggested choosing m = 3,…,7 69.The statistical complexity measure is defined with the product form as a function of the Bandt and Pompe probability distribution P associated with the time series. (Cleft[Pright]) is represented as33$$Cleft[Pright]=frac{J[P,U]}{{J}_{{max }}}{H}_{s}[P]$$
(5)
where ({H}_{s}left[Pright]=Hleft[Pright]/{{log }}m!) is the normalized permutation entropy. (J[P,U]) is the Jensen divergence$$Jleft[P,Uright]=left{Hleft[frac{P+U}{2}right]-frac{H[P]}{2}-frac{H[U]}{2}right}$$
(6)
which quantifies the difference between the uniform distributions U and P, and ({J}_{{max }})is the maximum possible value of (Jleft[P,Uright]) that is obtained from one of the components of P = 1, with all the other components being zero:$$Jleft[P,Uright]=-frac{1}{2}left[frac{m!+1}{m!}{{log }}left(m!+1right)-2{{log }}left(2m!right)+{{log }}(m!)right]$$
(7)
For each value of the normalized permutation entropy (0le {H}_{s}[P]le 1) there is a corresponding range of possible statistical complexity (Cleft[Pright]) values. Thus, the upper (({C}_{{max }})) and lower ((C_{{min }})) complexity bounds in the C-H plane are formed. The periodic sequences such as sine and series with increasing and decreasing (with ({H}_{s}[P]=0)) and completely random series such as white noise (for which (Jleft[P,Uright]=0) and ({H}_{s}[P]=1)) will have zero complexity. Furthermore, for each given value of the (0le {H}_{s}[P]le 1), there exists a range of possible values of the statistical complexity, ({C}_{{min }}le C[P]le {C}_{{max }}). The procedure for evaluating the complexity bounds ({C}_{{min }}) and ({C}_{{max }}) is given in Martin et al.72. We utilized the R package ‘statcomp’73 to evaluate the statistical complexity (C) and the permutation entropy (H) using the command global-complexity() for the order m = 6, and the command limit_curves(m, fun = ‘min/max’) was utilized to evaluate the complexity boundaries ({C}_{{min }}) and ({C}_{{max }}). In this study, we constructed two C-H planes: (1) C and H was computed for each hourly acoustic file during different seasons. The resulting lengths of C and H during spring, summer, and autumn-winter are similar to the number of hours in the particular season (Table 2). (2) C and H was computed every 4–5 days for the seasonal SPL. The resulting length of C and H obtained during spring was 9 points (each value of C and H for every 109 h), and in summer and autumn-winter was 12 points (each value of C and H for every 112 and 120 h).Determining predictability and dynamics (linear/nonlinear) using EDMIn this study, we utilized EDM to quantify the predictability (forecasting) and distinguish between the linear stochastic and nonlinear dynamics in the seasonal soundscape SPL. EDM involves phase-space reconstruction via delay coordinate embeddings to make forecasts and to determine the ‘predictability portrait’ of time series data40. The first step in EDM is to determine the optimal embedding dimension (E), and this is obtained using a method based on simplex projection41. The simplex projection is carried out by dividing the dataset into two equal parts, of which the first part is called the library and the other part is called the target. The library set is used to build a series of non-parametric models (known as predictors) for the one step ahead predictions for the E varying between 1 and 10. Then the model’s accuracies are determined when the model is applied to the target dataset and the prediction skill (⍴) for the actual and predicted datasets is measured. The embedding dimension with the highest predictive skill is the optimal E.For the appropriate optimal E chosen, the predictability profile of the time series data is obtained by evaluating ⍴ at Tp = 1, 2, 3, … steps ahead. The flat prediction profile of the ⍴–Tp curve indicates that the time series is purely random (low ⍴) or regularly oscillating (high ⍴). In contrast, a decreasing ⍴ as Tp increases may reject the possibility of an underlying uncorrelated stochastic process and indicate the presence of low-dimensional deterministic dynamics. However, the concern with the predictability profile is that it may exhibit predictability even if time series are purely stochastic (such as autocorrelated red noise). Hence, a nonlinear test can be performed by using S-maps (sequential locally weighted global linear maps) to distinguish between linear stochastic and nonlinear dynamics in the time series dataset by fitting a local linear map. S-maps similar to simplex projects provide the forecasts in phase-space by quantifying the degree to which points are weighted when fitting the local linear map, which is given by the nonlinear localization parameter θ. When θ = 0, the entire library set will exhibit equal weights regardless of the target prediction, which mathematically resembles the model of a linear autoregressive process. In contrast, if θ > 0, the forecasts of the library provided by the S-map depend on the local state of the target prediction, thus producing large weights, and the unique local fittings can vary in phase-space to incorporate nonlinear behavior. Consequently, if the (⍴–θ) dynamics profile shows the highest ⍴ at θ = 0 that is reduced as θ increases, it represents linear stochastic dynamics. If the ⍴ achieves the highest value at θ > 0, then the dynamics are represented by nonlinear dynamics.In this study, the R package “rEDM”74 was used to evaluate the optimal E, prediction profile (⍴–Tp), and dynamics profile (⍴–θ) for the seasonal SPL dataset. While evaluating these entities, the data points are equally into two parts, and sequentially the first half is chosen as the library set and the other as the target set. The length of the library and the target set for spring, summer, and autumn-winter are 492, 672, and 720. The command EmbedDimension() was used to determine the forecast skill for the E ranging from 1 to 10 and the optimal E with the highest forecast skill (Supplementary Fig. 2) was chosen. In this study, we found that for all seasons, the optimal E was 2. The (⍴–Tp) was evaluated for Tp varying between 1 and 100 using the command PredictInterval() and the (⍴–θ) was evaluated using the command PredictNonlinear() for θ = 0, 0.0001, 0.0003, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 0.5,0.75, 1.0, 1.5, 2, and 3 to 20.StatisticsThe nonparametric Kruskal–Wallis test, followed by post hoc Bonferroni’s multiple comparisons, was used to test differences in the seasonal H and C that were obtained directly from the hourly acoustic data during chorusing hours, as well as the H and C obtained for the seasonal SPL and the seasonal forecast skill. The statistical calculations were performed using the R package “agricolae” 75. More
