Collective thermoregulation in bee clusters

Samuel A. Ocko, L. Mahadevan


Swarming is an essential part of honeybee behaviour, wherein thousands of bees cling onto each other to form a dense cluster that may be exposed to the environment for several days. This cluster has the ability to maintain its core temperature actively without a central controller. We suggest that the swarm cluster is akin to an active porous structure whose functional requirement is to adjust to outside conditions by varying its porosity to control its core temperature. Using a continuum model that takes the form of a set of advection–diffusion equations for heat transfer in a mobile porous medium, we show that the equalization of an effective ‘behavioural pressure’, which propagates information about the ambient temperature through variations in density, leads to effective thermoregulation. Our model extends and generalizes previous models by focusing the question of mechanism on the form and role of the behavioural pressure, and allows us to explain the vertical asymmetry of the cluster (as a consequence of buoyancy-driven flows), the ability of the cluster to overpack at low ambient temperatures without breaking up at high ambient temperatures, and the relative insensitivity to large variations in the ambient temperature. Our theory also makes testable hypotheses for the response of the cluster to external temperature inhomogeneities and suggests strategies for biomimetic thermoregulation.

1. Introduction

Honeybees are masters of cooperative thermoregulation, and indeed need to be in order to survive during winter, raise their brood, cook predatory wasps, etc. They do this collectively by forming swarms wherein a fertilized queen leaves the colony with about 2000–20 000 bees, which cling onto each other in a swarm cluster, typically hanging on a tree branch, for times up to several days, while scouts search for a new hive location [1]. During this period, the swarm cluster regulates its temperature by forming a dense surface mantle that envelopes a more porous interior core. At low ambient temperatures, the cluster contracts and the mantle densifies to conserve heat and maintain its internal temperature, whereas at high ambient temperatures the cluster expands and the mantle becomes less dense to prevent overheating in the core. Over this period, when the swarm is in limbo before moving to a new home the cluster adjusts its shape and size to allow the bees to maintain and regulate the core temperature to within a few degrees of a homeostatic set point of 35°C over a wide range of ambient conditions.

The swarm cluster is able to perform this thermoregulatory task without a centralized controller to coordinate behaviour in the absence of any long-range communication between bees in different parts of the cluster [4]. Instead, this behaviour of a swarm cluster emerges from the collective behaviour of thousands of bees [5] who know only their local conditions. Early work on swarm clusters, and the related problem of winter clusters [6,7], used continuum models for variations in bee density, and temperature as determined by the diffusion of heat in a metabolically active material. Most of these models [811] assumed that the bees know their location and the size of the cluster, contrary to experimental evidence. However, a new class of models initiated by Myerscough [12] is based on local information, as experimentally observed. These models are qualitatively consistent with the presence of a core and mantle, but are unable to explain how a high core temperature persists at low ambient temperatures. Further refinements of these models that account for bee thermotaxis and also use only local information [13,14], yield the observed mantle–core formation and good thermoregulation properties at low ambient temperatures, while allowing for an increased core temperature at very low ambient temperatures, observed in some clusters [2,4,15,16]. However, the thermotactic mechanism that defines these models of bee behaviour causes the cluster to break up at moderate to high ambient temperatures, unlike what is observed.

Here, we present a model for swarm cluster thermoregulation1 that results from the collective behaviour of bees acting based on local information, yet propagates information about ambient temperature throughout the cluster. Our model yields good thermoregulation and is consistent with experiments at both high and low temperatures, with a cluster radius, temperature profile and density profile qualitatively similar to observations, without leading to cluster breakup. In §2, we outline the basic principles and assumptions behind our model. In §3, we formulate our model mathematically and characterize the parameters in it. In §4, we solve the governing equations using a combination of analysis and numerical simulation, and compare our results qualitatively with observations and experiments. We conclude with a discussion in §5, where we suggest a few experimental tests of our theory.

2. Model assumptions

Any viable mechanism for cluster thermoregulation consistent with experimental observations should have the following features: (i) the behaviour of a cluster must result from the collective behaviour of bees acting on local information, not through a centralized control mechanism. (ii) The bee density of a cluster must form a stable mantle–core profile. (iii) The cluster must expand (contract) at high (low) ambient temperatures to maintain the maximum interior ‘core’ temperature robustly over a range of ambient temperatures.

This suggests that any quantitative model of the behaviour of swarm clusters requires knowledge of the transfer of heat through the cluster, the movement of bees within the cluster, and how these fields couple to each other. The basic assumptions and principles behind our model are as follows:

  • — The only two independent fields are the packing fraction of bees in the cluster ρ (1-porosity) and the air temperature T. These then determine the bee body temperature and metabolism which can be written as a function of the local air temperature; convection of air in the form of upwards air currents depends entirely on the global bee packing fraction and temperature profiles.

  • — We treat the cluster as an active porous structure, with a packing fraction-dependent conductivity and permeability. Bees metabolically generate heat, which then diffuses away through conduction and is also drawn upwards through convection. The boundary of the cluster has an air temperature equal to the ambient temperature.

  • — Cold bees prefer to huddle densely, whereas hot bees dislike being packed densely. In addition, bees attempt to push their way to higher temperatures. The movement of bees is determined by a behavioural variable which we denote by ‘behavioural pressure’ Pb (ρ, T), which we use to characterize their response to environmental variables such as local packing fraction and temperature. In terms of this variable, we assume that the bees move from high to low behavioural pressure, a notion that is similar in spirit to that of ‘social forces’ used to model pedestrian movements [17]. Here, we must emphasize that behavioural pressure is not a physical pressure. For packing fraction to be steady, behavioural pressure must be constant throughout the cluster.

  • — We assume that the number of bees in the cluster is fixed and that the cluster is axisymmetric with spherical boundaries, whose radius R is not fixed (figure 1). At higher temperatures, the clusters become elongated and misshapen, so that the assumption is no longer accurate; however, this does not change our results for thermoregulation qualitatively. In general, to determine the shape of the cluster, we must account for both heat and force balance, but we leave this question aside in the current study.

    Figure 1.

    A cluster at an ambient temperature of (a) 11°C and (b) 27°C, approximately 12 cm across (Photo courtesy of [18]). Note that number of bees inside the cluster is nearly constant; change in cluster width is a result of changes in bee packing fraction [1]. (c) Schematic of interior (Ω), and boundary (δΩ), with mantle–core structure. Graphic are radial, vertical directions in polar coordinates, R is the cluster radius.

