Abstract
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 buoyancydriven 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 longrange 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 [8–11] 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 thermoregulation^{1} 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 ρ (1porosity) 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 fractiondependent 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’ P_{b} (ρ, 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.
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 largescale 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 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 that 3.1where k(ρ) is the packing fractiondependent conductivity (power/[distance × temperature]), 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(ρ) = k_{0}(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 that 3.2 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 fractiondependent permeability (distance^{2}), α_{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 used to describe the permeability of randomly packed spheres [19].
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 temperatureindependent, with 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 temperaturedependent 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 temperaturedependent 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 function 3.4where ρ_{m}(T) = min{ρ_{max}, ρ_{0} − α_{bee}T} 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 (∂P_{b}/∂ρ > 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 nearzero 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}(T_{a}).
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 .
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 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 variables and write the dimensionless form of equations (3.1)–(3.4) as 3.5 3.6 3.7 3.8 3.9 3.10for the packing fraction and temperature profiles of the cluster which depend on the dimensionless parameters ρ_{max}, ρ_{0}, α_{bee}, χ, T_{a}, k_{0}, κ_{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}(T_{a}). At equilibrium, (3.4) then implies that the behavioural pressure in the mantle and thus throughout the cluster will be −T_{a}χ. 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. P_{b}(ρ(T,T_{a}),T) = −T_{a}χ which we solve to find 3.11where we have made the substitutions c_{0} = α_{bee −} χ, c_{1} = χ. Intuitively, the c_{0} term characterizes the sensitivity of core temperature to ambient temperature, but the cluster cannot fully adapt at low T_{a} through this term alone; adaptation at lower ambient temperatures requires the c_{1} term, which is eventually responsible for overheating in the core at very low T_{a}. In figure 2, we graph the packing fraction ρ(T,T_{a}) obtained by tracing contours of equal P_{b} 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 longrange communication. Although we have mapped P_{b}(ρ,T) → ρ(T,T_{a}) for just one choice of behavioural pressure, our approach will work for any equation of state P_{b}(ρ,T) which uniquely defines a stable packing fraction, i.e. with ∂P_{b}/∂ρ > 0.
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 k_{0} = 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 c_{0}, c_{1} 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 T_{a} ∈ [−0.7,0.8] corresponding to a real ambient temperature T_{a} ∈ [0,30]°C. For the choice of parameters ρ_{0} = 0.85, c_{0} = 0.45 and c_{1} = 0.3, we find that within this wide range of T_{a} 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].
5. Discussion
We have shown that it is possible to provide a selforganized 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 twoway 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 pressuretaxis’ dynamical law (appendix F), we find no linear instability that leads to channellization without an ‘anemotaxis’ mechanism, where behavioural pressure increases with 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 beelevel 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 CO_{2} [3,20–22]. 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 nearinability to reproduce implies that the difference between the individual and the collective is nearly nonexistent, 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.
Acknowledgements
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 cluster
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 O_{2} 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 balance
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 k_{0}(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 find
A typical cluster has about 10 000 bees, which is about 1.4 kg, giving R_{0} = 7 cm. Plugging these values in (appendix A.2.), we find dimensionless conductivity and permeability to be
A.2. Dimensional analysis
Our conditions for heat balance imply that while our equation for behavioural pressure reads
We set our unit of length to be R_{0}, so that the volume constraint becomes ∫∫∫ρ dv = (4π/3)𝒱. We make the transformation ∇ → ∇/R_{0}. Then our heat balance equations read
On making the substitutions 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 becomes
We divide all terms in the heat equation by the base metabolism to yield
Making the substitution and the substitution leads to the set of equations Finally, making the substitutions for the coefficients leads to our full dimensionless set of equations for heat balance while the dimensionless behavioural pressure reads
Our model has seven parameters ρ_{max}, ρ_{0}, α_{bee}, χ, T_{a}, k_{0}, κ_{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, to be zero, and the temperature at which the packing fraction becomes zero, to be 1. We can then make slightly different substitutions for the coefficients
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 substitutions
We now find that there are now only five free parameters ρ_{max}, χ, T_{a}, k_{0}, κ_{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 where would mean that behavioural pressure has units of distance^{2}/time. A more complicated evaporation/condensation model would have the form 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/cm^{3} 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)/m_{bee} bees/cm^{3}.
The Myerscough model assumes that the bee density depends only on local temperature, and thus can be written as P_{b}(ρ,T) = ρ−ρ_{m}, where .
The Watmough–Camazine model defines a dynamical law for the bee density via the equations where μ(ρ) > 0 is a motility function and χ(T) is a thermotactic function. This may be written as where the behavioural pressure P_{b} is defined as
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:
repeat
Solve for at fixed temperature and density, then solve for T at fixed and
Expand or contract the bee packing fraction and temperature profiles to normalize total bee volume
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} = ρ, T_{new} = 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 http://web.mit.edu/socko/Public/PublishedCode/PRS2013BeePaperSimulations.zip).
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 s_{ij} and the vertical coordinate z_{ij}, where where n is the radius of the cluster in cells. All cells with are in the interior of the cluster, while all cells with are at the exterior of the cluster (figure 4), subject the the boundary conditions T_{ij} = T_{a}, P_{ij} = 0, ρ_{ij} = 0.
The volume of each cell with coordinates (i, j) is 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 The area shared by cell (i, j) and its vertical neighbour (i, j + 1) is 2πws_{ij}.
C.2.1. Heat equations
For heat balance, we discretize our (now dimensionless) heat equations as
The first term corresponds to metabolic heat generation, and the second term is heat conduction, where 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; The third term represents convective heat transfer, with u_{ij,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 function
C.2.2. Solving for buoyancydriven flow
The air flow from cell (i, j) to neighbouring cell (i′, j′) is given by 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. for all cells (i, j). This yields a set of linear equations that can be solved easily.
Appendix D. Role of temperaturedependent metabolic rate
To formulate our temperaturedependent 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). where the base metabolism, has units of power/volume.
The simulated cluster with temperaturedependent metabolic rate has the same qualitative features as when we use temperatureindependent metabolic rates. We set to be 0.00175 W cm^{−3}, half as high as when was set temperatureindependent, because now it represents the minimal metabolic rate, not the uniform one (figure 6). Using dimensional analysis, this results in the values of k_{0}, κ_{0} being doubled, giving k_{0} = 0.4, κ_{0} = 2.
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 equations where the behavioural pressure is now set by the temperature a distance of L_{bee} inside the cluster rather than by the ambient temperature, so that the equation for density is E1Assuming L_{bee} = 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 L_{bee} of about 0.14. We simulate the system for dimensionless total bee volumes of 0.5, 1 and 3, with the same parameters of c_{0} = 0.45, c_{1} = 0.3 and ρ_{max} = 0.8 as were used in the paper. In the first system, to test thermoregulation with this effect, we choose L_{bee} = 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 L_{bee} = 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 T_{a} = −0.7, T_{core} 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.
Appendix F. Linear stability of cluster
Solving for linear stability using simple ‘behavioural pressuretaxis’ dynamics (appendix F.1.), we find that all clusters simulated at a temperatureindependent metabolic rate are stable. However, at low ambient temperatures, clusters with a temperaturedependent (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 temperaturedependent metabolic rate, where an increased core temperature results in a greater net metabolic rate, aggravating the problem.
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 beelevel 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 pressuretaxis’ behaviour, where gives us the complete set of dynamical equations
Here we vary J_{0} 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 matrix 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_{ϕ}; These temperature and bee packing fraction perturbations will in turn give a change in airflow, pressure and behavioural pressure All of these will change the temperature and bee packing fraction time derivatives which will be proportional to . For each wavenumber k_{ϕ}, we construct the stability matrix and study its spectrum. To first order where ∇_{sz} is the gradient in the s and z directions, 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 by must have the form where is the sz component of . We solve for pressure using the condition At k_{ϕ} = 0, this condition simply becomes
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
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
The airflow is solved using the set of linear equations δP is set such that the divergence in the ϕ direction negates the divergence in the s,z directions, for all cells (i,j).
Footnotes
↵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.
 © 2013 The Author(s) Published by the Royal Society. All rights reserved.