Will biomimetic robots be able to change a hivemind to guide honeybees’ ecosystem services?

We study whether or not a group of biomimetic waggle dancing robots is able to significantly influence the swarm-intelligent decision making of a honeybee colony, e.g. to avoid foraging at dangerous food patches using a mathematical model. Our model was successfully validated against data from two empirical experiments: one examined the selection of foraging targets and the other cross inhibition between foraging targets. We found that such biomimetic robots have a significant effect on a honeybee colony’s foraging decision. This effect correlates with the number of applied robots up to several dozens of robots and then saturates quickly with higher robot numbers. These robots can reallocate the bees’ pollination service in a directed way towards desired locations or boost it at specific locations, without having a significant negative effect on the colony’s nectar economy. Additionally, we found that such robots may be able to lower the influx of toxic substances from potentially harmful foraging sites by guiding the bees to alternative places. These effects also depend on the saturation level of the colony’s nectar stores. The more nectar is already stored in the colony, the easier the bees are guided by the robots to alternative foraging targets. Our study shows that biomimetic and socially immersive biomimetic robots are a relevant future research target in order to support (a) the bees by guiding them to safe (pesticide free) places, (b) the ecosystem via boosted and directed pollination services and (c) human society by supporting agricultural crop pollination, thus increasing our food security this way.


Introduction
Biomimetic robots, that can naturally interact with living organisms, are a rather novel type of technology still requiring exploring their full potential. Such mechatronic agents could be a key mitigator to tackle the upcoming challenges for our ecosystems and societies: today's ecosystems are crumbling under multiple anthropogenic stressors, such as the ever-increasing land use [1], various forms of overexploitation and overharvesting of habitats as a consequence of increasing human populations [2,3], human-induced habitat loss and fragmentation [4], pollution [5], invasive species [6] and climate change [7]-to mention just a few external stressors that affect ecosystems negatively in a synergistic manner [8,9]. In their ultimate combined effect, anthropogenic stressors reduce the diversity and population sizes across the full spectrum of life, from microorganisms [10] to fungi [11] and wild plants [12], and across the wide range of animal life, as the cascading effects have been first reported in insects [13], followed by birds [14] and other vertebrates [15]. Ecological change can cascade through the layers of an ecosystem, as it was found in mathematical model analyses [16] as well as in meta-analyses of empirical datasets [17]. In consequence, a loss of diversity can cause extinctions, as high diversity and complex interspecific interaction networks are a key criterion for ecosystem stability [18]. Ultimately, a global tendency of diversity loss is observed, a phenomenon classified in several studies as being the beginning of the upcoming sixth mass annihilation or mass extinction in earth's history [19][20][21][22], marking the onset of a period of rapid loss of diversity. Such periods usually ramp up over several tens of thousands of years [21], suggesting the worst to come in the far future, while others predict severe civilization collapse effects already in a few decades [23].
Biomimetic technology can be a mitigation measure in this diversity crisis. Several studies suggest technology-based active ecosystem restructuring ('ecosystem hacking') as a mean to stabilize fragile ecosystems, e.g. by using biomimetic or biohybrid robotic technology [24,25]. A recent study suggests that biomimetic robots could 'rewire' ecosystem networks via ecologically immersive robots: in a 'proof of principle' experiment it was demonstrated that robots can mediate between a group of honeybees and a group of fish to enable a novel case of collective decision making on an interspecies level [26]. Ecosystem network restructuring may be a radical, but promising approach towards ecosystem stabilization, as several studies indicate the relevance of the network structure for ecosystem robustness and resilience [27][28][29]. Other studies suggest biomimetic and biohybrid technology for active ecosystem monitoring for conservation [30] or boosting specific important keystone species with technological support in order to strengthen ecosystems [31].
Here we focus on a potential application of biomimetic robots in order to support honeybee colonies' wellbeing in a biohybrid system consisting of living bees and mechatronic devices [31]. For providing foraging information artificially to honeybee foragers, several variants of biomimetic waggle dancing robots have been proposed [32][33][34][35], which imitate the dances that successful foragers perform to recruit other bees to their foraging targets [36,37]. For example, such robotic dancers could be used to keep bees from potentially dangerous foraging sites (contaminated with pesticides, fungicides, herbicides …) by guiding the bees to alternative food sources. Thus, such biomimetic technology can up-or downmodulate specific interactions within the ecological plant-bee network. In their extreme forms, such link modulations can lead to a rewiring of the interaction network by creating novel links, by advertising a yet unknown foraging target or by up-modulating one foraging target maximally to the cost of other foraging targets that get abandoned in consequence.
Network analysis of ecosystems revealed that interaction networks involving pollination services provided by animals are extremely fragile [38], thus the ongoing loss of pollinating insects is especially alarming. Among the crucial pollinators, honeybees are economically the most important, but endangered group [39]. 84% of European crops are (at least to some extent) dependent on animal pollination, 32% of the global crops species that are for direct human use are pollinated by wild bees and honeybees, another 11% can only profit from honeybee pollination, while a similar fraction of crop species profits only from wild bees [40]. Thus, this study puts the importance of the European honeybee (Apis mellifera L.) on the same level as the combined ecosystem pollination by the thousands of other wild pollinators together, making it the most important pollinator species for human food supply. However, given that the European honeybee forms massive colonies providing many thousand foragers per colony, the pollination service by honeybees cannot be overstated. Studies in the US have shown that five out of seven important crop species are limited by lack of pollination service in many locations across the country. Based on their findings, the authors predict the US-wide economic value of wild pollinators to be $1.5 billion for these seven crops, while the value for managed honeybees is estimated at $6.4 billion, underlining the economic importance of honeybees [41]. Across the most important 46 insect-pollinated crops a worldwide economic value of €625 billion was estimated [42], further stressing the economic value of pollination services. In addition, managed honeybees provide economic value beyond the pollination service, as they produce honey, wax and other goods. The global total honey market is estimated to be 2.4 million tons of honey worldwide [43], with a US market value for honey in 2019 and 2020 of around $300 million per year [44] and a global market of $7.84 billion per year in 2020 [45].
There are several reasons why honeybee colonies are an interesting location for biomimetic waggle dancing robots: first, honeybees as a keystone species establish their colonies on such a high complexity level that they are considered to be 'superorganisms' . Therefore, one or a few robots can centrally act in a hive and influence thousands of pollinators, in this way maximizing the ecosystem-service effect of the technologically augmented colony. Second, honeybees are the species providing the highest pollination value per species and also providing the highest direct economic value of all insects. Third, the fact that the technology never leaves the hive prevents potential environmental pollution. Fourth, eusocial insects like honeybees have evolved a rich set of social interaction mechanisms (like the waggle dance), a prerequisite for robots to engage in a meaningful animal-robot interaction. Fifth, honeybees are currently massively endangered, as the syndrome 'Colony Collapse Disorder' kills a high percentage of colonies every year [46]. Applying the robots to guide the bees to safe foraging targets could support their well-being.
In our study we use a mathematical model to predict the effect of such a biohybrid system consisting of bees and robots. We predict these effects with varying colony sizes, varying robot numbers and in varying environmental and colony conditions. Our main research questions are: In order to answer our research questions, we constructed a mathematical model that is capable of capturing the empirically observed dynamics of a full honeybee colony consisting of up to tens of thousands of adult bees. Our model is an ordinary differential equation (ODE) mean-field model describing the rates of change of 13 key system variables, that either represent quantities of bees (group sizes) that are in a specific behavioral state (e.g. foragers on a trip, foragers loaded with nectar in hive) or materials (nectar, honey) that are in a specific state of processing or storage, whereby our model focuses on the foraging of nectar and the choosing of nectar foraging targets. Table 1 shows these key system variables. The full model is described in the supplementary information file (SI), here we describe the most important aspects of our model.
The key structure of our model, as depicted in figure 1, is consisting of two submodels, one modeling the colony's foraging dynamics and one modeling the colony's nectar storage dynamics. An important regulatory feedback loop modulates the foraging according to the available storage capacity of the colony, while the collective foraging decisions of the bees are in parallel modulated by environmental conditions and by the biomimetic robots that are in the focus of our study. Three metrics predict the expected ultimate effects of the robots onto the vegetation, and onto the colony.
Our model follows the principles of system dynamics modeling [47] in which a stock and flow diagram (figure 2) graphically represents systems of connected ODEs in a way that the required conservation of mass within the system is always guaranteed. We model the foraging target selection system as a closed system (no bees are lost or hatch) and the nectar processing and storage system as an open system, as fresh nectar enters the colony through foraging while the foragers consume a fraction of it. The 'foraging dynamics' submodel handles the dynamics of task partitioning within each of the two foraging cycles that emerge when the modeled colony has to choose between two nectar sources (foraging targets A and B) that differ in their energetic quality concerning the distance of the foraging targets (energetic costs for flights) and the sugar concentration of the nectar offered there (energetic gains). Figure 2(i) shows the basic structure of this submodel, it follows the principle structure of an individual-based model's behavior state machine developed in a previous study [48,49].
The 'storage dynamics' submodel handles the storage dynamics inside of the hive, as nectar foragers do not store the collected nectar by themselves, but hand it over to another group of bees, the storer bees, via trophallactic (mouth-to-mouth) transfers. These storer bees then walk to the storage areas of the hive and search for appropriate storage cells. Figure 2(ii) shows the basic structure of this submodel, which basically follows the inspiration of a 'common stomach' model that was originally built to predict the regulation of nectar and pollen foraging of honeybee colonies [50].
For specific experimental analyses in this study, we have also built smaller additional submodels that are aimed at predicting (1) the energy accumulated by the colony, as the foraging provides energetic gains but also holds costs, (2) the pollination service provided to plants, and (3) the dynamics of toxic substances that might accumulate due to foraging on a specific foraging target that is treated with toxic substances. See sections A-I of the supplementary information file (SI) for detailed descriptions of those submodels.

