Localization Is Correlated with Scenting.

In the 2D experimental arena, worker bees placed at the opposite corner from the queen explore the space and eventually aggregate around the queen’s cage after $$ \sim 30$$ min (Fig. 1C). For each scenting bee $$ i$$ at time $$ t$$, we collected its position, $$ {\mathbf{s}}_{i,t}^p$$, and direction of scenting, $$ {\mathbf{s}}_{i,t}^d$$ (unit vector). We then correlated the scenting events with the spatiotemporal density of bees in the arena by treating $$ {\mathbf{s}}_{i,t}^p$$ and $$ {\mathbf{s}}_{i,t}^d$$ as a set of gradients that define a minimal surface of height $$ f\left(x,y,t\right)$$. This surface height $$ f\left(x,y,t\right)$$ corresponds to the probability that a randomly moving nonscenting bee will end up at position $$ \left(x,y\right)$$ by following the local scenting directions of scenting bees (see SI Appendix, section E for formal definitions and derivations).

To show that the attractive surface $$ f\left(x,y,t\right)$$ is correlated with the final concentration of bees $$ \rho \left(x,y,t^{\prime }\right)$$, we compute the normalized mutual information, $$ NMI\left({⟨f⟩}_t;⟨\rho \left(t_{end}\right)⟩\right)$$, between the attractive surface averaged over the entire experiment (Fig. 1D) and the density of the bees averaged over the last 5 min of the experiment. The mutual information measures the information that the two variables, $$ {⟨f⟩}_t$$ and $$ ⟨\rho \left(t_{end}\right)⟩$$, share, and determines how much knowing one variable reduces uncertainty about the other. The NMI is scaled between zero (no mutual information) and one (perfect correlation). We averaged the density of bees over the last 5 min of the experiment to capture the density of the entire group and avoid discrete peaks of density resulting from movement of individual bees (see SI Appendix, section E for definitions). For the experiment shown in Fig. 1 C–E, the NMI is 0.21. The top right region in the arena with maximal $$ {⟨f⟩}_t$$ corresponds to the location around the queen, illustrating the correlation between the scenting directions and the queen localization of the bees. We performed $$ N=14$$ experiments of various group sizes (rounded to the nearest 10), and report an average NMI of 0.11 ($$ {\sigma }^2=0.002$$). See SI Appendix, Table S1 for values for individual experiments, and see SI Appendix, Fig. S3 for pairs of the average frame and the average attractive surface for all experiments.

A Characteristic Distance between Scenting Bees.

The aggregation process is accompanied by a sharp increase in the average number of scenting bees, followed by a slow decay, until there are few scenting bees and the majority of the bees are aggregated around the queen (Fig. 1E, orange curve). See Movie S3 for the experiment shown in Fig. 1 C–E, with the attractive surface reconstruction and time series data. We also characterize the temporal dynamics of the distance between neighboring scenting bees, which are treated as adjacent points in the Voronoi diagram for each frame (see example diagram in SI Appendix, Fig. S2 E and F and section C for more details). Throughout the growth in the number of scenting bees, the distance between scenting bees decreases to a minimal value, then increases as the number of scenters decreases (Fig. 1E, green curve).

To show the reproducibility of the behavioral assay and assess the effect of group size in the aggregation process, we tested 14 group sizes ranging from approximately 180 to 1,000 bees. Three example experiments with $$ N_{\mathrm{b}\mathrm{e}\mathrm{e}\mathrm{s}}=320$$, 500, and 790 are shown in Fig. 2A and Movies S4–S6. Each example includes three snapshots showing the process of aggregation around the queen’s cage over 60 min. The average numbers of scenting bees over time for all experiments of various group sizes are shown in Fig. 2B. Generally, there are more scenting bees over time as the number of bees in the arena increases. Across densities, we observe a typical characteristic of a sharp increase in scenting bees at the beginning as the bees initially search for the queen, and a slow decrease as they slowly aggregate around the queen’s cage. Based on observations, we assume that bees in tight clusters, such as those at the very beginning, do not usually scent, as they are not sufficiently spread out to fan their wings.