A model based on these assumptions can be used to study both the equilibria and dynamics. Since a swarm cluster is a constantly changing network of attachments between bees, as bees grab onto and let go of nearby bees, and can also detach and reattach themselves at different points on the surface of the cluster, these ‘microscopic’ dynamical processes allow the cluster to quickly equilibrate to changes in ambient conditions [2,4]. Thus, while we will focus on the resulting equilibria, we will also briefly consider the slow dynamical modes that allow it to respond to large-scale weak forcing, as they are particularly important in the context of cluster stability.

3. Mathematical formulation

The two independent fields in our model are bee packing fraction Embedded Image and the air temperature T(r, t); while both of these are functions of space (r) and time (t), we will focus primarily on the equilibrium behaviour of these fields. For a static cluster of bees modelled as an active porous medium that generates heat and is permeable to air, the heat generated metabolically must balance the heat lost due to conduction and convection, so thatEmbedded Image 3.1where k(ρ) is the packing fraction-dependent conductivity (power/[distance × temperature]), Embedded Image is the metabolic heat production rate of the bees per unit volume and 𝒞 is the volumetric heat capacity of air (energy/[volume × temperature]). We model the conductivity of the cluster as arising from a superposition of random convection currents within the cluster which are suppressed at high bee packing fraction and the bare conductivity of the bees treated as a solid, and approximate this by a function k(ρ) = k0(1 − ρ)/ρ (table 1). Although this form diverges as ρ → 0 and random convection currents are unsuppressed, the ρ never vanishes in the interior of the cluster, so that this limitation is not a problem. Likewise, by bounding ρ from above, we prevent k from vanishing in the cluster. We further assume that the mean flux per unit area u (distance/time) is determined by Darcy's law for the flow of an incompressible buoyant fluid through a porous medium, so thatEmbedded Image 3.2 Embedded Image 3.3Here, equation (3.2) relates u to the effects of thermal buoyancy and the pressure gradient,2 while equation (3.3) is just the incompressibility condition, with κ(ρ) the packing fraction-dependent permeability (distance2), αair is the coefficient of thermal expansion (temperature–1), γair is the specific weight of air (pressure/distance), η is the viscosity of air (pressure × time) and P is air pressure. We assume that the permeability of the cluster may be approximated via the Carman–Kozeny equation Embedded Image used to describe the permeability of randomly packed spheres [19].

View this table:
Table 1.

Table of quantities and associated units.

In general, the bee metabolic activity is not a constant, and depends on a number of factors such as temperature, age, oxygen and carbon dioxide concentration, etc. [3,4,20]. To keep our model as simple as possible, we start by assuming that the metabolic rate is temperature-independent, with Embedded Image and show that this is sufficient to ensure robust thermoregulation, setting apart the details of the calculations that show that our model can also yield robust thermoregulation with a temperature-dependent metabolic rate (appendix D).

To close our set of equations, we still must relate ρ(r) to T(r), which we do by making a hypothesis that bees respond to local packing fraction and temperature changes by changing their packing fraction to equalize an effective behavioural variable which we call the behavioural pressure. Our formulation of this behavioural pressure relies on two assumptions that are based on observations

  • (1) In their clustered state, bees have a natural packing fraction which is a function of the local temperature. This natural packing fraction decreases with increasing temperature, and increases with lower temperature, until it reaches a maximum packing fraction ρmax. Effectively, cold bees prefer to be crowded, whereas hot bees dislike being crowded, consistent with a variety of experiments in the field and in the laboratory [4,7,15,18].

  • (2) In addition to having a temperature-dependent natural packing fraction, bees also like to push their way towards higher temperatures. This will cause areas of equal local temperature to pack more densely at low ambient temperatures, consistent with observations [4].

With these constraints in mind, a minimal model for bee behavioural pressure suggests the piecewise functionEmbedded Image 3.4where ρm(T) = min{ρmax, ρ0αbeeT} is the natural packing fraction. The constant ρ0 is dimensionless and represents the baseline for natural packing fraction; αbee has units of temperature−1 and describes how the natural packing fraction changes with temperature, while χ also has units of temperature−1 and describes how bees push their way towards higher temperatures. Behavioural pressure diverges as ρρmax to enforce the maximum packing density.3

