Self‐organizing cicada choruses respond to the local sound and light environment

Abstract Periodical cicadas exhibit an extraordinary capacity for self‐organizing spatially synchronous breeding behavior. The regular emergence of periodical cicada broods across the United States is a phenomenon of longstanding public and scientific interest, as the cicadas of each brood emerge in huge numbers and briefly dominate their ecosystem. During the emergence, the 17‐year periodical cicada species Magicicada cassini is found to form synchronized choruses, and we investigated their chorusing behavior from the standpoint of spatial synchrony. Cicada choruses were observed to form in trees, calling regularly every five seconds. In order to determine the limits of this self‐organizing behavior, we set out to quantify the spatial synchronization between cicada call choruses in different trees, and how and why this varies in space and time. We performed 20 simultaneous recordings in Clinton State Park, Kansas, in June 2015 (Brood IV), with a team of citizen‐science volunteers using consumer equipment (smartphones). We use a wavelet approach to show in detail how spatially synchronous, self‐organized chorusing varies across the forest. We show how conditions that increase the strength of audio interactions between cicadas also increase the spatial synchrony of their chorusing. Higher forest canopy light levels increase cicada activity, corresponding to faster and higher‐amplitude chorus cycling and to greater synchrony of cycles across space. We implemented a relaxation‐oscillator‐ensemble model of interacting cicadas, finding that a tendency to call more often, driven by light levels, results in all these effects. Results demonstrate how the capacity to self‐organize in ecology depends sensitively on environmental conditions. Spatially correlated modulation of cycling rate by an external driver can also promote self‐organization of phase synchrony.