Two submodels capture two important processes
These submodels interact with each other. For example, the largest two submodels depict the two important work-cycles (foraging and storing dynamics) that are involved in nectar foraging [51]. It is crucial for the efficiency of the colony that the workloads in these two work-cycles and the corresponding substance flows are in balance, as otherwise foragers would waste time by waiting for storer bees to unload, or vice versa. Thus, honeybees employ a set of behavioral feedback loops that consider queuing times and searching times in order to regulate the waggle dancing activity. These key regulatory mechanisms are captured in our model in the form of self-enhancing positive and regulatory negative feedback loops.

Forager recruitment through waggle dancing
Bees that return from their foraging trip can perform waggle dances that can recruit inactive bees to the previously visited foraging target. The longer they dance, and the more inactive bees are in the hive, the more recruitment will happen. In our model, this is implemented following mass-action-law mechanics [52], as for this recruitment process to happen a dancer must interact locally with an inactive bee that follows the dance and picks up the location of the foraging target that is encoded into the dance by the dancing bee [36]. This process creates a positive feedback loop in the system that drives the collective decision making.

Choosing among foraging targets decentrally
The number of dance rounds, and thus the duration of dancing of a returned foraging bee is determined by the energetic gain G and the energetic costs C of the foraging trip, both measured in Joules. It was observed in honeybees that the relative energetic net gain (G − C) /C is the best predictor for waggle dance lengths [37]. This way, energetically better foraging targets are advertised longer, thus the positive feedback loop that recruits novel foragers to that target is enhanced over others. Ultimately this allows the colony to always recruit the vast majority of foragers to the energetically optimal foraging target in the environment [53]. It was also found that waggle dance durations correlate with the overall nutritional status of the colony [37]. Figure 3(i) shows these dance-response curves, which we implemented in our model, and how they change with the nutritional colony status. Our model implementation follows a regression model of empirical data [37].

Abandoning foraging targets
Foragers abandon their foraging targets and become inactive again. These bees can then be either re-recruited to their previous foraging target or be recruited to another one. The rate of abandonment was found to depend again on the energetic quality of the foraging target on which the bees are currently foraging: the better this food source is, the less likely a bee is to abandon it. Thus, better foraging targets are abandoned slower than worse ones [53]. This way the colony keeps its ability to revert previously made foraging decisions. This process establishes a negative feedback loop in the system that allows the colony to react to dynamic patterns of floral nectar offerings. Locations of optimal nectar flows can change multiple times during a foraging day and nectar sources are locally limited, thus abandonment can prevent the colony from overexploitation and being trapped in local optima. Figure 3(ii) shows how we have modeled the abandonment of foraging targets depending on the bees' foraging flights' relative net energetic gains, following the implementation from a previous agent-based model [48] of the honeybee foraging system.

The unloading process
The returning foragers provide a nectar inflow to the hive that needs to be taken over by storer bees. If the nectar influx rate is higher than the unloading rate, a waiting queue emerges on the forager side, yielding delays in the foraging cycle. We model this process again following a mass-action law mechanism, as the loaded forager and the unloaded storer bee need to locally interact. We assume here a constant number of available storer bees and consider the crop load volumes of returning foragers, the number of foragers to unload and the capacity of the crops of storer bees, following empirical studies [51,54]. This is part of a negative feedback loop that regulates the foraging according to the colony's unloading capacity.

The storing process
The storer bees need to find suitable storage cells in the hive's comb before they can return to the interaction area with empty crops in order to unload the next loaded forager bee. The length of this storage cycle depends on the nutritional status of the colony, the fuller the combs already are, the longer it takes the storer bee to perform the task [51]. This is part of a negative feedback loop that regulates the foraging also according to the colony's nutritional status and can significantly modulate the colony's foraging strategy . Important functions of key microscopic mechanisms, which model the key behavioral decisions of honeybee foraging: (i) waggle dance rounds to perform, if a dance happens (see SI equations (10a) and (10b)), (ii) abandoning their foraging target (see SI equations (12a) and (10b) or (iii) deciding whether or not to dance for their foraging target at all (see SI equation (7)). Functions and their constant parameters are taken from various empirical studies and models [37,48,53,55,56]. A full list of parameters is given in table SG1 in the supplementary information file (SI). by adjusting its exploitation-to-exploration ratio. It may lead to having saturated colonies choose differently amongst foraging targets than hungry colonies, an ecophysiological aspect that we will also analyze in this study.