We emphasize that our model assumes that individual bees have no independent homeostatic set points for temperature or packing fraction; the behavioural pressure depends on a combination of ρ and T. Furthermore, we note that we only allow stable packing fractions (∂Pb/∂ρ > 0); thus, bee packing fraction in a cluster at equilibrium can never be lower than the natural packing fraction. In our formulation, a higher behavioural pressure at near-zero packing fractions represents the basic aggregation behaviour that creates the cluster. Additionally, the absolute behavioural pressure is of no consequence; only gradients and relative values are important.4 Further evidence for a behavioural pressure comes from experiments [6] which observed a relation between the ambient temperature and the core density, even when the core temperature remained approximately constant.

To complete the formulation of our problem, we need to specify boundary conditions for the temperature and packing fraction fields. As we have stated earlier, the surface bees will be at the ambient temperature; furthermore, as they can freely expand and contract to minimize their behavioural pressure, we assume that they will be at their natural packing fraction ρm(Ta).

Comparing our model with earlier models such as the Myerscough model [12] and the Watmough–Camazine model [13], we note that both fit into this general behavioural pressure framework (appendix B). However, our work differs from these studies in that we synthesize and generalize the implicit relation between bee behaviour and their environmental variables in terms of the behavioural pressure, with a form that combines elements from both models and is consistent with experimental observations. In addition, we account for the role of fluid flow and convection of heat in the cluster, which breaks the vertical symmetry, and is potentially relevant in heat transport at high Rayleigh number Embedded Image.

3.1. Dimensionless equations

To reduce the number of parameters in our model, we make our equations dimensionless. We note that the total number of bees in the cluster is constant; in our continuum model, this means that the total bee volume within the cluster ∫∫∫ρ dv is fixed. Rather than defining cluster size by number of bees, we define it in terms of the dimensionless total bee volume 𝒱 = ∫∫∫ρ dv/𝒱0, where 𝒱0 is the total bee volume of a typical cluster, i.e. the average volume of a bee times the average number of bees in a cluster (appendix A 1). Upon setting a characteristic length scale to be the radius of a sphere of volume Embedded Image we write the constraint on total bee volume as ∫∫∫ρ dv = (4π/3)𝒱. Scaling the temperature so that a typical ambient temperature of 15°C → 0, and the goal temperature of 35°C → 1, we use the dimensionless variablesEmbedded Image Embedded Image Embedded Image and write the dimensionless form of equations (3.1)–(3.4) asEmbedded Image 3.5Embedded Image 3.6Embedded Image 3.7Embedded Image 3.8Embedded Image 3.9Embedded Image 3.10for the packing fraction and temperature profiles of the cluster which depend on the dimensionless parameters ρmax, ρ0, αbee, χ, Ta, k0, κ0, 𝒱.5

3.2. Information transfer through equalization of behavioural pressure

On the outer boundary of the mantle, the air temperature is equal to ambient temperature, so that minimizing the behavioural pressure requires that the packing fraction at the mantle will be the natural packing fraction at the ambient temperature ρm(Ta). At equilibrium, (3.4) then implies that the behavioural pressure in the mantle and thus throughout the cluster will be −Taχ. This means that throughout the cluster, we may write the local bee packing fraction as a function of both the local temperature and the ambient temperature, i.e. Pb(ρ(T,Ta),T) = −Taχ which we solve to findEmbedded Image 3.11where we have made the substitutions c0 = αbee − χ, c1 = χ. Intuitively, the c0 term characterizes the sensitivity of core temperature to ambient temperature, but the cluster cannot fully adapt at low Ta through this term alone; adaptation at lower ambient temperatures requires the c1 term, which is eventually responsible for overheating in the core at very low Ta. In figure 2, we graph the packing fraction ρ(T,Ta) obtained by tracing contours of equal Pb which allows us to write the equilibrium local packing fraction everywhere in terms of the conditions at the boundary. Bees respond to their local conditions and move accordingly, and these variations in packing fraction propagate information about ambient temperature throughout the entire cluster without long-range communication. Although we have mapped Pb(ρ,T) → ρ(T,Ta) for just one choice of behavioural pressure, our approach will work for any equation of state Pb(ρ,T) which uniquely defines a stable packing fraction, i.e. with ∂Pb/∂ρ > 0.

Figure 2.

Graphical representation of equations (3.9)–(3.11), showing how local packing fraction is obtained through behavioural pressure. At low Ta, bees on the mantle will pack at ρ = ρmax (upper circle), and behavioural pressure is high throughout the cluster leading to higher interior packing fraction (upper solid line) and overpacking. At high Ta, bees on the mantle will pack more loosely (lower circle), and low behavioural pressure throughout the cluster leads to low interior packing fraction (lower solid line). The solid lines run along contour lines of equal Pb(ρ,T), as behavioural pressure must be uniform through the cluster. Coefficients used are ρ0 = 0.85, αbee = 0.75, ρmax = 0.8, χ = 0.3. (Online version in colour.)

Having obtained equation (3.11) from equations (3.9) and (3.10), from now on we work with equation (3.11) directly. The set of equations (3.5)–(3.8) and (3.11) with the boundary condition that the surface temperature of the cluster is the ambient temperature completes the formulation of the problem to determine ρ(r), T(r).

4. Simulations and results

While the boundary of the cluster is assumed to have spherical symmetry, the temperature and packing fraction fields inside do not need to have spherical symmetry owing to the effects of convection. However, the fields still have cylindrical symmetry, and therefore we can represent ρ(r), T(r) as ρ(s,z), T(s,z), where s is the distance from the central axis and z is the height. With this coordinate representation, we solve the governing equations (3.5)–(3.8) and (3.11) in a spherically bounded domain using a simple discretization scheme with 30 values of s and 60 values of z (appendix C.2.).