SSound data acquisition and processing
We began with 26 devices, mostly smartphones, optimised for voice recording using a variety of recording software with a typical sampling rate of 44100 Hz. Sound waveforms were either recorded as .wav files or converted to .wav prior to further analysis. While approaching the main recording site, we activated two recording devices ahead of time and placed one in a woodland patch discontiguous from the main recording site, as a long baseline reference. This off-site recording was made at a distance of 2960m east of location 8 (GPS coordinates 38.94244,, and its transform is presented in Fig. S25 for comparison with the other recording locations. To ensure our ability to compare the cycle phase of cicadas at two widely separate locations with exact timing, we activated the devices and ensured that both recorded a single preliminary reference event (a balloon pop) and recorded continuously thereafter. Recording devices were scattered evenly along forest paths and placed at ground level, then left undisturbed for the duration of measurement.
The recording of a secondary reference event on the devices after they had been collected again enabled us to verify later that all devices had recorded the same length of time between reference events and that their timings were thus comparable. One older device was rejected at this stage for 'losing time'. Four other devices were rejected for losing data due to equipment failure, logistical problems, or human error. Twenty useable locations were identified for further analysis.
We evaluated the sound spectral profile of the cicadas at each location we recorded, and applied adaptive non-causal filters to isolate the sound content associated with the cicadas, as follows. Audio data was divided into 0.1 second segments and D R A F T the sound power spectrum was found within each. We trimmed off data from outside the period when all recording devices were active and undisturbed. Each segment was Fourier transformed to identify the sound frequencies present. The total power inside a band (2 to 7 kHz) that matched the cicada spectrum was evaluated for each segment (this is called 'band-pass filtering'). We found the mean sound power spectrum (by Fast Fourier Transform) of all segments, and took the associated amplitude of each frequency in the range associated with the cicadas to make a filter profile. Mean sound power spectra are shown in fig.S1. Next, for every segment we multiplied the Fourier transform coefficients at these frequencies by this profile, prior to calculating its sound power. This multiplication exaggerates the contribution due to the cicadas at the expense of other intermittent sounds (this is called 'matched filtering'). Investigating cicada sound production, Nahirney et al. (2006) reported the characteristic frequency to be 7.1-7.4 kHz, higher than the spectral peak found in our cicada recordings. The resonant frequency of the resonant cavity (which amplifies the sound) may not match the tymbal frequency (generated sound), which might explain why the actual sound output of our cicadas does not precisely match these previously reported values.
The cicada sound volume values of the segments formed a time series which we then wavelet-transformed. Each time series of cicada sound volume was de-meaned and normalized to have unit variance, then wavelet transformed with a complex mother wavelet having a center frequency of f0 = 4 (Addison, 2002). Individual transforms of data before trimming, after trimming, after 2 to 7 kHz band-pass filtering, and after matched filtering are included for all recording locations at the end of this document, showing the effect of the cleaning process on the wavelet transform in each case (Figs S2-S21 and S25). A diagram showing the steps in data cleaning is presented in Fig. S26, and the wavelet transform results (shown in Figs S2-S21d) are presented in one combined figure in Fig. S27.
Relative amplitude and phase values of wavelet transforms for a particular point in time are plotted for each location in figure S22, and the evolution of these values over time is shown in forestclip.mp4. The mean sound power spectra of the locations are compared in figure S1. The mean wavelet energy plot over all locations was then obtained, and the average cicada cycling frequency in the forest was found at each point in time (the maximum in this plot for each time). For each location, at each point in time, we extracted a wavelet value giving the magnitude and phase at this variable frequency.

SCanopy light level
The sky was overcast and varied slowly in brightness over time, although the recording contains several sharp discontinuities associated with the video function automatically adjusting exposure. We took the mean pixel brightness of the recorded image as an indicator of sky brightness for the whole forest. Discontinuities were identified and removed by hand, with the segments on either side moved back into alignment. This produced a smoothly varying light curve, with some variability due to wind and other canopy disturbance. Although the field of view was partly obstructed by the canopy overhead at this observation location, the area of sky seen through the gaps in the canopy corresponds to a wide field of view -a patch of cloudy sky sufficient to cover the ground area of our observations. This approach was vindicated by the mean cicada activity being found to correlate with the sky brightness measured at this central location.
We tested the correlation coefficient between the light curve and the changing average cycling period of the whole forest. To evaluate statistical significance we generated artificial 'surrogate' light curves by rearranging the values from the light curve in random ways, so as to approximately preserve the spectral characteristics and autocorrelation function of the light curve, using an amplitude adjusted Fourier transform method (Schreiber & Schmitz, 2000), and then repeating the moving average (Methods) to get a curve of comparable autocorrelation. We also tested for a correlation with average amplitude of cycling in the same way.

STesting for 'overhearing'
To check that observed phase coherence was due to spatial synchrony between the cicadas in the neighborhood of each recording location, and was not attributable to microphones 'overhearing' cicadas from neighboring locations, we checked the most phase-coherent locations, 4 and 11, and examined temporal variability of cycling amplitude at each location along with temporal variability in phase agreement between these locations. We selected these sites as being close, with high apparent synchrony of sound cycling. As usual, we are limited by the use of cellphone recorder hardware and software that does not give us an absolute measure of sound energy, but can be used to follow relative changes, sufficient to track the 5 second cycling and changes in the amplitude of such cycles. In this check we suppose that most of the sound energy being recorded at each site originates 'locally', ie. in the immediate vicinity of the microphone, and some small fraction (perhaps zero) originates at the site ostensibly monitored by the other microphone. Both sites are found to vary over time in terms of the volume of sound and the amplitude of cycle oscillations which are recorded, and this variation is at least partially independent, indicating that at least some of the sound and the cycles originate locally at each site. When two independently varying (unsynchronised) signals are combined in this way, the correlation between the signal mixtures is maximised when the amplitude of the signal originating at one location is maximised and the other is minimised. This does not occur if the signals are not independent (synchronised). Thus in the case that apparent synchrony is due to microphone crosstalk with independent local sound sources, we would D R A F T expect the cross talk to produce maximum apparent synchrony when the independent local cycles are 'quiet' and the distant source 'loud', with the local sound being 'drowned out' by the sound from the alternate source. The limiting case would be local silence, and the resultant (relatively quiet) recording being 100 percent due to sound from the distant source, giving 100 percent synchrony. This effect is not observed at either site. When either of the local cycling powers drops below average, the synchrony level drops below the level found when both are above average. This is in accordance with our simulations in which high amplitude cycling more easily demonstrates (true) synchronisation.If phase coherence was attributable to 'overhearing', we would expect phase agreement to be strongest when cicadas at one location were producing large-amplitude volume fluctuations and those at the other location were not, resulting in recorded cycling in the second location being more dominated by the strong fluctuations at the first location. Defining phase agreement at time t as Re(e i(φ 4,t −φ 11,t ) ), where φ4,t is the phase of volume oscillations at location 4 and φ11,t is the phase at 11, we divided each set of transform values into the half for which the amplitude was above average and the half for which it was below average. We checked the phase agreement values obtained when: location 4 had above average cycling amplitude and location 11 had above average cycling amplitude (Fig. S23a); location 4 had above average cycling amplitude and location 11 had below average cycling amplitude (Fig. S23b); location 4 had below average cycling amplitude and location 11 had above average cycling amplitude (Fig. S23c); and location 4 had below average cycling amplitude and location 11 had below average cycling amplitude (Fig. S23d). The highest agreement was found for Fig. S23a, indicating that synchrony was greatest when cicada volume cycling at both locations was high. Agreement was lower for Fig.  S23b, indicating that the observed agreement was not due to the cicada sound from location 4 being recorded at location 11. Similarly, agreement was lower for

SAn entropy-based measure of phase clumping
We applied an entropy-based measure of non-uniformity in the whole-forest phase distribution. It is based on the approach of Kraskov et al. (2004) to calculating the entropy of a distribution (of, in our case, phase values) by ordering N observed values and working with their first differences. The first differences between the sorted phase values sum to 1 turn around the circle and are small inside clusters of similar phases, large outside. The entropy of the distribution of phases at a moment in time in the forest is [1] Approximating the probability density f (φ) between (sorted) observed values φn and φn+1 as 1 N (φ n+1 −φn) (equivalent to assuming a 1 in N probability of a value falling between φn and φn+1) we can approximate entropy as Here the summation is over n = 1, . . . , N and φN+1 is interpreted as referring to φ1, given the circular nature of phases. Phase differences are computed modulo 2π. In principle, this approach could fail as phase values of locations grow at different rates over time and sometimes overtake each other, implying infinite probability density when two values become equal. In practice, none of our phase values are exactly equal but the entropy estimate is 'spiky', including many transient peaks (Fig. S24). Another approach would be to base the entropy estimate on the spacing between k th nearest neighbors (Singh et al., 2003), but this is inappropriate for our weakly synchronous system where we see only pairs and small clusters of locations with high synchrony. Evaluating entropy based on the k th nearest neighbor would miss clusters of size less than k. We calculated the entropy estimate for each 0.1 seconds using the phase values drawn from the wavelet transforms for each 0.1 seconds, and then we smoothed by a moving average over 10 minutes.

SSpatially synchronous models
Our model cicadas were encoded by variables representing how far through the calling or not-calling period each cicada had reached at a point in time. When calling, each cicada would count down at a fixed rate until calling ceased, with a total calling time of 1.5 seconds. When not calling, each cicada would count up, at a variable rate, to a total that could depend on the volume of sound produced by other cicadas at that moment in time; then it would recommence calling. We modeled 25 trees, each with 100 cicadas. The cicadas were simulated over 10,000 seconds with a time step of tstep=0.05 seconds. The first state variable of cicada m, ξm(t), was non-zero only if the cicada was calling, and counted down from 1 to 0 over 1.5 seconds until the cicada stopped calling. The second variable, ψm(t), was a phase variable that accumulated from 0 to 1 when the cicada was not calling. This variable increased by a variable increment each time step. The increment was equal to a constant c plus the sum of a randomly determined fixed increment which was unique to each cicada and a randomly determined temporally varying increment: ψm(t + 0.01) = ψm(t) + c + 0.0002 m + 0.0075υm,t. [4]

D R A F T
Here c = 0.01 and and υ were both normally distributed random numbers with mean zero and variance 1. One value of m was produced, once and for all, for each cicada; this was done independently for different cicadas. The υm,t were independently produced for different values of m or t. When ψm reached a threshold of 1 full cycle, ψm = 1, it was reset to zero while the cicada called for 1.5 seconds. Thus an average cicada will call for 1.5 seconds exactly then rest for approximately 5 seconds, then repeat. Variability (and asynchrony) in calling rate between cicadas results from variable resting time.
With no interactions between the cicadas, we used a burn-in to reach a stable distribution of calling/not calling cicadas to use as a start point. For the first part of the simulation (0 to 2500 seconds), the cicadas called randomly.
We then introduced a series of interactions, which caused cicadas to call, and ψ to reset to zero, prior to reaching one full cycle. The sound of other cicadas was taken to reduce the calling threshold by an amount proportional to sound volume; cicadas could now hear other cicada calls, encouraging them to call too. In each tree, the sound volume level, Vi, was set equal to the fraction of cicadas that were calling in that tree, vi. This sound stimulus reduced the threshold (previously at ψ = 1) for a cicada in that tree to resume calling, by 0.2Vi. This change caused the cicadas in a tree to become synchronous. Thus the second quarter of the simulation (2500 to 5000 seconds) was dominated by periodic cycling of cicada call volume in each tree, with no relationship between trees because cicadas could not yet hear other cicadas in other trees.
For the third part of the simulation (5000 to 7500 seconds) we introduced interactions between locations. Trees were taken to be spaced in a 5-by-5 grid with one unit of spacing between nearest neighbors, so the maximum inter-tree distance was √ 32 units. In each tree the sound volume Vi was increased by a sum of terms j uij proportional to the fraction of cicadas calling in every other tree vj, subject to an inverse square law which produced a 32-fold reduction in sound between the closest and the most distant trees, and including a term to represent attenuation of sound by absorption within the forest, with a coefficient k = 0.1 (set to produce a further approximately 2 times reduction in sound between the most distant trees). Thus the addition to volume Vi at tree i due to cicadas in tree j was uij = 0.04vje −kd ij /d 2 ij where dij is the distance (in units as above) between trees i and j. The cycling in all the trees now synchronized, with the peaks and troughs occurring at the same time in all trees.
For the final part of the simulation (7500 to 10000 seconds), we slightly increased the mean cycling rate of all the cicadas to c = 0.0125, simulating their common response to an environmental change. This increase in rate produced an increase in mean volume (as less time was spent not calling), cycling amplitude and spatial synchrony between the trees.
To investigate the implications of a higher mean cycling rate for the amplitude and spatial synchrony of the cycles, we ran the simulation with a range of periods (tstep/c) for the 'not-calling' part of the cycle, from 3.5 to 5.5 seconds, with the result that faster calling produced systematically greater synchrony, as observed in our empirical data.
As described in methods, three alternate simulations were run in which the component of sound volume arising from other trees was adjusted at each tree to represent: 1. Extreme attenuation; the components of volume uij due to cicadas in trees more than 1 unit away were set to zero. 2. 'Realistic' signal delay; the components of volume due to cicadas in other trees were systematically delayed (by the use of past rather than current values), assuming trees were 75m apart and the speed of sound is 330m/s. 3. Extreme signal delay; the delay was taken to be 10 times the 'realistic' value. cassini) in 0.1s windows during the period all devices were in place in the forest. Volume between 2000 and 7000 Hz was obtained by band-pass filtering as described in SI Appendix, S2. d) Same as c, but volume between 2000 and 7000 Hz was obtained by matched filtering as described in SI Appendix, S2. cassini) in 0.1s windows during the period all devices were in place in the forest. Volume between 2000 and 7000 Hz was obtained by band-pass filtering as described in SI Appendix, S2. d) Same as c, but volume between 2000 and 7000 Hz was obtained by matched filtering as described in SI Appendix, S2. cassini) in 0.1s windows during the period all devices were in place in the forest. Volume between 2000 and 7000 Hz was obtained by band-pass filtering as described in SI Appendix, S2. d) Same as c, but volume between 2000 and 7000 Hz was obtained by matched filtering as described in SI Appendix, S2.     during the period all devices were in place in the forest. c) Log magnitude of the transform of sound volume between 2000 and 7000 Hz (the frequencies of sound emitted by M. cassini) in 0.1s windows during the period all devices were in place in the forest. Volume between 2000 and 7000 Hz was obtained by band-pass filtering as described in SI Appendix, S2. d) Same as c, but volume between 2000 and 7000 Hz was obtained by matched filtering as described in SI Appendix, S2. during the period all devices were in place in the forest. c) Log magnitude of the transform of sound volume between 2000 and 7000 Hz (the frequencies of sound emitted by M. cassini) in 0.1s windows during the period all devices were in place in the forest. Volume between 2000 and 7000 Hz was obtained by band-pass filtering as described in SI Appendix, S2. d) Same as c, but volume between 2000 and 7000 Hz was obtained by matched filtering as described in SI Appendix, S2. during the period all devices were in place in the forest. c) Log magnitude of the transform of sound volume between 2000 and 7000 Hz (the frequencies of sound emitted by M. cassini) in 0.1s windows during the period all devices were in place in the forest. Volume between 2000 and 7000 Hz was obtained by band-pass filtering as described in SI Appendix, S2. d) Same as c, but volume between 2000 and 7000 Hz was obtained by matched filtering as described in SI Appendix, S2. during the period all devices were in place in the forest. c) Log magnitude of the transform of sound volume between 2000 and 7000 Hz (the frequencies of sound emitted by M. cassini) in 0.1s windows during the period all devices were in place in the forest. Volume between 2000 and 7000 Hz was obtained by band-pass filtering as described in SI Appendix, S2. d) Same as c, but volume between 2000 and 7000 Hz was obtained by matched filtering as described in SI Appendix, S2. during the period all devices were in place in the forest. c) Log magnitude of the transform of sound volume between 2000 and 7000 Hz (the frequencies of sound emitted by M. cassini) in 0.1s windows during the period all devices were in place in the forest. Volume between 2000 and 7000 Hz was obtained by band-pass filtering as described in SI Appendix, S2. d) Same as c, but volume between 2000 and 7000 Hz was obtained by matched filtering as described in SI Appendix, S2.  is not due to 'overhearing'; it is high when both locations exhibit strong oscillations. a) Phase agreement when both locations were oscillating strongly. b) Phase agreement when location 4 was oscillating strongly and location 11 was not. c) Phase agreement when location 11 was oscillating strongly and location 4 was not. d) Phase agreement when neither location was oscillating strongly. Methods and details in SI Appendix, S4.  S24. Entropy of wavelet phases through time, plotted against 0.1, 1, 5, 50, 95, 99, 99.9 percentiles of entropy values obtained for random draws from a uniform phase distribution, with no smoothing. The percentiles are for the distribution of entropies of 20 independently uniformly distributed phases, corresponding a null hypothesis of independent cicada phases at the different recording locations. The empirical entropy is typically below the median (50th percentile, at approximately 1.3). Fig. S25. Wavelet transform showing strong cycling of volume for long baseline reference recording device, varying around a period of approx five seconds. Note this recording ends before the end of the usable period of the other devices presented here, so transforms of two filtered versions of the signal are presented for all available data. Colour indicates log magnitude of fluctuation with a particular periodicity as this varies over time. a) Log magnitude of transform of volume in 0.1s windows through whole signal, starting from beginning of recording. b) Log magnitude of transform of available data: average FFT magnitude from 2000 to 7000 Hz (band-pass) in 0.1s windows starting from main pop. c) Log magnitude of transform of available data: average FFT magnitude mulplied by spectral filter in 0.1s windows starting from main pop. d) is the log magnitude of average transform of sites at the forest location for comparison, 3km away. Delete segments outside useable time range

D R A F T
As a check, perform step X below: shown in Fig S2-S21a As a check, perform step X below: shown in Fig S2-S21b As a check, perform step X below: shown in Fig S2-S21c This plot of the wavelet transform of the trimmed and filtered cicada volume timeseries is shown in Fig S2-S21d Extract time-varying 'forest frequency' from the average local maximum shown in Fig. 1, then extract wavelet transform values from all the sites at this frequency to obtain inter-site phase synchrony results.