Decision whether to waggle dance or not
The longer a forager bee has to wait for unloading of its foraging load the less likely it will perform a waggle dance that will recruit additional foragers, as the colony is already out of balance in its unloading process [55]. This closes the large important negative feedback loop in the system, in which also the processes described in (4) and (5) are important components. This feedback loop also regulates the two work cycles of foraging and storing in a decentralized way. Figure 3(iii) shows how we modeled the dependency of mean waggle dance performance rate on the mean queuing delay of foragers. Figure 3(i) shows the variants of stimulus-response curves for dance length, which depend on the colony's nutritional status. The variable CFS (t) ('colony fill status') resides between 0.0 (empty storage combs) and 1.0 (full storage combs). It expresses the fraction of filled nectar or honey cells in combs within the space of all storage combs. The variable τ Queue is the average queuing delay a forager experiences upon return to the hive while waiting for unloading. The variables Q A (t) and Q B (t) are the average quality each forager experiences on their visited foraging targets. The given sugar concentrations of the foraging targets are based on sugar concentrations given in an empirical study [37,56] also considering the target distances to the hive and the foraging loads reported with artificial sugar water feeders [37]. The linear functions that describe the bees' dance durations are shifted to different x-axis intersections (dancing thresholds) while keeping their slope (figure 3(i)). Starved colonies show a very low dance-threshold at a sugar concentration of approx. 0.6 mol l −1 (corresponding to a quality of 12.86 J J −1 ), while well balanced colonies show the threshold at a sugar concentration of approx. 1.2 mol l −1 (quality of 24.23 J J −1 ), and massively saturated colonies show a dancing threshold at a sugar concentration of approx. 1.6 mol l −1 (quality of 30.15 J J −1 ) [37]. Our model's dance threshold functions (figure 3(i)) correspond closely for those reported in literature for starving, normal and filled colonies [37] with settings of CFS (t) ∈ {0.1, 0.6, 0.8} . The dancing threshold (intersection with the horizontal axis) is maximally 1.0 J J −1 away from the empirically observed data while the slopes of the curves fit exactly to the regression model derived from empirical data [37].

Adapting to various hive conditions
As a key component for our study, we model the effect of biomimetic waggle dancing robots on the foraging system. These robots are considered in the foraging submodel as additional waggle dancers for one of the two foraging targets. These robotic waggle dancers need not be considered in the storage submodel, as they do not carry scented nectar loads that need to be stored, although some robots might provide nectar samples to dance following bees, as an extra incentive. However, these nectar amounts are so small that they can be neglected for the storage dynamics.
The two submodels of our model are connected, so that a large negative feedback loop that regulates the overall system emerges. The storage submodel predicts for every point in time how many foragers can be currently unloaded. This is determined by the available storer bees that arrive empty at the dance floor area of the hive to unload the next forager. These empty available storers pose the overall unloading capacity (in µl), which ultimately determines how many foragers can be unloaded, as they bring a certain quantity (also in µl) of nectar to be stored. From the ratio of unloading capacity to required unloading, the is calculated. This unloading-driven transition represents the causal connection from the storage dynamic submodel to the foraging dynamics submodel. The two models are also connected causally in the other direction: the foraging dynamics submodel predicts for every point in time the number of foragers that come back to the hive and the nectar loads (in µl) that they bring back. This is the input flow that adds to the stock of still unloaded nectar N Foragers (t) in the storage dynamics submodel.
Overall, this creates a large negative feedback loop in the system, as it is depicted in figure 2(ii). This large negative loop counteracts the positive feedback loop of waggle dance recruitment. It ultimately determines the ratio between exploration and exploitation in the emergent foraging strategy. This ratio differs with the quality of offered foraging targets in the environment, with the ratio of forager population to receiver population in the hive, as well as with the nutritional state of the colony, which is its already accumulated nectar or honey stores. As our main objective is to investigate with the model how and when biomimetic waggle dancing robots may influence this collective decisionmaking process in honeybees, this large negative feedback loop is depicting the relevant key components of the honeybee colony self-regulation.

Numerical methods
In order to conduct our experimental simulation runs, our ODE model is solved with the RK4A method (4th-order Runge-Kutta method with automatically adjusted step size) and a time step of 1.0 s. We decided to use this numerical solver as we are modeling a continuous process and RK4A is used to solve differential equations. The short time step is necessary as we assume one round of a waggle dance to last for 2 s, as the duration of the waggle phase was found to scale linearly with the distance to the foraging target between 0.2 s and 0.65 s [57,58] for the most common foraging distances and can increase above 2.0 s for very far-away foraging targets. In this variant the chosen step-size is the maximum step size that is used for solving the equations, however, the algorithm automatically decreases the step size whenever needed based on the observed error. In addition to the waggle phase the bees then need to complete the dance round (make a sort of u-turn) before starting the next dance round. Thus, we assumed 2.0 s as an approximation of the duration of one waggle dance round for our mean field model.

Simulated environments and their dynamics
In total we simulated three distinct types of environmental dynamics (foraging target quality changes) in combination with three types of hive conditions to answer our research questions. The first type of simulation experiment confronted the simulated colony with a choice between two foraging targets of different energetic quality during the morning. At noon, the qualities of the foraging targets were suddenly altered and stayed then constant for the rest of the afternoon.
In our second type of simulation experiments, two foraging targets of equal energetic quality were available for the simulated colony during the morning. At noon, the sugar concentration at one foraging target suddenly increased, while the other target stayed unchanged. Both above-mentioned experimental settings were used to answer Q1 and Q2.
A third type of simulation experiment was employed to answer Q3 to Q8: it reflected a more natural and plausible foraging condition for agriculturally used honeybee colonies. Two available foraging targets that differ in their quality consistently throughout the whole runtime of the experiment were simulated. In these experiments, one foraging target always offered a better energetic quality than the other. At halftime of the experiment (noon), we assumed that a potentially harmful (toxic) substance was applied to the crops at the qualitatively better, but now unsafe, foraging target. Simultaneously, biomimetic waggle dancing robots were activated in our simulation. These robots advertised the safe, but less profitable, alternative foraging target, which the colony had not chosen for mass-exploitation previously. This experiment was performed with various energetic quality differences between the available foraging targets.
The metric used to answer Q1 and Q2 considers the predicted sizes of forager populations that were recruited to each foraging target in our simulation. Therefore, we observed the sum of all bees that resided in the blue stocks in our stock-andflow model shown in figure 2(i) for quantifying the currently active foragers on foraging target A, and the sum of all bees that resided in yellow stocks in figure 2(i) for quantifying the currently active foragers on foraging target B.
Additionally, the amount of total harvested energy gained from both foraging targets during the day was measured by first accumulating the returned crop loads of all foragers returning from either foraging target. Then, these amounts of nectar were converted into amounts of energy (in Joule) considering the sugar concentrations of those nectar loads by using the equation given in [37]. We plotted our data in 15 min wide intervals, analogously to validation data from empirical studies [53,59].
The metric used to answer Q3, Q4 and Q5 is given by the fraction of the population of all active foragers that were visiting either target in our simulation runs. The metric used to answer Q6 is given by the total number of pollination flights that were performed to foraging target A in the afternoon, which is the robotadvertised target. In order to answer Q7, the total number of tainted (toxic) loads that were returned by foragers from target B during the afternoon period in our simulation runs was considered, assuming that each returning forager transported one tainted load back to the colony. In order to answer Q8, the accumulated energy (Mio. Joule) from both targets during the day was observed, as described above.

