The origin of traveling waves in an emperor penguin huddle

Emperor penguins breed during the Antarctic winter and have to endure temperatures as low as −50 °C and wind speeds of up to 200 km h−1. To conserve energy, they form densely packed huddles with a triangular lattice structure. Video recordings from previous studies revealed coordinated movements in regular wave-like patterns within these huddles. It is thought that these waves are triggered by individual penguins that locally disturb the huddle structure, and that the traveling wave serves to remove the lattice defects and restore order. The mechanisms that govern wave propagation are currently unknown, however. Moreover, it is unknown if the waves are always triggered by the same penguin in a huddle. Here, we present a model in which the observed wave patterns emerge from simple rules involving only the interactions between directly neighboring individuals, similar to the interaction rules found in other jammed systems, e.g. between cars in a traffic jam. Our model predicts that a traveling wave can be triggered by a forward step of any individual penguin located within a densely packed huddle. This prediction is confirmed by optical flow velocimetry of the video recordings of emperor penguins in their natural habitat.

traveling wave can be triggered by a forward step of any individual penguin located within a densely packed huddle. This prediction is confirmed by optical flow velocimetry of the video recordings of emperor penguins in their natural habitat.

Introduction
The emperor penguin (Aptenodytes forsteri) is the only vertebrate species that breeds during the severe conditions of the Antarctic winter. To withstand this harsh environment and conserve energy, emperor penguins have developed striking morphological, physiological and behavioral adaptations.
The emperor penguin is the largest living penguin species [1], which together with its compact shape, short extremities and low thermal conductance (1.3 W m −2 • C −1 ; [2]) contributes to minimize heat loss [3]. As the emperor penguins have no nest for breeding and no individual territory, they can form dense clusters of thousands of individuals, so-called huddles, which provide them with an effective protection against cold temperatures and wind [4]. In a huddle, the body surface temperature of the penguins can rise within less than 2 h to 37 • C [5]. Huddling occurs most frequently during breeding in the midst of winter. Contrary to other penguin species, only male emperors incubate their single egg, by covering it in an abdominal pouch above their feet [1]. The breeding penguins can therefore perform small, careful steps [6].
The dynamics of huddling has previously been studied by analyzing the temperature and light intensity pattern recorded with sensors attached to individual penguins [5]. However, for ethical and economic reasons, this approach can only be applied to a small number of individuals within the huddle. To investigate the dynamics of a huddle as a whole, a recent study used timelapse video recordings and reported that small density waves travel through the huddle at regular time intervals of approximately 35-55 s [7]. Gradually, this leads to large-scale reorganization and movement of the entire huddle. However, the mechanism by which these density waves travel through the huddle, and how they are triggered, is currently unknown. A recent model in which the reorganization of the huddle structure is assumed to depend on the wind exposure of the individual penguins predicts that the huddle position gradually moves leewards [8]. According to this model, the huddle movements solely originate from penguins exposed to wind at the huddle boundary, whereas the penguins within a huddle remain stationary, which is in conflict with the observations [7].
In this study, we address two questions: what are the rules and the mechanisms that govern the motion of individual penguins within a densely packed huddle during a traveling wave? And how are the traveling waves triggered? To answer these questions, we developed a simple but biologically plausible model, based on a coupled system of differential equations for the positions and velocities of the individuals. In this model, we consider that the movement in a huddle does not arise at some higher organizational or hierarchical level, but is governed only by the interactions between neighboring penguins [9]. We further assume that all of the penguins are equal, that is, we do not impose a social or hierarchical structure to the huddle. Our approach therefore resembles other many-particle models used to describe collective behavior [10][11][12] in flocks of birds [13][14][15], schools of fish [16,17] or traffic jams [18,19].