Our choice of the dimensionless parameter k0 = 0.2 is constrained by experiments [3], while we estimate κ0 = 1 though a simple calculation assuming the bees to be randomly packed spheres, and ρmax = 0.8, slightly higher than the maximum packing fraction of spheres, as bees are more flexible (appendix A.1.). However, the parameters defining bee movement and behaviour, namely c0, c1 and ρ0 are experimentally unknown. Guided by the general observation that physiological performance is often improved by changing parameters while basic mechanisms remain unchanged, we optimize these parameters to achieve robust thermoregulation, i.e. the core temperature remains close to 1(35°C) over a range of scaled ambient temperatures Ta ∈ [−0.7,0.8] corresponding to a real ambient temperature Ta ∈ [0,30]°C. For the choice of parameters ρ0 = 0.85, c0 = 0.45 and c1 = 0.3, we find that within this wide range of Ta and a factor of three in 𝒱, the dimensionless core temperature stays in the range 0.7–1.3 corresponding to a real temperature range of 29–41°C, while the core temperature itself increases monotonically with ρ0 in an analytically solvable way (appendix A.2.).

Our simulations also capture the qualitative mantle–core structure of the cluster (figure 3) with a dense mantle surrounding a sparse core. We also find that at high ambient temperatures, the cluster expands and the mantle thins, and at low ambient temperatures the cluster contracts and the mantle thickens. Furthermore, we find that the core temperature, defined as the maximum temperature of the cluster, is higher at low ambient temperatures resulting from ‘overpacking’, consistent with experiments of [2,4,15,16], and predicted earlier [13]. Finally, we note that the temperature profile is vertically asymmetric due to convection, causing the point of maximum temperature to rise above the geometric centre of the cluster, as observed in experiments [2,4].

Figure 3.

Adaptation with temperature-independent metabolism. (a) Comparison of temperature and packing fraction profiles at high (0.8) and low (–0.7) ambient temperature, where the dimensionless total bee volume 𝒱 is 1. (b) Core temperature as a function of ambient temperature and total bee volume shows that as 𝒱 increases the core temperature increases but adaptation persists over a range of Ta (also plotted to guide the eye.) (c) Cluster radius as a function of ambient temperature and total bee volume showing how clusters swell with temperature, consistent with experiment. For bee packing fraction, we choose coefficients of ρ0 = 0.85, c0 = 0.45, c1 = 0.3, ρmax = 0.8. For heat transfer, we choose coefficients of k0 = 0.2, κ0 = 1. (Online version in colour.)

5. Discussion

We have shown that it is possible to provide a self-organized thermoregulation strategy in bee clusters over a range of observed ambient temperatures in terms of a few behavioural parameters. Our theory fits into a broader framework for understanding collective behaviour where the organism responds to the environment, but in doing so, changes the local environment and its behaviour until a common steady state is reached. Here, our model takes the form of a two-way coupling between bee behaviour and local temperature and packing fraction, quantified in terms of an effective behavioural pressure whose equalization suffices to regulate the core temperature of the cluster robustly. Although our choice of the form of the behavioural pressure is likely too simplistic, it is consistent qualitatively with experimental observations, and we think it provides the correct framework within which we can start to quantify collective behaviour. Furthermore, our strategy might be useful in biomimetic settings. Our formulation for behavioural pressure shows that in a cluster, bee packing fraction should depend only on local temperature and the temperature of bees at the boundary, which effectively control the surface packing fraction. These dependencies might be measured by applying different temperatures to the surface and interior of an artificial swarm cluster. This observation provides an immediately testable prediction: a cluster may be ‘tricked’ into overpacking and overheating its core by warming the bees just below the surface, while exposing the surface bees to a low temperature to increase behavioural pressure. As pointed out earlier, experiments and observations of bee core temperature and bee packing fraction [6] are consistent with these ideas, although a direct experiment of this type does not seem to have been carried out. Preliminary analysis of winter clusters shows that bee density can be written as a function of the local temperature (A. Stabentheiner 2013, personal communication).

Our model adjusts well in response to changes in ambient temperatures, but it does not have the same level of tolerance to different total bee volume that honeybees exhibit. We have used a continuum model where the bees on the surface are exposed to ambient temperature from all sides, but in reality the first layer of bees is hotter than the ambient temperature and supports a large temperature gradient driven by heat from interior bees. This means that the surface bees feel an average temperature higher than ambient because of the interior bees, and we predict this should give a large cluster a lower behavioural pressure than a small cluster at the same ambient temperature. This should reduce the sensitivity of core temperature to total bee volume, and we have confirmed this in simulations (appendix E).

We close with a description of some possible extensions of our study. Currently, our model ignores changes in shape of the cluster associated with force balance via the role of gravity, and the associated effects on thermoregulation. In reality, the cluster is a network of connections between bees which changes shape and size due to a combination of mechanical forces and heat balance, and a complete theory must couple these two effects as well.