Model validation
Before we used our model to investigate the effect of biomimetic robots, we first validated it by comparing specific model predictions with corresponding experimental results from literature. Our research question Q1 is aimed at this model validation step. We chose two seminal experiments for this validation step: the first one investigated collective choice (and re-choices) of honeybee colonies in response to environmental fluctuations [37]. The second one revealed a cross-inhibition effect between potential foraging targets [59].

Validating with collective foraging decisions
In order to confirm that our model is capable of producing plausible predictions of collective behaviors of honeybee colonies, we used an empirical experiment as a benchmark [53]. In this experiment two artificial ad-libitum feeders, placed equidistant (400 m) to the colony in opposite directions, offered sugar water to a colony's foragers and the number of foragers on both feeding sides was recorded. The colony starved before the experiment for many days due to weather conditions (cold, rain). In the experiment the first 4 h (morning phase) offered a high sugar concentration solution, thus a high energetic reward, on one feeder (2.5 mol l −1 ) and a low sugar concentration of the other feeder (1.0 mol l −1 ). At noon, the environmental situation was flipped, and the former suboptimal feeder became the new global optimum by now offering a high sugar concentration (2.5 mol l −1 ) while the former global optimum worsened to 0.75 mol l −1 . These experiments showed that the honeybee colony first mass-recruited to the global optimum in the morning and then, after the environmental change, reverted its foraging strategy to massively forage at the new global optimum. Figure 4(i) shows the empirically observed data (published in [53]) and figure 4(ii) shows the corresponding model prediction when we parameterized our colony to a similar setting: both foraging cycles were initialized with 12 bees (1 bee dancing, 0 bees on trip, 8 bees loaded in hive, 3 bees unloaded in hive), and the colony was set to have initially 110 inactive foragers. Our model's parameterization follows as closely as possible a previously published study [53]: this study contains an empirical experiment and an accompanying model. This model was specifically tailored to predict these very specific experimental conditions of that study and assumed a small forager population. Our model contains more stocks (that require initialization) than the previous study. However, we initialized our model in a way that the overall number of foragers in each cycle are identical. For measuring the flight activities on both foraging targets, we considered the sum of all bees in the four stocks that are associated with a specific foraging target A or B (blue or yellow stocks in figure 2(i)). Compared to the results of that model, we get very similar model predictions with our model, having the advantage that we will later be able to parameterize our model to many other foraging and colony-state conditions.
An experimental setup with artificial ad-libitum feeders is bringing side effects, as these sugar-water feeders are not reflecting very well the natural foraging conditions in which honeybees usually have to collect their nectar at blossoms. With artificial feeders the bees collect the fluids easily and, thus, have rather high crop loads when returning to the colony [53]. We represent this in our model parameterization by also assuming high forager loads by setting v Load = 48 µl/bee. High nectar influx into the colony will activate many storer bees. This is because a colony increases the number of storer bees [55] in those periods, as foragers that experience long queuing delays perform 'tremble dances' , which activate additional storer bees. This process is not part of our model. However, to reflect the outcome of this process, we set n Foragers = n Storers assuming enough storer bees in the colony. We used a low initial colony fill status CFS = (0) to model a 'starved' colony.
Our model can also predict the energetic dynamics of the colony, as it is shown in figure 4(iii). We find that foraging predictions are associated with a certain period of lowered energy influx. These 'energetic costs' are caused by environmental change. These predictions are consistent with predictions published based on (far more) complex agent-based models [48,49]. The fact that we also find very similar dynamics here is further validation of our model.

Validating with cross-inhibition effects
In a second empirical experiment, a feeder setting very similar to the setting used in section 3.1 was used. The main difference to the previous validation experiment was that in a 'morning phase' both equidistant ad-libitum feeders offered the same sugar concentration of 1.25 mol l −1 , and at noon, the sugar concentration at one feeder was increased to 2.5 mol l −1 for the afternoon, while keeping the other feeder unchanged.
In consequence, the bees foraged on both sources at a medium level in the morning. After the sugar concentration was increased at one feeder, the colony mass-recruited foragers to this feeder, while abandoning the other, which did not change in its quality. Clearly, a cross inhibition effect was observed [59], as shown in figure 5(i). We kept the model parameterization used before and only adjusted the sugar concentrations and the initial number of bees in each foraging cycle accordingly. Figure 5(i) shows the empirical data (published in [59]) and figure 5(ii) shows our model predictions with high agreement to this data. Again, we also investigated the energetic dynamics of the colony and, as expected, found that the improvement of the environmental offerings is also reflected in the energy accumulation in the colony, as is shown in figure 5(iii).
Overall, we find that our rather simple meanfield model is in good agreement with classical empirical studies on honeybee's collective foraging decision making. Based on these findings we can now explore our focal research questions via model predictions.

From artificial feeders to natural foraging
To answer our research question Q2, we developed our model further in order to better reflect the natural conditions that full-sized honeybee colonies encounter in the field, in contrast to experimental procedures conducted with small observation hives and artificial sugar water feeders. To address natural conditions, we made the following changes to our model parameterization and initialization: (1) Literature describes that foragers return to the hive with only partially filled crops when they have to fly to multiple nectar sources on their trips [54,60], in contrast to artificial ad-libitum feeders, where they tend to fill up the full crop volume quickly. We assumed here a crop volume of v Load = 21.8 µl/bee, as it is reported for periods of high natural nectar flows [54]. In times of low nectar flows in the environment, they return with even lower loads. (2) We increased the colony sizes to 30 000 bees, of which n Foragers = 7 024 are assumed to be foragers. We initialized our model with 7024 forager bees as a compromise of the data given from literature that varies vastly, attributing between 25% and 50% of the adult bees to being potential foragers [59,61,62]. Additionally, the number of frames in the hive was increased from 2 to 15 frames in order to reflect a full-sized honeybee colony. Finally, we considered 12 bees to be initially already foraging on each food source, to stay consistent with the previous foraging target choice experiments. (3) We further assumed that large colonies will have accumulated nectar reserves to supply them over the course of the seasonal weather changes, thus we assumed an initial colony fill status of CFS (0) = 0.6, as it is indicated as a medium colony fill status (slight nectar/honey gain per day) in figure 3(i).
Our simulation experiments show that changing these factors individually can have a large impact on the resulting prediction. Increasing only the colony size, see figure 6(i), keeps the overall dynamics same as for small scale observation hive colonies of the choice experiment (see figure 4). Reducing the nectar loads to the observed levels under natural foraging conditions seems to be an absolute game changer, as the colony does not redecide for a new foraging target in the afternoon period and keeps on foraging on the inferior foraging target for the rest of the day, see figure 6(ii). Adjusting only the colony fill status to less 'hungry' settings does not significantly change the predictions compared to the hungry colony in figure 4, it just reduces the number of predicted foragers, see figure 6(iii). However, applying all adaptation changes together yields a new view on the collective foraging decision making: the colony peaks in its foraging activity at approximately halftime of the morning phase, as saturation of storer bees kicks in at this point in time and this down-modulates the foraging cycle to fit to the storage cycle. The foraging approximately 'levels off ' at a certain forager number, which is significantly below the level of all available foragers. In the afternoon the colony reacts with a sharp decision making to the environmental fluctuation, see figure 6(iv).
We also validated our model with more natural parameterization and initialization against another simulation of the cross-inhibition experiment that is shown in figure 5 for a small colony in an observation hive foraging on artificial feeders. Again, we first adapted each coefficient and initialization individually, as we already did with the choice experiment in figure 6. Finally, we applied all changes of coefficients and initializations simultaneously. These simulation runs are shown in figure 7. In this experiment our model predicts comparable results to the choice experiment before: increasing only the colony size, see figure 7(i), keeps the predicted overall dynamics of the cross-inhibition experiment of small colonies in observation hives, see figure 5. Also, the model predictions strongly resemble the original experiment when reducing the nectar loads to the observed levels under natural foraging conditions (figure 7(ii)). Adjusting only the colony fill status to a 'less hungry' setting changes the whole foraging dynamics of the colony: the revision of the colony's initial foraging decision in the afternoon takes much longer than in the original experiment setting, comparing figures 5 and 7(iii). Applying all adaptation changes in parallel, our model predicts similar dynamics as in the original cross-inhibition experiment, just with higher foraging activities, see figure 7(iv). Compared to the original cross-inhibition experiment until approximately halftime of the morning phase the foraging activity rises at both targets. The foraging equilibrium for these environmental situations is reached by the colony. In the afternoon the colony reacts with a sharp decision making for the then better foraging target, at the expense of the forager population that frequents the other target, see figure 7(iv). Such a cross-inhibition is often a telltale-sign for the existence of a competition for a limited shared resource, which is here (1) the limited number of foragers that can be recruited via waggle dances, which then (2) compete for the limited number of storer bees to unload the nectar to, which then (3) compete for the limited suitable storage places for nectar in the combs.

