Emergent periodicity in the collective synchronous flashing of fireflies

Left alone, Photinus carolinus fireflies flash without an intrinsic period, making it uncertain when they may flash next. Yet when gathering at the mating lek in large swarms, these fireflies transition into predictability, synchronizing with their neighbors with a rhythmic periodicity. Here we propose a mechanism for emergence of synchrony and periodicity, and formulate the principle in a mathematical framework. Remarkably, with no fitting parameters, analytic predictions from this simple principle and framework agree strikingly well with data. Next, we add further sophistication to the framework using a computational model featuring groups of random oscillators via integrate-and-fire interactions controlled by a tunable parameter. This agent-based model of P. carolinus fireflies interacting in swarms of increasing density also shows quantitatively similar phenomenology and reduces to the analytic framework in the appropriate limit of the tunable parameter. We discuss our findings and note that the resulting dynamics follow the style of a decentralized follow-the-leader synchronization, where any of the randomly flashing individuals may take the role of the leader of any subsequent synchronized flash burst.


Introduction
Physical systems consisting of several interacting entities often exhibit large-scale properties which are distinct from the capabilities of each entity taken individually: this is the well-known concept of emergence. Emergence has been observed and studied in both inanimate and animate systems, including famously groups of animals (1; 2). Animal collective behavior broadly designates dynamical patterns that are unsupervised consequences of the accumulation of low-level interactions between neighboring individuals (3; 4; 5). One simple yet compelling manifestation of emergence in the natural world is in the form of firefly flash synchronization (6; 7; 8; 9; 10). For example, when su ciently many Photinus carolinus fireflies congregate into a mating swarm (lek), they start to align their flashes on the same tempo, creating a mesmerizing display which has captivated the curious mind of many. This serves to strengthen their species-specific signal and heighten the ability for conspecific males and females to identify one another (6; 11; 12). In addition to collective synchrony, a more careful examination of P. carolinus' flashing pattern further reveals another non-trivial signature: emergent periodicity. Indeed, in their natural habitat, these fireflies produce periodic bursts of flashes occurring with great regularity, with a period of roughly 12s (6; 11). Surprisingly, when put in isolation, a single firefly does not appear to show any regularity about when it emits its flash trains (8), where intervals between flash trains vary between a few seconds to a few minutes apart. How, then, can a multitude of interacting fireflies exhibit a specific frequency that does not appear to be encoded in any single one of them? This paper proposes an explanation for this emergent periodicity. Since firefly flashing is mostly guided by behavior, rather than physiology, one possible explanation is that fireflies indeed "know" their collective frequency, but simply do not show it in the absence of a cooperative and responsive crowd. While this hypothesis is hard to dismiss rigorously, we will show that there exists a simple stochastic mechanism that explains the convergence towards a common, well-defined period between flash bursts as the number N of fireflies increases.