Our heat balance equation and estimates (appendix A.1.) suggest that while convective terms are responsible for the asymmetry in the temperature profile, they do not play an important global role in thermoregulation. However, our calculation assumes a uniform packing which does not accurately represent the microscopic structure of the cluster, and thus we may have underestimated the Rayleigh number and the importance of convection. At high ambient temperatures, swarm clusters are observed to ‘channellize’, where channels open up to ventilate. Studying a simple ‘behavioural pressure-taxis’ dynamical law (appendix F), we find no linear instability that leads to channellization without an ‘anemotaxis’ mechanism, where behavioural pressure increases with Embedded Image instead we see only one kind of linear instability which comes from the mathematical necessity of fixing cluster radius. This raises the question of whether there are anemotaxis mechanisms. If not, how can channellization result from bee-level dynamics and mechanics? We have also neglected active cooling, which includes fanning, evaporative cooling (which can give up to approx. 50°C of cooling), diffusion of heat through diffusion of bees and effects of oxygen and CO2 [3,2022]. Finally, we have also neglected any implications of bee age distribution, despite knowledge of the fact that younger bees tend to prefer the core and produce less heat, while older bees prefer the mantle and can produce more heat [4,16,18]. Accounting for these additional effects will allow us to better characterize the ecological and possibly even evolutionary aspects of thermoregulation.

Thermoregulation is a necessity for a wide variety of organisms. When achieved collectively, individuals expend effort at a cost that accrues a collective benefit. The extreme relatedness of worker bees in a cluster and near-inability to reproduce implies that the difference between the individual and the collective is nearly non-existent, so that cost and benefit are equally shared. However, many other organisms are faced with the ‘huddler's dilemma’ [23]; expending individual metabolic effort is costly, and benefits a group that is only partially related. Because genetic relatedness, metabolic costs, individual temperature and spatial positions are all easily measurable [24], a collective thermoregulatory system is an ideal context in which to study the tangible evolution of cooperation and competition by building on our current framework both theoretically and experimentally.

Funding statement

For partial financial support, we thank the Henry W. Kendall physics fellowship (S.A.O.), the MacArthur Foundation (L.M.) and Human Frontiers Science Program grant RGP0066/2012.


We thank James Makinson, Madeline Beekman and Mary Myerscough for sharing videos of swarm clusters; Anton Stabentheiner for sharing data on winter clusters; all of them as well as Brian Lee, Zhiyan Wei and Ee Hou Yong for helpful suggestions; and an anonymous referee for very helpful and detailed suggestions on the manuscript.

Appendix A. Units, parameter estimation and dimensional analysis

A.1. Estimation of heat conductivity, metabolic rate and permeability

For the metabolic rate and heat conductivity, we use estimates from the experiments of Southwick [3], where a cluster of 4250 bees (608 g) is put into a set of roughly planar, parallel honeycombs, and the temperature profile and oxygen consumption are measured. The bees are roughly uniformly distributed with a bee packing fraction ρ of about 0.5, for which the metabolic rate is roughly uniform, and the temperature is well approximated with a parabolic profile, with a temperature of 34°C at the core, and 11°C at the edge 9.5 cm from the core, parallel to the combs. The combs insulate well, so heat transfer occurs primarily in the two directions parallel to the combs. Then, within the clusterEmbedded Image Embedded Image

The oxygen consumption rate is measured to be 6.5 ml min−1, which gives a volumetric metabolism of 0.0035 W cm–3, assuming a bee specific weight of 1 g cm–3, and an oxygen to energy conversion of 3.5 ml O2 kg min−1 ≡ 0.0012 W g−1. This metabolic rate agrees well with the experiments of Heinrich [4]. We now have all the pieces to calculate the conductivity using the conduction heat balanceEmbedded Image

At ρ = 0.8, the maximum packing fraction we allow, the conductivity becomes close to the value of 2.4 × 10−3 W (cm°C)–1 for fur and feathers [25] as suggested by Southwick which gives us some level of confidence in the functional form k0(1−ρ)/ρ we have chosen for conductivity.

To estimate κ0, we use the Carman–Kozeny equation [19]. The average bee weighs about 0.14 g, which corresponds to a sphere of diameter approximately 0.65 cm. In the absence of any detailed information about the bee structure in the cluster, we treat the cluster as a system of randomly packed spheres. It would be interesting to measure and better understand convective gas and heat transfer within swarm clusters. Using this diameter in the Carman–Kozeny equation, we findEmbedded Image

A typical cluster has about 10 000 bees, which is about 1.4 kg, giving R0 = 7 cm. Plugging these values in (appendix A.2.), we find dimensionless conductivity and permeability to beEmbedded Image

A.2. Dimensional analysis

Our conditions for heat balance imply thatEmbedded Image Embedded Image Embedded Image Embedded Image while our equation for behavioural pressure readsEmbedded Image Embedded Image

We set our unit of length to be R0, so that the volume constraint becomes ∫∫∫ρ dv = (4π/3)𝒱. We make the transformation ∇ → ∇/R0. Then our heat balance equations readEmbedded Image Embedded Image Embedded Image Embedded Image

On making the substitutionsEmbedded Image leads to the goal temperature of 35°C yielding a dimensionless temperature of unity, and a typical ambient temperature of 15°C corresponding to a dimensionless temperature that vanishes. Then our system of equations becomesEmbedded Image Embedded Image Embedded Image Embedded Image Embedded Image Embedded Image

We divide all terms in the heat equation by the base metabolism Embedded Image to yieldEmbedded Image Embedded Image Embedded Image Embedded Image

Making the substitutionEmbedded Image and the substitution Embedded Image leads to the set of equationsEmbedded Image Embedded Image Embedded Image Embedded Image Finally, making the substitutions for the coefficientsEmbedded Image Embedded Image Embedded Image leads to our full dimensionless set of equations for heat balanceEmbedded Image Embedded Image Embedded Image Embedded Image while the dimensionless behavioural pressure readsEmbedded Image Embedded Image