Predicted impacts of biomimetic dance robots
In order to answer our research questions Q3 and Q4 we developed an agriculturally inspired scenario, in which a potentially useful application for biomimetic waggle dancing robots can be investigated: • We assumed that there are two crop fields on which the modeled honeybee colony can forage. • These two foraging targets can differ in quality.
They were set (for simplicity) at equidistant locations from the hive but might differ in the nectar concentration they offer. • Foraging target B always offered a better quality, thus it was the preferred foraging target by the colony. • At halftime of the experiment a potentially harmful (toxic) substance was applied to foraging target B. • At this point in time the robots started to advertise foraging target A as a safe alternative.
As a metric for the observed effect, we observe the dynamics of the forager populations that are recruited to each foraging target and express these population sizes as fractions of the total number of available foragers in the simulated colony. Figure 8 shows the emerging foraging strategies predicted for three distinct experimental settings, in which the foraging conditions differ significantly. Figure 8(i) shows a simulation in an environment in which two foraging targets of rather poor, but slightly different quality (1.25 mol l −1 vs. 1.6 mol l −1 ) are available to the colony. Figure 8(ii) shows a similar setting but here we allowed the simulated honeybee colony to choose between two rather rich foraging targets of slightly differing qualities (2.2 mol l −1 vs. 2.5 mol l −1 ). Finally, figure 8(iii) shows a scenario in which one foraging target was energetically very rich while the other one was rather poor (2.5 mol l −1 vs. 1.25 mol l −1 ), in analogy to the choice experiments shown in figures 4 and 6. In all 3 environmental settings the colony was left for the morning to reach its equilibrium foraging strategy for the specific environmental situation. At noon, we assume that the preferred foraging target has become unsafe to the bees and thus a population of robots (number varied) actively starts to advertise the foraging target that is not chosen by the colony, in order to influence the colony's decision making.
Our results show that in all three environmental scenarios, a maximum of 86% of all potential foragers are predicted by our model to forage at the energetically better foraging target (figure 8) in the morning phase when no robots are active. In addition, the final activity of the foragers that collect nectar at target A increases with the number of robots that are actively advertising this foraging target during the afternoon period in all three scenarios.
Our model predicts that in both scenarios where the quality of the foraging targets is only slightly different (compare figures 8(i) and (ii)), that the waggle dancing robots will be able to significantly alter the honeybee colony's decision making towards foraging at target A, so that after approx. 2-3 h a majority of forager bees will fly to this safe foraging target that is advertised by the robots. However, in the scenario where the choice has to be made between two targets of very different qualities, see figure 8(iii), even 100 robots are predicted to be unable to divert the colony's initial foraging decision to a majority level.
If two foraging targets with rather poor quality figure 8(i) are to be chosen from, our model predicts that it will take 25 waggle dancing robots to draw the majority (68%) of all active foragers, resembling 46% of the total forager population, to target A. However, this effect is predicted to saturate quickly with increasing numbers of robots, as 100 robots can direct only 74% of all active forager bees, resembling 50% of the total forager population, finally foraging at target A, which is a marginal increase of effect for a four-fold increase of the number of robots. In contrast to that, already a single robot is predicted to guide 23% of all active foragers towards target A, resembling 15% of the total forager population.
In the scenario with two rather rich foraging targets (figure 8(ii)) to choose from, a population of 100 waggle dancing robots is predicted to be required to alter the colony's foraging decision significantly by leading the majority (69%) of all active foragers, resembling 47% of the total forager population, to gather nectar at target A. Under such circumstances, the effect is not as strongly saturating with increasing robot numbers as in the previously tested scenario, see figure 8(i), thus our simulation results suggest that more robots might attract an even larger forager population to the advertised safe target A. Also, differently to the previous scenario, one single robot is predicted to direct 12% of all active foragers foraging target A, resembling only 9% of the total forager population.
Comparing these three scenarios, we find that the revision of the colony's initial foraging decision is predicted to emerge earlier with two rather poor foraging targets than under the conditions with two rather rich foraging targets. We classify the prediction as a revision of the previous foraging strategy as soon as more foragers fly to target A than to target B. Additionally, we find that a larger robot population speeds up and strengthens this induced decision-reverting process. With two foraging targets of very different quality, such a reversion of the decision making is not predicted to occur, however with larger robot populations it can be assumed to happen late in the afternoon.