Fig. 2.

Distance between scenting bees is invariant to bee density. (A) Example experiments of various bee densities, each with three snapshots showing the process of aggregation around the queen’s cage over 60 min. (B) The average number of scenting bees over time for all experiments of various densities. There are more scenting bees with higher number of bees in the arena. Across densities, we observe a typical characteristic of a sharp increase in scenting bees at the beginning and decrease over time. (C) The average distance to the queen over time for all experiments of various densities. This distance gradually decreases as bees aggregate around the queen. The positive correlation between distance to queen and number of scenting bees is evidence for the functional importance of scenting events to the problem of queen finding. (D) The average distance between scenting bees over time for all experiments of various densities. Across densities, these distances are fairly constant at the beginning around the peaks in B, and increase over time. (E) The average distance between scenting bees at the peak in number of scenting bees as a function of density. Error bars represent the standard deviation over all of the distances for each density. The distance between scenting bees is approximately constant across densities, except for a low-density outlier. A linear model (green line) is fit to the data, excluding the outlier, with a slope of $$ m=0.00031\hspace{0.17em}\mathrm{c}\mathrm{m}/{\mathrm{n}}_{\mathrm{b}\mathrm{e}\mathrm{e}\mathrm{s}}$$.

The time evolution of the average distance from all bees in the arena to the queen is shown in Fig. 2C. There are no clear distinctions between the temporal dynamics of this distance as a function of density. For most densities, there is a sharp drop in this distance to the queen very early on ($$ \sim $$100 s to 200 s), corresponding to the early increase in number of scenting bees. The average distance to the queen then gradually falls as bees make their way to her vicinity.

We also measured the distance between neighboring scenting bees over time for all group sizes (Fig. 2D). This distance shows a characteristic gradual increase over time. At the beginning of the aggregation process at the peaks in scenting bees, this distance is fairly constant across all sizes, ranging from 5.00 cm to 7.35 cm, with an outlier from the lowest group size experiment at 9.06 cm (Fig. 2E). A linear model is fitted to the data (excluding the outlier) with the slope $$ m=0.00031\hspace{0.17em}\mathrm{c}\mathrm{m}/{\mathrm{n}}_{\mathrm{b}\mathrm{e}\mathrm{e}\mathrm{s}}$$. The presence of this constant distance between nearest-neighbor scenting bees during the initial stage of the aggregation suggests the possibility of a pheromone concentration threshold that turns on scenting for individual bees in this collective communication network. As more bees aggregate around the queen, the bees collectively “turn off” scenting. The higher distance between scenting neighbors later in the recordings suggest that the few remaining scenting events are more stochastic.

These experimental results, scenting behavior with wing fanning to direct pheromones, the threshold-dependent triggering of this behavior, and a characteristic distance between scenting bees, serve as core ingredients for our agent-based model. We will use the model as a proof of concept of our proposal of the mechanistic localization behavior, described in more detail in the next section. The goal of our modeling is to test hypotheses of the mechanisms behind the phenomenon and explore the possible emergent patterns that arise to assess the effect of the directional signaling strategy employed by the bees.

Model Predicts Optimal Signal Propagation within a Range of Behavioral Parameters.

We systematically explore a range of values for the directional bias $$ w_b$$, and the threshold $$ T$$, for various numbers of bees in the arena $$ N$$. Parameters $$ w_b$$ and $$ T$$ are potential behavioral parameters that bees could adjust based on input from the environment. The directional bias $$ w_b$$ represents the magnitude of the directional advection–diffusion of pheromone released by a bee (see Fig. 3B for a comparison of isotropic [$$ w_b=0$$] vs. advection–diffusion [$$ w_b=10$$]). The threshold $$ T$$ is the local concentration of pheromone above which a bee is activated from the random walk (i.e., a bee scents or walks up the concentration gradient when this threshold is met). Each combination of the three parameters ($$ w_b$$, $$ T$$, and $$ N$$) is a simulation repeated 20 times. All other parameters, including $$ C_0$$ (the initial pheromone concentration), $$ D$$ (diffusion coefficient), and $$ \gamma $$ (decay constant) from SI Appendix, Eq. 6, remain constant across all simulations.