Our model has seven parameters ρmax, ρ0, αbee, χ, Ta, k0, κ0, with an additional parameter for the total bee volume 𝒱, which varies from cluster to cluster.

A.2.1. Note on further dimensional analysis

We note that, when the metabolic rate is temperature independent, the goal temperature and the typical independent ambient temperature have no bearing on the actual behaviour of the model, only on whether it represents effective thermoregulation. Then, we may set the temperature where the packing becomes maximally dense, Embedded Image to be zero, and the temperature at which the packing fraction becomes zero, Embedded Image to be 1. We can then make slightly different substitutions for the coefficientsEmbedded Image Embedded Image Embedded Image

We may remove the parameter 𝒱 by setting the length scale to be the fully packed radius of this particular cluster rather than the fully packed radius of a typical cluster. Doing this makes the total bee volume constraint become ∫∫∫ρ dv = (4π/3) and requires the substitutionsEmbedded Image

We now find that there are now only five free parameters ρmax, χ, Ta, k0, κ0. We do not carry out this extra analysis in the main body of the paper because this causes us to lose sight of what the goal core temperature, typical ambient temperatures and typical cluster sizes are.

A.3. Units of behavioural pressure

Our model is based on behavioural pressure being uniform at equilibrium. The units and typical values of behavioural pressure are unknown: any set of dynamical equations for bee movement will result in the same equilibrium where behavioural pressure, whose units depend on the set of dynamical equations used, remains constant. For example, a simple taxis model is one whereEmbedded Image would mean that behavioural pressure has units of distance2/time. A more complicated evaporation/condensation model would have the formEmbedded Image where σ is the evaporation and condensation radius and 𝒯 is some transfer coefficient would give units of 1/time. A yet more complicated model involving mechanical compressibility or viscosity would have yet another set of dimensions for behavioural pressure. All of these models would, however, yield the same static solution.

Appendix B. Behavioural pressure formalism and its antecedents

To understand how our behavioural pressure formalism fits in with previous models, we compare them within this framework. We note that previous models have defined density in terms of bees/cm3 instead of packing fraction; to go between the two, bees may be assumed to have water density, and a packing fraction of 1 corresponds to (1 g)/mbee bees/cm3.

The Myerscough model assumes that the bee density depends only on local temperature, and thus can be written as Pb(ρ,T) = |ρρm|, where Embedded Image.

The Watmough–Camazine model defines a dynamical law for the bee density via the equationsEmbedded Image Embedded Image where μ(ρ) > 0 is a motility function and χ(T) is a thermotactic function. This may be written asEmbedded Image where the behavioural pressure Pb is defined asEmbedded Image

We note that as μ(ρ) > 0 for all ρ, the density component is minimized as ρ → 0, and the density at the surface must be fixed to prevent the cluster from falling apart, unlike in our behavioural pressure framework. By contrast, the Watmough–Camazine behavioural pressure allows the packing fraction at the mantle to become too high. Additionally, because the behavioural pressure can be divided into temperature and density components, the point of highest density will always be at the same temperature, regardless of size or ambient temperature, which is not observed experimentally. This is in contrast with our formalism.

Appendix C. Numerical methods

C.1. Method of solving for equilibrium

  • To reach equilibrium, we use an iterative scheme, described by the following pseudocode:

    • Embedded Image 

    • Embedded Image 

    • Embedded Image

  • repeat

    • Solve for Embedded Image at fixed temperature and density, then solve for T at fixed Embedded Image and Embedded Image

    • Embedded Image

    • Embedded Image Embedded Image 

    • Embedded Image 

    • Embedded Image 

    • Embedded Image 

    • Expand or contract the bee packing fraction and temperature profiles to normalize total bee volume

    • Embedded Image

    • Embedded Image

    • Embedded Image

  • Embedded Image

  • until converged