Predicted impacts of the colony's food status
In order to answer our research question Q5 we conducted the following specifically designed simulation experiment. For this experiment, the same three environmental scenarios were simulated that we used in the experiments to answer Q4. However, instead of varying the number of robots we varied the colony fill status while keeping the number of robots constant at 10 robots. The initial colony fill status was varied from CFS = 0.1, which represents a very hungry colony, to CFS = 0.9, which represents a very saturated colony. On the CFS scale, an initial value of CFS = 0.5 represents medium filled nectar stores in the colony.
In the first scenario, where the simulated colonies had to choose between two sources with rather poor qualities, figure 9(i), our model predicts the emergence of three rather distinct foraging behavior patterns. These patterns depend on the simulated colony's fill status: very hungry colonies (CFS (0) = 0.1) will strongly choose with 93% of the total forager population to collect nectar at the more profitable foraging target B in the morning. This behavior will not change in the afternoon, although ten robots advertise for the other foraging target A then. Medium saturated colonies (CFS (0) = 0.5) are predicted to resemble a foraging strategy similar to the one observed in the choice experiment before (see figure 4), despite the absence of any environmental changes: during the morning period the majority of the forager population (85%) will forage at the energetically better target B. When ten robots advertise for the inferior foraging target A, up to 32% of the total forager population will be reallocated to the advertised target A in the afternoon. In contrast to that, rather saturated colonies (CFS (0) ⩾ 0.7) are predicted to show a foraging strategy similar to the one that was observed in the cross-inhibition experiment (see figure 5): during the morning period there will be no clear collective decision made towards either foraging target.
Of the total forager population, 32% will forage at target A and 40% will forage at target B. However, in the afternoon, when ten robots advertise the energetically inferior foraging target A, up to 50% of the total forager population will forage at the advertised inferior target A and the foraging on the better source B will decrease to 14% of the total forager population.
In the second scenario, with two rather rich foraging targets to choose from, our model predicts two distinct foraging behavior patterns to emerge, see figure 9(ii). Medium filled or hungry colonies (CFS ⩽ 0.5) are predicted to make medium-strong foraging decisions in the morning by choosing the optimal target B decisively (between 69% and 71% of the total forager population) but by also allocating many foragers to the slightly inferior target A (between 15% and 29% of the total forager population). This foraging decision is predicted to be kept throughout the afternoon, despite ten robots advertising foraging target A. In contrast to that, colonies with a higher colony fill status (CFS ⩾ 0.7) are predicted to exhibit a strong decision making towards target B in the morning period (68% of the total forager population). Our model predicts that ten robots will be able to revert these colonies' initial foraging decisions in the afternoon: colonies with (CFS = 0.7) will reallocate 43% of the total forager population, Figure 9. Effect of the colony's fill status (CFS (t)) on the emerging foraging strategy in three environmental conditions. (i) Two poor foraging targets of slightly differing qualities. (ii) Two rich food sources of slightly differing qualities. (iii) A choice between two targets of very different qualities. In the afternoon period (brown background) ten robots are activated, advertising the foraging target, which the bees have not chosen previously. The dashed vertical line indicates when the foraging targets change in quality.
to the advertised and inferior target A, while colonies with (CFS = 0.9) will reallocate 49% of the total forager population, to target A.
In the third scenario, where the simulated colonies had to choose between two targets of very different qualities, our model predicts that the emerging foraging strategy will follow a stricter correlation with the colony fill status than in the before mentioned scenarios, see figure 9(iii). Very hungry simulated colonies (CFS = 0.1) will allocate nearly 100% of the total forager population foraging at the optimal target B in the morning while the inferior foraging target A will be almost neglected. This will not change in the afternoon period, even if ten robots advertise the inferior target A. However, the more a colony is saturated with nectar, the smaller is the number of foragers that it will allocate to target B in the morning. In the afternoon, when ten robots advertise the inferior foraging target A, our model predicts that this will further lower the number of foragers on target B significantly. With a medium filled colony (CFS = 0.5) for example, our model predicts that ten robots advertising the inferior target A will lower the total forager population from 86% (in the morning) to 66% (in the afternoon). For very nectar saturated colonies (CFS = 0.9) our model predicts that there will be the strongest reallocation to the advertised inferior foraging target A in the afternoon, raising the number of foragers there from 0.02% to 23% of the total forager population.
Comparing across these three scenarios, we find that the revision of the colony's initial foraging decision through the effect of biomimetic robots is predicted to work best in environmental conditions where the qualities of the available foraging targets are only slightly different from each other and with colonies that have a rather high nectar fill status.

