THE DYNAMICS AND STRUCTURE OF GROUPS: T W O C A S E STUDIES O F T H E C O M M O N H O N E Y B E E By James Watmough B.A.Sc (Engineering Physics) University of British Columbia M . Sc (Mathematics) University of British Columbia A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T T H E REQUIREMENTS OF FOR T H E D E G R E E OF D O C T O R OF PHILOSOPHY in T H E FACULTY OF G R A D U A T E STUDIES D E P A R T M E N T OF MATHEMATICS AND T H E INSTITUTE OF APPLIED MATHEMATICS We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA September 1996 © James Watmough, 1996 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of B r i t i s h C o l u m b i a , I agree that the L i b r a r y shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics and the Institute of A p p l i e d Mathematics The University of B r i t i s h C o l u m b i a 2075 Wesbrook Place Vancouver, C a n a d a V 6 T 1Z1 Date: abstract I study the problem of self-organized pattern formation by groups of organisms. "Self- organized" groups are coordinated by interactions between all individuals i n the group rather than by a leader or hierarchical organization. The information and stimuli necessary to coordinate individual activities are communicated through the group by these interactions. I begin by describing a modelling philosophy where as many details of the individual behaviours are included i n the model as are necessary for both biological and mathematical completeness, but no group behaviours are fed into the model. T w o specific phenomena are addressed using this approach. The first phenomenon is honey bee cluster thermoregulation. I address the question "can cluster thermoregulation be explained using only the observed behaviours of individual bees?" and show that a temperature dependent behaviour of the bees is sufficient to produce the observed global thermoregulation of the cluster. Previous models of cluster thermoregulation are discussed i n light of this modelling approach. The model also is able to make testable predictions about the density profiles of the cluster. The second phenomenon modelled is the transmission of a pheromone through a honey bee colony. The problem is of interest to beekeepers, since the pheromone is known to suppress swarming. I answer the question "How does congestion w i t h i n a honey bee colony affect the transmission of a pheromone through the hive?" T h i s question is central to a current hypothesis on the connection between the queen's ability to produce pheromones, colony size and congestion, and swarming. The cluster thermoregulation model involves a behavioural dependence on the cluster ii temperature. In contrast, the transmission model involves a direct exchange of a chemical pheromone between the bees. Thus, w i t h the thermoregulation model I examine interactions mediated through an environmental state, and w i t h the transmission model I examine direct interactions between individuals. I study these models using both Eulerian and Lagrangian models. The Lagrangian model is studied using cellular automata simulation. The Eulerian model consists of a system of partial integro-differential equations and is studied using a combination of numerical simulations and mathematical analysis. The results of the models suggest that the two interactions produce similar group properties. In the thermoregulation model, information about the global structure of the group diffuses through the cluster as thermal energy. In the transmission model, pheromone diffuses through the colony by the random motions of the bees. iii Table of Contents abstract ii List of Tables vii List of Figures viii acknowledgements ix 1 Introduction 1 2 Background 4 2.1 Group structure and group dynamics 4 2.2 A brief survey of self-organizing systems 7 2.3 M o d e l l i n g group dynamics 9 2.3.1 Lagrangian and Eulerian models 9 2.3.2 A simple random walk .' 12 2.3.3 Density dependent diffusion and chemotaxis 14 2.3.4 R a n d o m walks w i t h pheromone exchange 16 3 The Thermoregulation Model 18 3.1 . Introduction 18 3.2 Background 20 3.2.1 Swarm and winter cluster thermoregulation 20 3.2.2 M i n i m a l requirements for a satisfactory model 21 iv v 3.2.3 3.3 3.4 4 Previous models 23 The M o d e l 24 3.3.1 Assumptions of the model 24 3.3.2 Derivation of the M o d e l 26 Simulations 36 3.4.1 Non-dimensionalization 36 3.4.2 Implementation of the constraint 38 3.4.3 Results of the simulations 39 3.5 Analysis of the steady state 48 3.6 Discussion 52 The Pheromone Transmission Model 55 4.1 Introduction 55 4.2 Background 58 4.2.1 Reproductive swarming, supersedure, and absconding 58 4.2.2 The honey bee queen substance 59 4.3 4.4 4.5 Determination of the transmission rate constants 62 4.3.1 Summary of the transmission pathways 62 4.3.2 Calculation of the transmission rate constants 65 The Lagrangian model 71 4.4.1 Development of the model 72 4.4.2 The cellular automata simulation 80 4.4.3 Results of the cellular automata simulations 84 The Eulerian model 91 4.5.1 Derivation of the Eulerian model 93 4.5.2 Nondimensionalization of the system 98 v 4.5.3 4.6 4.7 5 Analysis •• 100 Predictions of the models 109 4.6.1 The Eulerian model 109 4.6.2 Comparison with the Lagrangian simulation 112 Discussion 113 Discussion 115 5.1 Self-organizing group dynamics . 116 5.2 Cellular automata simulations 117 5.3 Future work 118 List of Notation 129 vi List of Tables 3.1 Thermoregulation model parameters and their values 35 4.2 Transmission rates 66 4.3 Pheromone transmission model parameters 81 vii List of Figures 3.1 P r i m a r y features of a thermoregulating cluster 19 3.2 Thermotactic Velocity 28 3.3 Metabolic heat production per bee versus temperature 32 3.4 Cluster Temperatures 40 3.5 Cluster Radius 42 3.6 T o t a l Cluster Metabolic Output 43 3.7 Metabolic O u t p u t per Bee 45 3.8 Temperature and Density Profiles 46 3.9 Metabolic O u t p u t ' 47 3.10 Dependence of profiles on X /D 49 3.11 Plots oip{T,T ) 51 4.1 Pheromone transmission pathways 63 4.2 Pheromone D i s t r i b u t i o n : Uncongested Colony 85 4.3 Pheromone Distribution: Congested Colony 86 4.4 Pheromone distribution histograms: congested vs. uncongested colonies . 87 4.5 W a x Pheromone Distribution: Congested Colony 89 4.6 W a x Pheromone Distribution: Uncongested Colony 90 4.7 Initial spread of pheromone 92 4.8 Worker mean pheromone level 0 a R 108 acknowledgements I wish to thank Scott Camazine for his direction on the thermoregulation model, Mary Myerscough for many helpful comments, Mark Winston and Keith Slessor for their support of the Q M P model, and especially Leah Edelstein-Keshet for endless hours of guidance and careful editing. I also owe many thanks to Brian Wetton for his suggestions on the methods used for the numerical simulations of the thermoregulation model. Chapter 1 Introduction Ecology presents us w i t h a vast array of patterns in space, time, and virtually every other aspect. A l t h o u g h patterns such as groupings within a population are often a result of either random dispersion or inhomogeneities i n the environment, many organisms form groups even i n homogeneous environments. Further, many of these groups appear to be formed and maintained by communication and interaction between individuals, rather than through responses of individuals to the actions of a leader or through any form of hierarchical control. M a n y complex activities of groups can be explained using individual behaviours based on (spatially) local conditions. Organism grouping has very important ecological consequences. If populations are not well mixed, as is the case when organisms cluster into groups, then the average density of the entire population and the local density seen by an individual may be vastly different. These effects can be seen i n the distribution of species which are subject to predator-prey interactions (Turchin, 1989) or density dependent growth rates (Harada et a l . , 1995) and i n the spread of epidemics (Sato, 1994). In this thesis, I examine the connection between the defining characteristics of a group, such as its size, motion, and function, and the interactions and behaviours of the individuals that compose it. I focus on two specific examples of honey bee group organization. However, the modelling approach is applicable to grouping i n more general situations. The honey bee is chosen as a study organism since b o t h the individual bees and the group are easily observable. Further, the social behaviour of the honey bee 1 Chapter 1. Introduction 2 (unlike most other species) can be observed both in a laboratory setting and in the field. Although most of the experiments referred to in this thesis were conducted on Apis mellifera L . , a race of the European honey bee, the phenomena and results likely are relevant to other social insect species. The following phenomena are investigated by the models: 1. Honey bee cluster thermoregulation: A t cool temperatures, honey bees form into clusters. The bees are able to precisely regulate the temperatures within a cluster, apparently without communication. The bees respond to decreases in their local temperature by huddling and shivering. These behaviours give rise to a thermoregulation of the cluster. Thus, cluster thermoregulation is an aspect of the group that can be explained by the behaviour and interactions of individuals. 2. The transmission of honey bee queen substance through a colony: Healthy honey bee queens produce a substance that is transferred through the colony primarily by direct worker-worker exchange. Workers respond to a decrease in the substance by rearing new queens and beginning preparations for reproductive swarming. The transmission is thought to be hindered by colony congestion. Hence, colonies will rear young queens and swarm when they grow too large. The distribution of the queen substance depends on the behaviour of the individuals. Each of the above examples stresses a particular form of interaction. In the first example, the bees do not interact directly, but rather the behaviour of a bee depends on its local temperature, which is influenced by the behaviour of the bee's neighbours. The temperature at the centre of the cluster appears to be unknown to the bees on the cluster periphery. Nor is the ambient temperature surrounding the cluster detectable by the bees at the cluster centre. Despite the fact that the bees in the cluster respond solely to their local conditions, the cluster is regulated as a unit. The diffusion of thermal Chapter 1. Introduction 3 energy through the cluster allows information about the response of individual bees to be communicated to their neighbours. It also provides the feedback that results i n cluster thermoregulation. In the second example, the bees interact directly by exchanging the queen substance. In both examples, the behaviour of the individuals depends only on local conditions. There appear to be no central organizing structures or pre-patterns, nor is there a control hierarchy amongst the individuals. Each individual follows its own agenda. The purpose of the models is to predict the observed behaviour of a group from the known behaviours of the individuals. More specifically, I am interested i n predicting how changes i n the behavioural patterns of the individuals produce changes i n the behavioural patterns of the group. For this reason, all assumptions made are based on the behaviour of the individuals. T h e models are studied analytically and through numerical simulations. B o t h cellular automata simulations and discretizations of partial differential equations are used. T h e results of the studies show that the principle of self-organization is relevant to the group dynamics of honey bee colonies. This means that the structure of the group can be predicted from the dynamics of the individuals. In Chapter 2, I give a brief survey of groups and the models used to study their formation. T h e emphasis is placed on how groups are formed rather than their function. In b o t h the systems discussed i n this thesis I argue that the group is self-organized. Current research on this hypothesis is reviewed i n Section 2.2. In Section 2.3 I present the methods that have been used to model group dynamics. The details of the model for honey bee cluster thermoregulation are presented i n Chapter 3, and those for the model for the transmission of the honey bee queen substance through a colony i n Chapter 4. T h e thesis closes w i t h a discussion of implications to modelling group dynamics (Chapter 5), and a discussion of future work. Chapter 2 Background 2.1 G r o u p structure and group dynamics Groups range from purely spatial aggregations such as herds and swarms, to aggregations in b o t h space and direction of motion such as fish schools, and to societies which may show patterns i n the age, size, and task distributions of their members. There are many distinct mechanisms that can result i n the formation and maintenance of a group. These mechanisms can be divided into external organizing forces that shape the structure of the group and internal communications that allow the group to organize or to change its structure i n response to external conditions. Thus, groups should be characterized not only by their global structure, but also by the interactions and the communication between the individuals in the group. In this thesis I w i l l use the term g r o u p d y n a m i c s to denote the communication and interactions between the individuals in the group, and g r o u p s t r u c t u r e to denote the defining characteristics of the group. For example, the group structure of a fish school is the aggregation i n space and angle, whereas the group dynamics are the responses of the individuals to their neighbours. These group dynamics result i n a local alignment and spacing of the fish and are the mechanisms that maintain the global group structure. The group structure is best described at the scale of the group, whereas the group dynamics are described at the scale of the individual. T h i s distinction can be illustrated by the following examples of swarming, herding, 4 Chapter 2. Background 5 and schooling behaviour. Swarms of midges and honey bees both form because each i n dividual is attracted to the same point. Midges form swarms over bright objects (Okubo, 1980), and honey bees form swarms around their queen. However, not a l l groups are formed by external influences. Indeed, even groups that appear to be the result of an external organizing force may have additional dynamics that shape a n d maintain the group structure. For example, although chemicals released by the queen provide a stimulus for honey bee swarm formation, the workers release additional signals, providing a feedback mechanism that helps to maintain the swarm's structure, indeed it is the worker-released chemicals that initially attract the queen (Winston, 1987, pp. 133,142). N o such feedback is present i n the midge swarm. The elements of the group structure of a fish school are its size, shape, and velocity. (In this case, the velocity is the rate of change of the average position of all fish i n the school.) These elements are determined by the motions of the individual fish. Individuals are attracted to their neighbours if they are separated by medium or long distances, but repulsed from their neighbours if they are too close. Several authors ( H u t h and Wissel, 1992; N i w a , 1994; Reuter and Breckling, 1994) have shown that changes i n these individual behaviours result i n changes i n the elements of the group structure mentioned above. Further, these individual behaviours are sufficient for the formation of the group structure. N o external influences are necessary. Even i n cases where there appears to be an external organizing force responsible for the group structure, it may not be the determining factor i n the pattern. For example, a flock of sheep can be directed by a shepherd and a sheepdog. However, it is the change in behaviour of the sheep i n the presence of the sheepdog that leads to the formation of the flock. T h e sheep respond to the dog by moving closer to their neighbours, rather than scattering. T h e dog does not control the motion of the individual sheep, but rather influences their behaviour and their response to their neighbours. In contrast, sudden, Chapter 2. Background 6 external movements, such as a waving arm, can cause a midge swarm to disperse. T h e individual midges respond directly to the external event. Responses to their neighbours serve only to avoid collisions, not to maintain the group's structure. G r o u p dynamics can be classified into two broad categories. Individuals may interact directly, by exchanging information w i t h each other, or indirectly, by exchanging information w i t h their shared environment. A direct interaction can be defined as any action by one individual that elicits an immediate response from, or change i n the behaviour of another individual (i.e., a communication). For example, schooling fish respond to the motions of their neighbours. A n example of indirect interaction is given by the chemotaxis model of Keller and Segel (1971b). T h e bacteria interact with each other by changing the local food density. Since the behaviour of the bacteria depends on the local food density, this allows for interactions and feedback between individuals (communication). Similarly, the trail networks of foraging ants are created by individuals who respond to the local concentrations of chemical markers by changing their direction of motion, and by reinforcing the markers. The two models presented i n this thesis contrast these two forms of group d y n a m ics. In the thermoregulation problem, the behaviour of the bees depends on the local temperature. Since the thermal energy diffuses though the cluster, information about the core temperature is 'hidden' i n the thermal flux and temperature i n the mantle. B y basing their behaviour on the local temperature and flux, the mantle bees are indirectly responding to the conditions i n the core. In the pheromone transmission problem, the bees interact by a direct exchange of a (chemical) pheromone. However, since the bees are moving randomly through the colony, the pheromone (and therefore information about the queen) diffuses through the colony. Chapter 2. Background 2.2 A brief survey of self-organizing 7 systems The formation of groups often occurs i n systems containing many similar individuals, none of which appears to take the role of leader. In some cases it is clear that the driving forces for group formation are external. However, i n many cases, there is neither an external nor a central force driving the formation of the group. T h e group structure arises from the group dynamics. In other words, the structure is a result of the interactions between i n d i v i d u a l members of the group. In cases where the group structure is a result of local interactions between individuals and requires neither a centrally located nor external regulating mechanism, it is classified as S e l f - O r g a n i z e d . Equivalently, we may say that group structure emerges from self-organizing group dynamics. W i t h respect to organism grouping, the purpose of the term "self-organizing" is to distinguish between groups that arise as a result of an external attractant, such as a mate or a food source, and groups, such as fish schools, that are formed and maintained without these external organizing influences. For example, the motion of a swarming midge depends explicitly on its position relative to the attracting marker. T h e movement is also influenced by the positions of their nearest neighbours. However, if the marker is removed, the swarm will disperse. In contrast, the motion of a schooling fish depends only on its position relative to its nearest neighbours. Assuming the fish do not orient relative to the geometric center of the school, the group is self-organized. The term has also been applied i n developmental biology to distinguish between the extremes of development as the execution of a 'hard-wired' genetic code, and development as a series or set of complex interactions between the growing cells. Here, the purpose of the term is to emphasis the importance of the interactions between the components of the system, and to clarify the assumption that these interactions do not necessarily depend on the global aspects of the developing organism. The genetic code does not completely Chapter 2. Background 8 determine the development of the organism, but rather specifies and regulates these interactions. The first serious examination of a self-organizing system was presented i n a series of papers by Ilya Prigogine and his co-workers (Prigogine and Nicolis, 1967; Prigogine and Lefever, 1968; Lefever et al., 1967). These papers examined an oscillating chemical reaction discovered by Belousov and the model for morphogenesis proposed by T u r i n g (1952) 1 and compared them to the hydrodynamic instabilities giving rise to Benard convection (Chandrasekhar, 1961). These were all classified as dissipative structures and thought to have great potential for the study of biological structures which "originate i n a dissipative m e d i u m " and are "maintained by a continuous supply of energy" (Prigogine and Nicolis, 1967). T h e driving force for the pattern is buried partially i n the kinetics of the reaction and partially i n the diffusion of the reactants through the system. There is no external force creating the pattern, although energy is required to maintain i t . These are purely physical patterns where the group dynamics are replaced by reaction kinetics, and the individuals by molecules. Segel and Jackson (1972) speculate on the relevance of these dissipative structures to pattern formation i n predator-prey systems. (Yates (1987) has collected many more examples of self-organized physical systems.) The ideas of Prigogine have been applied to the organizational capabilities of social ants. M a n y species were shown to display self-organized trail formation (Pasteels et a l . , 1987; A r o n et al., 1989a; Beckers et al., 1990) or brood sorting (Deneubourg et a l . , 1991). Simulations have been used to show that patterns of brood, pollen, and nectar placement i n honey bee hives may arise from the behaviour of the individual bees (Camazine, 1991). It is also possible that many honey bee colony activities are regulated by trophallaxis 2 (Winston, 1987, pp 195-196). If forager behaviour is based on the brood nest worker Excellent reviews of oscillating reactions are given by Degn (1972) and Tyson (1976), and several demonstrations are described by Field (1972) and Lefelhocz (1972) The exchange of food between nestmates. 1 2 Chapter 2. Background 9 demand for pollen, nectar, and water and on the availability of food from the brood nest workers then there is a feedback between foragers and brood nest worker behaviour allowing foragers to respond to changes i n the conditions and requirements of the colony. More recently, the premise of self-organization has been applied i n several models of group movement. N i w a et al. (1994) have postulated equations (based on the conservation of momentum for the schooling of fish) which can be expressed as evolution equations for the position and momentum (velocity) of each individual fish. T h e changes i n the momentum are due to fluid forces and 'social' forces. These social forces are not based strictly on the behaviour of the individuals, but are designed to align the motion of each individual w i t h the mean motion of the group. Unfortunately, this makes it difficult, if not impossible, to estimate the form and magnitude of the forces from observations. Gueron and L e v i n (1993) have developed a model for herd movement which also is based on hypothetical 'social' forces. In their model, the social forces are functions of the relative positions of the nearest neighbours, and the group dynamics are based on local information only, not on global properties of the group structure. Hence, their model supports the hypothesis that herding is self-organized. However, their implementation of the social forces i n modelling the interactions brings the focus of the model away from the individual and again makes the determination of the parameters difficult. 2.3 2.3.1 Modelling group dynamics Lagrangian and Eulerian models In this thesis, I use models based on a Lagrangian description as well as models based on an Eulerian description. T h e Lagrangian description focuses on changes i n the state of each individual, whereas the Eulerian description focuses on changes i n the number of individuals at any given position. (See, for example, O k u b o , 1980, p. 3.) In terms Chapter 2. Background 10 of fluid mechanics, where the concepts originated, a Lagrangian description takes the perspective of a particle moving with the fluid, whereas an Eulerian description takes the perspective of the fluid motion relative to a rigid frame of reference. T h i s frame need not be fixed i n space. For example, a translating or rotating coordinate system is appropriate for studying travelling or rotating waves, an approach that has been used widely i n the mathematical biology literature i n the form of reaction diffusion equations. The difference between the Lagrangian and the Eulerian descriptions is nicely portrayed using the example of elastic materials. A n Eulerian description fixes the coordinate system relative to the observer, whereas a Lagrangian description uses a coordinate system that stretches w i t h the deformations of the material. A l t h o u g h originally intended to describe continuous fluids, Lagrangian descriptions are natural for studying group dynamics. Since we are interested i n basing the models on the behaviour of the individuals i n the group, it is natural to formulate equations describing the dynamics of the individuals. Eulerian models, that is, systems of equations written using an Eulerian description, are derived using the conservation of global descriptors, such as the mass of individuals. It is important to note, however, that when performing a mass balance it is necessary to make hypotheses about the dynamics of the individuals. These hypotheses may be difficult to present using an Eulerian description, and assumptions about the behaviour may not be clear. A l t h o u g h the two descriptions allow for different approaches to modelling group dynamics, they are related, and Eulerian models can be derived from Lagrangian models (Section 2.3.2). In another approach, G r i i n b a u m (1994) discusses how the parameters used i n Lagrangian models, which describe the behaviour of the individuals, can be used to develop Eulerian models, whose parameters describe the collective properties of the group. Since Lagrangian models often consist of large dynamical systems representing the evolution of the state of each individual, they present many difficulties for analysis and are Chapter 2. Background 11 usually studied using cellular automata ( C A ) simulations. C A simulations are numerical methods that use a coarse discretization of both the independent and dependent variables (Ermentrout and Edelstein-Keshet, 1993; Wolfram, 1984). Since Lagrangian models focus on the behaviour of individuals, C A simulations of these models follow a finite number of organisms as they move over a fixed lattice of spatial sites. The state variable representing the surrounding medium is discretized over' the same lattice. Eulerian models may be cast i n several forms. P a r t i a l differential equations can be used for predator-prey interactions or reaction-diffusion equations (Murray, 1989), where the spatial motions are small over the timescale of the interactions. P a r t i a l integro- differential equations can be used if there is pattern formation i n both space and aspect, and if the communication is non-local i n either space or aspect. One example of a physical process forming groups by such a mechanism is the aggregation of actin filaments (Civelekoglu and Edelstein-Keshet, 1994). T h e filaments diffuse i n both position and angle, and can interact and bind w i t h other filaments. T h e binding can take place w i t h filaments at any angle, w i t h a probability that depends on the relative angle. Thus, the binding can be described using an integral kernel. Finally, integro-difference equations can be used for systems such as plant dispersal (Neubert et a l . , 1995) where the interactions occur at discrete points i n time (generations) but globally i n space. Such systems are best modelled by a difference equation where the population of plants at the next generation depends on an integral function of the population i n the present generation and a kernel describing the global dispersion. In the following sections I present a method for developing Eulerian models. T h e models follow the examples of Keller and Segel (1971a; 1971b) and G u r n e y and Nisbet (1975) and consist of a partial differential equation describing the evolution of the density of individuals i n the state space. In each of these presentations the Eulerian models are developed by first proposing a Lagrangian model for the motion of the individuals, and Chapter 2. Background 12 then using the infinitesimal mean and variance of the individual motion i n developing an equation for the evolution of the density of individuals at each point i n space. I illustrate this technique using a simple random walk i n one dimension as the L a grangian model and show that this leads to the diffusion equation for an Eulerian model. I then make two simple extensions to this basic model. T h e first extension looks at systems i n which the individuals respond to an environmental state, such as a temperature or a chemical concentration. T h e second extension examines systems i n which the organisms have an additional state (such as a chemical concentration) that changes during contacts w i t h other individuals. T h e Eulerian models assume that the dynamics of the individuals depend on the distribution of the group rather than on exact individual positions. These 'social' forces must be carefully derived to allow reasonable parameter estimates to be made. 2.3.2 A simple random walk The simple random walk forms a good first approximation for a Lagrangian model, since its parameters can be easily estimated from observations of individual movements. Consider an organism that moves with steps of length 8 at the end of time intervals of length r , such that steps to the left or to the right occur w i t h equal probability. T h e position of the organism at time t is a random variable satisfying X(t) = X(t-r)±8. (2.1) For my purposes, the infinitesimal variance can be defined as the rate of increase of the mean squared displacement, and the infinitesimal mean as the drift, or rate of increase of the mean position. A s shown below, the mean squared displacement of the organism increases linearly with number of steps taken, and therefore w i t h time, and the mean position remains constant. T h e derivation follows that of Berg (1983). Chapter 2. Background 13 I use the notation E [ X ( t ) ] to denote the e x p e c t e d v a l u e of the random variable X at time t ( K a r l i n and Taylor, 1975). For this simple random walk E[X(t)]=E[X(t-r)}, (2.2) since the steps occur equally to the left and to the right. Thus, the infinitesimal mean is zero. T h e square of the position is X (t) = X {t -T)± 2 2 28X(t -r) + 8. 2 Since steps to the left and right occur with equal probability, the expected squared displacement is E[X (t)]=E[X (t-T)} 2 + 8. 2 2 This implies that the square of the displacement w i l l , on average, increase by 8 each step. 2 Since there are t/r steps i n a time interval of length t, the expected squared displacement of an organism starting at the origin is E[X (t)] = 2 ^t. Thus, the infinitesimal variance is a = - . r (2.3) 2 The central limit theorem ( K a r l i n and Taylor, 1975) implies that if 8 ,r —> 0 so that 2 a 2 = 8 /T remains constant, then the simple random walk converges to the d i f f u s i v e 2 p r o c e s s Y(t), where dY(t) is a normally distributed random variable w i t h mean zero and variance a . 2 T h e Fokker-Planck equation ( K a r l i n and Taylor, 1981) describes the evolution of the probability density function u(x,t) for the process Y(t) as (2.4) Chapter 2. Background 14 T h i s is similar to observing the system on a time scale much larger than the step time, r, and on a spatial scale much larger than the individual step size, 8. However, the diffusive process, Y(t), allows for arbitrarily large displacements i n any time interval and an apparent 'instantaneous' movement to infinity. This arises because, even though the expected displacement zero, the velocity of the individuals, <5/r, becomes large i n the limit as <5 ,r —>• 0. This implies that observations must also be restricted to a spatial 2 scale small enough to be quickly traversed by the individuals. In general, if an organism is performing a simple random walk w i t h an infinitesimal mean of p(x,t) and an infinitesimal variance of a (x,t), 2 then the probability density function for the position of the organism at time t is given by the Fokker-Planck equation (2.5) 2.3.3 D e n s i t y dependent diffusion a n d chemotaxis If many non-interacting organisms are released from the origin at an i n i t i a l time t = 0, then the distribution of the organisms will also evolve according to Equation (2.5). If the organisms are interacting, then this model no longer applies. However, Eulerian models can be developed which have the same form as the Fokker-Planck equation, and where the motility of the individuals is assumed to depend on the local population density (Gurney and Nisbet, 1975; G u r t i n and M a c C a m y , 1977). If the individual movement depends only on the positions of the nearest neighbours, we can make the approximation that the movement depends on the density of individuals i n a neighbourhood around the individual. For the case of fish schooling, where individuals align w i t h their nearest neighbours, this assumption appears reasonable (Turchin, 1980). It also holds if the behaviour of the individual is influenced by accumulated or averaged observations of neighbouring individuals over a time interval long enough so that it(x, t) Chapter 2. Background 15 is a reasonable approximation for the behaviour of a single neighbour averaged the time interval. T h e approach taken by (Okubo, 1980, p. 11,195), is to use Fick's diffusion model, which is a statement of conservation of mass, but assume a non-linear flux. Fick's diffusion equation for the density u(x, t) is du ¥ = -V-J. (2.6) where J is the flux of organisms past the point x at time t. For the case of pure (density dependent) diffusion with an infinitesimal variance that is independent of position, the flux can be written as J = -D{u)Vu. (2.7) T h e function D(u) is known as the motility of the organism. For organisms that enact a random walk biased down the population density gradient, Gurney and Nisbet (1975) arrive at the flux J = - — (1 + 4/U-u)Vii, (2.8) IT where p measures the magnitude of the bias. W h e n the drift is non-zero, this may be added as a second flux proportional to the density IA(X, t). For example, a random walk which is biased by a chemical concentration (Keller and Segel, 1971a; Sherratt, 1994) can be represented by the equation du = - V • (x(c)uVc - D(u)VuJ . ~dt (2.9) This is referred to as the chemotaxis equation. The variable c(x, t) is the concentration of a chemical i n or on the substrate. T h e coefficient of chemotaxis, X, may, i n general, depend on c. T h e parameters i n these equations describe the bulk properties of the motion of the organisms through state space. Chapter 2. Background 16 Unless significant simplifications can be made to the system before developing a simulation of these Eulerian models, it is often better to simulate the system using a L a grangian model. This Lagrangian model incorporates the full interacting particle system. T h a t is, the simulation should model the dependence of the individual behaviour on the configuration of individuals rather than the density as described by it(x, t). T h i s approach also leads to a discretization of the model that has a scale appropriate to the dynamics of the individuals. T h e cost of the time required to simulate a differential equation for each individual is balanced by the ability to model the interactions i n more detail, w i t h a focus on the dynamics of the individuals rather than on parameters describing their bulk properties. However, many simplifications can be made to Eulerian models which allow quick simulation and analysis of abstractions of the full system. One can only conclude that b o t h Eulerian and Lagrangian descriptions of the system should be examined. 2.3.4 R a n d o m walks with pheromone exchange A s another extension of the basic random walk, I associate a (chemical) pheromone level w i t h each individual. A s the individuals move randomly, they interact by exchanging the pheromone. If p is the pheromone level of individual i then this might be represented, l by the system of differential equations -jr = -IP at 1 + £K(p\p>)w(X* - X*). (2.10) j = 1 (Superscripts refer to the individual organisms.) T h e chemical level of each individual decays exponentially, and the chemical is exchanged with other individuals according to the kernels K and w. T h e first kernel is the rate of pheromone transmission from individual j to individual i, and the second kernel is the neighbourhood over which the encounters take place. Chapter 2. Background 17 The equations are coupled, since the exchange occurs between each pair of individuals. However, the equations can be decoupled by replacing the exchange rate w i t h the expected exchange rate given that the organisms have the density u(x,p,t). dp ~^ = ^P 1 l + J N 0 f°° f J K(p\p)w{x -x)u(x,p,t)dxdp, (2.11) t n where TV is the number of organisms. T h a t is One possible simplification is that w(x) is a delta function. T h a t is, the neighbourhood of the interactions is small compared with the observed scale, as would be the case for interactions requiring contact. W i t h this simplification, the previous relation can be written as l- = - d 7 ^ + TV | ° ° K(p\p)u(X\p,t)dp. (2.12) o If the organisms move randomly over Q, with a zero drift and a variance of cr , then 2 the Eulerian model is given by the Fokker-Planck equation for the p.d.f. u(x,p,t) with an added drift through pheromone levels. T h a t is, du <7 d / fc=Y dp~( ~ Jo 2 Au+ JP NU f°° ^^H^Q^dq). \ (2.13) This is a (non-linear) partial integro-differential equation describing the evolution of the density of the organisms over both (physical) space and pheromone level. T h i s equation can also be derived from a Boltzmann-like equation, and has been applied to the situations where individuals i n a colony rank themselves according to the outcomes of a series of bouts (Jager and Segel, 1992). The outcome of each bout changes an aspect of the participants referred to as their dominance. Thus, i n effect, the individuals exchange dominance during chance encounters. These extensions are used i n the following chapters to develop models for the cluster thermoregulation problem and the pheromone transmission problem. Chapter 3 The 3.1 Thermoregulation Model Introduction Honey bees will aggregate to' form clusters at cool temperatures. T h e y are able to regulate the temperature of the core of the cluster to within a few degrees of 35°C over a surprisingly wide range of ambient temperatures. A s the ambient temperature drops, the bees on the periphery of the cluster will form a dense mantle as shown i n Figure 3.1. T h i s reduces heat loss by reducing the surface area of the cluster and by l i m i t i n g air flow through the cluster. Surprisingly, thermoregulation is accomplished w i t h apparently no communication between the bees at the centre of the cluster and those i n the mantle. In this chapter I develop a model for honey bee cluster thermoregulation. T h i s model is based on the hypothesis that the behaviour of the individual bees is dependent on their local temperature, and that the thermoregulation of the cluster is a direct result of this behaviour. This chapter is based on a joint project recently published by W a t m o u g h and C a mazine (1995). The authors proposed a model for the self-organized thermoregulation of honey bee clusters that is based on the mechanics and energetics of the individual bees. Scott Camazine is a behavioural ecologist at Pennsylvania State University. He suggested the problem, the appropriate, biologically accurate assumptions, and the questions that the model should address. I was responsible for the derivation, analysis, and simulation of the model, and the interpretation of the results. 18 Chapter 3. The Thermoregulation Model 19 Core: • lower densities Figure 3.1: Schematic of a thermoregulating cluster showing the features that should be predicted by the model. Chapter 3. The Thermoregulation Model 20 T h e model addresses the following questions: 1. C a n the temperature and density profiles of a thermoregulating cluster be modelled on the basis of the experimentally observed individual behaviour? 2. C a n such a model account for observed responses of the cluster to changes i n the ambient temperature? T h e challenge i n addressing these points is to develop a model that uses o n l y the individual behaviour and local interactions of the bees as inputs, and does not use the global properties of the cluster as constraints or influences on the individual behaviour. This challenge has not been met by previous models of thermoregulation. T h e model consists of two partial differential equations for the temperature and density profiles and an integral constraint for the radius. Thus, solutions of the model specify the temperature and density profiles and the radius of the cluster as functions of ambient temperature and time. T h e assumptions of the model are all based on the behaviour and interactions of the individual bees. In Section 3.2,1 introduce the background necessary for the model and a brief review of previous models. In Section 3.3, I begin with a statement of the assumptions made about the individual behaviour of the bees, and then use these assumptions to develop the mathematical model. The full model is studied using a numerical simulation (Section 3.4). The steady state solutions are analysed i n Section 3.5. 3.2 3.2.1 Background S w a r m a n d winter cluster thermoregulation Precise thermoregulation by aggregations of honey bees occurs i n two distinct situations. T h e first is during r e p r o d u c t i v e s w a r m i n g . Swarms ranging from 2 000 to 20 000 bees Chapter 3. The Thermoregulation Model 21 may remain exposed on a tree branch for several days (Heinrich, 1981a; Heinrich, 1981b; Heinrich, 1985). A t ambient temperatures below 15°C, the swarm regulates its core temperatures w i t h i n a few degrees of 35°C and its mantle temperature near 15°C. For warmer ambient temperatures, the core temperature remains regulated near 35°C, whereas, the mantle temperature w i l l stay a few degrees above the ambient temperature. T h e second situation involves a w i n t e r c l u s t e r of 10 000 to 30 000 bees (Seeley, 1978; Seeley, 1983). Over an ambient temperature range of -15°C to 10°C, the core temperature of the cluster is maintained between 18°C and 32°C and the mantle temperature is maintained between 9°C a n d 15°C (Heinrich, 1985; Simpson, 1961; Southwick and Mugaas, 1971). In contrast to homeothermic animals, honey bee colonies do not possess a centralized physiological mechanism capable of monitoring and adjusting temperatures w i t h i n the cluster (Southwick, 1983). In fact, it has been demonstrated that the thermoregulation of the cluster does not rely on any form of acoustic communication or communication using volatile chemicals, or on any contact between bees on the mantle and those i n the core (Heinrich, 1981a). This evidence supports the assumption that thermoregulation is a result of the local collective actions of many individual bees. T h a t is, it is a result of each bee responding to the local temperature according to a simple behavioural algorithm. We may assume that the bees on the mantle are not aware of the conditions of the bees in the core, nor are the bees i n the core producing excess heat i n order to warm the bees on the mantle. T h e regulation of the core temperature results from each bee attempting to remain at a comfortable temperature within the confines of the behaviours permitted it by its own individual physiology. 3.2.2 M i n i m a l requirements for a satisfactory model Based on his observations, Camazine suggested the following criteria for a suitable model: Chapter 3. The Thermoregulation Model 22 • T h e model should predict the global properties of the cluster such as its radius, the core and mantle temperatures, and the core and mantle densities. In particular, the model should predict the formation of the denser mantle surrounding a core of warmer bees. • In keeping w i t h the premise that cluster thermoregulation is self-organized, the behaviour of the individual bees should be the only input to the model. • T h e model should not include a specific core or mantle temperature s e t p o i n t to be maintained. It is clear that the bees do not thermoregulate to keep the core temperature at a preset value. In fact, at cooler ambient temperatures, the core temperatures of large clusters can rise to levels high enough to k i l l the bees (Heinrich, 1981a). This is clearly not the result of regulation to a specific core temperature. • T h e behaviour of the bees should depend solely on local conditions (specifically the local temperature and density), and not on the bee's position w i t h i n the cluster. Most importantly, bees that are not i n the core should have no information about the core temperature, and behaviour should not depend on the ambient temperature. A l t h o u g h it is possible that behaviour depends on position, it is more likely that the behaviour depends on the local conditions, and that any apparent positional dependence of behaviour is a result of the positional variations i n temperature and density. • T h e model should predict a core temperature that increases w i t h decreasing a m bient temperatures, as observed i n several experiments (Fahrenholz et a l . , 1989; Heinrich, 1981a; Heinrich, 1981b; Southwick and Mugaas, 1971). Chapter 3. 3.2.3 The Thermoregulation Model Previous 23 models The models of Omholt and L0nvik (1986), Omholt (1987), and Lemke and Lamprecht (1990) show that the temperature profiles of the clusters are consistent with those of a heat producing sphere. They assume that heat loss is due to convection from the surface of the sphere and that heat transfer through the sphere is due to conduction only. Unfortunately, these models all begin with several assumptions about the global behaviour of the cluster that should, ideally, be predictions of the model. Specifically, the cluster radius is specified a priori, rather than requiring it to be a model prediction that can be tested against experimental data. Further, each of these models assumes that the metabolic heat production of a bee is a function of its position within the cluster, and that the behaviour of the individuals depends on both their position in the cluster and the ambient temperature. These assumptions require that individuals have knowledge of the global properties of the cluster and are therefore difficult to verify experimentally. For this reason, a model that does not assume a positional dependence on behaviour would be much more desirable. In short, both the global behaviour of the cluster and the local behaviour of the individuals in the cluster are used as inputs to a model intended to connect individual behaviour to group structure. A much improved model has been developed recently by Myerscough (1993). Her model uses the known metabolic rate of the individual bees to predict the global temperature profile within the cluster. Her model does predict cluster thermoregulation; however, it is unable to predict the observed increase in core temperature with decreasing ambient temperature. A n adjustment to her model (Myerscough, 1994), which includes a set point for the mantle temperature, does predict this cluster response. In her models, Myerscough assumes that the density at any point within the cluster can be determined a priori from the local temperature. Although this is plausible, it Chapter 3. The Thermoregulation Model 24 is not clear that bees i n the core of the cluster would be able to achieve this optimal density when the bees on the periphery of the cluster are forcing a closer packing by their inward motions. Further, this assumption specifies a target density for the bees rather than a mechanism for achieving that density. In the terminology suggested i n the previous chapter, it specifies the group structure and not the group dynamics. T h e model presented here draws heavily on Myerscough's model, but extends it by removing the specification of the dependence of the density on the temperature. I use the observed motions of the bees to predict the global density profile. 3.3 The Model The model presented i n this thesis is an attempt to identify the mechanisms responsible for thermoregulation. In contrast to previous models, this model incorporates the motions as well as the metabolism of the bees. In addition, I base the model solely on the responses of individual bees to their local environment (the group dynamics). T h e observed global features and responses of the cluster to changes i n the environment (the group structure) are not used as inputs, but are rather derived as outputs of the model. 3.3.1 Assumptions of the model The model is based on the following assumptions: A l : J E a c h bee bases her behaviour exclusively on the temperature she experiences i n her local environment. A 2 : | If a bee is too warm, she will move i n the direction of decreasing temperature. If she is too cold, she will move i n the direction of increasing temperature. A 3 : g T h e heat production of an individual bee is based on her metabolic rate. Below a threshold temperature, a bee increases her metabolism above the resting level by Chapter 3. The Thermoregulation Model 25 shivering with her flight muscles. Above this threshold temperature, she remains at her resting metabolism. The resting metabolism is an increasing function of temperature as given by Kammer and Heinrich (1974). | A4: j Heat is transferred through the cluster predominantly by conduction rather than by convection. | A5: | The cluster is spherically symmetric. Assumption A l is based on the experiments of Heinrich (1981a), which failed to show that the thermoregulation is coordinated by any form of chemical, physical, or acoustic communication among the bees. In contrast, each bee appears to rely on information taken from her immediate surroundings. The key mechanisms for thermoregulation are movement (Assumption A2) and heat production (Assumption A3). The assumed movements are based on the observation that the bees on the mantle will huddle closer together as the ambient temperature falls (Heinrich, 1981a; Ribbands, 1953; Southwick, 1983). Unfortunately, these observations were qualitative rather than quantitative. In contrast, the heat production data (Assumption A3) are based on several quantitative experiments. Assumption A5 is, in fact, a simplification. According to the observations of Heinrich (1981a), swarms are nearly spherical at lower ambient temperatures. A t ambient temperatures above 15°C, they will elongate and become more parabolic in shape. The top of the parabola will be deformed to the shape of the tree branch or other protrusion from which the swarm hangs. The simplification holds over wider range for winter clusters. Assumption A4 may be valid at lower ambient temperatures. Below ambient temperatures of 15°C, the bees on the mantle are observed to huddle close together, restricting the flow of air through the cluster. However, at higher ambient temperatures, visible Chapter 3. The Thermoregulation Model 26 channels through the cluster open and serve to rapidly convect heat from its centre (Heinrich, 1981b; Heinrich, 1981a). Thus, at ambient temperatures above 15°C, convection w i l l become more important than conduction as a heat transfer mechanism, and A s s u m p t i o n A 4 w i l l no longer be valid. 3.3.2 D e r i v a t i o n of the Model The model presented here is Eulerian. I develop equations for the evolution of the bee density using conservation of mass (Equation (2.6)). The flux is based on small random deviations from the motion along the temperature gradient (Assumption A 2 ) . The temperature profile evolves according to the classical heat equation w i t h a source and an added flux due to the movement of the bees. According to the criteria suggested by D r . Camazine (Section 3.2), the model should predict the cluster radius, and the temperature and density distributions w i t h i n the cluster. The cluster is three dimensional; however, w i t h the assumption of spherical symmetry (Assumption A 5 ) , positions i n the cluster can be defined by the radial distance from the centre of the cluster. Let r be the distance from the centre of the cluster, and let R(t) be the radius of the cluster at time t. Thus, 0 < r < R(t). Let T(r, t) and p(r, t) be the temperature and density of bees at position r and time t. The cluster radius and the density profile will be constrained by the condition (3.1) where N is the total number of bees i n the cluster. In the following paragraphs I derive the equations describing the evolution of the density and the temperature profiles of the cluster. The B e e D e n s i t y P r o f i l e will change due to the movements of the bees. As- sumption A 2 implies that the motion of the bees can be divided into two components: Chapter 3. The Thermoregulation Model 27 movements that are aligned w i t h the temperature gradient (thermotactic) and movements that are not (diffusive). I assume that the motion of the bees is a random walk w i t h a drift velocity proportional to the temperature gradient. B y A s s u m p t i o n A 2 , the drift velocity w i l l be i n the direction of increasing temperature if the bees are too cold, and reversed if they are too hot. This can be modelled by a random walk w i t h a drift dT velocity of X(T) — . T h e t h e r m o t a c t i c c o e f f i c i e n t , X(T), must be positive if the local temperature is below the h u d d l i n g t e m p e r a t u r e Th, and negative if the local temperature is above this level. This indicates that the bees move up the temperature gradient if they are too cold, and down the temperature gradient if they are too warm. The parameter Th can be measured experimentally, thus providing support for Assumption A 2 . However, these experiments have not yet been done. I assume X{T) to be a hyperbolic tangent. Thus, the behaviour can be characterized by a preferred huddling temperature Th, a l i m i t i n g t h e r m o t a c t i c coefficient X , 0 and a r e s p o n s e r a n g e T. r The thermotactic coefficient is X(T) = X 0 t a n h ( ^ ^ ) . (3.2) The function is plotted i n Figure 3.2. If the variance of the random walk is modelled using the motility D(p), then the flux of bees past the point r at time t is J { , T) = ~D{p)^ P P r + X(T)p^, (3.3) and the density profile evolves as B y introducing a dependence of the motility on the local bee density, it is possible to prevent overcrowding within the cluster. I assume that the motility is higher i n regions Chapter 3. The Thermoregulation Model 10 15 20 25 30 Temperature (° C) 28 35 40 Figure 3.2: For simplicity, the thermotactic coefficient is assumed to be the hyperbolic tangent X(T) = X t a n h ( ( T - T)/T ). This function is plotted above using X = 1, Th = 25 and T = 4, to show the shape of the curve. T h e essential properties of the thermotactic coefficient are that it be decreasing, and that it pass through zero at the preferred huddle temperature Th0 r h r 0 Chapter 3. The Thermoregulation Model 29 of very high density. This results i n bees spending less time i n regions where the density is too high, and lowers the probability of overcrowding. Heinrich (1981b) notes that the average density w i t h i n a cluster is roughly 2 bees/cm at higher ambient temperatures. A s 3 the temperature falls to 0°C the a v e r a g e density increases to a m a x i m u m of 8 b e e s / c m . 3 To obtain average densities i n this range, I assume an increase i n the m o t i l i t y if the density is above the m a x i m u m d e n s i t y p max = 10 bees/cm . For example, 3 max) (3.5) A l t h o u g h this is plausible, there is no experimental evidence that it is biologically accurate. It is done strictly as a convenient mechanism for preventing overcrowding i n the model. A l t h o u g h the parameters of the functions X(T) and D(p) have not yet been measured experimentally, realistic estimates can be made as follows. I assume that the variance arises from steps of less than 1 m m i n length at 1 second intervals, and that the thermotactic velocities are less than 1 m m / s i n a temperature gradient of 1 0 ° C / c m . This can be modelled as a step process (Section 2.3.1). A s shown i n Section 2.3.1, E q u a t i o n (2.3), this leads to maximums of 10~ c m / s for the base motility D 2 2 0 and 1 0 - 2 c m / s / K for 2 the thermotactic coefficient X . a T h e T e m p e r a t u r e D i s t r i b u t i o n within the cluster is modelled using a classical heat equation w i t h a source term. T h e rate of change of thermal energy a distance r from the centre of the cluster is due to three sources: heat transfer by movement of the bees (Assumption A 2 ) , heat transfer by air conduction (Assumption A 4 ) , and heat production by the bees (Assumption A 3 ) . I assume that the heat capacity of the bees is due to their water content. T h e average mass of a honey bee is 115 m g (Heinrich, 1981a). Hence, the heat capacity of a single Chapter 3. The Thermoregulation Model 30 bee is approximately 0.481 J / K . A t the observed densities of two to ten b e e s / c m , this 3 is two orders of magnitude higher than the heat capacity of dry air. I therefore assume that the fraction of the thermal energy of the cluster contained i n the air is negligible, so that this energy is contained mostly i n the bees. Thus, the thermal energy density of the cluster is cpT where c is the heat capacity of a single bee. The rate of change of the thermal energy of the cluster can be expressed as J ^ ) T {r M ,T)) = 2 P +pf{T). (3.6) The left hand side is the rate of increase of the thermal energy density of the cluster. T h e thermal flux J T is the combined energy flux due to bee motion and conduction. T h e final term is the heat production per unit volume of the cluster, where f(T) is the m e t a b o l i c output per bee. The energy flux through a shell of radius r is the sum of flux due to bee movement and flux due to air conduction. Recall that there is no heat transfer by air convection through the cluster (Assumption A 4 ) . T h e motion of the bees given by E q u a t i o n (3.3) implies a thermal energy flux cTJ . Combining this with the energy flux due to the p t h e r m a l c o n d u c t i v i t y of the cluster, k(p), gives the total thermal flux J (p,T) T = cTJ (p,T) p - k(p)^. (3.7) The conductivity of the cluster will decrease as the density increases. Southwick (1985) suggests that the conductivity may drop as low as 0.3 m W / m / K i n the cluster mantle since the bees are able to trap air i n their body hair, forming an insulating layer. This value of k is comparable to the conductivity of dry air. Conductivities may be slightly higher at lower densities due to increased air movement. Recall that the m a x i m u m density of the bees is p max by bees is p/p axm = 10 bees/cm . 3 Thus, the local fraction of the cluster occupied For simplicity, I assume conductivity decreases linearly w i t h density. Chapter 3. The Thermoregulation Model 31 T h a t is, k(p) = k - (k - ki)(p/p ) 0 0 (3.8) max where ka is the conductivity of a loosely packed cluster (p « p x) and k\ is conductivity ma of a densely packed cluster. T h e function f(T) represents the production of heat both passively, by the normal metabolic output of the bees, and actively, by shivering the flight muscles. T h e rate of metabolic heat production of a resting bee increases exponentially w i t h temperature ( K a m m e r and Heinrich, 1974). Lemke and Lamprecht (1990) suggest that the resting metabolic output of a single bee increases by a factor of 2.4 for every 10°C increase i n the local temperature. This factor is i n agreement w i t h the data cited by Seeley and Heinrich (1981), who also suggest a resting metabolism of / 1 8 = 1 m W at 18°C. T h i s is expressed by the function Q(T) = f 2.4 - ° ^ ° \ T 18 c 10 (3.9) c 18 A s the local temperature drops below a comfortable level, the bees begin to shiver. T h e bees are able to produce large amounts of heat by dislocating their wings and v i b r a t i n g their flight muscles. I assume that the bees begin shivering when the temperature falls below the s h i v e r t e m p e r a t u r e T , and that their metabolic heat production increases s exponentially to a m a x i m u m of f m can be estimated at f m at a local temperature of T . The m a x i m u m o u t p u t m = 50 m W / b e e , based on several experiments ( K a m m e r and Heinrich, 1974; Seeley and Heinrich, 1981; Heinrich, 1985). T h e values of T and T s m have not been measured. Thus, the dependence of the metabolic output of a bee on the Chapter 3. The Thermoregulation Model 32 M a x i m a l Shivering ( T ) m Temperature (° C ) Figure 3.3: T h e graph shows / ( T ) , the rate of metabolic heat production of a single bee as a function of the local temperature. Above T = 16°C, heat production doubles approximately every 10°C, and is / i = 1 m W at 18°C. Below T the bees w i l l increase their heat production by vibrating their flight muscles. I assume heat production increases exponentially to f = 50 m W as the temperature drops from T to T . For this figure, T = 13. s s 8 m s m m temperature is described by the following function: T < T, fr, m T -T s /CO Q{T ) S Q{T) fr, M s), T - T T m <T < T, s (3.10) T T > T,. T h i s function is sketched i n Figure 3.3. T h e m a x i m u m f m and the function Q(T) are known from experiments, and are joined to form a continuous function given by E q u a tion (3.10). Chapter 3. The Thermoregulation Model 33 E q u a t i o n (3.6) can be simplified using Equation (3.4) to yield T h i s equation, together with Equation (3.4), which is repeated here for convenience, I-?IK^ <"» make up a system of P D E s for the evolution of the temperature and the density of the cluster. The B o u n d a r y C o n d i t i o n s on the temperature and bee density at the centre of the cluster are dp f)<n\ dr r=0 = 0 dT f)T\ —— and lr=0 = 0, t>0. (3.13) These conditions arise.from the assumed spherical symmetry of the cluster a n d smoothness at the origin. A c c o r d i n g to N e w t o n ' s L a w o f C o o l i n g , The rate of heat loss by convection from the surface of the cluster is proportional to the difference between the cluster's surface temperature a n d the a m b i e n t t e m p e r a t u r e T . T h e rate of heat loss w i l l depend on a the properties of the cluster surface as well as on the moisture content a n d speed of the surrounding air, and will vary over the surface of the cluster and over time. However, these variations should be rapid, and should not result i n increased heat loss from any specific position on the surface. Heat loss can be approximated using the average rate of convection from the entire cluster. Balancing the heat energy generated by the bees w i t h the heat lost from the cluster surface leads to the boundary condition dT\ = h (r(R(t),t)-T«y ' r=R(t) c OO, (3.14) at the edge of the cluster. T h e coefficient h represents the average rate of heat loss from c the cluster surface due to convection. For heat transfer by free convection of air, the c o e f f i c i e n t o f c o n v e c t i v e h e a t t r a n s f e r (h ) ranges from 6 to 30 W / m / K c 2 (Kreith Chapter 3. The Thermoregulation Model 34 and B o n n , 1986, table 1.4). T h e temperature profiles given by Heinrich (1981a, Figure 1) can be used i n conjunction w i t h Equation (3.14) to estimate the ratio ( T _ T \ k(p(R(t),t) dT h. r=R(t) 1 to 3 c m . (3.15) Thus, the above range for h (0.6 to 3 m W / c m / K ) suggests that the conductivity near 2 c the edge of the cluster, k(p), is i n the range of 0.6 to 9 m W / c m / K . These values are comparable to heat conduction coefficients of 0.3 m W / c m / K for d r y air, 6 m W / c m / K for d r y sand, 10 m W / c m / K for moist sand, and 6 m W / c m / K for water. In addition, I impose the constraint of constant density • p(R(t),t)=p , (3.16) R at the edge of the cluster. To my knowledge, the actual conditions at the boundary of the cluster are not well characterized. It is likely that additional mechanisms are responsible for the cohesion of the cluster, for example, the density p R may depend on either the ambient temperature, T , or on the cluster radius, R. However, these conditions are not a known at this time. This boundary condition is chosen for its simplicity. For the results presented i n this thesis I use the value p R = 2. Further experimental a n d theoretical research is needed to improve this condition. I n s u m m a r y , honey bee cluster thermoregulation can be understood using the following scenario. Consider a cluster of bees at equilibrium at an ambient temperature above the huddling temperature T V A s the ambient temperature falls, the rate of heat loss from the surface of the cluster (see Equation (3.14)) will increase, causing the temperature on the surface of the cluster to decrease. A s the surface temperature falls below T V the bees on the surface will begin to move closer to their neighbours, increasing the local density of the cluster, and decreasing the radius R. If the ambient temperature continues Chapter 3. The Thermoregulation Model 35 Parameter Page 0.00062 c m / s / K M a x i m u m thermotactic coefficient. 27 T 25°C Huddling temperature. 27 T 4°C Temperature range for thermotactic motion. 27 D 0.0046 c m / s M o t i l i t y at comfortable densities. 29 c 0.481 Heat capacity of the bees. 30 Symbol Value x 2 0 h T 0 2 J/K/bee k 7.2 m W / c m / K Conductivity of loosely packed cluster. 31 h 3.6 m W / c m / K Conductivity of densely packed cluster. 31 fis 1 mW/bee Metabolic rate of resting bee at 18°C. 31 fm 50 m W / b e e M a x i m u m metabolic rate of shivering bee. 32 T 15°C Temperature below which bees shiver. 32 T 0 s Nm 13°C Temperature below which a l l bees shiver. 32 2 000 - 30 000 bees Number of bees. 26 h 1.8 m W / c m / K Coeficient of heat conduction. 33 T {-30,-25,...,20}°C Ambient temperature. 33 PR Pmax 2 bees/cm Density on cluster perifery. 34 M a x i m u m local density. 29 x c 2 3 10 bees/cm 3 Table 3.1: Thermoregulation model parameters and their values. Recall that K (Kelvins), J (Joules) and W (Watts) are SI units for temperature, energy, and power respectively. Chapter 3. The Thermoregulation Model 36 to fall, so that the surface temperature falls below the shiver temperature T , the bees s on the cluster surface will begin to shiver, increasing their metabolic heat production. 3.4 Simulations In this section, I present the numerical simulations of the system derived i n the previous section. The system is first non-dimensionalized to determine the appropriate time scale of the problem. I then use Baumgarte's (1972) technique to replace the constraint by a differential equation which has the constraint as a stable solution. 3.4.1 Non-dimensionalization Using the time dependent transformation *(<•,«> = j £ j , (3-17) the moving boundary can be removed from Equations (3.4) through (3.11). This is equivalent to measuring the spatial position relative to the size of the cluster. N o w x = 0 corresponds to the centre of the cluster, and x = 1 to the cluster edge. Equa- tions (3.1), (3.4) and (3.11) can be non-dimensionalized by introducing the following rescaled variables: <• = R*(t*) = t^rt, u(x,t*) = ^R(t), p(x,t*) = K ^T(r,t), PR The radius is scaled by the ratio of conductivity to surface convection. T h i s ratio is a measure of the spatial scale of the energy transfer. The time scale is a combination of the spatial scale and the motility of the bees. The temperature scale is determined from the ratio of thermotactic to diffusive motion. Thus, it is a scale appropriate to the responses of the bees. Finally, the density is scaled by the boundary density. These scalings are Chapter 3. The Thermoregulation Model 37 chosen to eliminate parameters from the equation for the evolution of the dimensionless density, and from the boundary conditions. Equivalently, we are considering the problem at a scale relevant to the movements of the bees. For convenience, I now drop the asterisks and use t a n d R(t) for the dimensionless time and radius. I also introduce the following dimensionless functions: KP) = D{p) X(u) = /CO X{T) k(p) /l8 and the following dimensionless parameters: kXp f 0 0 R 18 h D 2 0 o cD p 0 a k a R P = ' hlN 4Tvklp R The parameter e is a measure of the rate of metabolic heat production of the bees relative to the rate of heat loss from the cluster. T h e parameter a is a measure of the time scale of the thermal fluctuations relative to the time scale of the movements of the bees. T h e parameter 8 is a measure of the cluster size with respect to the spatial and density scales. The dimensionless equations are (after dividing Equation (3.11) through by p) 1 8 dp dt a du dt Rx 1 2 dx d 2 pR x dx 1 Liip) dp 2 > 2 +R 2 \ p dx , .du\ du dx j dx ,. . 0 < x < 1, x dRdu R dt dx (3.19) w i t h the boundary conditions dp dx = x=0 du dx 0, (3.20) = x=0 0, R(t)(u - u )\ a x=1 , Chapter 3. The Thermoregulation Model where u = T X /D a a 0 0 is the dimensionless ambient temperature. T h e constraint is g(p,t) = R 5 3.4.2 38 Implementation Jo x p(x, t) dx - p = 0. (3.21) 2 of the constraint I use Baumgarte's technique (Baumgarte, 1972) to implement the constraint i n the numerical simulation of the above equations. This technique involves replacing the constraint g(p,t) = 0 w i t h the differential equation (3.22) where 7 is a positive constant. T h e constraint g(p,t) = 0 is an asymptotically stable solution to this equation. Hence, any numerical errors introduced during the computations will remain small. Differentiating the constraint with respect to time leads to the second constraint ^-g(p,t) dt dp = 3R ^ f x pdx +R f dt Jo Jo 2 1 2 3 x ^dx. dt 1 (3.23) 2 Substituting — using Equation (3.19) into the above equation obtains d_ dt g(p,t) = 3i? ^ J 2 + R x pdx 2 1d U p o n integrating the last integral by parts this simplifies to d . ^dR . „ / , . du* . .dp (3.25) x=l Using E q u a t i o n (3.21) and Equation (3.25), Equation (3.22) becomes ~ dR „ / . <9p x . . du x p(x, t) dx-p^j = 0. + 7 2 (3.26) x=l Finally, dividing through by R gives the rate of change of R(t) as 2 dR ~dl dp -R\ cTx- dx1 / M X{u)p . . du 7 x=l ( \ f Jo 1 ^( '*) a;2 a: d a ; ~ (3.27) Chapter 3. 3.4.3 The Thermoregulation Model 39 Results of the simulations E q u a t i o n (3.19) and Equation (3.19) w i t h boundary conditions given by E q u a t i o n (3.20) and the constraint of Equation (3.27) were discretized using the method of lines. Second order differencing was used to approximate the spatial derivatives and Simpson's method was used for the integrals. Since oscillations were expected due to the thermotactic motion, a fourth order R u n g e - K u t t a algorithm was used for the temporal derivative. Core a n d mantle temperatures, a n d cluster radius: Figure 3.4 shows the steady state core and mantle temperature obtained from the numerical simulations as functions of the ambient temperature for two cluster sizes. These are comparable (at least qualitatively) to the experimental data of Heinrich (1981a, Figure 2) on the core and mantle temperatures of swarms. T h e mantle temperatures of Heinrich's swarms remained near 15°C for ambient temperatures below 15°C, whereas my results continue to decrease w i t h decreasing ambient temperatures. However, the mantle temperatures of winter clusters (Fahrenholz et al., 1989, Figure 2) were observed to increase linearly as ambient temperature increased from —10°C to 0°C. T h e core temperatures shown i n the figure for the 8000 bee cluster also match quite well w i t h the data of Fahrenholz (1989) i n the range - 1 0 < T < 0°C, and the data of Southwick and Mugaas (1988, a Figure 4) i n the range —15 < T < 5°C. T h e results show that the core temperatures do a not increase indefinitely as ambient temperatures decrease. T h e results also suggest that bees may be able to survive extremely cold temperatures as long as energy is available to the mantle bees for shivering. Chapter 3. The Thermoregulation Model 40 Ambient Temperature °C Figure 3.4: Results of numerical computations of Equations (3.19), (3.19), and (3.22) w i t h the boundary conditions given i n Equation (3.20). Note that the core temperature shows a slight negative correlation with the ambient temperature; whereas the mantle shows a positive correlation. B o t h the core and the mantle temperatures remain nearly constant for temperatures below 5°C. There is also a definite increase i n core temperature w i t h increased cluster size. The parameter values used are shown i n Table 3.1. T h e number of grid points was 20. Chapter 3. The Thermoregulation Model 41 These solutions were obtained by running the simulation at a constant ambient temperature until the solution appeared to converge. Convergence was assumed if the solutions d i d not change for several hours . 1 Convergence occurred over a large portion of the region of parameter space relevant to the problem. However, 'stable' oscillations and diverging solutions were also observed. T h e diverging solutions occurred when core temperatures rose faster than the movement of the bees could compensate. These solutions are not considered relevant to the problem, since the bees' responses to such dire situations have not been included i n the model. Similarly, although oscillating solutions have been observed (see discussion) these are thought to be linked to oscillations i n CO2 levels i n the cluster, and should therefore be examined using a more detailed model. Oscillations are a common feature of regulating systems, and the methods used by the bees to control or damp these oscillations deserve further study. Figure 3.5 shows the cluster radius for the same simulations. T h e larger cluster shows very little change in the cluster radius for ambient temperatures below 5°C. T h i s may be an indication of the conflict of interest between the bees i n the mantle and those in the core. Since the core is warmer than the huddle temperature Th, and the mantle is cooler than Th, the bees i n the core are moving outwards, while the bees on the mantle are moving inwards. T h e change i n radius is somewhat faster for higher ambient temperatures, presumably because the densities are lower and there is less resistance to the mantle bees moving closer together. Figure 3.6 shows the dependence of the total metabolic heat output of the cluster on b o t h the ambient temperature and the cluster size. For T < 0°C, these results are coma parable w i t h the experimental data of Heinrich (1981a, Figure 9) and Southwick (1983; 1988). For ambient temperatures above 0°C the data of Heinrich (1981a), which are taken from swarms, decrease linearly with increasing ambient temperature. 1 Times mentioned are relevant to the model and are not computation times. However, Chapter 3. The Thermoregulation Model 42 Ambient Temperature °C Figure 3.5: G r a p h of the steady state cluster radius versus the ambient temperature. T h e results are taken from the same simulations as those of Figure 3.4 Chapter 3. The Thermoregulation Model 43 800 700 600 500 400 300 200 100 10 15 20 Ambient Temperature °C Figure 3.6: Numerical computations of the metabolic output of the entire cluster for several cluster sizes. These computations are taken from the same simulations as the temperatures of Figure 3.4. Note that at temperatures below 0°C, the metabolic output of the cluster increases w i t h a decrease i n ambient temperature. A l s o , the metabolic output per bee is relatively independent of cluster size. Chapter 3. The Thermoregulation Model 44 Southwick's data (Southwick, 1988), which are taken from winter clusters, show a broad m i n i m u m of metabolic output at ambient temperatures near 10°C, similar to our model. Figure 3.7 shows the same data as the metabolic output per bee. Note that the output per bee does not change significantly between the two cluster sizes. A l t h o u g h , the metabolic output per bee is slightly lower for the larger cluster. Temperature a n d Density Profiles: T y p i c a l temperature and density profiles are shown i n Figure 3.8. There is a region of higher density corresponding to the mantle and a lower density in the core. However, this mantle does not extend to the edge of the cluster, as experimental observations suggest. Note that only a relatively small increase in the radius of the cluster is required for a significant reduction i n the thickness of the mantle. Since the cluster is spherical, the decrease i n the thickness of the mantle does not result i n a significant decrease i n the number of bees i n the mantle. In fact, the model predicts that the majority of the bees w i l l be found i n the mantle even at higher ambient temperatures. Figure 3.9 shows the metabolic output of the bees as a function of their position w i t h i n the swarm. T h e upper graph shows the cumulative fraction of the metabolic output of the swarm attributed to each shell of the numerical simulation. There were 20 such shells (grid points). Thus, for the two lowest temperatures shown, over 70% of the metabolic output of the cluster was produced by bees i n shells 19 and 20. These bees represent the outer 30% (by volume) of the cluster. A t 10°C, 50% of the heat is produced by the bees i n shells one through fourteen, which represent only 34% of the total volume of the cluster. T h e lower graph shows the local metabolic heat production per unit volume divided by the average metabolic heat production per unit volume of the entire cluster. This clearly shows that at high temperatures the heat is produced uniformly over the volume of the cluster, whereas at low ambient temperatures, the bees Chapter 3. 6 I 05 I ' -30 The Thermoregulation Model 1 i -25 1 i -20 45 1 1 1 i -15 i -10 i -5 1 1 1 r i 0 i 5 i 10 i 15 I 20 Ambient Temperature °C Figure 3.7: T h i s figure shows the same data as Figure 3.6, but as metabolic output per bee rather than as metabolic output of the entire cluster. Note that there is surprisingly little difference between the small and the large clusters. T h i s agrees w i t h the experiments of Heinrich (1981a). Chapter 3. The Thermoregulation Model 46 T a = 10 Ta = 0 T a = -10 Ta = -20 2 a, a 1 1 1 1 1 1—: / / co CD CP /// /.•••''/ /// CO CD _ / //' 5 r£2 Q T a = 10 Ta = 0 T a = -10 Ta - -20 \ — - _ \ 6 a o 1 '// X - \ \ ^ \ Radius c m Figure 3.8: Numerical computations of the temperature and density profiles for a cluster of 16000 bees. The thickness and density of the mantle decreases as the ambient temperature increases. T h e decrease i n thickness of the mantle does not imply that there are fewer 'mantle bees'. T h e number of bees i n the mantle will increase as the square of the cluster radius. T h e temperature profiles are i n a rough qualitative agreement w i t h the data of Heinrich (1981a). Note that these graphs are i n the dimensional variables. Thus, the size of the cluster is given by the endpoint of the curves. Chapter 3. The Thermoregulation Model ^ 47 3.5 - ft +J 2.5 - ft 2 O 2" 0 1 0 1 5 1 10 ' 15 1 20 G r i d point Figure 3.9: These data are taken from the same calculations as Figure 3.8. T h e upper graph displays the cumulative fraction of metabolic output produced by each shell. Note that the volume increases as the cube of the radius, so that even at higher ambient temperatures, the majority of the heat is produced by the outer t h i r d of the cluster. T h e lower graph shows the local heat production per unit volume divided by the total heat production per unit volume to correct for this increase i n the volume w i t h radius. Chapter 3. The Thermoregulation Model 48 on the mantle are the primary heat producers. Figure 3.10 shows the effect of changing the ratio of thermotactic motion to diffusive motion, X /D . 0 0 Recall that this ratio has units of temperature. A s the ratio is increased, the bee density profile becomes sharper. However, the core temperatures increase as the ratio is decreased. For a sufficiently low thermotactic coefficient, the bees no longer maintain the higher mantle density, and the core temperature drops below the huddling temperature (which was set to Th = 25 for this figure). 3.5 A n a l y s i s of the steady state I examine the steady state solutions of a simplified system i n which p(p) is a constant. T h i s allows the study of the dependence of the cluster responses on the thermotactic coefficient X(u). Assuming the motility to be constant allows for unrealistically high densities. Since p(p) is dimensionless, we can assume it to be identically one. To simplify the system, I first use the steady state solution of E q u a t i o n (3.19) to determine the local density of the cluster as a function of the local temperature and the mantle temperature u = u(l). R T h e steady state density will be a solution of 1 9 R?x dx 2 (x ( l-X{u)p ^]\=^ 2 \ w i t h the boundary conditions grated once to obtain d d \dx dp dx w r d x 0<x<l, (3.28) = 0, and p(l) = 1. E q u a t i o n (3.28) can be inte- x=l dp , x du T h e constant of integration must be zero since both . ^ and — are zero at the origin. dx dx D i v i d i n g E q u a t i o n (3.29) through by p and integrating from the edge of the cluster Chapter 3. The Thermoregulation Model 40 1 49 i i i i i I I j 1 X/o oD//OD D=o 9 ;iIX O / ;2i1//53 X o2 -= XO/DO=1/15 1 i — - """"""'] O o <D 30 H i 0 1 i 2 l I 0 1 3 i 2 4 l 3 5 I 4 6 l i 5 7 6 8 l 7 9 l 8 9 10 i I 10 Radius c m Figure 3.10: T h e temperature and density profiles resulting from variations i n the ratio X /D . increasing this ratio leads to a sharper peak i n the density profile, and a decrease i n the core temperatures. However, for sufficiently low values of the ratio, the mantle can no longer be maintained, and the core temperature decreases significantly. T h e densities are given i n b e e s / c m , and the temperatures i n degrees Celsius. 0 0 3 Chapter 3. The Thermoregulation Model 50 inwards, obtains the relation X(u)u = [ (3.30) dx, x (3.31) X{u) du. JUR Thus the density profile p can be expressed as a function of the temperature profile u and and the (dimensionless) mantle temperature UR p(u, u ) = exp ( / R \JUR (3.32) X(u) du Since X(u) is a monotonic decreasing function with a zero at Uh = T D /X , 0 the density are dimensionless temperatures h 0 w i l l have a m a x i m u m at Uh provided UR < Uh < u(0). ( uu -— u^ u \ h Using X(u) = — tanh I — - — ), where Uh and X r yields Xr / c o s h ( ^ ^ ^ P( ^R)= V cosh ( (3.33) -u U u h X V r In dimensional quantities, this relation is TX r cosh T T-T p(T, T ) = p R 0 ~D7 r R cosh (3.34) h T T Figure 3.11 shows two plots of this function. Figure 3.11 (a) shows the dependence of p(T, 10°C) on the ratio X /D . 0 0 Recall that this ratio was a measure of the temperature scale of the problem (page 36). Increasing this ratio has the effect of raising the peak in the density profile. This result is not surprising, since increasing the ratio represents an increased dominance of tactic, or directed motion, which should increase aggregation about the preferred huddling temperature T . h dependence of the profile on this ratio. W h a t is important is the exponential T h e densities change from relatively flat to Chapter 3. The Thermoregulation Model 51 Figure 3.11: These two plots show the dependence of the bee density on the local cluster temperature, and the mantle temperature. This dependence is predicted by the steady state model for constant conductivity and motility, as shown i n Equation (3.34). (a) T h e dependence of the profile on the thermotactic ratio X /D . (b) T h e dependence of the profile on the thermotactic range T . 0 r 0 Chapter 3. The Thermoregulation Model 52 unrealistically peaked over a narrow range of this ratio. Figure 3.11 (b) shows the effect of T on the temperature profile. This parameter affects the sharpness of the profile, but r does not have such a strong effect on the height of the maximum. 3.6 Discussion T h e range of the model is restricted by the assumption that convective heat transfer through the cluster is negligible (Assumption A 4 ) . This assumption is definitely false for ambient temperatures above 15°C, and possibly below this threshold as well. Unfortunately, a model incorporating convection would require a far more complicated geometry than the model presented here. Southwick and M o r i t z (1987) observe that the air w i t h i n the cluster is purged periodically by organized fanning. This is followed by a slower influx of fresh air. They speculate that the fanning activity of the bees is triggered by high levels of carbon dioxide. This suggests that convective heat exchange may also occur i n a periodic manner when the increased density of the mantle reduces the air flow through the cluster, and causes a buildup of carbon dioxide. If such a buildup of carbon dioxide coincides w i t h the increased metabolism of overheating cores, which is likely, then this may also provide a signal to bees on the mantle to allow air though the cluster to cool the swarm. Heinrich (1981a) performed several experiments to determine the effect of introducing CO2 into swarms. Experiments with individual bees (Heinrich, 1981a), conclusively showed that a burst of C 0 will stimulate a temporary rise i n metabolic 2 output. T h e results on whole swarms, however, were highly variable. H i s speculations i n light of the experiments are that low concentrations of C 0 core will result i n a temporary swarm core temperature. increase i n 2 injected into the swarm heat production and a possible rise of the Whereas, injections of high concentrations of CO2 into the swarm core "invariably result i n mantle opening, and heat loss." Thus, it appears that Chapter 3. The Thermoregulation Model 53 excessive amounts of CO2 i n the swarm core can trigger the opening of channels through the cluster which result in a rapid and significant heat loss. T h e formation of the mantle layer, and an inverse relationship between the core temperature and the ambient temperature results from the dynamics of the individual bees. These results have been observed i n several experiments, but have not been explained by any previous models. Myerscough's basic model (Myerscough, 1993) d i d not produce the observed inverse relationship between the core and ambient temperatures. However, this feature was produced w i t h her set point model (Myerscough, 1994). O u r model, which uses a similar metabolic rate per bee as Myerscough's basic model, is able to reproduce this feature. Thus, this feature may be explained as either a result of shivering to a set point, or of the thermotactic movements of the bees. Based on the available experimental evidence, both hypotheses have equal merit. A l l of these models produce broader temperature profiles i n the core region than those observed i n live clusters. Our model predicts the formation of the mantle and the core previously observed. These predictions are grounded on the motions of the individual bees. M u c h experimental work is now necessary to measure the parameters introduced i n the model, and to examine the density profiles of thermoregulating clusters i n more detail. Unfortunately, all observations of density profiles have been qualitative. These results support the speculation (Moritz and Southwick, 1992) that the active thermoregulatory behaviour is performed by the bees on the mantle. A s the ambient temperature falls below a critical value, these bees begin to shiver i n an attempt to keep their local temperature above a minimum set point. This results i n an increase i n the temperature at the cluster core. Indeed, this heating of the core has been noted to approach lethal temperatures i n some isolated cases (Heinrich, 1981a). T h e fact that the bees on the mantle w i l l attempt to keep warm even at the expense of overheating the core of the cluster strongly suggests that the thermoregulation is self-organized. Such an Chapter 3. The Thermoregulation Model 54 overheating should not arise if either the heating is performed by the core bees i n order to keep the bees on the mantle above a certain set point, or if the heating and huddling are performed by the bees on the mantle i n order to keep the bees i n the core at a set point. T h e overheating of the core suggests that the bees on the mantle are huddling and shivering i n an attempt to keep themselves within a certain temperature range, and that this desire is sometimes i n conflict with the attempts of the core bees to remain at the preferred temperature of 36°C. Chapter 4 The Pheromone Transmission M o d e l In the previous chapter I examined a system of bees whose interactions are mediated through an external field (the cluster temperature). The system studied i n this chap- ter also has an external field (a chemical concentration i n the wax comb of the hive). However, the dominant interaction is a direct exchange of chemical between the bees. 4.1 Introduction The behaviour of honey bees depends strongly on the presence or absence of a healthy queen. Colonies containing a healthy queen are referred to as q u e e n r i g h t (Winston and Slessor, 1992). In general, queenright colonies show a higher pollen intake and a greater honey production per worker than queenless colonies (Higo et al., 1992). Workers in a queenless colony may develop ovaries and attempt to raise new queens, whereas in queenright colonies, both queen rearing and worker ovary development are usually suppressed (Winston, 1987). One exception is that queen rearing always precedes colony swarming i n the early spring, even i n queenright colonies. Butler (1960b) showed the existence of a substance produced by the queen that suppresses queen rearing and worker ovary development, and Seeley (1979) demonstrated that this substance is transmitted to the workers by direct contact. However, it is only recently (Slessor et al., 1988; W i n s t o n et al., 1991) that the substance involved has been isolated and synthesized. It is now known to be a non-volatile pheromone produced 55 Chapter 4. The Pheromone Transmission Model 56 in the queen's mandibular glands and transmitted through the colony by direct contacts between workers and the queen, and by secondary contacts between workers. T h e pheromone is also deposited on the wax comb by the queen and then retrieved by the workers. N a u m a n n et al. (1991) showed that the pheromone is removed from the system by absorption by both the bees and the wax. W i n s t o n et al. (1991) found that the application of the pheromone to a colony can delay swarming, allowing colonies to grow to far larger sizes than normal. It is also known that supplying a congested colony with additional hive space can delay swarming (Winston, 1987). A n hypothesis linking these two observations is that colony congestion hinders transmission of the pheromone and decreases the amount of pheromone received by each bee, and that this decreased pheromone level is a trigger for swarming (Seeley, 1979; Seeley, 1983; W i n s t o n , 1987). In this chapter, I develop a model for the transmission of the queen pheromone through a honey bee colony. T h e following questions are addressed by the model: • W h a t is the effect of colony congestion on the transmission of the pheromone? • H o w important is the wax as a transfer mechanism? • H o w long does the pheromone remain in the system after the queen is removed? Since b o t h relieving colony congestion and applying the pheromone can delay swarming, it is important to determine the connection between colony congestion and pheromone transmission. T h e second question has not yet been answered by either experiment or theory. T h e model should be able to predict the level of pheromone on the wax. Finally, the pheromone signal must not persist for long times, since the colony must begin rearing new queens i n the event of the death or loss of the old queen. Experiments indicate a time lag of one half hour between removal of a queen and the loss of queenright behaviour Chapter 4. The Pheromone Transmission Model 57 (Winston and Slessor, 1992). T h e decay of the pheromone level i n the colony should reflect this. In Section 4.2, I introduce the background for the model. T h e model is based on the experiments of N a u m a n n et al. (1991), which isolated the mechanisms of pheromone transmission. In Section 4.3, I summarize these experiments and show how they can be used to determine the rates at which pheromone is produced and transferred though the colony. T h e model is developed i n two parts, a Lagrangian model (Section 4.4) and an Eulerian model (Section 4.5). A Lagrangian model is developed i n Section 4.4.1. T h e model consists of a system of (stochastic) differential equations describing the evolution of the pheromone on the bees and the wax. T h e system is described as a finite number of bees, a single queen, and a wax substrate. Each bee is characterized by a spatial position and a pheromone level. T h e wax is represented as a continuous distribution of pheromone. A cellular automata ( C A ) simulation is used to examine the Lagrangian model (Section 4.4.2). T h e results of the simulations are presented i n Section 4.4.3. In Section 4.5.1, I develop an Eulerian model for the evolution of the worker probability density function. T h e model consists of a system of partial integro-differential equations ( P I D E s ) . In order to make the transition from the Lagrangian model to the Eulerian model, I replace the behaviour of each individual bee w i t h the average or expected behaviour of the bees i n a small neighbourhood. This simplified model is used to estimate the expected pheromone level of each bee and the spatial distribution of pheromone on the bees and i n the wax. T h e model is similar to the previous model for cluster thermoregulation. However, the exchange of thermal energy has been replaced by the exchange of a chemical pheromone, and the bees have become pheromone sinks rather than energy sources. In its simplest form, this model reduces to a heat equation w i t h the queen represented by a randomly moving source. This is discussed i n Section 4.5.3. Chapter 4. The Pheromone Transmission Model 58 Using the model, I predict that a reduced motility of the bees w i l l lead to a larger variation i n the distribution of the pheromone. It is therefore plausible that colony congestion may constrict movement and hinder the transmission of the queen pheromone. T h e analysis also reveals that the wax does not play a significant role i n the pheromone transmission. A l t h o u g h large amounts of pheromone can be found on the wax, only a small fraction is picked up by the workers. Finally, the decay rate of the pheromone i n the entire system is consistent w i t h the time lag seen between the removal of the queen and the loss of queenright behaviour. These results are discussed i n detail in Section 4.7. 5 4.2 4.2.1 Background R e p r o d u c t i v e s w a r m i n g , supersedure, a n d absconding Honey bee colonies reproduce by swarming. In the early spring, when the colony is large, the workers w i l l begin rearing young queens. Soon thereafter, the old queen and a portion of the workers will leave the hive i n search of a new nest site. Swarm prevention is of major interest to beekeepers the world over, since swarming not only results i n a loss of over half the colony population, but also a significant amount of potential honey production that season. Queen rearing is a necessary precursor to swarming. T h e primary stimuli for queen rearing are thought to include "colony size, brood nest congestion, worker age distribution, and reduced transmission of queen substances." (Winston, 1987, p 191). However, each of these stimuli is influenced by resource abundance. A l l of these stimuli are necessary for swarming. For example, the addition of frames to a hive w i l l delay swarming. Y e t , colonies are commonly utilizing only 50% to 70% of the comb when they swarm (Winston and Taylor, 1980). Thus, relieving congestion can delay swarming, although uncongested colonies will swarm. Pheromones produced by the queen are thought to be Chapter 4. The Pheromone Transmission Model 59 a significant indicator of both the health and age of the queen and of colony congestion. There are two variations of swarming that are less common (Winston, 1987, p 197). Supersedure can occur if the old queen does not leave w i t h a swarm after the onset of queen rearing. The young queen emerges from her cell and kills the old queen. Supersedure is presumably triggered by a reduction in the efficiency of the old queen either due to age or injury. It is plausible that her pheromone production is no longer sufficient to suppress queen rearing, but the other triggers necessary for swarming are absent. However, this has never been proven. Absconding refers to the swarming of the entire colony and the abandonment of the old hive. Absconding may be the response of queenright colonies to a reduction of resources. 4.2.2 The honey bee queen substance Butler and Simpson (1958) showed that a specific substance produced in a mated honey bee queen's mandibular gland can suppress queen rearing. The substance is also present in virgin queens, but at much lower levels than in mated queens (Butler, 1960a). The substance has recently been identified as a blend of five chemicals (Winston and Slessor, 1992), and is now referred to as the queen mandibular pheromone (QMP). It is nonvolatile, and is transferred from the queen to the workers by direct contact and v i a the wax comb of the hive. The experiments of N a u m a n n et al. (1992) show that the pheromone is transmitted through the colony at a fixed composition. The chemicals do not separate during transfer. Q M P appears to be used by the bees both as a primer and as a releaser pheromone. Releaser pheromones elicit an immediate behavioural response when detected by a bee. T h i s effect is seen in the change in behaviour observed when bees come into contact w i t h the pheromone. Primer pheromones do not necessarily invoke an immediate response, but trigger a change in the behavioural repertoire or physiology of the bees. Chapter 4. The Pheromone Transmission Model 60 T h e preparation for swarming and the rearing of new queens are both possible effects of primer pheromones. If the queen is removed from the hive, the bees become visiblyagitated. However, if they are then sprayed w i t h a solution of synthetic Q M P , the behaviour of the bees returns to that of a queenright colony (Winston et al., 1991). T h e primer effects of the pheromone have been identified as suppression of worker ovary development, suppression of queen rearing, increased colony size, increased foraging and pollination, and increased per capita honey production (Winston et al., 1991). Since queen rearing and colony congestion are precursors to swarming, the pheromone also appears to delay swarming. However, the amount of the pheromone produced by the queen does not decrease prior to swarming (Seeley and Heinrich, 1981). One hypothesis is that if pheromone transmission is hindered by colony congestion (Winston and Taylor, 1980), then queen rearing would no longer be suppressed i n larger colonies. To discuss this hypothesis, we must first examine the mechanisms of Q M P transmission through the colony. T h e pheromone is produced in the queen's mandibular gland and transferred over the surface of her body during grooming (Butler and Simpson, 1958). Seeley (1979) and van der B l o m (1992) have made careful observations of the behaviour of bees that contact the queen. Bees i n the immediate vicinity of the queen are attracted to her and form a r e t i n u e around her. U p o n entering this retinue, bees begin removing pheromone onto their antennae. After a short time, a fraction of these retinue bees begin removing pheromone using their mouthparts. These two behaviours are quite distinct, and are termed a n t e n n a t i n g and l i c k i n g respectively. After leaving the queen's retinue, bees are more likely to contact other workers. They also perform fewer labour acts than other workers. This change i n behaviour lasts approximately 30 minutes after leaving the queen's retinue (Seeley, 1979). During this time these bees are referred to as m e s s e n g e r s . Seeley (1979) observed that the queen directly contacts only 35% of the workers over a ten hour period. The remainder of the workers receive pheromone either from the Chapter 4. The Pheromone Transmission Model 61 messenger bees or indirectly v i a the wax (see below). V a n der B l o m (1992) has shown that the messenger bees are not specialists. Retinue and messenger behaviour is displayed by all workers of the age normally found i n the brood area. Neither the probability of a bee joining the queen's retinue, nor the probability of a retinue bee becoming a licker depend on previous retinue encounters. However, P a n k i w (1994) discovered that the number of workers that display a retinue response when presented w i t h a glass lure coated with Q M P depends on the amount of pheromone present. She observed that the fraction of bees attracted to the lure increased w i t h the level of the pheromone up to a maximum, but then decreased as pheromone levels increased above this maximum. Since this maximum was above the level normally found on mated queens, it is likely that the decrease in the response represents a revulsion to unnaturally high concentrations of the pheromone. T h e experiments of Naumann et al. (1991) showed that pheromone on a bee's body surface is slowly absorbed through the cuticle. Pheromone is also deposited on the wax by b o t h the queen and the workers, and is acquired by workers as they inspect cells and move over the wax. T h e pheromone deposited on the wax is gradually absorbed into the wax. These experiments also demonstrated that pheromone is transferred during the licking and antennating behaviours of the workers, and that large quantities of pheromone are ingested during licking behaviour. In summary, Q M P is produced by the queen, removed by absorption through the bees' cuticles and into the wax, and ingested by the workers. The steps i n this process are shown i n the schematic of Figure 4.1. Since the pheromone is non-volatile and is transferred though the colony by contact, it is possible that colony congestion hinders movement and therefore slows pheromone transmission through the colony. T h e movement of the queen is composed of small scale local motions and stationary behaviour, punctuated by larger scale movements from one Chapter 4. The Pheromone Transmission Model 62 area of the hive to another. The queen spends most of her time either laying eggs, grooming, or exchanging food (Seeley, 1979; Fergusson and Free, 1980). D u r i n g each of these activities the queen is either stationary or moving short distances. Five percent of the queen's time is spent moving from one area of the brood comb to another. These larger scale movements consist of fairly straight runs of several centimetres, which often end on a different frame. The workers also spend much of their time engaged i n stationary activities such as grooming, exchanging food, or inspecting cells. Their movements consist of small displacements of 3 to 6 cm approximately once per minute. The displacements are concentrated i n straight walks of 1 to 10 seconds duration. Messenger workers tend to move faster, and hence show larger displacements (Seeley, 1979). P r e l i m i n a r y observations (Jeff Pettis, unpublished) suggest that the queen's movements are not hindered by crowding, whereas the worker displacements are reduced by crowding. 4.3 Determination of the transmission rate constants 4.3.1 Summary of the transmission pathways T h r o u g h many careful experiments, Slessor et al. (1988) isolated and synthesized Q M P . Using the synthesized pheromone, N a u m a n n et al. (1991; 1992) were able to determine the rates of pheromone transfer between the queen, the workers, and the wax. Figure 4.1 shows a schematic of the pathway for pheromone transmission through the colony. N a u m a n n et al. (1991) conducted several experiments to determine the rates of pheromone transfer at each step. The notation used for the rate constants is the same as that used i n their paper. N a u m a n n et al. (1991; 1992) assume that all pheromone transfers follow first order kinetics. This implies that the amount of pheromone transferred, absorbed, or ingested at each stage is proportional to the pheromone level of the source. The rate constants used by N a u m a n n et al. (1991) are defined as follows: Chapter 4. The Pheromone Transmission Model 63 Figure 4.1: The pathways for Q M P transmission were determined by N a u m a n n et al. (1991) by a series of experiments. The pheromone is produced by the queen, and transmitted to the workers through direct contact, indirectly through the wax, and by secondary and tertiary worker-worker exchanges. Naumann's experiments show that a large fraction of pheromone removed from the queen is immediately ingested. Pheromone is removed from the system by ingestion, and by absorption by the bees and the wax. The rate constants ki are defined on page 62-64 k 0 is the r a t e at which pheromone is t r a n s f e r r e d to the queen's body s u r f a c e . A large portion of the pheromone produced i n the mandibular gland is ingested during this process. This ingestion is included i n the rate constant ko. Once on the queen's body, the pheromone is available to the workers. k\ is the r a t e o f p h e r o m o n e a b s o r p t i o n t h r o u g h t h e q u e e n ' s c u t i c l e . Over time, pheromone becomes tightly bound to the cuticle or absorbed into the body and can no longer be transferred to the workers. This is treated as a pheromone sink i n the model. &2 is the rate of pheromone absorption through the worker's cuticle. Pheromone is absorbed into the workers i n the same manner as it is absorbed into the queen. Once pheromone is absorbed through the cuticle, I assume that it is removed from the system. Thus, I assume that there is no reverse Chapter 4. The Pheromone Transmission Model 64 process of ingested or absorbed pheromone returning to the body surface. &3 is the r a t e at w h i c h the queen transfers p h e r o m o n e to the workers. The amount of pheromone transferred is i n direct proportion to the amount of pheromone present on the body surface of the queen as well as the number of workers i n her retinue. This rate constant is not first order and can be separated into the two rates k and kg. 8 /c is the 4 r a t e at which workers acquire pheromone f r o m the wax. This is the amount of pheromone removed by a single bee from a wax surface per unit time. Note that this rate does not depend on the area of the wax, since the bee only occupies a portion of wax during each instant. Thus, the units are still those of a first order rate constant. k is the r a t e o f d e p o s i t i o n o f p h e r o m o n e o n t h e w a x by the queen. k§ is the r a t e o f d e p o s i t i o n o f p h e r o m o n e o n t h e w a x by the workers. 5 &)7 is the rate of absorption of p h e r o m o n e into the wax. A l t h o u g h the phero- mone absorbed into the wax cannot be picked up by workers, it can be detected by them. In fact, observations suggest that pheromone i n the wax is used as a cue by the workers (Free, 1987). However, since this cue is not directly related to queenright behaviour, I model this mechanism as a sink. ks is the r a t e o f p h e r o m o n e r e m o v a l f r o m q u e e n o r m e s s e n g e r b y a s i n g l e antennating bee. k$ is the r a t e o f p h e r o m o n e r e m o v a l f r o m q u e e n o r m e s s e n g e r b y a s i n g l e licking bee. Pheromone is transferred from the queen to the workers during worker licking and worker antennation. Pheromone transfer during each of these activities occurs at two Chapter 4. The Pheromone Transmission Model different rates. 65 T h e rate constants k$ and kg are the specific removal rates by single antennating and licking bees respectively. T h e rate of pheromone transferred from the queen directly to the workers (fc ), depends on the numbers of workers engaged i n each 3 activity. Thus, while k is not a first order rate constant, it is assumed that k$ and k 3 9 are. 4.3.2 Calculation of the transmission rate constants The rate constants were determined i n a series of experiments that isolated each link in Figure 4.1. Since only relatively high concentrations of the natural pheromone can be detected, the only rates that can be determined using current assay methods are the rates of production and absorption of pheromone by the queen. Slessor (1988) has developed a synthetic pheromone that can be labelled w i t h t r i t i u m . T h e radioactivity of the sample can be used to accurately measure quantities of the tritiated pheromone at much lower levels than the natural pheromone. Thus, the tritiated pheromone can be used to determine the rates of pheromone transfer between workers and between the workers and the wax. In the remainder of this section, I summarize the experiments of N a u m a n n et al. (1991) and calculate the values of above rate constants. T h e results are summarized i n Table 4.2. T h i s section is included to show how the rate constants were determined. T h e mechanisms responsible for pheromone transmission are not fully understood, nor do the experiments point to a unique model. Hence, I feel it is important to present the complete detail of how the transmission rate constants were determined. A b s o r p t i o n i n t o w a x (A/7): To determine the rate of pheromone absorption into the wax, N a u m a n n et a l . (1991) coated several sections of wax w i t h tritiated pheromone and left them undisturbed for Chapter 4. The Pheromone Transmission Model notation description 66 range of values production by queen 1 ng/s absorption by queen (8.0±2.4)xlfJ- S"-1 absorption by workers 6.8xlCr s " h k mean transfer from queen to worker h deposition on wax by queen (3.5 ± 1.4)xl0- s--1 (6.2 ± 2.8)xl0" S"-1 (8.0 ± 2.4)xl0" S"-1 h deposition on wax by workers k k$ absorption into wax ( 1 . 3 ± 0 . 1 ) x l 0 " " s--1 ( 1 . 0 ± 0 . 2 ) x l 0 " " s--1 acquisition by antennators (3.5 ± 1.2)xl0" k acquisition by lickers ( 3 . 2 ± 0 . 6 ) x l 0 " " S"-1 h 4 7 9 transfer from wax to workers t 4 4 1 4 5 5 4 3 5 S "-1 3 Table 4.2: A summary of the pheromone transfer rates for the pathways shown i n F i g ure 4.1. T h e values and notations are taken from N a u m a n n et al. (1992). f values vary from 50 pg to 10 ng/s. (page 68) varying times. To determine the measure of the biological availability of the pheromone on the wax, they then introduced workers to the wax for five minutes. After this time, the workers were dissected, and the pheromone they acquired was measured. T h e data of N a u m a n n et a l . (1991, Figure 8) indicate that the amount of pheromone retrieved from the wax by the workers decreases from 25 ng to 4 ng i n 30 minutes (1800 s). A s s u m i n g that both the rate of pheromone absorption into the wax and the rate of pheromone pickup from the wax by the workers are first order, the rate constant k-j satisfies 4 = 2 5 e ~ Thus, k = l . O x l O 7 - 3 s - 1 1 8 0 0 f c 7 . . T h e error i n the data is approximately 20%. T h e data (Naumann et a l , 1991, Figure 8) do not in fact show continued first order absorption into the wax, but rather a steady state of 2 pg. This may indicate that there is an upper limit to the amount of pheromone that can be absorbed into the wax, or that there may be a reverse process whereby pheromone previously absorbed becomes available to the workers again. However, the amount of pheromone returned to the Chapter 4. The Pheromone Transmission Model 67 workers through this mechanism would be negligible. T r a n s f e r f r o m q u e e n t o w a x (k ) 5 and absorption by queen (ki): To determine the rate of pheromone transfer from the queen to the wax, N a u m a n n et al. (1991) placed queens coated with tritiated pheromone on a section of wax comb. After varying times, the queens were removed, and both the queens and the comb were analysed to measure the amount of labelled pheromone transferred. If we assume that the transfer of pheromone from the queen to the wax, and the absorption of pheromone by the queen's body both follow first order kinetics, then Q = -hQ-hQ, (4.1) W = k Q, (4.2) 5 where Q is the pheromone present on the body surface of the queen, and W is the pheromone present on the wax. Pheromone is removed from the queen by absorption through her cuticle (ki) and by transfer to the wax (k ). 5 Note that i n the experiment, the total pheromone (both pheromone on the wax and pheromone absorbed into the wax) was measured. Hence, the absorption of pheromone into the wax {k-j) does not affect W. G i v e n that Q(0) = 14 ng and W ( 0 ) = 0, the solution of Equations 4.1 and 4.2 is Q(t) = Q(0)e- , (4.3) W(t) = J^(Q(0)-Q(t)). (4.4) (fcl+fc5)t The d a t a collected after 5 minutes gave an average of Q(300s) = 10.7 ng, and W(300s) = 0.25 ng. Combining these values w i t h the value for kj determined i n the previous paragraph gives k = 7 x l 0 5 data. - 5 s , and ki = 9 x l 0 - 1 - 4 s . - 1 There are no errors given for the Chapter 4. The Pheromone Transmission Model Production by queen 68 (ko): In order to determine the rates of queen pheromone production, N a u m a n n et al. (1991) removed several queens from their hives and isolated them i n petri dishes for specified times. A f t e r their isolation, individual queens were removed from the dishes, washed w i t h methanol to measure the pheromone levels on their body surface, and dissected to determine the amount of pheromone i n their mandibular glands and other body parts. (These experiments measured natural, not synthetic pheromone.) T h e amount of pheromone that could be removed from the body surface with methanol appeared to reach a steady state of 1.0 ± 0.3 /ig after two hours. Assuming that the pheromone is transferred to the body surface from the mandibular gland at a constant rate (ko), and that pheromone is removed from the body surface by absorption (ki) and transferred to the petri dish (A? ), a steady state balance implies that k — (ki + k )Q = 0, where Q is the 5 0 5 pheromone on the body surface of the queen. This gives a value of 0.9 ± 0.3 ng/s for ko, the rate of pheromone transfer from the mandibular gland to the body surface of the queen. T h e variations quoted do not account for the possible error i n the estimates of k\ and k . 5 The pheromone found i n the glands increased w i t h the isolation time. Naumann et a l . (1991) assume that this increase is due to a steady production of pheromone, and calculated an average production rate of 204 pg per day, or 2.4 ng/s. However, the production rates observed for individual queens varied between 12 and 400 pg per day. The reasons for this large variation are not known, but may be due to differences i n age, species, or time of day. T h e steady increase of pheromone i n the gland is puzzling, especially since the amount on the body surface appears to have reached a steady state. In this thesis I take ko to be the r a t e o f p h e r o m o n e t r a n s f e r t o t h e b o d y surface o f t h e q u e e n , and assume it to have an order of magnitude of 1 n g / s . T h i s avoids a l l Chapter 4. The Pheromone Transmission Model 69 concern about ingestion of pheromone as it is transferred from the gland to the body surface. T r a n s f e r b y l i c k i n g (k$) a n d a n t e n n a t i n g (/eg): To determine the rate of pheromone transfer from the queen to workers and the rate of transfer from workers to workers, N a u m a n n et a l . (1991) allowed workers to contact a glass lure coated with tritiated pheromone. In five seconds of antennating contact w i t h a lure coated w i t h 8.5 ng of pheromone, workers acquired 1.5 ± 0.5 pg of pheromone each. Thus, k = (3.5 ± 1 . 2 ) x l 0 ~ s . D u r i n g 30 seconds of licking contact w i t h a lure 5 s - 1 coated w i t h 250 ng of pheromone, workers acquired 23 ± 4 ng of pheromone each. Thus, kg = (3.2 ± 0 . 6 ) x l 0 -3 s . However, approximately 50% of the pheromone transferred _ 1 during licking contacts was immediately ingested by the workers. M e a n transfer from queen to workers (k ): 3 The ratio of lickers to antennators i n the queen's retinue depends on the activity of the queen, but is usually near 1:9 (Seeley, 1979; N a u m a n n et al., 1991; W i n s t o n and Slessor, 1992). T h i s ratio decreases when the queen is moving, and increases slightly when she is stationary. If we assume that the ratio of lickers to antennators is fixed at 1:9 then the two transfer rates for lickers and antennators calculated above can be replaced by the average rate per worker given by k = (9k + k )/10. Since 50% of the pheromone is ingested 3 8 9 during licking encounters, the fraction of the average ingested is I = (0 x 9 + 50)/10 = 5%. This assumption is used for the Eulerian model of Section 4.5. Chapter 4. The Pheromone Transmission Model 70 A b s o r p t i o n b y w o r k e r s (A^): To determine the rate of absorption of pheromone by the worker cuticle, N a u m a n n et a l . (1991) treated workers w i t h 14 ng of tritiated pheromone and isolated them for various times. Analysis of the proportion of pheromone that could be washed from the worker's body using methanol after these times suggested a half-life of 13 minutes (k = 8 . 9 x l 0 ~ s 4 2 - 1 ) for live workers and 17 minutes (k = 6 . 8 x l 0 ~ s 2 4 - 1 ) for dead work- ers. These two results are not significantly different, nor is there a significant difference between these values and those for the adsorption by the queen (fci). T r a n s f e r b e t w e e n w o r k e r s a n d w a x (A; , ki): 6 To determine transfer from worker to worker v i a the wax, N a u m a n n et a l . (1991) first allowed workers to contact a pheromone-coated lure for 30 seconds, and then placed them on a section of wax. After a further five minutes, these workers were removed and analyzed. U p o n removal from the wax, the 'messenger' workers were replaced w i t h 15 other workers for a further five minutes. analyzed. Finally, these workers and the wax were In the first five minutes, individual workers each transferred 130 ± 9 pg of pheromone to the wax. Analysis showed these workers to have 3.3 ng of pheromone each. Since the analysis counted the total pheromone on the workers, including that bound to the cuticle, it follows that k = (1.3 ± 0 . 1 ) x l 0 ~ 6 4 s _ 1 . D u r i n g the second five minutes, each worker recovered 2.2 ± 0.9 pg of the 130 ± 9 pg of pheromone available on the wax. T h i s gives a value of (6.2 ± 2 . 8 ) x l 0 " 5 s " for k . 1 4 A l t h o u g h this value is not the same as that determined by N a u m a n n et a l . (1991), it does agree w i t h the one calculated i n N a u m a n n et al. (1992). There are i n fact, two sets of experiments presented i n N a u m a n n et al. (1991) that can be used to determine A; , and there is a difference of an order of magnitude between the two results. T h e 4 Chapter 4. The Pheromone Transmission Model 71 value for fc calculated i n N a u m a n n et al. (1991) was based on an experiment which 4 involved an application of the tritiated pheromone directly to a section of wax. The value presented here was based on an experiment i n which the pheromone was removed from a lure by a worker and then deposited on the wax by that worker. A similar experiment was performed by N a u m a n n et al. (1992). Since the second setup is a more realistic representation of the physical process, the values presented here are expected to be reliable. The difference i n retrieval by the workers i n the two situations is difficult to explain. However, it does appear that the absorption of pheromone into the wax is not a first order mechanism. 4.4 The Lagrangian model In this section, I develop a system of stochastic differential equations (SDEs) that describes the rate of change of the positions and pheromone levels of the individual bees and the level of pheromone on the wax. The equations are based on the transmission pathways and rate constants calculated in the previous sections. T h i s section is divided into three parts. In Section 4.4.1, I list the assumptions made about the behaviour of the bees and derive the S D E model. In Section 4.4.2, I examine the system using a C A simulation. Finally, i n Section 4.4.3 I discuss the results of the simulation. The Lagrangian model is also a staging ground for the development of the E u l e r i a n model presented i n Section 4.5. The motivation for developing the stochastic model rather than going directly to the Eulerian model is twofold. First, this approach provides a link between the Eulerian model and the C A simulation. Second, the Lagrangian description provides a better framework for outlining the assumptions of the model and basing the model on the observations of individual behaviour introduced i n Section 4.2. Chapter 4. 4.4.1 The Pheromone Transmission Model 72 Development of the model I begin this section w i t h the definition of the variables used for the model and a description of the domain representing "the hive. This is followed by a listing and discussion of the assumptions of the model. The. development of the model is divided according to the three m a i n types of pheromone transfer occurring i n the system, absorption, exchange between bees and wax, and exchange between bees. A b s o r p t i o n can be modelled by a simple exponential decay; however the remaining exchange rates require a more detailed implementation. F i r s t , I discuss the rates for pheromone exchange between the bees and the wax and use these to develop an equation for the pheromone level of the wax. Second, I develop the notation used for pheromone exchange between bees. Finally, I combine these processes to develop the equations for the pheromone level of the queen and then the workers. T h e hive is represented as a square with sides of length R and is designated by Q = (0, i?)x(0, R). T h e typical artificial nest consists of a hive box w i t h five frames, each 9 inches by 12 inches. Bees can move from frame to frame v i a holes i n the frames and spaces around the frame edges. Since the bees can access both sides of a frame, it seems reasonable to choose R — 83 c m , which would give an area similar to that of the total area of five frames. Swarming usually occurs soon after the colony reaches 20 000 bees ( M a r k W i n s t o n , personal communication). Hence, I am interested i n colonies of that size. (Note that this is roughly 3 bees/cm .) 2 The complete system, consisting of the wax frame, N bees, and a single queen, is represented by the following state variables: Chapter 4. The Pheromone Transmission Model w(x,t) level of pheromone per unit area on the wax substrate at (x,t), x (i) position of the queen at time t, Q(t) pheromone level of the queen at time t, Xi(t) position of bee i at time t, Pi(t) pheromone level of bee i at time t. q 73 The positions x, Xi, and x q are vectors i n Q, and the pheromone levels Pi, Q, and w are non-negative real numbers. (The index i runs from 1 to N.) T h e wax pheromone level, w(x,t), represents only the pheromone on the wax that is available for pickup by the worker bees. Pheromone that has been absorbed into the wax is assumed to have been removed from the system. T h e pheromone level of the queen, Q, is the amount of pheromone on the body surface of the queen. Assumptions: The following three assumptions are used for the model presented i n this subsection: | A l : | Worker and queen movements are random walks w i t h zero drift. I assume that the movements of the bees are not aligned w i t h any of the state variables in the system. T h i s does not imply that the movement of the bees is completely random. For example, the queen's movement would presumably be correlated w i t h the positions of filled and empty brood cells. T h e assumption states that the queen's movements are not correlated w i t h either the distribution of pheromone i n the wax or the positions of the workers. Similarly, the worker movements are not correlated w i t h the pheromone levels or the position of the queen. (Although the workers respond to the queen's presence, there is no indication that they actively seek her out.) Chapter 4. The Pheromone Transmission Model 74 | A2: | Individual bees occupy a finite area, and deposit pheromone evenly over this area. Worker and queen interactions occur only if these areas overlap. T h i s assumption affects both the pheromone exchange between the bees and the wax, and the pheromone exchange between the bees. The pheromone is presumably deposited on the wax by 'footprints' left by walking bees and by contact between the bee and the wax during cell inspections and other activities. | A3: | Worker retinue formation is a releaser effect of Q M P (see page 59). Therefore, workers w i l l form retinues around any object coated w i t h sufficient quantities of the pheromone, be it queen, worker, or lure. Strictly speaking, the 'retinue' refers to the group of bees surrounding the queen at any given time. However, for the purpose of studying the transmission of the phero- mone, there appears to be no reason to distinguish between the retinue response of the workers to the queen and their response to messenger workers. Experiments indicate that workers will form retinues about glass lures coated w i t h the pheromone. O n this basis I assume that retinue formation and licking and antennating behaviour depend solely on the amount of pheromone on the queen. Further, I assume that pheromone acquisition from the queen and pheromone acquisition from messenger bees use similar mechanisms. Thus, for the purpose of this model, the sole difference between the queen and the messengers is the fact that the queen is the source of pheromone, whereas the workers are a sink. To simplify the presentation of the model, bees acquiring pheromone from messengers are also referred to as retinue bees, simply because they are displaying retinue behaviour. To my knowledge, there have been no observations of multiple worker retinues formed around messenger bees. In fact, Seeley (1979) records only antennations and inspections Chapter 4. The Pheromone Transmission Model 75 of messengers by single bees. I suggest that this may be due to the lower pheromone level of the messenger bees relative to that of the queen, rather than any further distinction between the messengers and the queen. Hence, the model should predict that multiple bee retinues rarely form around messenger bees. Worker a n d queen movement: Under Assumption A l , a bee's position is a random variable w i t h zero drift and a variance of a . I assume that the bees cannot enter or leave the hive. This can be modelled by 2 imposing reflecting boundary conditions. N o other positional dependence of the motion is assumed. Recall that the worker movements can be approximated by displacements of between 0.5 and 6 c m each minute "(Section 4.2). I assume that the bees are either stationary or moving w i t h a fixed velocity. T h e bees switch between these behaviours randomly and spend on average T stop = 54 s stationary and T move = 6 s moving. T h e movement speed is always 1 c m / s . T h e movement direction is a random variable w i t h a uniform distribution that changes with each transition from stationary to moving. A s discussed i n Section 2.3.1 this is a jump process with a step length i n the range 0.5 < 8 < 6 c m at time intervals of r « 60 s. Such a process has a variance i n the range 0.008 < o /2 < 0.3 c m / s . 2 2 Pheromone exchange with the wax: Assumption A 2 states that each worker occupies a specific region. I take this to be a disc of area A ^ . T h e queen is assumed to occupy a slightly larger disc of area A . T h e bees q are not rigid discs w i t h well defined edges. Thus, the regions occupied by neighbouring bees may overlap. T h e area Ab is the area over which the bee interacts w i t h the wax and other bees. In order to allow for a similar approach for both the C A model and the Eulerian model which follow, I define the neighbourhood function Chapter 4. The Pheromone Transmission Model *y 6 76 1 if \x — y\ < r, 0 otherwise. (4.5) = If b is the 'radius' of a bee then r S b xy dx = irb = A . 2 (4.6) b The typical radius of a worker is b = 0.5 c m . Queens are roughly twice the size of the average worker. Thus, Ab = 6 7r = 2 0.257T c m , and Aq = 2A . 2 b A section of wax at position x receives pheromone from each bee that overlaps position x. I assume that these bees deposit pheromone evenly over their area of contact at the rate k . Thus, the rate of pheromone deposition at position x is 6 ( - ) 4 " 7 1=1 units pheromone per unit area per unit time. T h e rate of pheromone removal from the wax at position x is proportional to the number of bees overlapping position x, and the amount of pheromone on the wax. This is expressed as N kw £ 6 b 4 . xx i=l units of pheromone removed per unit time per unit area of the wax. A s shown i n Figure 4.1, the rate of change of the pheromone level of the wax is due to absorption (fc ), deposition by workers (k ), loss to workers (fc ), and deposition by 6 7 the queen (fc ). 5 4 In each case, let the pheromone transfer be first order. T h e level of pheromone on the wax satisfies ^ = -kinj + f > ^ 6 - k w)6 4 b XXt + h^5l . Xq (4.8) Chapter 4. The Pheromone Transmission Model 77 Pheromone is absorbed into the wax (thereby becoming unavailable for pickup by workers) at the rate k . T h e first term i n the summation is the total pheromone on all workers 7 that overlap position x, divided by the total area over which this pheromone is deposited (Ab). T h e second term i n the summation is simply the number of workers w i t h i n one bee radius of x. T h e final term is the pheromone deposited by the queen evenly over her area of contact w i t h the wax. Worker-queen and worker-worker pheromone transfer: To model the retinue behaviour, I introduce the matrix of exchange rates for the pheromone transfer between the bees, (e^ is the rate of pheromone removal from bee i by bee j.) r k if bee j is antennating bee i, k if bee j is licking bee i, 0 otherwise. 8 9 (4-9) Since a l l exchanges are treated identically, exchange with the queen is represented by appending the row i = q to the matrix. T h e entries are random variables w i t h values in the three states {k , fc , 0}. The times spent i n each state depend on the pheromone levels 8 9 of the bees i n the encounter radius. Pankiw (1994) performed experiments to determine the fraction of bees that contact a messenger as a function of the pheromone level of the messenger. C o m b i n i n g her results with Seeley's (1979) observations of the mean time spent licking and antennating, we can make assumptions about the times bees spent i n the retinues of both the messengers and the queen. Pheromone exchange can occur between bees whose encounter discs overlap. Thus, the probability of entering a bee's retinue depends on the pheromone level of that bee and its location. Pankiw's data on retinue response to various amounts of pheromone Chapter 4. The Pheromone Transmission Model 78 indicate a quadratic dependence on the logarithm of the pheromone level (Pankiw et al., 1994). However, the m a x i m u m of this quadratic is above the normal pheromone level of the queen. I approximate the relevant section of her results w i t h a linear function of the logarithm of the pheromone level. Thus, the probability of entering a retinue increases linearly with the logarithm of the queen (or messenger) pheromone level from a probability of zero at a t h r e s h o l d l e v e l T t o a m a x i m u m at the s a t u r a t i o n l e v e l S. If the times spent near the queen (or messenger) before entering her retinue are exponentially distributed with a mean of Tenter then At p | i enters retinue of j y C - enter At r i n time interval At T enter ns< , Pj logfe/T) log(5/r) c i f r < p . < 5 > (4.10) " * XiX 0 if P j < T. I have assumed that bees exchange pheromone only if their separation is less than twice their radius. Hence, the encounter radius of a bee is 26 and the encounter area A = e (26) 7T = 7r c m . 2 2 Retinue bees begin as antennators and become lickers after spending an average time of T ant antennating. Once bees become lickers, they remain lickers until they leave the retinue. The average time spent i n a retinue is denoted T . Recall that the mean times stay spent i n each activity were observed by Seeley (1979). If bee j is i n the retinue of bee i then P r { i leaves retinue of j i n time interval At} = (At/T )8f. . stay x + (1 - 6? ) (4.11) Further constraints are that workers cannot be i n more than one retinue at a time, and that messengers cannot be i n a retinue. In summary, the above probabilities define the transitions of the entries of the matrix {eij} between the values {kg, kg, 0}. Chapter 4. The Pheromone Transmission Model 79 T h e level of pheromone o n the queen: Referring to Figure 4.1, the level of pheromone on the queen's body surface increases due to transfer from the mandibular gland (fc ) and decreases due to absorption (k\), 0 deposition onto the wax (k ), and transmission to the workers (e j). These losses are 5 q proportional to the level of pheromone on the queen's body surface. Hence, the pheromone level of the queen evolves according to the relation (4.12) T h e level of pheromone o n the workers: The pheromone level of an individual bee changes due to absorption, loss to wax, gain from wax, and exchange w i t h other workers or the queen. T o compute the rate of change of the pheromone level of the workers, I first develop an expression for the rate of pheromone exchange with other workers, and then for the rate of pheromone transmission from the wax to the worker. Pheromone absorption into the cuticle and transmission to the wax are simply proportional to the pheromone level of the worker. The rate of pheromone removal from messenger bee i by workers i n her retinue is SjLi ijPi e ( p g / ) - However, the rate at which pheromone is acquired by retinue bee j is (1 — Ij)eijPj, s where Ij is 0 or 0.5 if retinue bee j is antennating or licking respectively. The difference between the two (j2iLi Ij ijPi) e 1S the amount of pheromone ingested by the retinue bees surrounding messenger i. T h e rate of increase of the pheromone level of worker i is N (4.13) The rate of pheromone transfer from the wax to bee i is found by integrating over Chapter 4. The Pheromone Transmission Model 80 the region occupied by this bee, j^w(y,t)6 dy^A w(x ,t). b x%y b (4.14) i The approximation assumes that w varies slowly on the length scale of an individual bee, so that 8 b xy f» A 8(x — y). b C o m b i n i n g a l l of the above effects, the pheromone level of worker i evolves as follows: -jl = -(k 2 + k )p + k A w(x ,t) 6 l 4 b i + g + e (l - h)Q{t). l qi (4.15) The first term is the rate at which pheromone is absorbed into the cuticle and lost to the wax. T h e second term is the rate of pheromone gain from the wax by contacts over the encounter area A . This gain is proportional to the level of pheromone on the wax under b the bee. T h e final two terms represent transfer between workers and transfer from the queen. In summary, Equations (4.8), (4.12), and (4.15) describe the evolution of the pheromone levels on the wax, the queen, and the individual workers. These are studied using a C A simulation i n the next section and are used as a basis for an Eulerian model i n Section 4.5. 4.4.2 T h e cellular automata simulation To study the S D E s developed i n the previous section, I have developed a cellular automata ( C A ) simulation. T h i s is a numerical simulation of Equations (4.8), (4.12) and (4.15), using a forward Euler discretization. In this section, I describe the C A simulation. T h e results are presented i n Section 4.4.3. T h e main considerations are the discretization of the hive, the approximation of the movement of the bees, and the description of the retinue formation in the discrete setting. T h e d o m a i n is approximated by a regular hexagonal grid. Each point on the grid has a state value indicating the local pheromone level i n the wax. I arrange the grid Chapter 4. Symbol The Pheromone Transmission Model Value Parameter 0.25TT c m 0.57T c m A e T 7T cm 2 6s - top move 54 s T T stay 30 s 1 s T 2 2 81 Page A r e a occupied by a single worker 76 A r e a occupied by the queen 76 Encounter area of workers 78 Mean time spent moving by worker 75 Mean time spent stationary by worker 75 Mean time spent i n retinue by worker 78 ant T ter 10 s Mean time of antennating behaviour 78 6s Mean time spent near queen before entering her retinue 78 Saturation pheromone level for retinue behaviour 78 T 1 Pg 100 pg Threshold pheromone level for retinue behaviour 78 M 4 bees M a x i m u m number of bees per grid point 80 N 20 000 Number of bees 72 R 83 c m Length scale of hive 72 x en s Table 4.3: Summary of parameters, and the values used i n the C A simulations. A d d i tional parameters are shown i n Table 4.2 Chapter 4. The Pheromone Transmission Model 82 so that the area of a single hexagon is the area occupied by a single worker (A ). b This requires a grid spacing of A x = 1 c m . Recall that the areas occupied by neighbouring bees may overlap (page 76). This implies that several bees may occupy a single grid point. The movement of the bees is broken into a speed and a direction of motion. T o simplify the simulation, bees move one grid spacing per unit time. T h e average speed of the bees is roughly 5 m m / s , which implies a time step of At — 2 s. T h e simulation uses reflecting boundary conditions. Thus, bees attempting to move off the grid are reflected off the boundary. Congestion is modelled using the following assumption: | A 4 : | Movement onto a grid point becomes less probable as the number of bees on that point increases. T h e probability reduces to zero for a finite number of bees (denoted M). T h i s assumption states that there is an upper limit to the number of bees that can occupy the same grid point. Further, restrictions to movement begin before this m a x i m u m is reached. Bees that are unable to move onto a grid point have their direction reversed. This condition does not hold for the queen who may move onto congested grid points. Recall that colony congestion does not hinder the queen's movements (Section 4.2). D a t a s t r u c t u r e s : T h e data for the simulation is contained i n an array of states for the wax and a vector of states for the bees. T h e values for p were taken as unsigned long integers (32 bits), and the values for x were taken from a hexagonal grid of R rows and R diagonals. T h e wax elements consist of a pointer to a linked list of bees on the grid point and an integer value for a pheromone level. T h e bee state is defined by a pheromone level, a number of retinue bees, and a pointer to a messenger bee. T h e use of the retinue counter serves to eliminate the possibility of messenger-messenger chains. There is no Chapter 4. The Pheromone Transmission Model 83 limit to the size of the retinues, and the encounter probabilities dd not depend on the retinues. T h e sole constraint is that a messenger bee cannot be i n a retinue. A t each time step, the worker state is updated according to one of the following three possibilities: 1. If the worker is a messenger, it will move randomly. 2. If the worker is not a messenger and is not in a retinue, it will search the neighbouring grid points plus its own grid point and (with the probability defined by E q u a t i o n (4.10)) enter the retinue of the bee with the highest pheromone level i n the neighbourhood. If the worker does not engage with the host, it will move as described above. 3. A worker i n a retinue will remain an average time Tgtay of leaving is At/Tstay Workers spend a average time T Thus the probability ant antennating before switching to licking behaviour. Thus, if the worker remains i n the retinue, it will begin licking w i t h probability At/T . ant If the worker is not on the same grid point as the messenger, it will move there with probability PFOIIOW T h i s allows for the retinue bees to remain w i t h the queen as she moves. T h e queen simply produces pheromone at a constant rate. Since she cannot enter a retinue, she just moves randomly, oblivious to the retinue of workers surrounding her. C o o r d i n a t i o n o f t h e t i m e i n t e r v a l s Each time step, the bees are processed i n the sequence stored i n memory. (This sequence remains the same throughout the simulation.) T h i s is followed by the processing of each grid point, also i n sequence. There has been some concern in the literature with the use of synchronous updates of cellular automata versus an asynchronous update (Huberman and Glance, 1993). Since the spatial patterns formed i n this simulation agree with the spatial patterns formed from the Eulerian formulation (Section 4.7), I do not expect this to be a problem. Chapter 4. The Pheromone Transmission Model 84 T h e purpose of the simulation is to determine the effect of the details of b o t h the movement of the queen and the discrete nature of the retinue on the pheromone distribution. E v e n though the number of bees i n the simulation may seem quite large, the number of messenger-worker encounters is not, especially when we consider the density of encounters per unit area of wax. For this reason, it is important to run the simulation w i t h a large number of workers. In addition, when one considers that the bees are distributed not just i n space, but i n pheromone level as well, the densities of bees are quite small. 4.4.3 Results of the cellular automata simulations The results of the simulations support the hypothesis that colony congestion hinders pheromone transmission. Figure 4.2 and Figure 4.3 show the results of two nearly identical simulation runs after one hour. T h e simulation run shown i n Figure 4.3 has 2000 bees, while the simulation run shown in Figure 4.2 has only 1500. E a c h shaded circle i n the figures represents a single bee. T h e lightest shading indicates a pheromone level below 1 pg, and the darkest shading indicates a pheromone level above 100 ng. Intermediate shadings increase w i t h the logarithm of the pheromone level. W h i t e spaces indicate a space without a bee. B o t h simulations use a 25 x 25 domain, w i t h a m a x i m u m of four bees per grid point ( M = 4). Therefore, the output is represented by a 50 x 50 hexagonal grid of bees and spaces. T h e open circles in the lower left had corner of the figure represent the position of the underlying grid points. Note that w i t h 2000 workers distributed over 625 spatial points, there are on average 3.2 bees per point, whereas 1500 bees distributed over the same number of points gives an average of only 2.4 bees per point. T h e size of the domain and the number of bees correspond to those found on a single frame i n the typical artificial hive used by bee keepers. However, under normal conditions, the queen is not restricted to a single frame Chapter 4. The Pheromone Transmission Model 85 Figure 4.2: Pheromone distribution in a uncongested colony. Each disc represents a single bee, w i t h the shading of the disc proportional to the logarithm of the bee's pheromone level. T h e lightest shade corresponding to bees with less than 1 pg of pheromone, and the darkest shade to bees w i t h more than 100 ng. T h e simulation was run w i t h 1500 bees on a 25 x 25 grid, and a maximum of four bees per grid point. Hence, the hive is represented by a 50 x 50 hexagonal array of bees and spaces. T h e open circles i n the lower left corner show the positioning of the underlying grid. Figure 4.3: Pheromone distribution in an congested colony. The simulation was run with 2000 bees and the same parameters as Figure 4.2. Each disc again represents a single bee. T h e shadings have the same significance as in Figure 4.2. Chapter 4. The Pheromone Transmission Model 87 Congested colony Mi 550 500 500 450 450 400 350 300 400 Mi SH 350 300 HH 250 CD CO SH CD WO CO S-H CD Uncongested colony O SH CD 200 250 200 150 150 um B 100 50 100 50 0 0 0 10 20 30 40 50 60 Pheromone level of worker 70 10 20 30 40 50 60 70 Pheromone level of worker Figure 4.4: These histograms show the distribution of pheromone on a l l workers i n the colony. T h e data are from the same r u n as those of Figure 4.3 and Figure 4.2. Thus, there are 2000 workers i n the congested colony and 1500 workers i n the uncongested colony. T h e height of each bar represents the number of workers w i t h pheromone levels in the given range. However, the rightmost bar i n each graph (centered above 64 pg) represents the number of bees with pheromone levels greater than 60 pg. Note that the majority of the 500 additional workers i n the congested colony have pheromone levels below 4 p g . Chapter 4. The Pheromone Transmission Model 88 as she is i n these simulations. Thus, the pheromone distribution should be qualitatively similar to that i n a natural colony, but the pheromone levels i n the natural colony should be lower. Methods for scaling these values are discussed i n Section 4.5. The figure shows two main differences between the congested and uncongested colonies. • T h e congested colony shows a larger portion of bees w i t h less than 1 p g of pheromone. • T h e average pheromone level per bee drops off more noticeably w i t h distance from the queen i n the congested colony. T h i s first observation is shown more clearly by the histogram of Figure 4.4. T h e number of workers i n each range of pheromone levels is similar except for the number w i t h pheromone levels below 1 pg. In the uncongested colony, the pheromone distribution shows less of a decay w i t h distance from the queen. Another interesting effect of congestion is the memory of the previous position of the queen shown i n Figure 4.3. This is another indication of the poor transmission of the pheromone away from the queen. The results clearly show that the transmission of the pheromone is hindered by restricted worker movement, even though the queen's movement is not restricted. Note however, that no trace of the queen's passage can be seen i n the pheromone distribution of the workers. This may be due to either a lower retinue size while the queen is moving resulting i n a slower transfer of pheromone, or a low transfer of wax from the wax to the workers. This will be discussed further i n the next section. The pheromone levels of the wax (from the same simulation runs) are shown i n F i g ures 4.5 and 4.6. T h e two figures do not differ significantly. T h e previous positions of the queen are retained longer i n the pheromone distribution of the wax than they are i n the pheromone distribution of the workers. Note that the shadings used i n the figures for the Chapter 4. The Pheromone Transmission Model 89 Figure 4.5: W a x pheromone distribution i n a congested colony. T h e data are from the same simulation r u n as that used i n Figure 4.3. Each disc now represents a single grid point. T h e lightest shading represents a pheromone level below 1 pg, and the darkest shading a pheromone level above 10 ng. A logarithmic scale was used for shadings between these levels. Chapter 4. The Pheromone Transmission Model 90 Figure 4.6: W a x pheromone distribution i n an uncongested colony. The data are from the same simulation run as that used i n Figure 4.2. T h e shadings have the same significance as i n Figure 4.5 Chapter 4. The Pheromone Transmission Model 91 wax have a different scaling than those used for the pheromone levels of the bees. T h e lightest shading of Figure 4.6 still represents a pheromone level below 1 pg. However, the darkest shading now represents a pheromone level above 10 ng. T h e path of the queen can be clearly seen i n both cases. Figure 4.7 shows the initial spread of pheromone from the queen for b o t h congested colonies and uncongested colonies. Although the differences are not great, the messenger bees spread faster i n the uncongested colony. In the next section I show that the half-life of the pheromone level on an isolated worker is roughly 15 minutes. Thus, the spread of pheromone on this time scale is significant. (Recall also that messenger behaviour typically lasts up to one half hour after leaving the queen's retinue.) In summary, colony congestion slows the spread of pheromone from the queen. T h i s results i n higher pheromone levels near the queen and lower pheromone levels farther away, and i n a larger number of workers without significant amounts of pheromone i n more congested colonies. T h e pheromone level of the wax does not change significantly w i t h colony congestion. 4.5 The Eulerian model In this section, I develop and study an Eulerian model for the pheromone transmission problem. T h e model examines the probability density of workers as a functions of their pheromone level and position i n the hive, and therefore consists of a system of equations for the probability density function (p.d.f.) describing the distribution of workers (Section 4.5.1). T h e equations contain a non-linear integral term that makes analysis difficult. A s a simplification, I examine the moments of the p.d.f. w i t h respect to the pheromone level of the workers (Section 4.5.3). T h e equations for these moments can be reduced to a pair of linear ordinary differential equations for the evolution of the expected Chapter 4. The Pheromone Transmission Model 92 Figure 4.7: A time sequence showing the initial spread of the pheromone i n a congested and an uncongested colony. The left column shows a sequence for a congested colony of 2400 bees on a 25 x 25 grid w i t h a maximum of 4 bees per grid point. T h e right hand column shows the same sequence for an uncongested colony of 1500 bees on a 25 x 25 grid w i t h a maximum of 4 bees per grid point. Note that the messenger bees (indicated by intermediate grey levels) have spread farther from the queen in the uncongested colony. Chapter 4. The Pheromone Transmission Model 93 pheromone levels of the wax and the workers as function of position. T h e solutions of these equations clarify the effect of congestion on the transmission of the pheromone. 4.5.1 Derivation of the Eulerian model Several assumptions are made i n deriving the simplified model of this section. I introduce these assumptions i n three stages: the distributions of the workers and the pheromone level of the wax, the retinue of the queen, and the retinue of the messengers. These assumptions are used to simplify the equations for the evolution of the pheromone levels on the wax, queen, and workers given by Equations (4.8), (4.12), and (4.15) respectively. For convenience, these equations are repeated here. rim — dQ dt -g N (1 n- = -krw + Zfa^-hw^ = k- = -(ki + k^pi + kiAbwixuQ + gi + egiil-IJQtt). 0 ^kx + ks + J^e^j + h^S^. Q- (4.16) (- ) 4 17 (4.18) The purpose of this section is to develop a system of equations (using an Eulerian description) to approximate the above system. This is accomplished by introducing a p.d.f. for the distribution of the workers over space and pheromone level. T h i s p.d.f. w i l l evolve according to the Fokker-Planck equation with an additional flux through pheromone space as described i n Section 2.3.1 (Equation (2.13)). The first set of assumptions introduce functions for the probability density of the workers and pheromone density on the wax. | A5: | T h e probability of finding a given worker at pheromone level p and position x at time t can be expressed by the probability density function u(x,p,t). A6:1 u(x, p, t) and w(x, t) do not vary significantly over the length scale of a single bee. Chapter 4. The Pheromone Transmission Model 94 According to A s s u m p t i o n A 5 , the probability that bee i is i n region S and has a pheromone level between a and b at time t is P r j x i G S, a < pi < b} = f f u(x,p,t) dpdx. Js Ja (4-19) Since every bee is located i n the hive and has a non-negative pheromone level, r poo / Jo Ja u(x,p,t)dxdp = l. (4.20) Further, since the probability given by Equation (4.19) must be non-negative for a l l values of S, a, and b, u{x, p, t) > 0, for a l l x e £1, p > 0. (4.21) A s s u m p t i o n A 6 allows the approximation S b XXi « A 8(x - ), b (4.22) Xi where 6(x) is the two dimensional delta function and A is the area occupied by a single b bee. T h i s implies that the expectation of the sum of an arbitrary function f(xi,pi,t) of a bee's state variables over all bees within one bee radius of position x can be expressed as an integral over a l l pheromone levels TV E 52f(xi,Pi,t)6 . .i=i b = NE[f(xi,p t)8 h b ], XiX x x = N f f Ja Jo f(xi,pi,t)A 8(xi = NA f(x,p,t)u(x,p,t)dp. roo b Jo b - x)u(xi,pi,t) dpidxi, (4.23) E q u a t i o n (4.23) is used to remove the dependence of the pheromone level on the wax on the exact distribution of the workers. T h a t is, the summation over all bees (the second term of E q u a t i o n (4.16)) is replaced by an integration over a l l pheromone levels. The second set of assumptions simplifies the structure of the retinues formed around the queen and the rate of pheromone transfer from the queen to the workers. Chapter 4. The Pheromone Transmission Model 95 A 7 : 1 T h e ratio of lickers to antennators in the queen's retinue is constant. | A8: | A l l bees in the vicinity, of the queen will enter her retinue. Recall, from the discussion in Section 4.3, that the ratio of lickers to antennators i n the queen's retinue depends on the activity of the queen, but is usually near 1:9. If we assume this ratio to be constant, then the two rates for licking and antennating can be combined into a single rate constant representing the rate of pheromone removal from the queen per worker in her retinue, k$. A s s u m p t i o n A 7 allows the replacement of the licking rate, kg, and the antennating rate, kg, w i t h the average transfer rate, k . Combining this w i t h A s s u m p t i o n A 8 , which 3 states that all bees within one encounter radius of the queen will be in her retinue, gives the approximation e i ~ k^S^ g , which, when combined w i t h E q u a t i o n (4.23), this can be used to replace the summation in the final term of Equation (4.17) w i t h an integration over all pheromone levels. It can also used to simplify the final term of Equation (4.18). T h e final two assumptions simplify the messenger-worker interactions. These simplifications are used to replace the exchange rate ^ w i t h an integration of the p.d.f. u(x,p, t) over all pheromone levels. A9:1 Messenger encounters always involve antennating and never involve licking. | A10:1 T h e fraction of the bees with a pheromone level of q and w i t h i n one encounter radius of messenger bee i that are in retinue of messenger bee i can be expressed as a function of Pi and q. Recall from Section 4.2 that only antennations have been reported during messenger encounters (Seeley, 1979). Hence, I assume that the transmission rates to kg, for all encounters. can be set A s discussed in Section 4.2, this also implies that ingestion Chapter 4. The Pheromone Transmission Model 96 is negligible during messenger-worker encounters. These assumptions were not made for the C A simulation of the previous section, since they should, ideally, result from the behaviours of the bees. They are made here i n order to develop a simpler Eulerian model. Let K(p, q) to be the probability of a bee w i t h pheromone level q being i n the retinue of a bee w i t h pheromone level p given that the distance between the two bees is less than one encounter radius. Using this notation and Assumption A 9 , the rate of pheromone transfer from messenger i to worker j is E[ ] = k K(pi, )8 eij 8 (4.24) e Pj Under A s s u m p t i o n A 9 , all messenger-worker exchanges are antennations. Recall that the ingestion fraction Jj is negligible for antennations (page 59). Summing this expression over the entire population gives the expected rate of pheromone transferred to bee i. Thus, the expected value of gi from Equation (4.13) is N E jiPj e j=0 'N E = (4.25) ijPi e (hK{Pj,Pi)Pj k NA J^ 8 E ksKipup^p^S^. - (K(p,pi)p- K(pi,p)p^u(xi,p,t)dp. (4.26) (4.27) The first term i n the integral is the expected rate of pheromone transfer to worker i from the potential hosts i n her encounter neighbourhood. T h e second term is the expected rate of pheromone transfer to workers i n her retinue. T h e expected rate of pheromone exchange is considerably smoother than the actual rate. However, the model should still give some indication of how the various parameters affect the transmission. Thus, given a worker with pheromone level p at position x, the expected rate of change of the pheromone level of that worker due to pheromone exchange w i t h other workers is E[si|Pi —p,Xi = x] = NA k G(p, e 8 u), Chapter 4. The Pheromone Transmission Model 97 where roo G(p,u)=l Jo (K(q,p)q- K(p,q)p)u(x,q,t)dq. (4.28) Thus Equation (4.15), the rate of change of pheromone on bee i is governed by the equation Pi = ~(k + h) 2 Pi + k A w(xi, t) + fc (l - I)A Q8(xi - x ) + E[g ]. 4 b q 3 q (4.29) t Using E q u a t i o n (4.29) and the derivation of the Fokker-Planck Equation i n Section 2.3.1 leads to a P D E for the p.d.f. u(x,p,t). W i t h the above approximations, the system of Equations (4.16), (4.17), and (4.18) becomes rod dw dQ — = —kyiu + N J • (k p — k Abw)u(x,p,t) dp + k Q6(x — x ), = k - (hi + h + Nk A Q Q = 3 ~ dp {~ ^ + + 4 q J o ' 5 u(x ,p,t)dp^Q, q p u + q (4.30) (4.31) ^ k AbWU I)A Q6(x - x )u + NA k G(p, q q e 8 u)), (4.32) where A is the Laplacian i n two spatial dimensions (i.e., the domain Q ) . (Note that NA q / °° u(x ,p, t) dt is the expected number of bees i n the neighbourhood of the queen.) 0 q In summary, Equations (4.30), (4.31) and (4.32) represent the evolution of the expected pheromone levels of the wax, the queen, and the workers respectively. I have approximated the rate of pheromone exchange between messengers and workers based on the distribution of worker pheromone levels at each point. T h i s simplification results in a continuous pheromone transfer between workers, and does not properly represent the fact that the workers move a finite distance between contacts. Chapter 4. 4.5.2 The Pheromone Transmission Model 98 Nondimensionalization of the system Equations (4.30), (4.31), and (4.32) can be non-dimensionalized using the following rescalings: independent variables dependent variables =R x u = (h + k )(k + 5 fc (l k 3 _ - I)k A fc (l = 0 3 (h + h)(k 1 t = — k + k e 2 - I)'koA ^ h + k _ k Nk (l w = k^kx + k^fa (4.33) 5 6 , 6 h) 0 + k )R 2 2 q 2 I)k A 3 0 q + k^R* For example, t* = it is the time i n dimensionless units. For convenience and readability I do not introduce new notation for the dimensionless variables. T h e reader may assume that all variables and functions are dimensionless from this point forward, unless clearly noted otherwise. Using the values summarized i n Tables (4.2) and (4.3) gives the following values for the rescalings: x = 83 cm, i = 21 min, u = Q = 1 pg, p = 34 pg, w = 4 x 10~ p g c m ~ , 6 _ 1 2 (4.34) 13pgcm~ . 2 T h e spatial scale x simply reflects the size of the domain. Recall that x is a position in two dimensions, and therefore has two components. B o t h components are scaled by x. T h e temporal scale is related to the rate of decay of pheromone on the workers. Recall that k is the rate of pheromone absorption into the workers cuticle, and that k is the 2 & rate of pheromone transmission from the workers to the wax. Hence, /Fin 2 fa 14 minutes is the half-life of the pheromone on an isolated worker. T h e pheromone level of the queen is rescaled by the ratio of the production rate k to the rate of pheromone removal from 0 the queen by absorption and transmission to the wax (ki + k ). Therefore, Q is the steady 5 state pheromone level of a solitary queen isolated on a section of wax. T h e rescaling of the Chapter 4. The Pheromone Transmission Model 99 pheromone level of the workers is the product of Q, the rate of pheromone removal from the queen to the workers (fc ), and the rate of removal of pheromone from the workers 3 (k + k ). Thus p is the expected steady state pheromone level of an isolated worker 2 6 in contact w i t h a pheromone source of pheromone level Q. Finally, the rescaling of the pheromone level of the wax is the product of p and the ratio of k (the rate of transmission e of pheromone from worker to wax) to kj (the rate of absorption of pheromone into the wax). Therefore, w is the steady state pheromone level of the wax i n a colony whose workers a l l have a pheromone level of p. These rescalings reflect the transmission of pheromone from the queen to the workers and finally to the wax. T h e scaling of the p.d.f. (u) preserves the value of the integral of Equation (4.20). These rescalings lead to the following dimensionless system of P I D E s : — = dt dw T —= dt dQ W T— = Q pAu — — (—pu + op f°° —w — /3w / Jo 1-Q-aQj^ f°° udp+ Jo f 00 + euG(p, u) + Q8(x — x„)u), TJWU (4.35) pu dp + jQSlx — x„), (4.36) u(x ,p,t)dp. (4.37) q The dimensionless parameters are <7 k^keAbN 2R?(k + hy k R*(k + k y 2 2 7 . k + k 2 ' r 7 q = k + k —, h + h' 2 6 a = 7 e 6 k R (k + 2 5 ' kR 2 8 2 b p kAN R?(k + k y 6 kjNA 6 k Tw 2 ' 7 ~ 2 fc ) 6 ^ fc (i-i)MrV 3 kAN z q R {k + k )' 2 1 5 Using the values summarized i n Tables (4.2) and (4.3) gives the following values for Chapter 4. The Pheromone Transmission Model 100 the parameters: A 4 1 T = 0.02, = 0.8, = 0.9, V = 0.007, e = 0.1, P = 0.04, 7 = 1, (4.39) a — 0.6. Note that the experimental error i n these values is of the same order as the values themselves- Hence, they are order of magnitude estimates only, and have been rounded to the most significant digit. T h e motility, / i , is a measure of the scale of worker movements relative to the dimensions of the hive and the decay rate of the pheromone on the workers. T h e half-lives of the pheromone on the wax and the queen i n relation to t are given by T and r respecq W tively. T h e transmission of pheromone from the wax to the workers is represented by the two parameters r\ and 8, which measure the transmission rates from wax to workers relative to the scalings of the workers and the wax respectively. T h e fact that r\ is small suggests that the pheromone transmission from the wax to the workers has little impact on the pheromone level of each worker. The rate of increase of the level of pheromone on the wax due to deposition by the queen is given by 7. T h e final two parameters, e and a, represent the rate of pheromone transmission between workers and from the queen to the workers respectively. T h e parameters measure the transmission rates per worker (k and k ) relative to the rates of absorption and deposition (fc + fc and k + k ) and 3 8 x 5 2 6 to the number of workers per unit area. Therefore, these parameters increase as colony congestion increases. 4.5.3 Analysis T h e fact that the integral term i n the system is quadratic in u makes analysis difficult. However, many properties of the system can be obtained from a study of the moments Chapter 4. The Pheromone Transmission Model 101 of the p.d.f. u(x,p,t) over p. These moments are defined as n(x, t) m(x, t) v(x, t) o u(p, x, t) dp, (4.40) pu(p,x,t)dp, (4.41) / Jo p u(p,x,t)dp 2 2 (4.42) m . E v o l u t i o n equations for these moments are derived by integrating Equation (4.35) over all pheromone levels. This requires several integrations by parts, which are simplified by the following two assumptions: k = 0,1,2. These assumptions ensure that the boundary terms arising from integration by parts vanish. Since there should be no bees with pheromone levels larger than that of the queen, A s s u m p t i o n A 1 2 is reasonable. T h e reasoning behind Assumption A l l is that, although the pheromone levels of the workers may be arbitrarily small, there should be no workers w i t h zero pheromone levels. This w i l l be violated if some workers initially have pheromone levels of zero. E v o l u t i o n equations for the moments are obtained by multiplying E q u a t i o n (4.35) through by p (k E {0,1, 2}) and integrating once. This leads to the following system of k evolution equations for the mean worker pheromone level, the wax pheromone level, and the queen pheromone level: Chapter 4. The Pheromone Transmission Model 102 n t = pAn, m t = / x A m — m + ^wra + QnS(x — x ), t = ^ A i ; — 2v + V m • V m v (4.43) (4.44) q /•oo + m Q ( l - n)5(x - x ) + 2e / puG(p,u)dp, q T"to^t T ^ q = —w — Pwn + m + 'jQ6(x — x ), = 1 - (l + a)Qn(x ,t). (4.45) (4.46) q (4.47) q Zero-flux boundary conditions are assumed for n(x, t), m(x, t), and v(x, t) at the edge of the domain. T h e analysis of these equations proceeds as follows: F i r s t , Equation (4.43) and E q u a tion (4.47) do not involve m , v, or w and can be solved separately. Second, E q u a - tions (4.44) and (4.46) are decoupled from Equation (4.45) and form a linear subsystem. These equations would be coupled to Equation (4.45) through the integral of G(p, u) that would appear i n Equation (4.44). However, this integral equates to zero under the assumption that there is no ingestion during pheromone exchange between workers (see A s s u m p t i o n A 9 and Equation (4.28)). T h i r d , since 77 < < 1, Equation (4.44) can be completely decoupled from Equation (4.46) by assuming r\ = 0. T h i s leaves E q u a t i o n (4.44) as a linear diffusion equation w i t h a randomly moving source. T h e speed of the source is compared to the motility of the bees to determine the effect of reducing bee motility on the distribution of the pheromone. Finally, the implications of the solution of E q u a t i o n (4.45) when e = 0 are discussed. E q u a t i o n (4.43) implies that the worker density evolves towards the globally stable steady state solution n(x,t) = 1. (Recall that n must satisfy the condition / n(x, t) dx = 1 Chapter 4. The Pheromone Transmission Model 103 where fl is the non-dimensional hive area, which is the unit square.) Since n(x,t) is a constant, Equation (4.47) is decoupled from the rest of the system. T h i s implies that the pheromone level of the queen also approaches a steady state of 1 Qs 1 + (4.48) a' T h e above observations imply that, after a time, the spatial distribution of workers becomes homogeneous and the pheromone level of the queen approaches Q . W i t h these s simplifications, Equations (4.44) and (4.46) become m t = pAm — m + rjw + Q 6(x — x ), (4.49) Tw t = —w — (3w + m + jQ 6(x — x ). (4.50) w s q s q Since these equations are linear, the solutions can be written as a series i n the eigenfunctions of the L a p l a c i a n w i t h zero-flux boundary conditions. The normalized eigenfunctions for the Laplacian on the unit square w i t h zero-flux boundary conditions are 4>jk{x) = 2 cos(j7rx) cos(/c7ry), where x = (x,y). T h e mean pheromone level of the workers and the pheromone level of the wax can be represented by the series m{x, t) V =£ w(x, t) j,fc=0 O-jk (4.51) <j)jk{x). \ 3k ) h Substituting this expansion into Equation (4.49), the coefficients of the series satisfy the relation { 11 -1 - L jk T) jk V l/r -(1 w + + / ? ) / T , J \ 3k b V/ Qs(f>jk(Xq), (4.52) 7 J where Ljk = ( j + k )ix p. T h e eigenvalues of this system are 2 jk A 2 2 --(l 2 +L jk 1+ + ± T,, 2V ( 1 + L * + r. / V 4 + H 'Tin (4.53) Chapter 4. The Pheromone Transmission Model 104 and the eigenvectors are 1 + 3 + V \ (4.54) 1 T h e time scale for the pheromone decay from the system is given by the largest eigenvalue which is \Q 0 — 1. This implies that upon removal of the queen, the slowest mode to decay is A o- Note that the time scaling is i = l/(k 2 0 + k ), which represents 6 pheromone absorption into the worker cuticle. I now examine three cases. First, to isolate the effect of worker movements from the effect of queen movements, I examine the case where the queen is fixed i n position. Second, i n order to estimate the average pheromone level per worker over the entire colony, I examine a case i n which both the queen and the workers have unrestricted longer range movements. This assumes that the bees move from one region of the hive to another on a faster timescale than the pheromone decay. A l t h o u g h the bees and the queen do not move this fast under natural colony conditions, it is physically possible. Finally, I examine the equations with a moving queen using numerical simulations. This allows a comparison to be made w i t h the results of the Lagrangian simulation presented in the previous section. Case I: fixed queen T h i s situation can be tested experimentally by placing the queen i n a wire cage, thus fixing her position, but allowing the workers to contact her. Workers from various positions w i t h i n the hive can then be removed and their pheromone levels analyzed. These pheromone levels can be compared with the solution described below as a test of the model. Unfortunately, the amount of pheromone typically found on the workers is below the m i n i m u m level detectable using current techniques. Once these techniques are Chapter 4. The Pheromone Transmission Model 105 refined, the model can be rigorously tested. If the position of the queen, x , is constant, then the solution of E q u a t i o n (4.52) is q ( O-jk c e V *) t ]k + d e 3k t,k+ ]k , 1 + L k ) , 1 + b where the constants Cj and dj k 1 + B + rpy v l + j(l j 3 ) +L) jk (4.55) j are determined from initial conditions using E q u a - k tion (4.51). Since 77 and B are small, the eigenvalues are both negative. Hence, the steady state solution is ( 1 = V Ws lim t—>oo m{x, t) J (4.56) \ w(x,t) j ~ ^{l ( n + B) + L ){l ]k \ <i> {x). jk 1 + B + 777 l (l + L ) (4.57) Since we are interested i n the case where 77 < < 1, we can approximate the steady K + 7 jk ) state solutions by m s = f; QsMx^jx) (4.58) + ^ - 6 + lQ S(X - Z g ) S 1 + /? + 0(77). (4.59) jk, (4.60) However, note that / ? < < ! , hence !> \-fc - 1 L and —rjjk v/ 1 \ 1 ; (4.61) Chapter 4. The Pheromone Transmission Model 106 where ~ indicates equality to at least 0(^Jrj + y/]3). T h i s indicates that, upon removal of the queen, the pattern indicating the queen's position remains i n the wax as the pheromone is absorbed, however, the movement of the bees causes the worker distribution to smooth out as the pheromone is absorbed. C a s e II: a v e r a g e d p h e r o m o n e l e v e l s Recall that the queen's movement consists of periods of roughly stationary behaviour while laying eggs and grooming, punctuated by long range movements (Section 4.2). Under the simplifying assumption that this movement makes the queen equally accessible to all bees, regardless of position, we can integrate Equations 4.49 and (4.50) over space and compute the steady state average pheromone levels. These are Qs = T~~ J 1+ a (4.62) m s « Q, (4.63) ». » s (4.64) Using the scalings of Equation (4.34) the dimensional values are Qs = S~ 1+ a ~ 720 ng, P 1 ' W 22 pg, + a w(l + y) + + a) * (4.65) 2 * 1 6 p g / C m - It is worthwhile to note that m can be expressed as s fc (l 0 m s - I)/(k + 2 R2( + k )/(k A ) kl 5 3 q h) + N- l 4 - D D j Chapter 4. The Pheromone Transmission Model 107 T h e numerator is the ratio of the queen's production rate to the workers decay rate. T h e denominator is the sum of the dimensionless area of the comb and the absolute colony size N. Thus, the pheromone level of the queen depends on a , which is a measure of colony congestion, whereas the mean pheromone level of the workers depends on the absolute colony size. C a s e III: moving queen Setting TJ = 0 i n Equation (4.49) and introducing the scaled pheromone level z = m / ( l + a) leads to the equation z = pAz — z + 6(x — x ). t (4-67) q Figure 4.8 shows the solution to Equation (4.67) for two values of p. T h e position of the queen (x ) is a random variable whose values are uniformly distributed over the q unit square and change at ten minute intervals. This reproduces the general features of the queens movement (page 62). B o t h simulations use the same sequence for x . T h e q graph uses a logarithmic scale for the pheromone level; thus, a value of zero indicates a pheromone level of m = 1/(1 + a). For smaller values of p, the pheromone level decays by several orders of magnitude as we move away from the queen. Workers i n the queen's retinue have pheromone levels of several hundred picograms, while those far away from regions where the queen has recently visited have pheromone levels well below one picogram. The second moment A s a final stage of the analysis, I examine the second moment of the worker pheromone level (Equation (4.42)). Recall that the integral i n Equation (4.45) arises from the exchange of pheromone between workers as defined by Equation (4.28). If the pheromone Chapter 4. The Pheromone Transmission Model 1 - 0.5 108 i i • •0.5 - - I / / .1 ..... 1 j ! i i / / i / / / /' /' \ \ / i j i i \ \ ; ; i i i i i i i / ; ; ; \ \ / ' . i / / \ V s \ 1 / ' s Figure 4.8: T h e two plots show solutions to Equation (4.67) for p = 0.1 (top) and H = 0.01 (bottom). T h e domain is the unit square. T h e contours show the mean worker pheromone level using a logarithmic scale, since the variation over the length of the hive is large. A s w i t h the C A results (Figure 4.3), the previous positions of the queen can be seen at the lower motility. However, the pheromone level drops far lower further away from the queen. T h e higher motility leads to a more uniform distribution of the pheromone. Chapter 4. The Pheromone Transmission Model 109 levels of the messengers are always higher than the pheromone levels of the workers i n their retinues, then K(p,q) = 0 for p < q (see Equation (4.28)) and the integral of E q u a t i o n (4.45) is always negative. Therefore, v < pAv + 2 / i V m • V m — 2v. t If e = 0, this relationship becomes an equality whose steady state solution is v = s ( 1 + / i J ,; • 2 ) 2 (4-68) Thus, i n the absence of exchange, increasing the motility of the bees increases the local variation of the pheromone distribution. This is interpreted as an increased long-range correlation i n worker pheromone levels, but a decreased short-range correlation i n worker pheromone levels. 4.6 Predictions of the models 4.6.1 The Eulerian model Using the analysis of the Eulerian model the following predictions are possible: • T h e level of pheromone on the queen is a decreasing function of colony congestion (Equation (4.62) and Equation (4.38)). • T h e mean pheromone level on the workers is a decreasing function of the number of bees (Equation (4.66)). T h e worker pheromone levels also decrease w i t h distance from the queen. This decrease is amplified by a reduced motility of the workers. • T h e level of pheromone per square centimetre of the wax comb is comparable to the level of pheromone on the workers (Equation (4.65)). However, the rate of pheromone retrieval from the.wax by the workers is much slower than the rate of pheromone removal from the workers (since r\ « 1). Chapter 4. The Pheromone Transmission Model 110 • U p o n removal of the queen, the disappearance of pheromone from the workers is predominantly due to absorption through the cuticle. Recall that the time scale of the system has been non-dimensionalized by the rate k + k (the sum of absorption 2 6 through the cuticle and transfer to the wax), which corresponds to a half-life of approximately 14 minutes. Hence, and eigenvalue of -1 indicates a decay at this rate. • T h e variance increases w i t h increasing p, (worker motility) under very broad assumptions on the retinue formation function K(p,q). T h e steady state value of the level of pheromone on the queen, Q , is inversely pros portional to the parameter a, and therefore decreases with increased colony congestion. Recall that a is proportional to the average worker density i n the colony (see E q u a tion (4.38)). Thus, the level of pheromone on the queen does not depend on the number of bees, but rather on the number of bees per unit area. Note that this dependence is of the form 1/(1 + a) (see Equation (4.62)), and that for the 'typical' colony preparing to swarm, a = 0.6 (see Equation (4.39)). Thus, as the number of workers doubles, raising a from 0.3 to 0.6 for example, there is only a 20% decrease i n Q (see Equation (4.62)). T h e s importance of this dependence is further decreased by the observation that the retinue response leads to crowding near the queen regardless of the average density of the colony. Hence, the average density of workers as seen by the queen is independent of the average density of the colony. In contrast, the average level of pheromone on a worker, m , does not depend on s colony congestion, but is inversely proportional to the s u m of the area of the hive and the t o t a l number of workers. This result is not surprising since pheromone produced by the queen is divided between the wax comb and the workers. Hence, the amount on any given worker decreases w i t h either an increased number of workers or an increased area Chapter 4. The Pheromone Transmission Model 111 of comb. T h e analysis leading to Equation (4.66) shows how the pheromone is divided up between the a r e a of wax and the n u m b e r of workers. T h i s effect is referred to as the dilution effect (Naumann et a l . , 1993). Note that p does not depend on N. Hence, the dependence of the average worker pheromone level on changing colony size, N, will, again, be due to changes i n the quantity 1/(1 + a). Thus, using the same calculations as i n the previous paragraph, the increase i n workers from 10,000 to 20,000 i n a 'typical' hive (see Table 4.3) will result i n a decrease i n the mean worker pheromone level by 20%. T h e motility of the workers is significant i n determining the distribution of pheromone through the colony, even if the movement of the queen is not hindered. Decreasing the motility increases the variation i n worker pheromone levels w i t h distance from the queen. T h i s is consistent w i t h the hypothesis that increased congestion of the colony restricts the movement of the workers and leads to an imbalance i n the distribution of pheromone. Restricted movement is represented in the model by a smaller value for the motility p. Decreasing p yields a slower decay for the higher spatial modes (smaller magnitude of the eigenvalues), which results i n a persistent spatial pattern to the pheromone distribution. T h e converse argument states that increasing the motility causes a faster decay of the higher spatial modes, which implies a faster spread of pheromone through the system and a more uniform distribution. In contrast to the 20% changes i n worker pheromone levels that arise from the dilution effect discussed in the previous two paragraphs, changes i n worker motility can lead to much larger changes i n pheromone levels (see Figure 4.8). A l t h o u g h the pheromone is divided in nearly equal proportions between the workers and the wax, the results indicate that the rate at which pheromone is removed from the wax by the workers is much smaller than the rate at which pheromone is absorbed into the worker cuticle. This is indicated i n the model by relation rj « 1, which implies that the rate of transfer of pheromone from the wax to the workers is very slow relative to the rate of pheromone removal from the workers. This suggests that normal pheromone levels Chapter 4. The Pheromone Transmission Model 112 in the wax will not present a significant signal to the workers. However, the pheromone level on tracks laid by the queen are significantly higher than the average, and the amount of pheromone acquired by workers moving over these queen tracks is be similar to the amount removed from messengers. E q u a t i o n (4.68) predicts that increasing the motility of the workers increases the local variance of the worker pheromone distribution. T h i s implies a trade-off between a distribution of pheromone over a large area but not evenly between the workers, versus a more even distribution over a smaller area. 4.6.2 Comparison with the Lagrangian simulation The results of the Lagrangian simulation suggest the following predictions • T h e path taken through the hive by the queen is not reflected i n the spatial pattern of the worker pheromone levels. • T h e path of the queen is reflected i n the pheromone levels on the wax. • In b o t h congested and uncongested colonies the pheromones levels of the workers are significantly higher near the queen. T h e effects of congestion can be seen i n the number of workers with high pheromone levels further from the queen. • In the congested colonies, pockets of higher worker pheromone levels may remain longer around previous positions of the queen. T h e first point can be explained using the fact that the rate of retrieval of pheromone from the wax by the workers is much slower than the rate of decay of pheromone on the workers. However, it is interesting to note the even when the queen leaves a strong trail on the wax, as was the case for the congested colony simulation shown i n Figure 4.3 and Figure 4.5, this does not appear to affect the worker pheromone distribution. T h e final Chapter 4. The Pheromone Transmission Model 113 two results were also obtained for the simplified model of Equation (4.67), as shown i n Figure 4.8. Thus, this prediction will stand so long as the exchange between workers conserves pheromone. 4.7 Discussion Colony congestion can affect the transmission of pheromone i n two ways. F i r s t , increasing the number of bees i n the colony while keeping the same production rate (k ), 0 will decrease the amount of pheromone received by each bee. T h i s is shown by the inverse relationship between steady state pheromone levels on workers (m ) and the number of s workers (JV) i n Equation (4.66). Seeley and Fell (1981) speculate that this effect should not be significant i n suppressing queen rearing, based on earlier results that congested, unpopulous colonies are likely to rear queens, yet populous, uncongested colonies are not (Seeley, 1979). M y results explain these experiments by determining how pheromone is split between the wax and the workers and showing that the average worker pheromone level would change by only 20% between the two experiments. Hence, the congestion of the colony measured as bees per unit area of u t i l i z e d comb is more significant than the absolute number of bees. T h e emphasis on utilized comb is necessary since colonies often do not use all of the available cavity space (Winston and Taylor, 1980). T h e second effect of colony congestion is the possible decrease i n bee motility w i t h increased crowding. B o t h the simulation of the full model and analysis of the linear model show that a decreased motility of the bees leads to a higher variance i n the pheromone distribution. Thus, if we assume that increased bee densities lead to decreased bee motility, then congestion will cause a significant reduction i n pheromone transmission. N a u m a n n et al (1993) conducted an experiment to determine the effect of colony congestion on pheromone transmission through the colony. Their results showed that Chapter 4. The Pheromone Transmission Model 114 the average pheromone levels on workers with detectable amounts of pheromone d i d not depend on colony size, but that the percentage of bees w i t h detectable amounts of pheromone was lower i n more populous colonies. The results of the model agree w i t h each of these points. Since the pheromone is divided between wax and workers, a large increase in worker population will produce only a small change in average worker pheromone level. Thus, i n the experiments, the average pheromone level was not noticeably different between the populous and unpopulous colonies. T h e result that a larger number of bees i n the populous colony receive no detectable levels of pheromone agree w i t h the results of the C A simulation. One interesting result of the experiments was that the pheromone was removed from the lure at a much slower rate that determined from the earlier experiments of N a u m a n n et al. (1991). In fact, the value suggested for k was 3 a full order of magnitude lower in the colony congestion experiments of N a u m a n n et al (1993) than i n the laboratory experiments of N a u m a n n et al. (1991). T h i s would reduce the magnitude of a and hence further reduce the effect of colony size on the average pheromone level per worker. N a u m a n n et al (1993) noted that no spatial pattern to the worker pheromone distribution was observed (i.e., a spatially homogeneous distribution observed) and concluded that worker motility was not important. However, their experiments d i d not have a continuous source of pheromone (i.e., ko = 0). In this case, my model predicts that the bulk of the pheromone is removed from the lure during the first hour of the experiment, during which, only one sample was taken. In order to see the spatial effects caused by congestion samples would have to be taken over a 15-30 minute time frame rather than a 24 hour time frame. Chapter 5 Discussion In this thesis I have developed models for two phenomena of self-organization i n the honey bee A. mellifera. T h e models assume that the behaviour of the individual bees depends only on local conditions and that the global structure of the group arises from the interactions between neighbouring bees. To emphasize these two aspects of the group, I have proposed the terms group structure to refer to global characteristics of the group, such as energy or pheromone distributions, and the term group dynamics to refer to the behaviour of the individuals which gives rise to the group structure. T h e goal of my research is to provide an analytical link between group dynamics and group structure. The two models presented i n this thesis contrast two broad classes of group dynamics. In the thermoregulation problem, the behaviour of the individuals depends on the local temperature a n d the individuals can influence the local temperature. T h e bees interact w i t h each other only through the temperature profile. Thus, the actions of the bees on the core are 'communicated' to the bees on the mantle by thermal diffusion. In the pheromone transmission problem, the bees interact directly by exchanging pheromone. Diffusion again plays an important role, but i n this case the diffusion of the pheromone results from the random motion of the bees. T h e similarity of the two problems is seen when one looks at the average pheromone level of the bees as a function of position. T h i s quantity is transmitted through the colony i n an fashion analogous to thermal diffusion. However, the thermoregulation model implicitly assumes that the temperature is smooth. T h a t is, the temperature difference between nearby bees is small. In contrast, 115 Chapter 5. Discussion 116 the pheromone transmission model examines the distribution of bees over both space and pheromone level. This leads to the interesting result that increasing the motility of the bees increases the local variance i n pheromone levels. 5.1 Self-organizing group dynamics M a n y of the collective activities of honeybees, and other social insects (Pasteels et al., 1987; A r o n et a l . , 1989b; Beckers et al., 1990; Deneubourg and Goss, 1989; Deneubourg et a l . , 1990; Deneubourg et al., 1991), can be referred to as Self-Organizing. T h i s term indicates that the group structure can be explained using the group dynamics, without requiring any prepattern, central control, or external organizing force. A s A s h b y (1962) points out, this term is very general and can be applied to nearly all systems. However, the most interesting feature of self-organizing groups is their ability to respond en masse to external stimuli. Some groups also show abilities to change their group structure as external conditions change. For example, the foraging patterns of many ant species may change when food sources change (Bernstein, 1975; Franks et a l . , 1991). A simulation showing how this switching of group structure may result from a slight change i n the group dynamics was studied by Watmough and Edelstein-Keshet (1995). T h i s feature may be used to further refine the definition of a self-organizing system. For example, comb construction by honeybees presents a beautiful regular pattern, but it is not clear that this should be termed as self-organizing, as the pattern does not depend on external conditions, but is invariably within a narrow construction specifications. T h i s may be an example of construction according to a prepattern. The prepattern being the similar dimensions of each worker bee. The hypothesis of self-organizing group dynamics has also been successfully applied to fish schooling. H u t h and Wissel (1992) have developed a simulation i n which the Chapter 5. 117 Discussion motion of each individual fish is determined solely from the position and velocity of its nearest neighbours. Their results show that self-organized schooling occurs for certain values of the model's parameters. Experimental observations of both individual and group behaviour for fish schooling are also available and have been compared w i t h the simulation results ( H u t h and Wissel, 1994). 5.2 Cellular automata simulations Organisms are discrete entities occupying a position i n a continuous spatial domain. Changes i n this position are usually continuous over time. Neither P D E models nor C A models are perfectly suited to this situation. P D E models, which usually i m p l y an Eulerian description, must assume that the number of individuals per unit area is large enough to work w i t h densities of individuals. This is often fine for spatial distributions; however, when one considers aspects such as pheromone level, age, or size, the number of individuals per unit volume of state space is usually small, and random fluctuations from the average density may be important. E x a m i n a t i o n of these systems must be done using a Lagrangian description. The chief difficulty with Lagrangian models is the large number of equations that must be solved. T h i s difficulty persists even though careful implementation of the C A models can reduce the number of computations dramatically. C A simulations of Eulerian descriptions are also common and can be solved much faster than their Lagrangian counterparts, allowing for more detail at m i n i m a l expense (of computer time and user patience). These simulations may model individual movement through patches or may lump similar individuals together to form 'Super-Individuals' (Scheffer et a l . , 1995). Despite the optimistic name, these simulations are Eulerian and do n o t focus on the individual (they simply use a coarse discretization of the state). Chapter 5. 5.3 Discussion 118 Future work T h e pheromone transmission model has similarities to the model for information collection and transmission by ant colonies presented by Adler and G o r d o n (1992). In their model, the ants are either aware of or ignorant of a food source and can share this i n formation by contacting other ants. M y model could be used to extend the knowledge of the food source to include the accuracy of the knowledge, thus, for example, allowing the accuracy (or possibly preference for a particular source) to diminish w i t h time. T h i s model could also be extended to collective decision making by social insects. The model used for cluster thermoregulation also has applications to wound healing and tumour growth. Each of these examples involves the growth of a spherical cluster of cells i n a medium which is being depleted locally as the cells grow (Chaplain et al., 1995). Cells move w i t h i n the cluster in response to the local conditions of the medium. This model is a three dimensional analogue of Keller and Segel's (1971a; 1971b) model for chemotaxis. The thermoregulation model does not take into account air convection through the clusters, which is probably the dominant mode of heat transfer for ambient temperatures above 15°C. Not only is this a challenging system to model, but there are no experimental results to suggest a proper starting point. The system is thermal convection through a porous medium and presents several serious challenges to modelling. First, the local porosity of the medium changes i n response to the local air temperature. T h i s change is neither uniform nor symmetric, since the bees create channels through the cluster which rapidly convect heat from the core. Second, the medium is a heat source whose local heat production increases with increasing air temperature. Finally, the medium has a free boundary w i t h unknown boundary conditions. Based on my present model, very few simplifications can be made without losing essential features of the system. Chapter 5. Discussion 119 The C A simulation developed i n this thesis for the pheromone transmission problem would provide a suitable platform for experimentation w i t h this system. The state variables of the wax would be replaced w i t h the local air temperature and velocity. The simulation could be used to examine the consequences of the various possible rules for the movement of the bees and to determine which movements allow for the opening of air channels through the swarm. The pheromone transmission model is a basis for many other systems. The code is written i n C + + , an object oriented language. This makes it easy to alter the dynamics of b o t h the individuals and the substrate. It also makes it a simple matter to study systems w i t h an inhomogeneous substrate. The pheromone simulation presented i n this thesis could be adapted to model both herding and schooling behaviour w i t h little difficulty. The results of these simulations would build on those of H u t h and Wissel (1992) and suggest approaches for mathematical models. Bibliography A d l e r , F . R . a n d G o r d o n , D . M . (1992), Information collection and spread by networks of patrolling ants, The American Naturalist 140(3):373-400 A r o n , S., Deneubourg, J . - L . , Goss, S., and Pasteels, J . M . (1989a), Functional selforganization illustrated by inter-nest traffic i n ants: the case of the Argentine ant, i n W . A l t and G . Hoffmann (eds.), Biological Motion, Proceedings, Konigswinter, (1989), V o l . 89 of Lecture Notes in Biomathematics, pp 533-547, Springer Verlag, B e r l i n A r o n , S., Pasteels, J . M . , and Deneubourg, J . - L . (1989b), Trail-laying behavior during exploratory recruitment i n the argentine ant, Iridomyrmex humilis ( M a y r ) , Biology of Behaviour 14(3):207-217 Ashby, W . R . (1962), Principles of the self-organizing system, i n H . von Foerster and J . George W . Zopf (eds.), Principles of self-organization, Pergamon Press, O x f o r d Baumgarte, J . (1972), Stabilization of constraints and integrals of motion i n dynamical systems, Comp. Math. Appl. Mech. Eng. 1:1-16 Beckers, R . , Deneubourg, J . - L . , Goss, S., and Pasteels, J . M . (1990), Collective decision making through food recruitment, Insectes Sociaux (Paris) 37(3):258-267 Berg, H . C . (1983), Random Walks in Biology, Princeton University Press, New Jersey Bernstein, R . A . (1975), Collective decision making through food recruitment: Foraging strategies of ants i n response to variable food density, Ecology 56:213-219 Butler, C . (1960a), Queen substance production by virgin queen honeybees (Apis mellifera L . ) , Proceedings of the Royal Entomological Society of London (A) 35:170-171 Butler, C . (1960b), T h e significance of queen substance i n swarming and supersedure i n honeybee (Apis mellifera L . ) colonies, Proceedings of the Royal Entomological Society 120 BIBLIOGRAPHY 121 of London (A) 35:129-132 Butler, C . and Simpson, J . (1958), T h e source of the queen substance of the honeybee (Apis mellifera L.), Proceedings of the Royal Entomological Society of London (A) 33:120-122 Camazine, S. (1991), Self-organizing pattern formation on the combs of honeybee colonies, Behavioral Ecology and Sociobiology 28:61-76 Chandrasekhar, S. (1961), Hydrodynamic and hydromagnetic stability, Oxford University Press, Oxford C h a p l a i n , M . A . J . , Giles, S. M . , Sleeman, B . D., and Jarvis, R. J . (1995), A mathematical analysis of a model for tumor angiogenesis, Journal of Mathematical Biology 33:744770 Civelekoglu, G . and Edelstein-Keshet, L. (1994), Modelling the dynamics of F-actin i n the cell, Bulletin of Mathematical Biology 56(4):587-616 Degn, H . (1972), Oscillating chemical reactions in homogeneous phase, Journal of Chemical Education 49(5):302-307 Deneubourg, J . - L . , A r o n , S., Goss, S., and Pasteels, J . M . (1990), T h e self-organizing exploratory pattern of the Argentine ant, Journal of Insect Behavior 3(2):159-168 Deneubourg, J . - L . and Goss, S. (1989), Collective patterns and decision making, Ethology, Ecology and Evolution 1:295-311 Deneubourg, J . - L . , Goss, S., Franks, N . R., Sendova-Franks, A . , Detrain, C , and Cretian, L. (1991), T h e dynamics of collective sorting; robot-like ants and ant-like robots, in J . Meyer and S. W i l s o n (eds.), Simulation of Animal Behaviour; From Animals to Animats, pp 356-365, M I T Press, Cambridge Ermentrout, B . and Edelstein-Keshet, L. (1993), Cellular automata approaches to biological modeling, Journal of Theoretical Biology 160:97-133 Fahrenholz, L., Lamprecht, I., and Schricker, B . (1989), T h e r m a l investigations of a BIBLIOGRAPHY 122 honey bee colony: thermoregulation of the hive during summer and winter and heat production of members of different bee castes, Journal of Comparative Physiology B 159:551-560 Fergusson, A . W . and Free, J . B. (1980), Queen pheromone transfer w i t h i n honeybee colonies, Physiol. Entomol. 5:359-366 F i e l d , R. J . (1972), A reaction periodic i n time and space, Journal of Chemical Education 49(5):308-311 Franks, N . R., Gomez, N . , Goss, S., and Deneubourg, J . - L . (1991), T h e blind leading the blind i n army ant raid patterns: Testing a model of self-organization (Ffymenoptera: Formicidae), Journal of Insect Behavior 4(5):583-607 Free, J . B . (1987), Pheromones of Social Bees, Comstock Publishing Assoc., Ithaca, New York G r i i n b a u m , D. (1994), Translating stochastic density-dependent individual behaviour w i t h sensory constraints to an Eulerian model of animal swarming, Journal of Mathematical Biology 33:139-161 Gueron, S. and L e v i n , S. A . (1993), Self-organization of front patterns i n large Wildebeest herds, Journal of Theoretical Biology 165:541-552 Gurney, W . S. C . and Nisbet, R. M . (1975), T h e regulation of inhomogeneous populations, Journal of Theoretical Biology 52:441-457 G u r t i n , M . E . and M a c C a m y , R. C . (1977), O n the diffusion of biological populations, Mathematical Biosciences 33:35-49 Harada, Y . , Ezoe, H . , Iwasa, Y . , M a t s u d a , H . , and Sato, K . (1995), Population persistence and spatially limited social interactions, Theoretical Population Biology 48:6591 Heinrich, B . (1981a), T h e mechanisms and energetics of honeybee swarm temperature regulation, Journal of Experimental Biology 91:25-55 BIBLIOGRAPHY 123 Heinrich, B . (1981b), Energetics of honeybee swarm thermoregulation, Science 212:565566 Heinrich, B . (1985), T h e social physiology of temperature regulation i n honeybees, i n B . Holldobler and M . Lindauer (eds.), Experimental Behavioural Ecology and Sociobiology, p p 393-406, Sinauer Associates, Inc., Sunderland Higo, H . A . , Colley, S. J . , and W i n s t o n , M . L . (1992), Effects of honeybee Apis mellifera L . queen mandibular gland pheromone on foraging and brood rearing, Canadian Entomologist 124:409-418 H u b e r m a n , B . A . and Glance, N . S. (1993), Evolutionary games and computer simulations, Proc. Natl. Acad. Sci. USA 90:7716-7718 H u t h , A . and Wissel, C . (1992), T h e simulation of the movement o f f i s h schools, Journal of Theoretical Biology 156:365-385 H u t h , A . a n d Wissel, C . (1994), T h e simulation of the movement of fish schools i n comparison w i t h experimental data, Ecological Modelling 7 5 / 7 6 : 1 3 5 - 1 4 5 Jager, E . a n d Segel, L . A . (1992), O n the distribution of dominance i n populations of social organisms, SIAM Journal of Applied Math 52(5):1442-1468 K a m m e r , A . a n d Heinrich, B . (1974), bumblebees, Metabolic rates related to muscle activity i n Journal of Experimental Biology 61:219-227 K a r l i n , S. and Taylor, H . M . (1975), A first course in stochastic processes, Academic Press, New Y o r k K a r l i n , S. and Taylor, H . M . (1981), A second course in stochastic processes, Academic Press, New Y o r k Keller, E . F . and Segel, L . A . (1971a), M o d e l for chemotaxis, Journal of Theoretical Biology 30:225-234 Keller, E . F . and Segel, L . A . (1971b), Travelling bands of chemotactic theoretical analysis, Journal of Theoretical Biology 30:235-248 bacteria: a BIBLIOGRAPHY 124 K r e i t h , F . and B o h n , M . S. (1986), Principles of Heat Transfer, Harper & R o w , New York, 4 th edition Lefelhocz, J . F . (1972), T h e color blind traffic light, Journal of Chemical Education 49(5):312-314 Lefever, R . , Nicolis, G . , and Prigogine, I. (1967), O n the occurence of oscillations around the steady state i n systems of chemical reactions far from equilibrium, Journal of Chemical Physics 47(3): 1045-1047 Lemke, M . and Lamprecht, I. (1990), A model for heat production and thermoregulation in winter clusters of honey bees using differential heat conduction equations, Journal of Theoretical Biology 142:261-273 M o r i t z , R . F . and Southwick, E . E . (1992), Bees as superorganisms: an evolutionary reality, Springer-Verlag, Berlin Murray, J . (1989), Mathematical Biology, Springer Verlag, B e r l i n Myerscough, M . (1993), A simple model for temperature regulation i n honeybee swarms, Journal of Theoretical Biology 162:381-393 Myerscough, M . (1994), T h e effects of age structure and low temperature setpoints on temperature regulation i n stationary honeybee swarms. N a u m a n n , K . , Slessor, K . N . , Prestwich, G . D . , and L a t l i , B . (1992), Intra-nest transmission of aromatic honey bee queen mandibular gland pheromone components: move- ment as a unit, Canadian Entomologist 124:917-934 N a u m a n n , K . , W i n s t o n , M . L . , and Slessor, K . N . (1993), Movement of honey bee (Apis mellifera 1.) queen mandibular gland pheromone i n populous and unpopulous colonies, Journal of Insect Behavior 6(2):211-223 N a u m a n n , K . , W i n s t o n , M . L . , Slessor, K . N . , Prestwich, G . D . , and Webster, F . X . (1991), P r o d u c t i o n and transmission of honey bee queen (Apis mellifera L . ) mandibular gland pheromone, Behavioral Ecology and Sociobiology 29:321-332 BIBLIOGRAPHY 125 Neubert, M . G . , K o t , M . , and Lewis, M . A . (1995), Dispersal and pattern formation i n a discrete time predator-prey model, Theoretical Population Biology 48:7-43 N i w a , H . - S . (1994), Self-organizing dynamic model offish schooling, Journal of Theoretical Biology 171:123-136 O k u b o , A . (1980), Diffusion and Ecological Problems: Mathematical Models, V o l . 10 of Biomathematics, Springer Verlag, Berlin O m h o l t , S. (1987), Thermoregulation i n the winter cluster of the honeybee, Apis mellifera, Journal of Theoretical Biology 128:219-231 O m h o l t , S. and L0nvik, K . (1986), Heat production i n the winter cluster of the honeybee, Apis mellifera. A theoretical study, Journal of Theoretical Biology 120:447-456 P a n k i w , T . , W i n s t o n , M . L . , and Slessor, K . N . (1994), Variation i n worker response to honey bee Apis mellifera L . queen mandibular pheromone (Hymenoptera: A p i d a e ) , Journal of Insect Behavior 7:1-15 Pasteels, J . M . , Deneubourg, J . - L . , and Goss, S. (1987), Self-organization mechanisms i n ant societies (I) : T r a i l recruitment to newly discovered food sources, i n J . M . Pasteels and J . - L . Deneubourg (eds.), From Individual to Collective Behavior in Social Insects, pp 155-175, Birhauser, Basel Prigogine, I. and Lefever, R . (1968), Symmetry-breaking instabilities i n dissipative systems. II, Journal of Chemical Physics 48(4):1695-1700 Prigogine, I. and Nicolis, G . (1967), O n symmetry-breaking instabilities i n dissipative systems, Journal of Chemical Physics 46(9):3542-3550 Reuter, H . and Breckling, B . (1994), Self-organization o f f i s h schools: an object-oriented model., Ecological Modelling 75/76:147-159 Ribbands, C . (1953), The behavior and social life of honeybees, Bee Research Association, London Sato, K . (1994), Pathogen invasion and host extinction i n lattice structured populations, BIBLIOGRAPHY 126 Journal of Mathematical Biology 32:251-268 Scheffer, M . , Baveco, J . M . , DeAngelis, D . L . , Rose, K . A . , and van Nes, E . H . (1995), Super-individuals: a simple solution for modelling large populations on an individual basis, Ecological Modelling 80:161-170 Seeley, T . D . (1978), Life history strategy of the honey bee, Apis mellifera, Oecologia 32:109-118 Seeley, T . D . (1979), Queen substance dispersal by messenger workers i n honeybee colonies, Behavioral Ecology and Sociobiology 5:391-415 Seeley, T . D . (1983), Honeybee Ecology, Princeton University Press, Princeton Seeley, T . D . and Fell, R . D . (1981), Queen substance production i n the honeybee (Apis mellifera), Journal of the Kansas Entomological Society 54(1):192-196 Seeley, T . D . and Heinrich, B . (1981), Regulation of temperature i n the nests of social insects, i n B . Heinrich (ed.), Insect Thermoregulation, pp 159-234, John W i l e y & Sons, New Y o r k Segel, L . A . and Jackson, J . L . (1972), Dissipative structure: an explanation and an ecological example, Journal of Theoretical Biology 37:545-559 Sherratt, J . (1994), Chemotaxis and Chemokinesis i n Eukaryotic cells. T h e Keller-Segel equations as an approximation to a detailed model, Bulletin of Mathematical Biology 56(1):129-146 Simpson, J . (1961), Nest climate regulation i n honey bee colonies, Science 133:1327-1333 Slessor, K . N . , K a m i n s k i , L . - A . , K i n g , G . G . S., Borden, J . H . , and W i n s t o n , M . L . (1988), Semiochemical basis of the retinue response to queen honey bees, Nature 332:354-356 Southwick, E . E . (1983), T h e honeybee cluster as a homeothermic superorganism, Comparative Biochemistry and Physiology 75A:641-645 Southwick, E . E . (1985), Bee hair structure and the effect of hair on metabolism at low BIBLIOGRAPHY temperature, 127 Journal of Apicultural Research 24:144-149 Southwick, E . E . (1988), Thermoregulation i n honeybee colonies, i n G . R . Needham, R. E . P. Jr., M . Delfmdado-Baker, and C . E . Bowman (eds.), Africanized Honeybees and Beemites, E l l i s Horwood L i m i t e d , Chichester, England Southwick, E . E . and M o r i t z , R . F . A . (1987), Social control of air ventilation i n colonies of honeybees, Apis mellifera, Journal of Insect Physiology 33(9):623-626 Southwick, E . E . and Mugaas, J . (1971), A hypothetical homeotherm: the honeybee hive, Comparative Biochemistry and Physiology 40A:935-944 Turchin, P. (1989), Population consequences of aggregative movement, Journal of Animal Ecology 58:75-100 Turing, A . (1952), T h e chemical basis of morphogenesis, Philosophical Transactions of the Royal Society of London, Series B 237:37-72 Tyson, J . J . (1976), The Belousov-Zhabotinsky Reaction, V o l . 10 of Lecture Notes in Biomathematics, Springer-Verlag, Berlin van der B l o m , J . (1992), Individual involvement i n queen-attending of worker honeybees, Insectes Sociaux (Paris) 39:237-249 W a t m o u g h , J . and Camazine, S. (1995), Self-organized thermoregulation of honeybee clusters, Journal of Theoretical Biology 176(3):391-402 W a t m o u g h , J . and Edelstein-Keshet, L . (1995), Modelling the formation of trail networks by foraging ants, Journal of Theoretical Biology 176(3) :357-372 W i n s t o n , M . L . (1987), The Biology of the Honeybee, Harvard University Press, C a m bridge W i n s t o n , M . L . , Higo, H . A . , Colley, S. J . , P a n k i w , T . , and Slessor, K . N . (1991), T h e role of queen mandibular pheromone and colony congestion i n honey bee (Apis mellifera L . ) reproductive swarming (Hymenoptera: Apidae), Journal of Insect Behavior 4(5):649-660 BIBLIOGRAPHY W i n s t o n , M . L . and Slessor, K . N . (1992), 128 T h e essence of royalty: honey bee queen pheromone, American Scientist 80:374-385 W i n s t o n , M . L . and Taylor, O . R . (1980), Factors preceeding queen rearing i n the Africanized honeybee (Apis mellifera), Insectes Sociaux 27(4):289-304 W o l f r a m , S. (1984), Cellular automata as models of complexity, Nature 311:419-424 Yates, F . E . (ed.) (1987), Self-Organizing Systems and The Emergence of Order, P l e n u m Press, New Y o r k List of Symbol Notation Description Section 2.3.1, page Background: t Time TV Number of organisms X(t) a random variables x Position vector 13 u(x, t) Density of individuals at position x 13 8 Spatial step size for random walk 12 T Temporal step size for random walk 12 J F l u x of organisms past position x 15 X Coefficient of chemotaxis 15 D(u) Diffusion coefficient or motility of organisms 15 E[A|73] Expectation of event A given event B 13 p Pheromone level of the i 16 i f (x) K e r n e l describing neighbourhood of interactions 1 12 -. 9 th at time t individual K(p ,p>) K e r n e l describing pheromone exchange l 7 Pheromone decay rate c(x, t) Chemoattractant concentration C h a p t e r 3, T h e r m o r e g u l a t i o n 16 16 16 16 .15 model: p(r, t) Cluster density at position (r,t) 26 T(r, t) Cluster temperature at position (r,t) 26 129 List of Notation 130 r radial position i n cluster 26 t time 26 J (p, T) p F l u x of bees past radial position (r, t) with density p and temperature T . . . 27 JT(P, T) T h e r m a l flux past radial position (r, t) with density p and temperature T .. 30 D(p) M o t i l i t y of bees at density p 27 X(T) Thermotactic index of bees at a Temperature T 27 X Q M a x i m u m thermotactic index at extreme temperatures 27 T h H u d d l i n g temperature: X(T ) 27 T Temperature range over which the thermotactic index increses 27 D M o t i l i t y of bees at comfortable densities 29 Pmax M a x i m u m physically possible density of cluster 29 c Heat Capacity of a single bee 30 k(p) C o n d u c t i v i t y of bees at density p 31 k C o n d u c t i v i t y of loosely packed cluster 31 &i C o n d u c t i v i t y of tightly packed cluster of bees 31 f(T) Metabolic heat production of a single bee at a temperature T 31 /is Metabolic heat production of a single bee at a temperature of 18° C 31 f M a x i m u m metabolic rate of shivering bee 32 Q(T) Metabolic heat production of a resting bee at temperature T 31 T Temperature at which bees begin to shiver 32 T Temperature at which bees shiver at m a x i m u m capacity 32 R(t) Radius of cluster at time t 26 N N u m b e r of bees i n cluster 26 h coefficient of convective heat transfer from the surface of the cluster 33 A m b i e n t temperature 33 r Q 0 m s m c T Q h = 0 List of Notation 131 PR Density at the edge of the cluster 34 x Dimensionless spatial coordinate within cluster 36 u(x, t) Dimensionless temperature at position (x,t) 36 p(x,t) Dimensionless density at position (x, t) 36 e Dimensionless measure of metabolic output 37 a Dimensionless measure of heat capacity per bee 37 3 Dimensionless measure of cluster size 37 ji(p) Dimensionless motility of the bees 37 X(u) Dimensionless thermotactic velocity of the bees 37 4>(u) Dimensionless metabolic output of a single bee 37 A Dimensionless cluster thermal conductivity 37 g(p, t) Constraint on cluster size 38 u Dimensionless ambient temperature 37 u Dimensionless temperature of cluster periphery 50 Uh Dimensionless huddle temperature 50 u Dimensionless thermotactic range 50 a R r C h a p t e r 4: P h e r o m o n e t r a n s m i s s i o n model R Dimensions of hive 72 Q T h e hive area 72 w(x,t) Pheromone level of the wax 72 x (t) Position of the queen 72 Q(t) Pheromone level of the queen 72 Xi(t) Position of worker i 72 Pi(t) Pheromone level of worker i 72 t time 72 q List of Notation a 132 variance of worker movements 75 8 The neighbourhood function 76 A A r e a of wax covered by queen 76 Ab A r e a of wax covered by worker 76 A Encounter area 78 ./V Number of bees in hive 72 2 r q e k ... kg Rates of pheromone transfer 62 0 6ij The encounter matrix 77 L fraction of pheromone ingested by bee i 79 Qi Rate of pheromone transfer to bee i 79 At Size of time steps 82 Ax Size of space steps 82 Tmove M e a n time spent moving by worker 75 T top M e a n time spent stationary by worker 75 T tay M e a n time spent i n retinue by worker 78 T M e a n time of antennating behaviour 78 T M e a n time spent i n vicinity of queen before entering her retinue 78 S Saturation pheromone level (on a single bee) for retinue behaviour 78 T Threshold pheromone level (on a single bee) for retinue behaviour 78 M M a x i m u m number of bees per grid point s s ant enter 80 u(x,p,t) worker probability density function 93 x Position i n the hive 93 p pheromone level 93 K(p, q) fraction of bees with Pi = p which enter retinues of messengers w i t h G(p, u) measure of pheromone transferred to a worker w i t h Pi — p P i = q .97 97 List of Notation 133 H dimensionless motility of workers 99 77 dimensionless measure of pheromone aquired from wax by workers 99 e dimensionless measure of pheromone transfer from worker to worker 99 T dimensionless timescale for the wax pheromone level 99 6 dimensionless measure of pheromone lost from wax to workers 99 7 dimensionless measure of pheromone transferred to wax from queen 99 a dimensionless measure of pheromone transferred from queen to workers . . . . 99 W r dimensionless timescale for the queen's pheromone level n(x, t) worker density at (x,t) 101 m(x,t) mean pheromone level of workers at (x,t) 101 v(x, t) variance of the pheromone level of workers at (x, t) 101 q 99 4>ij{x) Eigenfunctions of the Laplacian in two spatial dimensions 103 aij(t) coefficients of the eigenfunction expansion of m(x,t) 103 Xij eigenvalues 103 Q Steady state pheromone level of the queen 103 m Steady state mean worker pheromone level 105 w Steady state pheromone level of the wax 105 v Steady state variance of worker pheromone level 109 z Pheromone level relative to m 107 7 Scaling for dimensionless variable s s s s s 98
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The dynamics and structure of groups : two case studies...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The dynamics and structure of groups : two case studies of the common honeybee Watmough, James 1996
pdf
Page Metadata
Item Metadata
Title | The dynamics and structure of groups : two case studies of the common honeybee |
Creator |
Watmough, James |
Date Issued | 1996 |
Description | I study the problem of self-organized pattern formation by groups of organisms. "Self-organized" groups are coordinated by interactions between all individuals in the group rather than by a leader or hierarchical organization. The information and stimuli necessary to coordinate individual activities are communicated through the group by these interactions. I begin by describing a modelling philosophy where as many details of the individual behaviours are included in the model as are necessary for both biological and mathematical completeness, but no group behaviours are fed into the model. Two specific phenomena are addressed using this approach. The first phenomenon is honey bee cluster thermoregulation. I address the question "can cluster thermoregulation be explained using only the observed behaviours of individual bees?" and show that a temperature dependent behaviour of the bees is sufficient to produce the observed global thermoregulation of the cluster. Previous models of cluster thermoregulation are discussed in light of this modelling approach. The model also is able to make testable predictions about the density profiles of the cluster. The second phenomenon modelled is the transmission of a pheromone through a honey bee colony. The problem is of interest to beekeepers, since the pheromone is known to suppress swarming. I answer the question "How does congestion within a honey bee colony affect the transmission of a pheromone through the hive?" This question is central to a current hypothesis on the connection between the queen's ability to produce pheromones, colony size and congestion, and swarming. The cluster thermoregulation model involves a behavioural dependence on the cluster temperature. In contrast, the transmission model involves a direct exchange of a chemical pheromone between the bees. Thus, with the thermoregulation model I examine interactions mediated through an environmental state, and with the transmission model I examine direct interactions between individuals. I study these models using both Eulerian and Lagrangian models. The Lagrangian model is studied using cellular automata simulation. The Eulerian model consists of a system of partial integro-differential equations and is studied using a combination of numerical simulations and mathematical analysis. The results of the models suggest that the two interactions produce similar group properties. In the thermoregulation model, information about the global structure of the group diffuses through the cluster as thermal energy. In the transmission model, pheromone diffuses through the colony by the random motions of the bees. |
Extent | 10795434 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-03-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0079783 |
URI | http://hdl.handle.net/2429/6182 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1996-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1996-148513.pdf [ 10.3MB ]
- Metadata
- JSON: 831-1.0079783.json
- JSON-LD: 831-1.0079783-ld.json
- RDF/XML (Pretty): 831-1.0079783-rdf.xml
- RDF/JSON: 831-1.0079783-rdf.json
- Turtle: 831-1.0079783-turtle.txt
- N-Triples: 831-1.0079783-rdf-ntriples.txt
- Original Record: 831-1.0079783-source.json
- Full Text
- 831-1.0079783-fulltext.txt
- Citation
- 831-1.0079783.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0079783/manifest