Behavioral experiments
A P. carolinus lek in its natural habitat contains several thousands of fireflies which display a robust collective flash pattern. They flash over the course of periodic bursts separated by a few seconds of total darkness (Fig. 1A, over a few seconds). Collective bursts in the swarm have a well-defined period (peak-to-peak) of about 12s (8). One could think, then, that each individual firefly also emits flash trains with about the same time period, and that the e↵ect of visual interactions is to align these individual trains on the same tempo. In other words, the swarm could be a set of coupled oscillators converging to a common phase, as has been described in previous models (13; 14; 15; 16; 17). Crucially, however, when a single firefly is taken out of the lek and placed in a large (2m 3 ) enclosing volume visually insulated from the rest of the group, all periodicity in the occurrence of flash trains is lost. The single firefly continues to emit sporadic bursts(Fig. 1B,C), but the time between successive flash bursts varies between a few seconds and a few minutes (Fig 1B  (8)). This suggests that individual inter-burst intervals (IBI) occur at random, and may thus depend on a variety of behavioral factors. When collecting measurements from 10 di↵erent fireflies recorded for several minutes under the same conditions, we are able to outline the distribution of inter-burst intervals for a single firefly in isolation (Fig. 1D, purple). (The underlying assumption here is that all fireflies have the same distribution of inter-burst intervals). Interestingly, as the number of fireflies within the enclosing volume is increased, a regularity in the time between bursts starts to emerge. At about N =15, the distribution of inter-burst intervals becomes very similar to that observed in the natural habitat. And for N =20, it is clear that there is a very strong collective periodicity separating flash bursts of about 12s, similar to that of the undisturbed swarm flashing just outside the tent (Fig 1D,E).
3 Proposed principle of emergent synchrony and periodicity, and its analytic formulation Here we propose the following paradigm, derive its mathematical formulation, and validate its predictions against experimental data: (1) Each time a firefly has finished flashing, it waits a random time t, drawn from a distribution b(t) which is identical for all individuals, before flashing again. (2) Upon flashing, a firefly stimulates all other fireflies to also flash. (3) After flashing, each firefly resets its internal waiting time to another random t. We denote by T b the collective interburst interval, that is the time between any two successive flashes produced in the swarm. The probability Q that T b is larger than some time T ? equals the probability that no firefly has flashed before T ? , i.e.
where c(t) is the cumulative probability corresponding to b(t), since all fireflies are independent until the moment they flash. The corresponding probability density function follows directly as: Thus we have set up a mathematical framework which takes as its input the experimentally observed interburst distribution, and makes specific predictions with no fine-tuning fitting parameters. If there exists a minimal time T 0 such that b(t) = 0 for all t < T 0 , and b(T 0 + ✏) > 0 for arbitrarily small positive values of ✏, then P N (T b < T 0 ) = 0 and P N (T b > T 0 ) vanishes for large N (since the integral is smaller than 1). Therefore, the distribution converges towards a single possible value, T b = T 0 for large N .
These predictions are consistent with the intuitive result that the shortest possible interburst interval is the only one that occurs in large, fully-connected, and instantaneously-stimulated groups of fireflies. We expect such a threshold minimum time to exist owing to physiological constraints, which prevent the fireflies from flashing continuously forever without pause. Intuitively, as the number of fireflies increases, there is a greater probability that at least one of those fireflies will flash at an interval close to the minimum.
Conceptually, in the idealization that at N ! 1 this distribution converges to a Dirac delta function, which tends to make the flashing patterns perfectly periodic with no variation (see SI). However, for a finite number N of fireflies, the distribution peaks at a value greater than T 0 , and has a specific non-zero width with decreases with increasing N (see SI). These specific predictions are spectacularly borne out by the experimental data. With no fine-tuning fitting parameter, and the experimentally observed single firefly distribution (Fig. 2F) as the only input to the mathematical framework, we see an excellent match between the N -dependent experimentally observed interburst distributions and the corresponding prediction from analytic theory (Figs. 2A-D). Moreover, the corresponding sharpening of the peak of the distribution (resulting in decreasing noise) with increasing N also quantitatively matches with the trend predicted by theory -see the plot of standard deviation vs. N in panel Fig. 2E). Through these compelling matches between predictions from the theory, sans fitting parameters, and the experimental observations, we establish the validity of the proposed principle for emergent synchrony and periodicity.
Furthermore, using the analytic framework the following rigorous results can be generally proved to hold for any input single firefly distribution: As the number N of fireflies increases, along with the variance, all the moments the interburst distribution monotonically decrease. In addition, the left-most mode shifts further towards the left with increasing N until it reaches T 0 . Taken together, what these predictions show is that for any input distribution shape, we are guaranteed to get emergent periodicity and synchrony through the proposed mechanism. We have provided detailed derivations of these predictions in the supplementary section.