Effects on important performance metrics
In order to answer our research questions Q6, Q7 and Q8 the same three environmental scenarios were simulated that we used in the experiments to answer Q4. The predicted effect of different robot numbers on the colony's pollination service, on the pesticide influx from a foraging target and on the colony's energy gain was measured in the three different settings investigated with our model. Figure 10(i) shows the predicted impact of the robots on the total number of pollination flights towards A, which is advertised by the robots in the afternoon phase of the simulation experiment. Our results show that in all experimental settings the total number of pollination flights to target A increases with a higher number of robots.
With two targets of rather poor quality to choose from (blue bars in figure 10(i)), our model predicts 14 841 pollination flights to target A, without robots advertising this target. There are 1.9 times more pollination flights (from 14 841 to 27 821 flights) predicted towards target A with one robot advertising target A, than without robots advertising target A. If the number of active robots is quintupled from one to five robots, the number of predicted pollination flights is 1.9 times as high as with one robot. If the number of active robots is further quintupled from 5 to 25 robots, the predicted number of pollination flights is 0.9 times as high as with five robots. If the number of active robots is further quadrupled from 25 to 100 robots, the number of predicted pollination flights is 0.3 times as high as with 25 robots. Joule] at the end of the day. It is assumed that the toxic substance (e.g. pesticide) was applied to foraging target B at noon and that also the robots (if any) were activated at this point in time to advertise the untainted foraging target. The total length of the simulation experiment was from 8:00 AM to 4:00 PM.
In the scenario with two rather rich foraging targets to choose from (green bars in figure 10(i)), the total number of pollination flights predicted towards target A is lower than in the scenario with two poor food targets to choose from independent of the robot number. Our model predicts 12 441 pollination flights to target A, without robots advertising this target. With one robot advertising target A the number of predicted pollination flights is 1.4 times as high as without robots advertising target A (from 12 441 to 17 650 flights). If the number of active robots is quintupled from one to five robots, the number of predicted pollination flights is 2.7 times as high as with one robot. If the number of active robots is further quintupled from 5 to 25 robots, the predicted number of pollination flights is 1.8 times as high as with five robots. If the number of active robots is further quadrupled from 25 to 100 robots, the number of predicted pollination flights is 0.7 times as high as with 25 robots.
In the scenario where the choice has to be made between two targets of very different qualities (orange bars in figure 10(i)), the total number of pollination flights predicted towards target A is lower than in the previous two scenarios independent of the robot number. Our model predicts 39 pollination flights to target A, without robots advertising this target. With one robot advertising target A the predicted number of pollination flights is 88.5 times as high as the number of pollination flights (from 39 to 3452 flights) without advertising robots. If the number of active robots is quintupled from one to five robots, the number of predicted pollination flights is 2.4 times as high as with one robot. If the number of active robots is further quintupled from 5 to 25 robots, the predicted number of pollination flights is 2.1 times as high as with five robots. If the number of active robots is further quadrupled from 25 to 100 robots, the number of predicted pollination flights is 1.1 times as high as with 25 robots.
In all three scenarios from figure 10(i) our model predicts a saturation process with larger robot populations concerning the predicted pollination enhancing effect. When comparing all three scenarios from figure 10(i), our model predicts that with two rather poor foraging targets the most pollination flights towards target A will be performed with and without robots, followed by the scenario with two rather rich foraging targets and the scenario with two foraging targets of very different quality. Also, a larger robot population will lead to more foraging flights in all three scenarios. However, in all three environmental conditions our model predicts that the effect of a higher number of pollination flights towards target A will saturate with larger robot populations. Figure 10(ii) shows the predicted impact of biomimetic robots on the total influx of pesticides gathered from target B. Our model predicts that in all experimental settings the total influx of pesticides will decrease with a higher robot number.
In addition, our model predicts that if two targets with poorer quality are to be chosen from (blue bars in figure 10(ii)), the total number of tainted loads carried into the hive by bees from target B is 155 293 with no robots advertising foraging target A. One robot advertising target A, in order to reduce the pesticide influx, will be able to decrease the number of tainted loads carried into the hive by 8%, while 25 robots will be able to reduce the toxic influx by 38%, and 100 robots can reduce this influx by 43% compared to the initial value without active robots.
In the scenario with two rather rich foraging targets to choose from (green bars in figure 10(ii)), the influx of tainted loads is predicted to be higher than in the scenario with two rather poor targets to choose from, independent of the robot number. Here the total number of tainted loads carried into the hive from target B, without robots advertising target A, is 157 911.
The pesticide influx will be reduced by 3% from the initial value without active robots if one robot is advertising target A. 25 robots will reduce the toxic influx by 27%, while 100 robots will reduce the toxic influx by 38% of the initial value without active robots. In consequence, our model predicts that the same fraction of the potential toxic influx (38%) will be reduced by applying 25 waggle dancing robots in the scenario with two rather poor foraging targets offered, whereas in the scenario with two rather rich foraging targets, 100 waggle dancing robots will be needed.
In the scenario where the choice has to be made between two targets of very different qualities (orange bars in figure 10(ii)), the influx of toxic loads is predicted to be the highest amongst all tested scenarios.
Here the total number of tainted loads carried into the hive from target B without robots advertising target A is 170 516.
The toxic influx will be reduced by 2% of the initial value without active robots, when one robot is advertising target A, while 25 robots will be able to reduce the toxic influx by 16%, and 100 robots will be able to reduce the toxic influx by 27% of the initial value without active robots.
Comparing across all three scenarios shown in figure 10(ii), we find that the higher the qualitydifference of the two foraging targets is, the higher is also the influx of toxic substances into the hive. This tendency is independent of the number of applied waggle dancing robots. Our model predicts that in such a scenario far more robots will be needed to decrease the influx of tainted loads, than in scenarios with two rather poor or two rather rich foraging targets to choose from. Figure 10(iii) shows the predicted impact of biomimetic robots on the colony's energy gain in the afternoon, when the robots are active. In all experimental settings the total accumulated energy will remain almost the same except for the most extreme setting (orange bars in figure 10(iii)).
Our model predicts that if two targets with rather poor quality are available (blue bars in figure 10(iii)), the total amount of accumulated energy at noon will be 32.2 Mio. Joule with no robots advertising target A. This is 39% less than when two rather rich foraging targets (51.7 Mio. Joule) (green bars in figure 10(iii)), or two targets of very different qualities (52.2 Mio. Joule) (orange bars in figure 10(iii)) are available and no robot is advertising target A.
In conclusion, our model predicts that the foraged nectar and therefore the energetic gain of the colony will remain almost the same, independent of the number of robots advertising the always energetically less profitable target A. In the worst tested case, where two targets of very different qualities are to be chosen from and 100 robots are advertising the energetic less profitable target A (orange bars in figure 10(ii)), our model predicts that the colony's energy gain will be reduced by 5%, from 52.2 Mio. Joule to 49.5 Mio. Joule.