Fig. 3.

Scenting model. (A) Our $$ L\times L$$ 2D simulation box is discretized into $$ l\times l$$ sized pixels. (B) Scenting bees can produce directional bias (example on Right with $$ w_b=10$$), compared to no directional bias on Left. (C) When a bee detects a local pheromone concentration that is above the activation threshold ($$ C\left(x,y,t\right)>T$$), the bee calculates the concentration gradient around it (using the nearest nine pixels, highlighted in different shades of green) and walks up gradient toward higher pheromone concentrations. (D) Example of a simulation showing a series of snapshots. The queen is shown as a red filled circle, and worker bees are shown as a red filled circle colored by their internal state: scenting (green), activated (orange), and nonactivated (gray). Activated bees perform a directed walk up the pheromone gradient, and nonactivated bees perform a random walk. The instantaneous pheromone concentration $$ C\left(x,y,t\right)$$ corresponds to the green color scale. Simulation parameters are $$ N=100$$, $$ w_b=30$$, $$ T=0.01\hspace{0.17em}{C}_0=0.0575$$, $$ D=0.6$$, and $$ \gamma =108$$.

We quantify aggregation around the queen through the average distance of worker bees to the queen, $$ ⟨d\left(t\right)⟩=1/N{\sum }_{i=1}^N\sqrt{{\left(x_i\left(t\right)-x_q\right)}^2+{\left(y_i\left(t\right)-y_q\right)}^2}$$, where $$ x_i\left(t\right)$$ and $$ y_i\left(t\right)$$ are the x and y positions of worker bee $$ i$$, respectively, and $$ x_q$$ and $$ y_q$$ are the x and y positions of the stationary queen bee. We also extract various other properties of the aggregation processes in the model: the queen’s cluster size and the number of clusters as additional measures of efficiency, the distance to the queen from the farthest active bee (i.e., a bee that is scenting or performing the directed walk up the pheromone gradient) to assess how far signals as a function of $$ w_b$$ and $$ T$$ propagate, and the number of scenting bees as a measure of energy expenditure. Our model predicts four distinct phases for all densities, which are determined by $$ \left(w_b,T\right)$$. Note that “phase” does not refer to a period in a time sequence but to the dynamics and outcome observed from different combinations of the parameters.

Phase 1.

For low values of both $$ w_b$$ and $$ T$$, the bees aggregate into small clusters that are homogeneously spread throughout the simulation box (Fig. 4E). This is reflected by a sharp initial decrease and then gradual decrease of $$ ⟨d\left(t\right)⟩$$ throughout the simulation (Fig. 4A), a consistently high number of scenting bees (Fig. 4B), a consistently small cluster of bees in the queen’s vicinity (Fig. 4C), and a consistently high distance of the farthest active bee from the queen (Fig. 4D). A simulation is provided in Movie S7.

Fig. 4.

Directional bias is associated with optimal search and aggregation in the scenting model. (A) The average distance of the worker bees to the queen as a function of time steps, for examples of the four different phases. (B) The average number of scenting bees as a function of time steps, for examples of the four different phases. (C) The average queen’s cluster size as a function of time steps, for examples of the four different phases. (D) The average distance of the farthest active bee to the queen as a function of time steps, for examples of the four different phases. (E–H) Example of a simulation from the four different phases, showing a temporal series of snapshots. The queen is shown as a red filled circle, and worker bees are shown as a red filled circle colored by their internal state: scenting (green), performing a directed walk up the pheromone gradient (orange), and performing a random walk (gray). The instantaneous pheromone concentration C(x,y,t) corresponds to the green color scale. Simulation parameters are $$ N=300$$, $$ C_0=0.0575$$, $$ D=0.6$$, and $$ \gamma =108$$. For phase 1, $$ w_b=0$$, $$ T=0.01$$ (E); for phase 2, $$ w_b=10$$, $$ T=0.5$$ (F); for phase 3, $$ w_b=30$$, $$ T=0.05$$ (G); and, for phase 4, $$ w_b=0$$, $$ T=1.0$$ (H). Phase 3, which is associated with the fastest aggregations around the queen, is highlighted in red.