Computational Model
The theoretical formulation in the preceding section is built on the assumption that all fireflies immediately start flashing upon seeing any other one flash, but in practice there could be some time delay or imperfect information transfer, which could be made shorter if the firefly sees additional fireflies flashing too. The rate at which this delay is shortened in proportion to the number of flashing fireflies is given by the behavioral coupling between the fireflies, labeled . When ! 1, this limit represents the idealization derived in the theory section: the strongly correlated limit, wherein a single firefly's flashing is su cient to immediately stimulate all others to also start flashing, while = 0 represents completely non-interacting fireflies.
Peaks with a height less than y ✏ = 0.01 are replaced with y ✏ to avoid zero-divide and divide-by-zero problems. T b 1 is defined as any peak with a value less than 3.5s seconds, which can otherwise be interpreted as the "noisy" regime. T b 2 is defined as any peak with a value greater than 3.5s seconds, otherwise known as the "periodic" regime. The transition between the two regimes generally appears around = 1.0 and stabilizes around = 2.7 (see Fig S6).

Formulation
We propose a simple numerical simulations based on the mechanism previously described. Following previous computational models (18; 19; 20), we implement a group of N fireflies whose flashing dynamics is governed by charging and discharging processes which represent the time between two flashes and the duration of a flash, respectively. These processes are determined by both an agent's internal characteristics and its interactions with the group. Specifically, the internal state of firefly i is characterized by variables V and ✏ whose evolution follows: which is a standard equation for the integrate-and-fire (IF) scheme. The firefly flashes and sets ✏ i = 1 when V = 1 for a duration T di , and is dark while charging (✏ i = 0) for time T si = T bi T di , where T bi represents the start-to-start inter-flash interval for firefly i, drawn directly from the input distribution envelope in Fig. 2F. However, the firefly may be "pulled" faster towards flashing if neighboring fireflies are flashing as well, which is represented by the third term, where ij = {0, 1} represents connectivity and ij is the coupling strength.
For simplicity, here we use all-to-all connectivity ( ij = 1, 8(i, j)) and vary the common interaction ij = N . The crucial di↵erence with prior IF implementations is that T b is a random variable whose value is drawn from our experimental distributions of interburst intervals (see Fig. S4) and resets, for each agent, every time they individually switch state. In this stochastic IF model the variability between flashes is accounted for, while maintaining the overall structure of the IF model.