Discussion and conclusion
The ODE model that we created for this study showed to be useful for investigating the potential effects of biomimetic waggle dancing robots in honeybee colonies. Based on the findings of our simulation experiments, we were able to answer all our research questions as follows.
We found that a simple mathematical ODE model is capable of predicting valid and plausible foraging strategies as they are empirically observed in honeybees, see figures 4 and 5 (Q1). This model can capture the behavior of full-sized colonies under natural foraging conditions, as it is shown in figure 6 (Q2). Our model predicts that one waggle dancing robot will be able to change the collective decision making slightly, see figure 8 (Q3). With a higher number of robots, this effect will strengthen, as is shown in figure 8, but this effect will also saturate if robot numbers get very high (Q4). In addition, we found that the colony's nectar fill status has an important influence on the colony's foraging decision making independent of the use and the number of robots, as is shown in figure 9 (Q5). Our model predicts that the effect of these robots on the colony's nectar economy will be negative but only very small even with a high number of robots, as is shown in figure 10(i) (Q6). However, the waggle dancing robots will be able to lower the influx of toxic loads by advertising a safer foraging target significantly, as is shown in figure 10(ii) (Q7). The robots will also be able to significantly raise the number of pollination flights towards the advertised target, as is shown in figure 10(iii) (Q8).
A sensitivity analysis of the model shows that the general qualitative pattern of model predictions emerges also in other variations of parameter values (for more information see the supplementary information document in section G). Furthermore, we analyzed the influence of stochasticity in the forager recruitment process by adding a stochastic variable ε (t) to the probability of a bee to encounter a dance to either target A or B apart from the chance to be proportional to the numerical ratio of performed waggle dances for those two targets. The higher this noise factor gets, the more undirected the colony's collective foraging decision will be, as can be seen in section H of the supplementary information. Also, the effect of the robots on collective decision making gets lower with increasing levels of noise (for more information see supplementary information figure SH1(c)). These analyses show that the model is structurally robust against changes of parameter values and dynamically robust against noise. However, the sharpness of the collective decision making and the effectiveness of the robots in guiding the bees will be logically decreased more and more with increasing noise strength.
A prototype of such a biomimetic waggle dancing robot is currently built by a research group in Germany [26]. However, no empirical data regarding the effects of these biomimetic waggle dancing robots onto the colony-level behavior is published yet. Thus, our mathematical model is built in a first-principles way, mechanistically capturing the known involved microscopic mechanisms of biomimetically recruiting bees to foraging sites. For the basic colony's behavior without such robots, there exists empirical data [50,53], which was used to validate our base model. Our predicted foraging dynamics with the expected effects of the bio-mimetic robots can be used as a valuable piece of information in the design process of such robots and to also inform the design of the empirical experiments with such robots, as soon as those are physically available.
Therefore, the overall picture emerging from these answers to our research questions leads us to drawing the following conclusions: introducing biomimetic robots into a honeybee colony can have positive effects on the hive's safety and on the economic use of honeybees in agriculture and ecosystem support. A strong effect of the robots on the pollination service can be expected, as robots advertising specific targets will raise the total number of pollination flights at these locations significantly (figure 10(i)). Such waggle dancing robots can have a significant effect on lowering the influx of toxic loads, e.g. from areas where pesticides have been sprayed. This is achieved by advertising alternative safer foraging targets ( figure 10(ii)). However, the robots can only compete so much with the dances that the bees perform, which keeps this desired effect significant but limited. Alternative methods, such as technologymediated directed disturbance of waggle dances of bees that are advertising harmful areas to their nestmates, could have a stronger effect for this purpose.
Our model shows that there are energetic costs in the nectar economy of the colony that emerge from diverting the bees to inferior foraging targets. However, these costs are marginal and likely neglectable for the final honey yield at the end of the season, see figure 10(iii). These energetic deficits induced by robots are most neglectable in saturated colonies, which are also those where the robots showed the strongest positive effects ( figure 9).
How does the behavior of the robots differ from the behavior of animals and how is this depicted in the model? Honeybees exhibit a natural cost-versus-gainrelated dance-length correlation, which the robots do not necessarily need to adhere to. In addition, waggle dancing robots do not bring any nectar to the hive, therefore have no interaction with storage bees and do not need to be unloaded. Thus, they need to be handled by our model very differently to the living bees. In consequence, biomimetic robots can purposefully advertise lower-attractive targets with a higher dance length. Thus, they can influence the collective decision making in ways that living honeybees cannot.
In summary, our key findings indicate that not only the effectiveness of a single robot to imitate a waggle dance and the number of such robots will be important for the ultimate effect onto the colony. Our results show, that even if such a robot perfectly imitates waggle-dances, thus if it is not discriminable from a natural dancer for the bees, other factors will play important key roles: our model predicts that the colony's hunger status and the availability of alternative nectar foraging sites in the colony's surrounding environment will have a significant impact on the effectiveness of such robots.
Our research suggests the following policies for applying biomimetic waggle dancing robots: (a) the advertised alternative foraging targets should have an as-high as possible energetic value, (b) the more saturated the colonies are the better the robots will perform and (c) it will be more advantageous to use fewer robots across many hives than many robots within a few hives, as the predicted saturation effect can be prevented this way, allowing to maximize the overall effect of the available biomimetic robot population that is available (figure 8).
We chose here an ODE modeling approach on purpose. Following the parsimony principle [63], we wanted to develop a model that is flexible concerning spatial placement and dynamic qualities of food sources, similar to the existing multi-agent models HoFoSim [48], HoFoReSim [49] and FlowerSim [64], which are computationally much heavier. In parallel, we wanted to incorporate the storage processes in honeybee colonies, thus we had to further expand a previously published ODE-based model from literature [53]. As an initial step, we made a first waggle dance robot analysis [65] based on the original literature model, and then decided that this specific foraging model needs to be extended in four directions: (1) free placement of foraging targets in different distances from the hive, (2) free dynamics of sugar concentration offered at these foraging targets, (3) modeling also the storage process, the queuing delays that arise from there and the influence of these queuing delays on the waggle dancing probability, as well as (4) the ability to also represent natural foraging conditions of full-sized hives that forage on crop fields instead of artificial feeders. The model presented here was found to be in the sweet spot of model-building complexity, as it is still rather simple in formulation and computational efforts, but reliable enough to produce valid predictions (see figures 4 and 5). In the environments that we simulate with our model, there is always one clear attractor that draws the system into one single equilibrium state, which depends on the pre-given coefficients of every scenario. These equilibria are governed by the major feedback loops in the system, as they are summarized in figure 1 and, in more detail, in figure SA1 in the supplementary information. The goal of our study was to investigate, based on these feedbacks, how strongly the robots will be able to shift this attractor and to alter the equilibrium state of the system.
Mathematical modeling is important to approach these questions, as honeybee colonies are complex adaptive systems (CAS), as they are defined by [66]. This CAS of the colony is embedded into an ecosystem, which is a CAS in itself [67]. Thus, the overall system can be expected to be very complex concerning its structure, but also concerning its emerging dynamics. Thus, mathematical modeling and computer simulation are key tools for any form of ecosystem intervention. We can only analyze a certain set of foraging conditions with our simple ODE model. If a higher granularity of space is needed, e.g. to incorporate the exact size or the specific location of crop fields, a spatially finer resolved modeling approach should be taken, like a cellular model or an agent-based model. However, the model we present here would be a rich source of inspiration and starting point for such a computationally more expensive model building approach. In our opinion, the most promising candidate as a starting point and substrate for such a finer resolved analysis will be the BEEHAVE model [68] and its descendants [69].
The concept of guiding pollination service needs also to be seen from an ethics perspective. There is a lot of (well founded) criticism against the idea of replacing key organisms, like honeybee foragers, with masses of free flying micro-robots. For example, the project 'RoboBees' [70] triggered a very critical discussion of such concepts [71,72]. Such flying microrobots could hardly effectively pollinate but potentially heavily pollute the environment while raising a set of legal problems in parallel [73]. The concept of waggle dancing robots in the hive is not replacing any forager bees at all, it just adds an additional information flow to the existing ones in the colony's hive mind, allowing humans to guide the foraging targets of normal living honeybees. The reasons could be to keep them from areas that could be harmful for them, to keep natural reserve areas for wild pollinators without competition from massive honeybee foraging, to pollinate specific floral patches or crops that lack natural pollination service, to ensure certain food security aspects or many other reasons. The concept of waggle dancing robots does not contradict our (hopefully) ongoing efforts to preserve and significantly re-naturalize the environment, it rather can act as an additional tool or contingency, which may be capable of slowing down the ecological decay and buying valuable time for restoration measures this way. Such robots can also not pollute the environment, as they will always stay in the hive under human control, as they can only dance and neither crawl nor fly. Even if the whole bee colony leaves the hive, the technology will stay behind under human operator control. These robots cannot reproduce, thus also do not carry the known dangers of genetically modified organisms [74] or imported species [75] to act as surrogate pollinators. Overall, such biomimetic waggle dancing robots seem to pose little risk and might offer high potential gains.
The guiding of foraging of honeybees is not purely restricted to positive recruitment by introducing additional waggle dancing only. It could also be achieved by suppressing specific waggle dance behavior, as it is a form of motion, and recent studies have shown that motion behaviors can be suppressed by emitting vibrational stimuli with accelerometers embedded into the dance floor comb [76]. Several motion-modulating methodologies can be also combined with each other, as is discussed in a recent article about 'Ecosystem hacking' for various types of physical stimuli on combination [24] on freshly emerged bees, while other forms are discussed for the queen and the comb bees in another model-driven analysis [25]. Addressing the application of such technologies across various places in the hive and addressing several task groups and age groups of bees is also described in [31]. The effects of waggle dance robots have only been tested to a limited extent in rather unnatural conditions (small observation hives, artificial feeders [32][33][34][35]). Our model-based analysis predicts the feasibility of guiding full-sized honeybee colonies with only a few robots (but more than one), assuming that they are as attractive as a living bee performing waggle dances. The future will show through biomimetic empirical studies if robots can reach this level of attractiveness for bees and if they can be autonomous and small enough to operate in the required numbers and with the required operational times in full sized hives. This is a technological challenge, but we assume it is a crackable one, thus we think that biomimetic robots are an interesting and promising approach to achieve a long-lasting and large-scale ecosystem level effect. Given that one single colony can have up to ten thousand foragers, which are foraging within a radius of 13 km around the hive for nectar and up to 6 km for pollen [62], thus cover an area between 113 and 530 square kilometers per colony, such robotic technology can have a significant effect on the spatial and the temporal distribution of the pollination service that is provided by a colony to its surrounding ecosystem.

Data availability statement
The data that supports the findings of this study are openly available online (https://zenodo.org/record/ 7590235).