The Spatial–Temporal Shapes of Clusters in Experiments and Simulations.

To connect the results of the model and the experiments, we analyze the clusters that change in shape and size over time in the experiments (for densities ranging from $$ N=250$$ to 400) and simulations (for density $$ N=100$$ where the percolation networks are best observed as elongated clusters). See SI Appendix, section I and Fig. S9 for how we isolate and define the clusters for this analysis.

To quantitatively compare the experiments to the different phases in the model, we measure farthest active distance, cluster area, cluster circularity, and cluster angle to the queen and show the time series of these attributes in Fig. 5. We exclude phase 4 simulations from these analyses, as the virtual bees are never activated, and no clusters form. First, in Fig. 5A, we show the distance to the queen from the farthest active bee (in experiments: scenting bees; in model: scenting or bee walking up the gradient) as a measure that differentiates the phases based on how far the signals propagate. In phase 1, this distance plateaus at a relatively high value when bees remain in small clusters throughout the arena. In phase 2, this distance plateaus at a relatively low value as the bees must be closer to the queen to be activated, and signals do not propagate as far as in phase 3, where we see a high initial distance as the percolation network forms, and a gradual decrease as bees cluster around the queen. The temporal dynamics of the experimental curve (green) are similar to the phase 3 simulation. Note that the oscillations in the experimental series are likely due to the spontaneous scenting events that may come from bees far away and not yet in the queen’s cluster. For each cluster attribute, we also measure the Pearson correlation between the experiments and the simulations in each phase, presented in SI Appendix, Table S4. For the farthest active distance, the coefficients between the experiments and the simulations in phases 1, 2, and 3 are −0.075, −0.662, and 0.759, respectively. For all attributes, the number of samples used to compute the correlation is $$ N=\mathrm{1,975}$$ (simulation data are down-sampled to match the number of data points in the experimental data).

Fig. 5.

Cluster properties of the experiments and simulations. For all panels, the green curve shows the average in experiments, and labeled gray curves show the average in different phases of the simulations. (A) The distance of the farthest active bee to the queen as a measure of how far the signals propagate. (B) The cluster area over time. (C) The cluster circularity over time. (D) The cluster angle relative to the queen over time.

Second, we quantify the cluster area for both the experiments and the simulations (Fig. 5B). In phase 1, the cluster around the queen is very small, as most bees are stuck in local equilibria. In phase 2, the cluster gradually increases and plateaus at a relatively intermediate size. In phase 3 and in the experiments, the cluster area quickly increases as the percolation network activates, and gradually decreases as bees approach and swarm around the queen. Shown in SI Appendix, Table S4, the coefficients between the experiments and the simulations in phases 1, 2, and 3 are 0.527, −0.784, and 0.762, respectively.

Third, we quantify the circularity of the clusters over time, in Fig. 5C, defined as $$ \left(4*Area*\pi \right)/\left(Perimete{r}^2\right)$$. A perfect circle has a circularity value of one. In phase 1, the cluster around the queen forms quickly and stays very small and relatively round, with a high circularity. Phase 2 simulations tend to have irregular clusters that form via bees finding the queen via the random walk, resulting in a relatively low to intermediate circularity values. In phase 3, we observe the clusters starting at relatively low circularity as the percolation networks activate, and becoming more circular over time as the compact cluster around the queen forms. The experimental clusters gradually increase in circularity similar to phases 2 and 3. Shown in SI Appendix, Table S4, the coefficients between the experiments and the simulations in phases 1, 2, and 3 are 0.684, 0.716, and 0.790, respectively.