Transition to Periodicity
This model exhibits a transition to group periodicity as interactions between agents are increased. We define the group interburst interval as the time between one flash and the next flash produced by any other firefly in the swarm. For example, consider the case of N = 20 (Fig. 3D). When = 0 each firefly behaves purely individually and interburst intervals tend to aggregate towards small values due to the random unsynchronized flashing of the N fireflies each with a flashing behavior typical of isolated individuals. This remains the case until the coupling strength, , becomes large enough that there is enough of collective entrainment to align the flashes of the group. This phenomenon starts emerging around = 1, and strengthens until a plateau at stable periodicity where ⇡ 2.7. In these regimes, when one firefly flashes it quickly triggers all others. All agents then reset their charging time at roughly the same moment, and the smallest T b defines the duration between this flash and the next. As a consequence, interburst intervals of the collective, T b , shift to a larger value corresponding to the smallest time between flashes for an individual firefly (t b0 ). We quantify this transition by examining the characteristic peaks in the T b distribution: a peak that occurs before the transition marked at T 1 b with corresponding probability P (T 1 b ), and similarity, a peak occurring after the transition marked at T 2 b with corresponding probability P (T 2 b ). The ratio of these probabilities, P (T 2 b )/P (T 1 b ) is low before the transition, and experiences a sharp increase at = 1 (Fig.3E). See SI Section 8.3.2 for full definitions. The high-coupling peak is also naturally sharper at increasing N : at larger N , the probability that some T c,i approaches the minimum possible T b is higher, resulting in more regularity the collective flashing pattern (Fig. 2F, and Fig. S5).
As our model has a single open parameter, namely the coupling strength , we conduct a detailed comparison of the model and experimental data to infer the most likely value of for the firefly system. A systematic parameter sweep over the values of and N provide a set of T b interval distributions. We use the slope of the relationship between P (T b 2 and P (T b 1 to identify the beta value at which the system stabilizes to emergent periodicity ( Fig. 3E and Fig. S5). These chosen values are indicated in red on Fig. 3A-D, and their corresponding T b distributions are shown in Fig. 2A-D.

Discussion and concluding remarks
In this work, we have proposed a synchronization mechanism that produces emergent periodicity and demonstrated its remarkable quantitative applicability to the synchronous periodic flashing of fireflies as observed in natural settings. In systems following our principle, individuals may behave erratically without any periodicity in their behavior, yet when brought together as a collective, their behavioral patterns become highly synchronized and periodic. Moreover, this e↵ect increases with the number of fireflies present through a simple and intuitive behavioral pattern. Using this principle, we successfully predict the qualitative sharpening of the peak of the distribution of interval between flashes by simply using the interval between flashes of isolated individual fireflies and without requiring any fitting parameter. Further, our computational model quantitatively builds on the predictions of the theory by letting the strength of coupling between fireflies vary and provides added insights.
Specifically, we have shown that the simple behavioral model presented in this paper successfully reproduces the experimental distributions of inter-burst intervals for groups of N fireflies (Figure 2, A-D). All the input parameters for the model come directly from experimental results in Ref. (8) and subsequent unpublished field season results from the Great Smoky Mountains: the wide distribution of inter-burst intervals for single isolated fireflies, the two model timescales of charging time and discharging time are both data-driven from Ref. (8). To demonstrate the simplicity of the model dynamics, we simulate bursts of only one flash in length. The only fitted parameter for the model is the coupling strength , which demonstrates a transition in the dynamics of the model where > 1.0 (Fig. 3).
If the number of fireflies increases indefinitely, or if there are visual obstacles in the environment, the assumption that each firefly can practically immediately perceive when another firefly starts flashing will no longer hold. In this case, a finite time delay in perceiving the onset of the flashing could lead to an inter-burst interval that is greater than what is expected for the ideal case. The resulting inter-burst interval distribution will consequently be shifted to the right compared to the distribution given by Eq. (2).
While the general ideas underlying the theory framework will continue to hold, the mathematical formulation will need more sophistication to take these subtler e↵ects into account. Existing mathematical models designed for emergent synchronization of individual oscillators could be extended to account for individual variability in the period of individuals. For example, the emergent periodicity predicted for the Kuramoto model converge on the mean frequency of isolated agents (21), and the dynamical quorum sensing converge on the low frequencies of the isolated agents (22). However, the fireflies converge on the highest frequency in the repertoire of isolated individual fireflies. While individual behavior may appear as extremely complex, collective behavior based on simple and credible behavioral rules converges towards a simple emergent phenomenon as we have demonstrated. This wait-and-start phenomenon might be observable in di↵erent biological systems as well.
The mathematical implementation of the proposed paradigm results in an interburst interval distribution which converges towards a unique possible value corresponding to the lower bound of the individual IB distribution, at increasing N . That means that in the limit of an infinitely large and entirely connected swarm, the smallest IBI always occurs. This is at odds with two empirical observations: 1) while most of the smallest IBI from an isolated firefly peak at 12s and more, there are some residual value between 5s and 12s; 2) natural swarms comprising thousands of fireflies do not exhibit a 5s period. We propose some explanation to reconcile these two facts.
First, fireflies are known to produce annex flash patterns, for instance for alarm, in addition to the primary courtship phrase. It is possible that isolated fireflies in a confining volume switch to di↵erent behavioral modes that produce atypical flash train with intervals less than what they would typically do in an unobstructed environment with responding peers. Secondly, it is possible that the swarm bu↵ers against unusual perturbations. More than finite-size e↵ects, the main caveat here is that the swarm is not all-to-all connected, as we showed previously (9). Therefore, the e↵ective N would be much smaller than the swarm size, and correspond more appropriately to the number of active fireflies in each other's field-of-view. As a consequence, the mode of the distribution would remain centered around 12s, rather than smaller value which would only occur at very large e↵ective N .
It is easy to imagine extensions of this work that leverage the spatial positions of individuals in the model using distance-or sight-dependent coupling to modify the adjacency matrix and add further complexity to the system, and this framework makes implementation of this idea ripe for a future endeavor. To provide direct evidence for the underlying mechanistic principles, further experiments are needed. A promising avenue consists of artificially and controllably tuning the interactions within the group, for example, artificial flash entrainment with an LED should be able to decrease the inter burst interval.