Materials and methods
To study the positional reorganization processes in a penguin huddle, we observed two different emperor penguin colonies at Pointe Géologie near the French research base Dumont d'Urville, and at Atka Bay near the German research station Neumayer III.
Recordings at the Atka Bay colony were made from an elevated (12 m), distant (115 m) position, using a Canon D400 camera with a 300 mm lens. High-resolution time lapse images were recorded every 1.3 s for a total of 4 h and analyzed off-line to detect and track penguin positions (appendix B and [7]). The air temperatures during the recordings (August 2008) varied from −33 to −43 • C, with wind speeds around 14 m s −1 .
The recordings at Pointe Géologie were made in June and July 2005. The video recordings were stored on miniDV tapes with a spatial resolution of 720 × 576 pixel at 25 frames s −1 . These recordings were of insufficient resolution for tracking of single penguins, and thus were analyzed using optical flow velocimetry from images that were 5 s apart. Each image was computed as the median of ten consecutive frames (0.4 s) to reduce the image noise and the digitization errors. The optical flow field was computed using a combined Lucas-Kanade and Horn-Schmuck algorithm [20][21][22].

Model
The spontaneous formation of regular particle clusters is a well-known phenomenon in manyparticle systems. For example, spherical colloidal particles with Lennard-Jones-like interaction potentials, characterized by attractive forces at large distances and steric repulsion forces at small distances, can relax into a low-energy configuration in the form of a densely packed cluster with a triangular lattice structure [23]. Systems of non-spherical particles, possibly combined with more complicated interaction potentials, such as dipole-dipole interactions, can form clusters where not only the particle centers are arranged periodically but also the axes of the particles are aligned [24]. Such systems may be used as simple models for the formation of clusters in more complex systems, such as huddles of emperor penguins.
In the case of colloidal or other inert systems, the particles are passively pulled along the vector sum of all of the interaction forces exerted by their neighbors, including thermal forces from the solvent. However, animals differ from these inert many-particle systems in that they are self-driven. They actively perceive the positions of their neighbors, process this information and then activate their locomotory system to move into an appropriate direction. Therefore, the movements of self-driven agents are more complex than those of inert particles. Moreover, since the animals use internal energy reservoirs to generate the locomotive forces, describing the motion of the animals as a many-particle system with suitably chosen interaction potentials would violate the energy conservation law [12].
Self-driven agent models have been successfully applied to understand human traffic dynamics [19]. A key concept in traffic models is the so-called desired velocity v des i of each agent i, which can depend in an arbitrary way on the internal state of i and on the positions and velocities of its neighbors [18]: The actual velocity v i (t) of an agent does not immediately follow the desired velocity v des i , but converges toward it with a characteristic relaxation time constant τ . This time lag accounts for the information processing time of the agent, the finite response times of its locomotory systems and the inertia effects. It is essential for the formation of stop-and-go waves in traffic jams [19] and, as we will demonstrate below, also for traveling waves in penguin huddles.