The intermediate steps can be solved as a system of linear equations. Note that this method does not add or remove cells to vary cluster radius and conserve the total number of bees; it grows and shrinks a fixed number of cells. The solution is considered to be converged when ρnew = ρ, Tnew = T, = 1 to within 10−10, which takes about 100–200 iterations, about a minute on a laptop. All simulations were done using MATLAB (Source code at

C.2. Discretization of space

To solve for the temperature and density profiles, we must first discretize the system. While spherical symmetry is broken due to convection, the system retains rotational symmetry about its axis. We therefore use cylindrical coordinates, where each cell is given indices (i, j), and has coordinates which represent the distance from the central axis sij and the vertical coordinate zij, whereEmbedded Image where n is the radius of the cluster in cells. All cells with Embedded Image are in the interior of the cluster, while all cells with Embedded Image are at the exterior of the cluster (figure 4), subject the the boundary conditions Tij = Ta, Pij = 0, ρij = 0.

Figure 4.

Cell layout of a cluster 30 cells in radius. Interior cells are coloured white, exterior cells are coloured grey.

The volume of each cell with coordinates (i, j) is Embedded Image where w = R/n is the width of each cell. Each cell (i, j) neighbours four other cells, (i, j + 1), (i, j−1), (i + 1, j), with the exception of cells which border the axis (i = 0), which only have three neighbours. The area shared by cell (i, j) and its outside horizontal neighbour (i + 1, j) is Embedded Image The area shared by cell (i, j) and its vertical neighbour (i, j + 1) is 2πwsij.

C.2.1. Heat equations

For heat balance, we discretize our (now dimensionless) heat equations as

Embedded Image

The first term corresponds to metabolic heat generation, and the second term is heat conduction, where Embedded Image is the sum of all (i′, j′) neighbouring (i, j), with the heat conductance between two neighbouring cells depending on the harmonic mean of the conductance of each cell; Embedded Image The third term represents convective heat transfer, with uij,i'j’ the air flow from cell (i, j) to (i′, j′) (Units of dimensionless volume/time due to discretization).

We have chosen an ‘upwinding’ scheme [26]; when air flows out of a cell, the outwards heat flux is determined by the temperature of that cell. When air flows into a cell from a neighbouring cell, the inwards heat flux is determined by the temperature of the neighbouring cell where the air originates. This scheme uses the Heaviside step functionEmbedded Image

C.2.2. Solving for buoyancy-driven flow

The air flow from cell (i, j) to neighbouring cell (i′, j′) is given byEmbedded Image where the air conductance between two cells again depends on the harmonic mean of the permeability of each cell, and the pressure is set so that air flow is conserved in every cell i.e. Embedded Image for all cells (i, j). This yields a set of linear equations that can be solved easily.

Appendix D. Role of temperature-dependent metabolic rate

To formulate our temperature-dependent metabolic rate, we note that bees have a higher base metabolic rate at high temperatures than at low temperatures. Our metabolic rate for high temperatures comes from experiments involving oxygen consumption in swarm clusters [2]. At moderate temperatures, bees on the mantle keep their body temperatures approximately 3°C above ambient temperature, and at below 15°C, they will ‘shiver’ to keep their body temperature at 18°C [2]. Assuming a constant coefficient of thermal transfer between a bee and the surrounding air, this leads to a formulation of metabolism by shivering at air temperatures below 15°C, and gives us the full piecewise function for metabolic rate (figure 5).Embedded Image where the base metabolism, Embedded Image has units of power/volume.

Figure 5.

Metabolic rate as a function of temperature.

The simulated cluster with temperature-dependent metabolic rate has the same qualitative features as when we use temperature-independent metabolic rates. We set Embedded Image to be 0.00175 W cm−3, half as high as when Embedded Image was set temperature-independent, because now it represents the minimal metabolic rate, not the uniform one (figure 6). Using dimensional analysis, this results in the values of k0, κ0 being doubled, giving k0 = 0.4, κ0 = 2.

Figure 6.

Adaptation with temperature-dependent metabolic rate. Note that cluster behaviour is nearly the same as for constant metabolism (cf. figure 3). (a) Comparison of temperature and packing fraction profiles at high (0.8) and low (–0.7) ambient temperature, where the dimensionless total bee volume 𝒱 is 1. (b) Adaptation is nearly the same as for temperature-independent metabolic rate. Ta is also plotted to guide the eye. (c) Cluster radius is also nearly the same as for temperature-independent metabolic rate. For packing fraction, we choose coefficients of ρ0 = 0.85, c0 = 0.5, c1 = 0.25, ρmax = 0.8. (Online version in colour.)

Appendix E. Role of finite bee size on thermoregulation

In this paper, we have defined the temperature at the boundary to be the ambient temperature, and we assumed the surface bees feel the ambient temperature, which sets the behavioural pressure accordingly. In reality, bees point their heads inwards, and feel a temperature gradient driven by the heat produced by interior bees. Therefore, it may be more realistic that the behavioural pressure is set by the temperature a slight distance inwards from the surface. If we include the effects of convection, this implies that spherical symmetry must be broken, and the temperature becomes not just a function of distance, but also dependent on angle. Therefore, to close our set of equations without having to delve deeper into the question of cluster shape, we must neglect upwards convection and only consider conduction. This gives us the system of equationsEmbedded Image Embedded Image where the behavioural pressure is now set by the temperature a distance of Lbee inside the cluster rather than by the ambient temperature, so that the equation for density isEmbedded Image E1Assuming Lbee = 1 cm, slightly shorter than the body length of a worker bee, for an average cluster radius of 7 cm we find that the dimensionless Lbee of about 0.14. We simulate the system for dimensionless total bee volumes of 0.5, 1 and 3, with the same parameters of c0 = 0.45, c1 = 0.3 and ρmax = 0.8 as were used in the paper. In the first system, to test thermoregulation with this effect, we choose Lbee = 0.14, and compensate for the slightly lower average behavioural pressure by increasing ρ0 to 0.95. In the second system, we test thermoregulation without this effect, and so we choose Lbee = 0, ρ0 = 0.85, as we did in the main body. We find that for large clusters, the temperature gradient created by the interior bees lowers behavioural pressure and loosens the cluster (figure 7). This mitigates overheating in the core and sensitivity to total bee volume, e.g. at Ta = −0.7, Tcore varies approximately 50% more with cluster size when behavioural pressure is set by ambient temperature rather than the temperature beneath the surface. We note that in the models of Myerscough [12] and Watmough–Camazine [13], a similar shielding of surface bees from ambient air and reduction of sensitivity to cluster size was achieved through use of a heat transfer coefficient, with units of power/(area temperature) between the cluster surface and ambient air.

Figure 7.

Comparison of core temperatures. (a) Behavioural pressure is set by the temperature slightly below the surface of the cluster, and sensitivity to total bee volume is mitigated. (b) Behavioural pressure is set by ambient temperature, and sensitivity of core temperature to total bee volume is considerably higher. Ta is also plotted in each to guide the eye. Note that (b) is very similar to figure 3b, as it results from the same equations except for convection. (Online version in colour.)

Appendix F. Linear stability of cluster

Solving for linear stability using simple ‘behavioural pressure-taxis’ dynamics (appendix F.1.), we find that all clusters simulated at a temperature-independent metabolic rate are stable. However, at low ambient temperatures, clusters with a temperature-dependent (appendix D) metabolic rate can be linearly unstable via an overheating instability (figure 8), which we believe to be a relic of fixing the boundaries of the cluster. In this instability, bees from the core move to the mantle, increasing the mantle thickness and insulation. This causes the core to heat up, increasing its behavioural pressure causing even more bees to move from the core to the mantle, leading to eventual runaway. We only see this instability when using a temperature-dependent metabolic rate, where an increased core temperature results in a greater net metabolic rate, aggravating the problem.

Figure 8.

(a) Cluster profile. (b) Contour plots of unstable eigenvector of linear response matrix. Circle represents boundaries of the cluster, and a temperature-dependent metabolic rate is used. J0 = 1, 𝒱 = 1.5, Ta = −0.7. ρ0 = 0.85, c0 = 0.5, c1 = 0.25, ρmax = 0.8. For heat transfer, we choose coefficients of k0 = 0.4, κ0 = 2. Note that the packing fraction very close to the boundary is fixed because ρ = ρmax. (Online version in colour.)

Dynamical behaviour that allows the cluster radius to vary would suppress this instability, as bees moving from the core to the mantle would result in an expansion of the cluster, increasing the surface area and cooling the core. However, a dynamical model which allows the boundaries of the cluster to change requires a better understanding of the bee-level structure and the mechanics within a swarm cluster, and we leave this aside here.

F.1. Methods for calculating linear stability

Having solved for the equilibrium state, we want to find if this state is stable or unstable. To do so, we must define some dynamical laws for bee movement. Choosing a simple ‘behavioural pressure-taxis’ behaviour, where Embedded Image gives us the complete set of dynamical equationsEmbedded Image Embedded Image Embedded Image

Here we vary J0 over a wide range to reflect the large variations in bee movement and temperature timescales, and define F, G to be the functionals that determine the dynamics of the system. We also emphasize some issues about these choices: (i) our form assumes a substrate that the bees move on. In reality, in a swarm cluster the bees are the substrate, and changes on one side of the cluster propagate mechanically to other parts of the cluster at a rate faster than the taxis rate. (ii) Bee movement is not necessarily a local movement down pressure gradients. Bees can disconnect from the cluster and reattach at a different point, which does not fit into the local gradient picture. (iii) There is a discontinuity in behavioural pressure and packing fraction at the boundary of the cluster, so this taxis model does not explain how the boundary of the cluster can change. Within these limits then, this choice of behaviour gives us a full set of differential equations. To determine whether the equilibrium state of the system is linearly stable, we must determine whether the linear response matrixEmbedded Image has positive eigenvalues. We note that the equilibrium temperature and bee packing fraction profile we solved for is symmetric with respect to rotation about the central axis (ϕ direction), and therefore we may partition the space of perturbations into subspaces defined by wavenumber kϕ; Embedded Image These temperature and bee packing fraction perturbations will in turn give a change in airflow, pressure and behavioural pressure Embedded Image Embedded Image All of these will change the temperature and bee packing fraction time derivatives which will be proportional to Embedded Image. For each wavenumber kϕ, we construct the stability matrix and study its spectrum. To first orderEmbedded Image Embedded Image Embedded Image where ∇sz is the gradient in the s and z directions, Embedded Image Regions where ρ = ρmax are locked at ρmax and not allowed to vary in bee packing fraction.

F.1.1. Solving for Δu, ΔP

The above equations depend on Δu, which is determined byEmbedded Image Embedded Image must have the formEmbedded Image where Embedded Image is the sz component of Embedded Image. We solve for pressure Embedded Image using the condition Embedded Image Embedded Image At kϕ = 0, Embedded Image this condition simply becomes Embedded Image 

Therefore, at a given ΔT, Δρ, we can solve Δu, ΔP from this set of linear equations.

F.1.2. Numerical computation of linear response matrix

To solve for stability, we discretize the system in the s, z directions in the same way that we did when solving for the equilibrium state. Because we are only solving for linear stability and the system starts off uniform in the ϕ direction, we do not need to discretize in the ϕ direction. For perturbations at a certain wavenumber kϕ, the temperature derivative is, to first order

Embedded Image

Note the slight modification of the upwinding terms, where we have also included a ϕ component to represent the influx or outflux of air in the ϕ direction.

The density derivative is, to first order

Embedded Image

The airflow is solved using the set of linear equationsEmbedded Image Embedded Image δP is set such that the divergence in the ϕ direction negates the divergence in the s,z directions,Embedded Image for all cells (i,j).


  • 1 While we do not consider the related problem of winter clusters here, a behavioural pressure formalism should also be applicable in this case.

  • 2 We note that in Darcy's law, a buoyant ‘body force’ may be added to the gradient of pressure in the same way that is done for the full Navier–Stokes equations.

  • 3 This discontinuity is not a problem in practice as our simulations do not work with behavioural pressure directly.

  • 4 The magnitude of these gradients and differences are not important, and allow us to omit a prefactor from |ρρm(T)|.

  • 5 Note this can be further reduced to five dimensionless parameters as we do in the appendix, but this causes us to lose sight of what the goal temperature, typical ambient temperatures and typical cluster sizes are.

  • Received November 6, 2013.
  • Accepted November 21, 2013.


View Abstract