Experimental Data
The individual and collective flashing of P. carolinus fireflies was recorded during 10 nights of field experiments in June 2020 in Great Smoky Mountains National Park (Tennessee, USA). The experimental protocol had been developed and implemented the previous year (8). In the natural swarm with hundreds to thousands of interacting fireflies, collective flashing consists of synchronous flashes every T f ' 0.5s, during periodic bursts T b ' 12s (Fig. 1C). However, is has been observed previously that individual fireflies in visual isolation do not exhibit burst periodicity. To characterize the onset of burst flashing, we performed experiments in a controlled environment. Fireflies were gently collected using insect nets, then placed individually in small plastic boxes, where species and sex were verified. Males were subsequently introduced into a secluded cuboid tent (approximately 1.5 ⇥ 2 ⇥ 1.5m 3 ) made of breathable black fabric and covered by a black plastic tarp to insure optimal visual isolation from fireflies on the outside. A GoPro Fusion 360-degree camera placed inside the tent recorded the entire volume at 30 or 60 frames-per-second (fps). Flashes were detected in video processing by intensity thresholding. Burst were identified as (temporal) connected components of flashes less than 2s apart. Interburst intervals ⌧ b were calculated as the duration between the first times of successive bursts. Tent experiments allow to observe the collective behavior of a small and known number of fireflies in interaction, while providing enough space for them to fly, hence reducing experimental artifacts from excessive confinement. We observed the flashing behavior of both individual fireflies in isolation and groups of 5, 10, 15 and 20 fireflies. We observed 7 individual fireflies alone in the tent, over durations between 30min and 90min. We observed that although these fireflies produced flash trains at a frequency of about 2Hz, the delay between successive trains was apparently randomly distributed, from a few seconds to tens of minutes. Then, we carried out 3 sets of experiments where the number of fireflies was increased to 5, then 10, then 15, then 20, each condition being maintained between 15min and 30min. As previously reported, collective burst flashing only appears at about 15 fireflies.

Derivation of mathematical formulation
The probability distribution P N (T ) of the inter-burst interval T of a group of N fireflies can be calculated as the probability distribution that one of the N flies starts flashing intrinsically at time T , while the rest (N 1) fllies haven't flashed until then. This is given by Eq (2), which is replicated below for convenience: The three successive terms on the right are, from left to right: the number of ways to choose the leader firefly: N ; the probability that the remaining N 1 fireflies have not flashed till time T , i.e., that they will flash later: ; and the probability density that leader firefly flashes at time T : b(T ).

Numerical demonstration
We use numerical calculations to demonstrate how synchronised periodicity arises in an arbitrary system which follows the extreme-value statistics used in our theory. Here, we take an arbitrary probability distribution (given by N = 1 label in Fig. 4) and plot the distribution of the minimum of N samples obtained from the N = 1 distribution. The distributions for arbitrary N are described by Eq. 4 as derived previously. As N increases, these distributions become sharply peaked with maximum probability peaked at a value larger than the minimum of the N = 1 distribution. For a system in which these quantities represent the interval between events, for large N , those events would become highly periodic as the width of the distribution narrows.