Wave model
To model the penguins as self-driven agents, we have adopted and modified the above equations of motion. We assume as a starting condition a triangular lattice structure for the penguin positions. If they are in this ideal, densely packed configuration, no spontaneous movements take place, except when a penguin triggers a wave. Furthermore, we assume that all of the penguins are oriented at all times in the same direction (x) and never rotate. Although in the video recordings curved or circular huddle structures accompanied with correlated rotations of individual penguins were observed, the minimal model presented here only considers the linear wave motions in the x-direction without rotations. Consequently, we consider only the velocities in the x-direction and leave the y-positions of the penguins constant.
This simplification allows for a direct application of the one-dimensional equation of motion used in the traffic systems [18], which describes the adjustment of the current (x-) velocity v i (t) to the desired velocity v des i of an agent i with relaxation time τ : Here, the desired velocity can assume one of Here, (x) is the Heaviside step function, and v step as well as d th are the free model parameters.
For defining x des i , we assume that there exists an ideal distance d 0 between neighboring penguins. In a perturbed configuration, the distances d i j to the neighbors j cannot all be d 0 . Thus, the desired position x des i is the position with the minimal overall deviation from d 0 , i.e. at the minimum of j (d i j − d 0 ) 2 . Unlike the traffic systems where only the car in front matters, in the case of the penguins, all of the six direct neighbors (numerically determined by the Delaunay triangulation [25]) are taken into account, except at the borders of the huddle, where the numbers of direct neighbors can be smaller (see appendix C for details).
The ideal distance d 0 between two penguins i and j can be expressed as the sum of two radii d i j 0 = r i + r j , which describe the desired free space of each animal individually. Note that the radius r i need not correspond to the physical radius of the animal. Generally, we assign the same radius to every penguin, however, we also explore a distribution of penguin radii with a Gaussian distribution centered around d lattice /2 and a standard deviation σ of 0.1 d lattice , which leads to an imperfect triangular structure.
The configuration of the huddle remains unchanged unless the triangular structure of the huddle becomes locally disturbed by a supercritical amount. Any penguin within the huddle can cause such a disturbance by performing a step forward, i.e. setting its desired velocity v des i to v step for a duration t 0 . We assume a constant rate R tr for such triggering events for any penguin, so that the triggering events are a Poisson point process with an average waiting time T tr = 1/(N R tr ), with N equal to the number of individuals of the huddle.
Thus, the parameters of the model are τ , d th , v step , d 0 , t 0 , R tr and N . However, the number of parameters can be reduced. The relaxation time τ is chosen to be τ = d th v step , so that a decelerating penguin stops at x des . Slight deviations from this ideal relaxation time τ lead to no fundamental change of the system dynamics but only to a small constant distortion of the huddles structure. Furthermore, the duration of a step t 0 as well as the lattice distance d 0 can be set to unity. As we will demonstrate below, the huddle size N has no effect on the wave speed and the overall dynamical properties of the system. Moreover, variation of the trigger rate R tr does not show any qualitatively new effects as long as 1/R tr for an individual penguin is larger than the time t 0 (at smaller rates the huddle would be in constant motion). Therefore, only d th and v step remain as effective parameters of the model.

Dynamical properties of the system
Any penguin within the huddle occasionally performs a step. This locally disturbs the triangular configuration of the huddle and triggers each of the neighboring penguins to also perform-after a small time-delay-a single step (figure 1). This delay depends on the actual distances between the penguins, the threshold distance and the speed of the movements. The resulting cascade of the reorganizations spreads out from the initiating penguin as a wave with a fixed speed and a rectangular-shaped front ( figure 1(d)). Because all of the motions are fueled by the internal energy reservoirs of the individual penguins, no attenuation of the wave speed and the amplitude occurs. After the wave has traveled through the whole huddle, the conformation of the huddle is the same as before the wave, but the whole huddle has now moved one step forward. This is in qualitative agreement with the video recordings [7], where a wave leads to no obvious distortion in the huddle conformation but moves the entire huddle forward. Both the video recordings and the simulations display intermittent dynamics, where the resting phases are interrupted by short periods of movements (figures 2(a) and (b)). It is not known at present, however, what causes an individual penguin to initiate a wave, or what keeps the penguins in the huddle to remain motionless for more than 30 s in between periods of movements. We next explore the behavior in the huddle when a traveling wave was triggered simultaneously or in close succession at two distant sites. As the frequency 1/T = N R tr of the waves increases linearly with huddle size N , the resulting global velocity of the huddle could be naively expected to also grow linearly with the huddle size if each triggering event would cause the huddle to advance by one step. For larger huddles, the intermittent behavior would eventually vanish and the penguins would be in constant motion, a feature not observed in the video recordings [7] or in our simulations. Instead, we observe in our model that the fronts of two simultaneous waves do not pass through each other but merge. After the merged wave has traveled through the entire huddle, each individual penguin has performed only a single step, despite multiple overlapping triggering events. This interesting emergent effect is possible because wave propagation in a huddle is a nonlinear phenomenon so that the superposition principle (mutual cancelation or amplification of two waves) does not hold.