Finally, in Fig. 5D, we show the angular deviation of the cluster orientations from the orientation leading to the queen. This relative angle plateaus for phases 1 and 2, in which the process of clustering around the queen is less dynamic and the cluster’s final shape is determined earlier. In phase 3 and the experiments, the cluster gradually orients toward the queen over time. We compute the angular correlation between the relative angle time series using $$ ⟨\cos \left({\theta }_{experiment}-{\theta }_{simulation}\right)⟩$$. Shown in SI Appendix, Table S4, the coefficient between the experiments and the simulations in phases 1, 2, and 3 are 0.49, 0.94, and 0.96, respectively.

Overall, we compare the structural phases we observe in the simulations and the spatiotemporal structures formed by the bees in the experiments. By extracting the four properties of the clustering dynamics and comparing the experimental data to the simulation data from phases 1 through 3, we show that the clustering of the bees in our experiments tend to behave inconsistently with virtual bees in phases 1 and 4. Hence, the real-world system likely exists within phase 2 and/or phase 3 of the model.

The Effect of Bee Density on the Phase Boundaries.

We use our model to construct phase diagrams for three different densities, low ($$ L^2/50$$), medium ($$ L^2/100$$), and high ($$ L^2/300$$). Across densities, low $$ T$$ and low $$ w_b$$ result in phase 1 with multiple clusters (Fig. 6). In this phase, the range of $$ T$$ increases, while the range of $$ w_b$$ decreases as density increases. On the other hand, phase 4 is the result of high $$ T$$ and high $$ w_b$$ values. With more bees in the arena, other phases expand, while the range of phase 4 gets smaller.

Fig. 6.

The effect of bee density on phase boundaries. Phase diagrams constructed from scenting model dynamics using summary heat maps in SI Appendix, Fig. S5 as a function of $$ T$$ and $$ w_b$$ for all three densities, $$ N=50$$ (A), $$ N=100$$ (B), and $$ N=300$$ (C). The black rectangles indicate the example simulations we showed in detail in Fig. 4. The area of the optimal phase (phase 3) grows with $$ N$$. Phase 3 in $$ N=50$$ is $$ \mathrm{103,\; 051}$$ pixels, or 15.9% of the total diagram; phase 3 in $$ N=100$$ is $$ \mathrm{137,\; 282}$$ and 21.2%, respectively; and phase 3 in $$ N=300$$ is $$ \mathrm{321,\; 807}$$ and 49.7%, respectively.

Phase 2 is typically the result of intermediate $$ T$$ values and low to intermediate $$ w_b$$ values. These higher $$ T$$ values prevent the percolation seen in phase 3, as the pheromone signals are not as far-reaching when bees are less sensitive to sense them to become propagators. Phase 3 with the successful aggregation via percolation is achieved with intermediate to high $$ w_b$$ and low $$ T$$ values. As density increases, the range of $$ w_b$$ and $$ T$$ for phase 3 is greater, as there are more bees to create and sustain the communication network in the arena. Phase 3 in $$ N=50$$ is $$ \mathrm{103,051}$$ pixels, or 15.9% of the total diagram; phase 3 in $$ N=100$$ is $$ \mathrm{137,282}$$ and 21.2%, respectively; and phase 3 in $$ N=300$$ is $$ \mathrm{321,807}$$ and 49.7%, respectively. When there is a higher density of bees, the bees are able to aggregate around the queen while not needing to be as highly sensitive (i.e., lower $$ T$$ value) to the signals. Across all densities, phase 3 never occurs when there is no directional bias (i.e., $$ w_b=0$$). This illustrates the importance of the directed signaling strategy to create an effective communication network within the model. The existence of a network of transmitters and receivers of pheromones that percolates across the entire computational simulation arena is crucial to achieve fast and successful aggregation around the queen.