Approach to Delta Function
We show that as the number of fireflies (N ) increases, the variance and all the moments of the interburst interval distribution decrease and the distribution eventually converges to a Dirac Delta function. From Eq. 4, the m th moment for N fireflies is, Let the function be defined as, thus, We expect the distribution of inter-burst intervals to terminate at some large value and not go on to infinity (at most, they are limited by the finite lifespan of the fireflies), thus, Now, at any given value of t, N (t)  N 1 (t). This inequality is strict whenever 0 < (t) < 1. Such a region exists unless b(t) is a Dirac Delta function. If b(t) is a Dirac Delta function, then P N ( ) hT m Thus, all moments strictly decrease as N increases. From Eq. 8, the variance for N fireflies is, Writing the second term initially as a multiple integral over the entire t, t 0 > 0 plane, In the preceding step we have used the symmetry of the integrand under t $ t 0 . The second term of Eq. 11 can be similarly written down: Combining, Thus, The two functions in the above integrand satisfy: 0  (t)  (t 0 )  1, using the properties of the cumulant function. Thus, Rearranged, this tells us that the integrand in Eq. 15 is non-positive (i.e.,  0) everywhere. Thus, we have proved that V N +1  V N . In other words, the variance of the flashing distribution monotonically decreases with increasing number of fireflies. Further, as N ! 1, N (t) ! 0 for all t above T 0 (which is the maximum value of t below which b(t) is 0). For values of t below T 0 , N (t) = 1 irrespective of N . Thus, from Eq. 8, which represents moments of the Dirac Delta function P N !1 (T ) = (T T 0 ). Thus, as the number of fireflies tends to infinity, the distribution of interburst intervals tends to a Dirac Delta function peaked at T 0 .

Behavior of Mode
For a single firefly interburst interval distribution b(t) that is continuous for t T 0 and di↵erentiable for t > T 0 (where T 0 is the maximum value of t below which b(t) is 0), we show that the left-most mode shifts to the left as the number of fireflies (N ) increases, unless it reaches T 0 , in which case it stays at T 0 on increasing N .
The mode would be the local maximum of distribution P N . Di↵erentiating Eq. 4, Let the left-most mode of P N be located at t = t ⇤ N . If t ⇤ N = T 0 , we have Now, on increasing the number of fireflies by 1, we still have Thus, the mode stays at T 0 . On the other hand, if t ⇤ N > T 0 , we have, Now, on increasing the number of fireflies by 1, we get Thus, P N +1 increases towards the left of t ⇤ N , i.e., T 0  t ⇤ N +1 < t ⇤ N . Thus, the left-most mode shifts to the left as the number of fireflies (N ) increases, unless it reaches T 0 , in which case it stays at T 0 .

Preparing input for the model
The input distribution for the model's inter-burst interval T b is sampled directly from envelope distributions that encapsulate observations of one firefly's inter-burst interval. These envelope distributions were generated using an interpolating -spline between bin centers of the histogram of the distribution, normalized so that the area underneath the envelope sums to 1. The protocol for generating this envelope distribution is as follows: 1. Read and clean the data (d) instantiate each of N agents with one value from tbs2 and run simulation Figure 5: Slope of the relationship between T b2 and T b1 as a function of beta. The stable beta is found for each N by choosing the last value for which the absolute value of the slope exceeds 1.0.

Peak inter-burst intervals
The peak values of the inter-burst intervals are also shown here for additional clarity. Note that as N increases, the peak T 2 b is slowly decreasing. As N continues to increase to infinity, we expect to see this value trend downwards until reaching T 0 b , as discuss in Fig 1.C. T 1 b is defined as the value of T b for which the probability density function is maximal at a given , where the value of T b < 0.5max(T b ). T 2 b is similarly defined as the value of T b for which the probability density function is maximal at a given , where the value of T b > 0.5 ⇤ max(T b ). This splits the distribution into two pieces, the upper and lower regimes. We detect the peaks using the scipy.signal.find peaks algorithm with a height = y ✏ = 0.01, a prominence or y-delta of 0.005, and an x-distance of 50. If the probability is less than y ✏ , we exclude the peak and replace its height with y ✏ .

Simulation parameters
All experiments carried out with this agent-based model were conducted via simulation. Code for the simulations and the data processing was written in Python and can be found at this link. The simulation outputs a time series of flashes and their positions. For each set of parameters, we ran simulations for one hundred trials of 30,000 timesteps each. Parameters can be varied run-by-run via command-line arguments, which made a grid search parameter sweep over coupling strength and number of fireflies N easily parallelizable. All other values required for the synchronization dynamics are instantiated from experimental observations as explained in the main text.