Wave speed
The wave speed through the huddle can be derived analytically (see appendix A) as the distance between two penguins, divided by the time t (d th ) a penguin needs to move the threshold distance W 0 is the Lambert W -function to take the acceleration of the forward step into account (see appendix A). Note that the acceleration and the threshold distance are not independent but are coupled, as explained above. The maximum velocity v step , the wave speed v wave and the distance between the penguins d 0 can be experimentally obtained (figure 2(a)), and from these values the threshold distance can be calculated according to equation (4). From the video recordings published in [7], we estimate a wave speed of around 12 cm s −1 , a maximum velocity during the steps of 1-2 cm s −1 and a distance d 0 of 34 cm. From these numbers we estimate a threshold distance d th of around 2 cm. Independently, a rough estimate of d th can also be directly obtained from the penguin trajectories as the amplitude of the distance fluctuations between two neighboring penguins during a step (see appendix B). From this, we find a value of d th = 1.84 ± 0.99 cm, which is similar to the value calculated from equation (4). Interestingly, this distance is comparable to twice the thickness of the compressive feather layer of around 1.2 cm [5]. This suggests that the penguins touch each other only slightly when standing in a huddle, without compressing the feather layer so as to maximize the huddle density without compromising their own insulation. More importantly, the analytical and the numerical results for the wave speed as a function of the threshold distance are nearly identical ( figure 3).

Features of the model
The traveling waves have been reported to have three main effects: global motion of the huddle, increasing the density of the huddle and allowing separate huddles to merge [7]. In the following, we explore whether our model also shows these features.

Maximizing density
High density requires a high degree of order in the huddle, ideally a strictly triangular lattice structure. Here, we test numerically how the huddle structure recovers from a disordered configuration. We distort a perfect triangular lattice with a lattice constant of d 0 by randomly moving the individual penguins by a Gaussian distributed distance d in the x-direction, which is the direction the penguins are facing. We find that the disorder in the huddle, expressed as the standard deviation σ = stdev(d i j ), disappears quickly within t < 2t 0 , where t 0 is the typical duration of a single step. This healing or annealing process takes place without the initiation of a traveling wave, and is driven only by the lattice disorder that causes the penguins to move according to equation (3). The residual degree of the disorder remains below the threshold distance ( figure 4). Imperfections of σ > 0.25, however, do not heal spontaneously in our simulation for two reasons. Firstly, the rearrangements are so large that before some of the penguins have had a chance to find their ideal resting position, they have collided with another penguin and thus have triggered another wave of rearrangements. Secondly, the ideal positions, determined from the distance to the surrounding penguins, move at the same time and with a similar speed as the movement of each individual penguin. Thus, the huddle never comes to a rest. In real penguin huddles, we occasionally observe a complete break-up of the ordered huddle structure. However, it is at present not clear if this phenomenon is triggered by the large and ongoing movements of one or several penguins in the huddle.

Merging of the huddles
In order to investigate how the model behaves when huddles merge, we start with two separate huddles some distance apart, both facing the same direction. If the waves are only initiated in the rear huddle, it approaches the huddle in front with every passing wave. Eventually, the rear huddle has covered the distance to the huddle in front, and as the two huddles touch, they merge seamlessly. All of the subsequent waves now pass through the united huddle (figure 5).

Initiation of the waves
In the simulations, the penguins move only forward in the x-direction. This is similar to the situation in a traffic jam. Indeed, the reorganization waves in the penguin huddles closely resemble the stop-and-go waves and the intermittent behavior in jammed traffic. There are some fundamental differences, however: while in traffic jams the waves always start at the front of the queue and travel upstream from there, the waves in a penguin huddle can travel in both directions and also spread sideways in a rectangular pattern. Moreover, the waves in a penguin huddle can be triggered from any penguin regardless of its position. In particular, there is no need for a leading 'pacemaker' penguin in our model.
To test whether these results recapitulate the experimental findings, we perform an optical flow analysis of the video recordings of the penguin huddles. For this, we select a circular-shaped huddle in which the penguins face and move in the counter-clockwise direction around a stationary center ( figure 6(a)). In this case, the lattice structure exhibits dislocation boundaries with a spiral pattern around the center, however over a shorter range, an approximately triangular order prevails. The optical flow is averaged in the y-direction (columnwise) over the lower half of the huddle to collapse the flow information onto a line (in the x-direction). The flow in the x-direction is then plotted versus time in a kymograph representation (figure 6(c)). Waves traveling from the left to the right side of the lower half of the huddle (in the direction of the penguin movements) are visible in the kymograph as lines with a negative slope. In the analogy, waves traveling from the right to the left against the direction of the penguin movements are visible in the kymograph as lines with a positive slope. The occurrence of the positive and the negative slopes shows that the waves can travel forward and backwards through the huddle ( figure 6(c)). Note that due to the larger lattice disorder in the case of a round huddle, the traveling waves may not always spread through the entire huddle but can die out over shorter distances (figure 6(c)1). The optical flow analysis also reveals that the waves start at different points in the huddle (corresponding to different ypositions in the kymograph). This leads to the conclusion that indeed the traveling waves can be triggered from different positions and therefore by different penguins in the huddle, which confirms the prediction of the model.

Conclusion
We present a simple model for the dynamical reorganization of emperor penguin huddles, based on an extension of the established one-dimensional models for traffic systems. Our model exhibits similarities with jammed traffic systems such as the stop-and-go waves and intermittent behavior, but also differences, such as that the waves can travel in all directions through a penguin huddle, which is confirmed here from the video recordings of the huddling penguins.
The model considers only the interactions between neighboring penguins, whereby a penguin moves to achieve or restore an optimal distance to the neighbors if this distance exceeds or falls below a threshold. Because each penguin is caged by surrounding penguins, a movement necessarily causes a propagation of the disturbances from the distance optimum and therefore leads to a traveling wave. The speed of the traveling wave depends on the optimal distance and the walking speed of the individual penguin, which both can be measured, and on the threshold distance, which can be inferred from the model. From our data, we estimate a threshold distance of 2 cm. We independently confirm this value from the fluctuations in the distances between the neighboring penguins during a traveling wave. Furthermore, the model confirms that the traveling waves help to achieve an optimally dense huddle structure, and describes how huddles can merge. More importantly, the model predicts that a traveling wave can be triggered from any penguin, regardless of its position in the huddle. This prediction is confirmed by an optical flow analysis of the video recordings from the huddling emperor penguins. threshold distance has been exceeded. This is in analogy to the saltatory propagation of the action potentials along the myelinated axons from one node of the Ranvier to the next node. In figure A.1 on the left side, a configuration of an undisturbed hexagonal packing is depicted, with a wave arriving from the left that has reached the rear penguins (cyan). These penguins now start moving to the right by a distance b, at which point they have pushed the desired position of the center penguin (red) beyond a critical distance d th . This situation is depicted in figure A.1 right. At this instance, the wave front jumps to the center penguin, which in turn starts moving to the right, and so forth. The effective wave speed is thus the time t (b) needed for the rear penguins (cyan) to move to the distance b, divided by the total distance a = √ 3 2 d 0 that the wave front has traveled during that time. Note that the rear penguins (cyan) continue to move beyond position b, but this has no further effect on the wave front.
To calculate the distance b, we first compute the positional error E(x i ) as the sum of the squared neighbor distance errors The desired position is the position of the minimum of E(x i ) We consider the origin of the coordinate system to be at the position of the central penguin.  (mean ± std) pixels, corresponding to 1.84 ± 0.99 cm, should be regarded only as a rough estimate for this parameter.

Appendix C. Assigning neighbors
As discussed above in section 3.1, the choice of neighbors is important for calculating the desired position x des . In the model, each penguin is assigned six direct neighbors as determined by the Delaunay triangulation [25]. However, care has to be taken with the penguins at the border of the huddle. Here, a simple Delaunay triangulation would link the penguins that are in no 'direct' neighborhood with each other (figure C.1(a)). Instead, we perform a Delaunay triangulation after including additional points (virtual penguins) outside the border of the huddle. To construct these points, a circle with radius 2r i is drawn around every border penguins i, and the points where two such circles intersect are selected. If a point falls within the desired region of an existing penguin, i.e. it is less than r j away from a penguin j, it is rejected ( figure C.1(b)). As the huddles can merge and thus the boundary of a huddle can change, these additional points are recalculated for every frame.