Mapping the stereotyped behaviour of freely-moving fruit flies

T o b e publishe d in The Journal of the R oyal So ciety Interfac e (2014) Mapping the stereot yp ed b eha viour of freely-mo ving fruit flies Gordon J. Berman, 1 Daniel M. Choi, 2 William Bialek, 1 and Josh ua W. Shaevitz 1 , ∗ 1 Joseph Henry L ab oratories of Physics and L ewis-Sigler Institute for Inte gr ative Genomics 2 Dep artment of Mole cular Biolo gy Princ eton University Princ eton, NJ 08544 (Dated: No v ember 11, 2018) A frequen t assumption in b eha vioural science is that most of an animal’s activities can b e described in terms of a small set of stereotyped motifs. Here we introduce a metho d for mapping an animal’s actions, relying only upon the underlying structure of p ostural mo v emen t data to organise and classify behaviours. Applying this metho d to the ground-based b ehaviour of the fruit fly , Dr osophila melano gaster , w e find that flies perform stereotyped actions roughly 50% of the time, disco vering ov er 100 distinguishable, stereot yp ed b ehavioural states. These include multiple mo des of lo comotion and grooming. W e use the resulting measurements as the basis for iden tifying subtle sex-sp ecific b eha vioural differences and revealing the low-dimensional nature of animal motions. I. INTR ODUCTION The concept of stereot yp y–that an organism’s be- ha viours can b e decomp osed into discrete, repro ducible elemen ts–has influenced the study of ethology , b e- ha vioural genetics, and neuroscience for decades [1, 2]. Animals p ossess the ability to mov e in a v ast contin uum of w ays, theoretically constrained only b y the biome- c hanical limits of their own morphology . Despite this, the range of behavioural actions typically p erformed by an animal is thought to b e muc h smaller, constructed largely of stereotyped actions that are consistent across time, individuals, and, in some cases, ev en sp ecies [3, 4]. A discrete b ehavioural rep ertoire can p otentially arise via a n um b er of mec hanisms, including mechanical limits of gait con trol, habit formation, and selectiv e pressure to generate robust or optimal actions. In many instances, the searc h for an individual behavioural neural circuit or gene b egins with the assumption that a particular ac- tion of in terest is stereot yp ed across time and individuals [5, 6]. Despite the cen trality of this concept, with few ex- ceptions [7–11], the existence of stereotyp y has not b een prob ed exp erimentally . This is largely due to the lack of a comprehensive and compelling mathematical frame- w ork for b eha vioural analysis. Here, we introduce a new metho d for quan tifying postural dynamics that retains an animal’s full behavioural complexity , using the fruit fly Dr osophila melano gaster as a mo del organism to dis- co v er and map stereot yp ed motions. Most prior metho ds for quantifying animal b ehaviour lie in one of t wo regimes. One of these is the use of coarse metrics such as a gross activity level (e.g. mean veloc- it y or n umber of times the organism crosses a barrier) or counting the relative frequencies of particular even ts engrained into the exp erimental set-up (e.g. turning left or right in a maze). While this approac h allows for high- throughput analysis of v arious organisms, strains, and ∗ Corresponding author: shaevitz@princeton.edu sp ecies, only the most gross asp ects of b ehaviour can b e ascertained, p otentially ov erlo oking the often sub- tle effects of the manipulations of in terest that are only apparen t at a finer descriptive level. The other com- mon approac h for b ehavioural quantification is to use a set of user-defined b ehavioural categories. These cate- gories, such as walking, grooming, or fighting, are co di- fied heuristically and scored either by hand or, more re- cen tly , via sup ervised machine-learning techniques [12– 16]. While the latter approac h allo ws for higher through- put and more consis ten t lab eling, it remains prone to h uman bias and an throp omorphism and often precludes ob jectiv e comparisons b etw een data sets due to the re- liance on sub jective definitions of b ehaviour. F urther- more, these analyses assume, a priori , that stereotyped classes of behaviour exist without first showing, from the data, that an organism’s actions can b e meaning- fully categorised in a discrete manner. Ideally , a b ehavioural description should manifest it- self directly from the data, based up on clearly-stated assumptions, each with testable consequences. The ba- sis of our approac h is to view b ehaviour as a tra jec- tory through a high-dimensional space of p ostural dy- namics. In this space, discrete b ehaviours correspond to epo chs in which the tra jectory exhibits pauses, cor- resp onding to a temp orally-extended b out of a particu- lar set of motions. Epo chs that pause near particular, rep eatable p ositions represent stereotyped b ehaviours. Moreo v er, momen ts in time in which the tra jectory is not stationary , but instead mov es rapidly , corresp ond to non-stereot yp ed actions. In this pap er, w e construct a b ehavioural space for freely-mo ving fruit flies. W e observ e that the flies ex- hibit appro ximately 100 stereotyped b ehaviours that are in tersp ersed with frequen t b outs of non-stereotyped b e- ha viours. These stereotyped b ehaviours manifest them- selv es as distinguishable p eaks in the b ehavioural space and corresp ond to recognizably distinct b ehaviours such as walking, running, head grooming, wing grooming, etc. Using this framework, we b egin to address biological questions ab out the underlying postural dynamics that generate behaviour, op ening the door for a wide range 2 FIG. 1. Schematic of the imaging apparatus. of other inquiries into the dynamics, neurobiology , and ev olution of b ehaviour. I I. EXPERIMENTS W e prob ed the sp ontaneous b ehaviours of ground- based flies ( Dr osophila melano gaster ) in a largely fea- tureless circular arena (Fig 1). Under these condi- tions, flies displa y a multitude of complex, non-aerial b eha viours such as lo comotion and gro oming, t ypically in v olving m ultiple parts of their b o dies. T o capture dy- namic rearrangements of the fly’s p osture, w e recorded video of individual b ehaving animals with sufficient spa- tiotemp oral resolution to resolv e moving b o dy parts suc h as the legs, wings, and prob oscis. W e designed our arena based on previous w ork which sho w ed that a thin cham b er with gently sloping sides pre- v en ts flies from flying, jumping, and climbing the walls [17]. T o keep the flies in the fo cal plane of our camera, w e in v erted the previous design. Our arena consists of a custom-made v acuum-formed, clear PETG plastic dome 100mm in diameter and 2mm in height with sloping sides at the edge clamp ed to a flat glass plate. The edges of the plastic co ver are slop ed to prev ent the flies from b eing o ccluded and to limit their ability to clim b upside-do wn on the co ver. The underside of the dome is coated with a rep ellent silane comp ound (heptane and 1,7-dichloro- 1,1,3,3,5,5,7,7-o ctameth ylte-trasilo xane) to prev en t the flies from adhering to the surface. In practice, we find that this set-up results in no b outs of upside-down w alk- ing. Ov er the course of these exp eriments, we studied the b eha viour of 59 male and 51 female D. melano gaster (Oregon-R strain). Eac h animal was imaged using a high-sp eed camera (100 Hz, 1088 × 1088 pixels). A prop ortional-integral-deriv ative (PID) feedback algo- rithm is used to keep the mo ving fly inside the camera frame by controlling the p osition of the X-Y stage based on the camera image in real time. In eac h frame we fo cus our analysis on a 200 × 200 pixel square contain- ing the fly . W e imaged each of the flies for one hour, yielding 3 . 6 × 10 5 mo vie frames p er individual, or ap- pro ximately 4 × 10 7 frames in total. All asp ects of the instrumen tation are controlled by a single computer us- ing a custom-written LabView graphical user in terface. Eac h of these flies was isolated within 4 hours of eclo- sion and imaging o ccurred 1-14 days after that. Flies w ere placed in to the arena via aspiration and were sub- sequen tly allow ed 5 minutes for adaptation b efore data collection (Fig S1). All recording o ccurred b etw een the hours of 9:00 AM and 1:00 PM, thus reducing the ef- fect of circadian rhythms, and the temp erature during all recordings was 25 o ± 1 o C . I I I. BEHA VIOURAL ANAL YSIS The general framework of our analysis is describ ed in Fig 2. Images are first segmen ted and registered in order to isolate the fly from the background and en- force translational and rotational inv ariance. After this, they are decomp osed into postural time series and con- v erted in to wa velet spectrograms, thus creating a spatio- temp oral represen tation for the fl y’s dynamics within the images. These sp ectrograms are used to construct sp ec- tral feature vectors that we em b ed into t w o dimensions using t-Distributed Sto chastic Neighbor Embedding [18]. Lastly , we estimate the probabilit y distribution ov er this t w o dimensional space and iden tify resolv able p eaks in the distribution. W e confirm that sustained pauses near these p eaks corresp ond to discrete b ehavioural states. A. Image segmentation and registration Giv en a sequence of images, w e wish to build a spatio- temp oral representation for the fly’s p ostural dynamics. W e start by isolating the fly within each frame, follow ed b y rotational and translational registration to pro duce a sequence of images in the co ordinate frame of the insect. Details of these pro cedures are listed in App endix A. In brief, we apply Canny’s method for edge detection [19], morphological dilation, and erosion to create a binary mask for the fly . After applying this mask, we rotation- ally align the images via polar cross-correlation with a template image, similar to previously developed methods [20 – 22]. W e then w e use a sub-pixel cross-correlation to translationally align the images [23]. Lastly , every image is re-sized so that, on a verage, eac h fly’s b o dy co vers the same n umber of pixels. An example segmentation and alignmen t is shown in Supplemen tary Movie S1. B. P ostural decomp osition As the fly b o dy is made up of relativ ely inflexible seg- men ts connected b y mobile joints, the num b er of p ostu- 3 FIG. 2. Ov erview of the data analysis pip eline. Ra w images of the Dr osophila melano gaster are segmented from the back- ground, rescaled to a reference size, and then aligned, creating a stack of images in the co-moving and co-rotating frame of the fly . These images are then decomp osed via PCA into a relatively lo w-dimensional set of time series. A Morlet wa velet transform is subsequen tly applied to these time series, creating a sp ectrogram for eac h p ostural mode separately . After nor- malization, eac h p oin t in time is mapp ed into a 2-dimensional plane via t-SNE [18]. Lastly , a watershed transform is applied to a gaussian-smo othed density ov er these p oints, isolating individual p eaks from one another. ral degrees of freedom is relatively small when compared to the 40,000 pixels in each image. Accordingly , a natu- ral represen tation for the fly’s posture would b e to enu- merate the relativ e angles of each of the fly’s app endages as a function of time [24 – 26]. Extracting these v ariables directly from the images, how ever, is prohibitiv ely diffi- cult due to occlusions and the complex fly lim b and wing geometry . As an alternative strategy , we find that nearly all of the v ariance in the 4 × 10 4 pixel images can b e explained b y pro jecting the observed pixel v alues onto a Euclidean space of just 50 dimensions. W e apply Principal Compo- nen t Analysis (PCA) to radon transforms of the images. PCA is a frequently-used metho d for con verting a set of correlated v ariables into a set of v alues of linearly un- correlated eigenmodes. Results from this analysis can be describ ed as the space spanned b y the eigen vectors of the data co v ariance matrix, C , corresp onding to the largest m eigen v alues out of the total laten t dimensionalit y of the data. While, in general, there is no rigorous man- ner to choose m , here, we will keep all mo des containing correlations larger than the finite sampling error within our data set. According to this heuristic, w e set m = 50 (Fig 3(c)), a n um b er of mo des explaining approximately 93% of the observ ed v ariation (Fig 3(d)). Details of this computation can b e found in App endix B. W e refer to these directions of correlated v ariation as p ostural mo des. As seen in Fig 3(b), these mo des are fly-lik e in app earance, but do not lend themselv es to in tuitiv e in terpretation. How ev er, pro jecting individual images onto these axes, we can conv ert a mo vie of fly b eha viour into a 50-dimensional time series, Y ≡ { y 1 ( t ) , y 2 ( t ) , · · · , y 50 ( t ) } , (1) as exemplified in Fig 3(e). C. Sp ectrogram generation The instantaneous v alues of the p ostural modes do not pro vide a complete description of b ehaviour, as our defi- nition of stereotyp y is inheren tly dynamical. Previously published studies of b ehaviour hav e searched for motifs – rep eated subsequences of finite length – within a b e- ha vioural time series [11, 27]. How ever, this paradigm is often confounded b y problems of temp oral alignment and relativ e phasing betw een the component time series. Ad- ditionally , certain b eha viours (for example, wing gro om- ing in Dr osophila ) in volv e m ultiple appendages mo ving at different time scales, thus complicating the choice of motif length. As an alternative to this approac h, we use a sp ectro- gram representation for the p ostural dynamics, measur- ing the p ow er at frequency f asso ciated with each p os- tural mo de, y k ( t ), in a window surrounding a moment in time, S (k , f ; t ). More sp ecifically , we compute the am- plitudes of the Morlet contin uous wa velet transform for eac h p ostural mo de [28]. Although similar to a F ourier sp ectrogram, wa velets possess a multi-resolution time- frequency trade-off, allo wing for a more complete de- scription of p ostural dynamics o ccurring at several time scales [29]. In particular, the Morlet wa velet is adept at isolating c hirps of p erio dic motion, similar to what w e observ e in our data set. By measuring only the ampli- tudes of the transform, we eliminate the need for precise temp oral alignment that is required in any motif-based analysis. Details of these calculations are shown in Ap- p endix C, and an example spectrogram is display ed in Fig 3(f ). F or the results presen ted here, w e look at 25 frequency channels, dyadically spaced b etw een 1 Hz and 50 Hz, the larger of which b eing the Nyquist frequency of the system. D. Spatial embedding S (k , f ; t ) is comprised of 25 frequency channels for eac h of the 50 eigenmo des, making each p oint in time represen ted b y a 1,250-dimensional feature vector en- co ding the p ostural dynamics. As correlations, often strong, exist b etw een the v arious mo de-frequency chan- nels, we exp ect that the dimensionality of the manifold con taining the observed v alues of S (k , f ; t ) should be v astly smaller. As such, we would like to find a low- 4 FIG. 3. Generation of sp ectral feature v ectors. (a) Raw image of a fly in the arena. (b) Pictorial representation of the first 5 p ostural mo des, ˆ x 1 − 5 , after in verse Radon transform. Black and white regions represent highlighted areas of each mo de (with opp osite sign). (c) First 1,000 eigenv alues of the data matrix (black) and shuffled data (red). (d) F raction of cumulativ e v ariation explained as a function of num b er of modes included. (e) Typical time series of the pro jection along p ostural mo de 6 and (f ) its corresp onding wa velet transform. dimensional represen tation that captures the important features of the data set. Our strategy for dimensional reduction of the feature v ectors is to construct a space, B , suc h that tra jecto- ries within it pause near a repeatable position whenev er a particular stereotyped b ehaviour is observed. This means that our em bedding should minimise any lo cal distortions. How ev er, we do not require preserv ation of structure on longer length scales. Hence, we chose an em b edding that reduces dimensionalit y b y altering the distances b etw een more distant p oin ts on the manifold. Most common dimensionality reduction metho ds, in- cluding PCA, Multi-dimensional scaling, and Isomap, do precisely the opposite, sacrificing lo cal verisimilitude in service of larger-scale accuracy [30 – 32]. One metho d that do es p ossess this prop erty is t-Distributed Sto chas- tic Neighbor Em b edding (t-SNE) [18]. Like other em- b edding algorithms, t-SNE aims to tak e data from a high-dimensional space and embed it into a space of m uc h smaller dimensionality , preserving some set of in- v ariants as b est as p ossible. F or t-SNE, the conserved in v ariants are related to the Marko v transition proba- bilities if a random walk is p erformed on the data set. Sp ecifically , w e define the transition probability from time p oint t i to time p oint t j , p j | i , to be proportional to a Gaussian k ernel of the distance (as of y et, unde- fined) b etw e en them: p j | i = exp  − d ( t i , t j ) 2 / 2 σ 2 i  P k 6 = i exp  − d ( t i , t k ) 2 / 2 σ 2 i  . (2) All self-transitions (i.e. p i | i ) are assumed to be zero. Eac h of the σ i are set such that all p oin ts hav e the same transition en trop y , H i = P j p j | i log p j | i = 5. This can b e interpreted as restricting transitions to roughly 32 neigh b ors. The t-SNE algorithm then embeds the data points in the smaller space while keeping the new set of transition probabilities, q j | i , as similar to the p j | i as possible. The q j | i are defined similarly to the larger-space transition probabilities, but are now, for tec hnical reasons, prop or- tional to a Cauc hy (or Student-t) k ernel of the p oints’ Euclidean distances in the em b edded space. This algo- rithm results in an em b edding that minimises lo cal dis- tortions [18]. If p j | i is initially v ery small or zero, it will place little to no constraint on the relative p ositions of the t w o points, but if the original transition probabilit y is large, it will factor significantly into the cost function. This metho d’s primary drawbac k, how ev er, is its p o or memory complexity scaling ( ∝ N 2 ). T o incorporate our en tire data set in to the em bedding, we use an imp ortance sampling tec hnique to select a training set of 35,000 data p oin ts, build the space from these data, and then re- em b ed the remaining p oin ts into the space as b est as p ossible (see App endix D for implementation details). Lastly , we need to define a distance function, d ( t i , t j ), b et w een the feature vectors. W e desire this function to accurately measure ho w different the shap es of t wo mo de-frequency sp ectra are, ignoring the ov erall m ul- tiplicativ e scaling that o ccurs at the b eginning and the end of b ehavioural b outs due to the finite nature of the wa velet transform. Simply measuring the Eu- 5 clidean norm betw een t wo spectra will b e greatly af- fected by such amplitude mo dulations. Ho w ev er, b e- cause S (k , f ; t ) is comp osed of a set of wa velet ampli- tudes, it m ust therefore be p ositive semi-definite. As suc h, if we define ˆ S ( k , f ; t ) ≡ S ( k , f ; t ) P k 0 ,f 0 S ( k 0 , f 0 ; t ) , (3) then we can treat this normalised feature vector as a probabilit y distribution ov er all mo de-frequency chan- nels at a given p oint in time. Hence, a reasonable dis- tance function is the Kullback-Leibler (KL) divergence [33] b etw een tw o feature v ectors: d ( t 1 , t 2 ) = D K L ( t 1 || t 2 ) ≡ X f ,k ˆ S ( k , f ; t 1 ) log 2 " ˆ S ( k , f ; t 1 ) ˆ S ( k , f ; t 2 ) # . (4) IV. RESUL TS A. Em b edded space dynamics Figure 4 shows the embedding of our spectral feature v ectors into tw o dimensions, the space ( z 1 , z 2 ), for all of the 59 individual male flies. W e first note that nearb y p oin ts hav e similar p ow er ( P k,f S (k , f ; t )), ev en though the embedding algorithm normalises-out v ariations in the total p o w er of the p ostural motions. Embedding the same data into three dimensions yields a very similar structure with less than 2% reduction of the embedding cost function (Eq. D1, Fig S3). W e generated an estimate of the probability density , b ( z ) b y con volving eac h p oint in the em b edded map with a gaussian of relativ ely small width ( σ = 1 . 5, Fig 4(b)). F ar from b eing uniformly distributed across this space, b ( z ) contains a large num b er of resolv ed local maxima. The lo cations of these p eaks provide a potential repre- sen tation for the stereotyped behaviours that the flies p erform. As exp ected, w e find that individuals displa y significan tly less intra- than inter-individual v ariation when their b ehavioural maps are compared (Fig S4). This space not only contains peaks, but the tra jectory through it also pauses at rep eatable lo cations. Through n umerical differen tiation of z 1 ( t ) and z 2 ( t ), we observe a t w o-state “pause-mo ve” pattern of dynamics. Typical time traces of z 1 ( t ) and z 2 ( t ) sho w this t yp e of tra jec- tory , with long stationary p erio ds in terspe rsed b y quic k b outs of mov emen t (Fig 5(a)). More quantitativ ely , w e find that the distribution of velocities within the embed- ded space is w ell-represen ted by a tw o-comp onent log- normal mixture mo del in whic h the the tw o peaks are separated by almost t w o orders of magnitude (Fig 5(b)). The distribution of p oints in the lo w-v elo cit y case (ap- pro ximately 45% of all time p oints) is highly lo calized with distinguishable p eaks (Fig 6). The high-v elo cit y p oin ts, in contrast, are more uniformly distributed. B. Beha vioural states The em b edded space is comprised of p eaks surrounded b y v alleys. Finding connected areas in the z 1 , z 2 plane suc h that clim bing up the gradien t of probabilit y densit y alw a ys leads to the same local maximum, often referred to as a w atershed transform [34], we delineate 122 re- gions of the em b edded space. Eac h of these con tains a single lo cal maximum of probability density (Fig 7(a)). When the tra jectory , z ( t ) pauses at one of these p eaks, w e find that each of these ep o chs corresp ond to the fly p erforming a particular stereotyped b eha viour. These pauses last an ywhere from .05 s up to nearly 25 s (Fig 8(a)). Observing segmen ts of the original movies corresp ond- ing to pauses in one of the regions, w e consistently ob- serv e the flies performing a distinct action that corre- sp onds to a recognizable b ehaviour when view ed by eye (Supplemen tary Movies S2-11). Man y of the mov ements w e detect are similar to familiar, intuitiv ely defined b e- ha vioural classifications such as walking, running, front leg grooming, and proboscis extension, but here, the seg- men tation of the movies into b ehavioural categories has emerged from the data itself, not through a priori def- initions. Moreo ver, we see that near-by regions of our b eha vioural space correspond to similar, yet distinct, be- ha viours (Fig 7(c)). This classification is consistent across individuals (Fig- ures 8-9, Supplemen tary Movies S3-11). The v ast ma jor- it y of these regions are visited by almost all of the flies at some p oint (Fig 8(b)). 104 of the 122 regions were visited b y ov er 50 (of 59 total) flies, and the remaining b eha viours w ere all low-probabilit y even ts, containing, in total, less than 3% of the ov erall activity . C. Beha vioural states as p erio dic orbits P erio dic orbits in p ostural mo v emen ts are suggestive of underlying lo w-dimensional dynamic attractors that pro duce stable b ehavioural templates [35]. These types of motifs hav e b een h ypothesized to form the basis for neural and mec hanical con trol of legged lo comotion at fast time scales [36]. Because our b ehavioural mapping algorithm is based upon similarities betw een p ostural frequencies exhibited at differen t times, a p otential hy- p othesis is that pauses in b ehavioural space corresp ond to perio dic tra jectories in the space of p ostural mo ve- men ts (Eqn 1). In our data, a fast running gait (the b ottom-most region of Fig 10(h)) corresponds to peri- o dic oscillations of the p ostural time series with a clear p eak at 12.9 Hz in the pow er sp ectral densit y (Fig 10(a)- (b)). This frequency is in go o d agreement with previous measuremen ts of the fly walking gait [37, 38]. T o systematically inv estigate the p erio dicit y of the p ostural dynamics, for eac h b eha vioural b out w e map time on to a phase v ariable, a cyclic co ordinate defined on the unit circle. This pro cess is usually referred to as phase reconstruction. The metho d we use, Phaser [39], 6 FIG. 4. Em b edding of feature vectors. (a) T raining set points embedded in to t wo dimensions via t-SNE. Color co ding is prop ortional to the logarithm of the normalization factor P k,f S ( k , f ; t ). (b) Probabilit y density function generated from em b edding all data p oints and conv olving with a gaussian ( σ = 1 . 5). p erforms Hilb ert transforms to construct phase estima- tions from several time series separately , then combines these estimates via a maxim um lik eliho o d estimate that uses F ourier-series based corrections. Here, we apply Phaser to the postural mo de time series, y k ( t ), treating the correlated motions along all 50 postural eigenmo des as synchronized oscillators. W e p erformed this recon- struction for each m ulti-cycle b ehavioural b out. After reconstructing the phases for all of the 5,483 b outs of fast running observed in male flies, w e observ e a clear pe- rio dic pattern across sev eral of the p ostural mo des (Fig 10(c)-(f )). This t yp e of analysis also brings additional insight in to the subtle distinctions b et w een our observed b ehavioural states. If we construct phase-av eraged orbits for sev en of the running b ehaviours, we observ e many differences in the gait dynamics (see Appendix E, Fig 10(g)). F or instance, we observ e an increase in many mo de ampli- tudes as the gait frequency increases (e.g. in modes 3, 12, and 13), as noted in previous work [40]. In addition, w e also see subtle changes in phase (e.g. in mode 4), as w ell as a near-elimination of a p erio d-doubled tra jectory (seen in mo de 14). This type of observ ation could allo w for a more thorough understanding of sp eed control and gait transitions in hexapo d lo comotion. W e also find oscillatory p ostural dynamics for other stereot yp ed b eha viours, with man y b eha viours resulting in a p erio dic orbit in p ostural space (Fig 10(i)). These b eha viours are found in many regions of b ehavioural space, suggesting that muc h of b eha viour is indeed con- fined to low-dimensional postural dynamics. It is impor- tan t to note that p erio dic tra jectories emerge directly from our analysis, even though the wa velet transform used to define our feature v ectors do es not preserv e phase information. D. Differences in b ehaviour b et w een males and females T o demonstrate the p ow er of this metho d to detect subtle differences in b ehaviour, we compared the b e- ha vioural spaces of male and female fruit flies by em- b edding the p ostural time series data from females into the behavioural space deriv ed from the male flies (Fig 4). Figure 11(a) displa ys the male and female behavioural probabilit y densities. W e find a striking difference b e- t w een the t wo sexes, with lo comotory behaviours greatly enhanced but resting and slow motions largely sup- pressed in females when compared to males. This is in agreemen t with previous results, showing that young fe- males are more active than their male coun terparts [41]. W e then sough t to isolate subtle b ehavioural differ- ences b et w een the sexes that are eviden t in the fine-scale structure of these maps. An example of this can b e seen in the “Wing Mov emen ts” p ortion of the behavioural space (the lo wer left corner of the map). First, w e ob- tained b oth male and female region-normalised (R-N) PDFs (Fig, 11(c)), where the in tegral of the b ehavioural space density within the ”Wing Mo vemen ts” region in- 7 FIG. 5. Dynamics within behavioural space. (a) Typical tra jectory segment through behavioural space, z 1 ( t ) (blue) and z 2 ( t ) (red). (b) Histogram of velocities in the embedded space fit to a t wo-component log-gaussian mixture mo del. The blue bar c hart represents the measured probabilit y dis- tribution, the red line is the fitted model, and the cyan and green lines are the mixture comp onents of the fitted mo del. tegrates to one. Within the space of wing mo vemen ts, w e identified regions that show statistically significan t differences b etw een the tw o sexes using a Wilco xon rank sum test [42] at eac h point in b eha vioural space. This test determines the lo cations of significant difference b e- t w een the median male PDF v alue and the median fe- male PDF v alue ( p -v alue < . 01). Regions where signifi- can t differences w ere found are indicated b y the dashed lines in Figure 11(d). P articular b ehaviours, suc h as left-wing gro oming, are sexually dimorphic (11(d) solid b ox, Mo vies S12-S13). Male-preferred gro oming includes a kic k of the middle leg on the left side of the b o dy that clears the profile of the wing and mo v es an teriorly b efore pushing bac k to w ards the p osterior. F emale-preferred gro oming lacks this additional leg mov ement. W e verified this differ- ences by isolating the mean p ostural-space orbits asso- ciated with each of these regions (Figures 11(f ), S6). Imp ortan tly , while these orbits are statistically differ- en t, the a v erage frequencies for the b ehaviours are not ( f male = 3 . 49 ± . 15 Hz v ersus f female = 3 . 28 ± . 08 Hz). W e note that these results are consistent across a large range of the b ehavioural-map smo othing parameter σ (Fig S5), such that fine-tuning of the spatial structure of the b ehavioural map is not necessary to obtain the results seen here. It should b e noted that future study is necessary to determine the ethological relev ance of these findings and to understand ho w muc h of the v ariance we observ e is related to the sp ecifics of our exp erimen tal paradigm. Ho w ev er, the fact that these distinctions are found with- out sp ecifically lo oking for an y of them – emerging only from underlying statistics of the b ehavioural map – pro- vides quantitativ e verification that the classifications we mak e are meaningful. Inherent in an y unsup ervised clas- sification metho d is the question of how to v alidate its accuracy . Here, there is no ground truth with which to compare, since a significant aim of our w ork is to disp ense with a priori b ehavioural definitions. How- ev er, by showing that meaningful distinctions and ag- glomerations can b e made b etw een different b ehavioural instances, we pro vide evidence that the approach intro- duced here can b ecome the basis undergirding a wide range of exp erimen tal in vestigations in to the b eha viour of animals. V. CONCLUSIONS The abilit y to map and compare the b ehavioural rep ertoire of individuals and p opulations of animals has applications b ey ond the study of terrestrial dynamics in fruit flies. Combined with to ols for genetic manipulation, DNA sequencing, neural imaging, and electrophysiology , the identification of subtle b ehavioural distinctions and patterns b etw een groups of individuals will impact deep questions related to the in teractions betw een genes, neu- rons, b ehaviour, and ev olution. In this initial study , we prob ed the motion of individuals in a largely feature- less en vironment. Extensions to more complicated situ- ations, e.g. where sensory inputs are measured and/or con trolled, genes are manipulated, or multiple individu- als are present, are readily implemented. Finally , we note that the only Dr osophila -specific step in our analysis pipeline is the generation of the p ostu- ral eigenmodes. Given mo vies of sufficient qualit y and length from differen t organisms, sp ectral feature vectors and b ehavioural spaces can b e similarly generated, al- lo wing for p otential applications from w orms to mice to h umans and a greater understanding of how animals b e- ha v e. A CKNO WLEDGEMENTS W e thank Yi Deng, Kieren James-Lubin, Kelsi Lind- blad, and Ugne Klibaite for assistance in data collection and analysis, as well as Jessica Cande, David Stern, and Da vid Sc h wab for discussions and suggestions. JWS and GJB also ackno wledge the How ard Hughes Medical Insti- tute Janelia F arm Visitor Program and the Aspen Center for Physics, where many ideas for this w ork were formu- 8 FIG. 6. Concentration of b ehavioural space during stereot yp ed mov ements. Comparison b et ween the densities generated during stereotyped (a) and non-stereotyped (b) ep o chs. FIG. 7. Segmen tation into b ehavioural regions. (a) Boundary lines obtained from p erforming a watershed transform on the PDF from Fig 4(b). (b) Integrated probabilities within each of the regions. (c) The organisation of b ehavioural space into regions of similar mov emen t types. Definition of regions is p erformed through visual assessmen t of movies. lated. This work was funded through a wards from the National Institutes of Health (GM098090, GM071508), The National Science F oundation (PHY–0957573, PHY– 1066293), the Pew Charitable T rusts, the Swartz F oun- dation, and the Alfred P . Sloan F oundation. D A T A ACCESSIBILITY Example co de for running the algorithms can b e found at h ttps://gith ub.com/gordon b erman/MotionMapp er. F or access to raw image data, please email JWS at shae- vitz@princeton.edu. APPENDICES App endix A: Image pro cessing T o isolate the fly from the background, w e apply Cann y’s method for edge detection [19], resulting in a bi- nary image containing the edge p ositions. W e then mor- phologically dilate this binary image by a 3 × 3 square in order to fill an y spurious holes in the edges and pro ceed to fill all closed curv es. This filled image is then morpho- logically eroded by a square of the same size, resulting in a mask. After applying this mask to the original image, w e now hav e our segmen ted image. 9 FIG. 8. Behavioural state dynamics. (a) Distribution of o c- cupancy times in all b ehaviours. (b) Num ber of individuals (out of 59 p ossible) that visit each b ehaviour at some point during observ ation. While our tracking algorithm ensures that the fly re- mains within the image b oundaries, the centre of the fly and the orien tation within the frame v ary o ver time. Ha ving obtained a sequence of isolated fly images, w e next register them b oth translationally and rotationally with resp ect to a template image. The template image is generated by taking a typical image of a fly and then man ually ablating the wings and legs digitally . F or our first step, we rotationally align. This is ac hiev ed through finding the angle that maximises the cross-correlation b etw een the magnitudes of the 2D p o- lar F ourier transforms for each image and the template. Because all translation information app ears in the phase of the 2D F ourier transform, this rotational alignmen t, based only up on the magnitude of the transform, is inde- p enden t of an y initial translations b etw een the images. Accordingly , once rotational alignment is ac hiev ed, we can subsequently register the images translationally via a cross-correlation. FIG. 9. Beha vioural space p eaks corresp ond to specific stereot yp ed b ehaviours. Selected regions within b ehavioural space are shown and are lab eled via the colour-coded legend on the right. Instances of dwells within each of these regions can be seen in Supplementary Movies 3-11. The examples displa yed in these movies are randomly selected and contain clips from many different flies, showing that the b ehavioural space provides a representation of b ehaviour that is consis- ten t across individuals. App endix B: Postural decomp osition from images The aim of the p ostural decomp osition is to tak e our set of 200 × 200 aligned images and create a low er- dimensional representation that can b e made into time series. Naively , one would simply perform PCA on the the images, using each pixel v alue as a separate dimen- sion. The fly images, how ever, contain too many pixels to analyse due to memory limitations. T o make this problem more tractable, we analyse only the subset of these pixels which ha ve non-negligible v ari- ance. Man y pixels within the fly image are either al- w a ys zero or alw a ys saturated, th us containing almost no dynamical information. Accordingly , we w ould like to use only a subsample of these measurements. The most obvious manner to go ab out this is to find the pix- els containing the highest v ariance and keep only those ab o v e a certain threshold. The primary difficulty here, ho w ev er, is that there is not an ob vious truncation p oint (Fig S2A). This is most likely the result of the fact that the fly legs can p otentially o ccupy the ma jorit y of the pixels in the image but only are present in a relatively small num b er in any given frame. Hence, many of these p eriphery pixels all ha ve similarly mo derate standard de- viations, making them difficult to differentiate. A more compact sc heme is to represen t the images in Radon-transform space, which more sparsely param- eterises lines such as legs or wing veins. After Radon transformation, the probabilit y densit y function of pixel- v alue standard deviations has a clear minimum and we k eep pixels whose standard deviation is larger than this v alue (Fig S2B). This results in keeping 6,763 pixels out of 18,090, retaining approximately 95% of the total v ari- ation in the images. If there are N images in our sample, w e can represent our data set, X as an N × 6 , 763-element matrix. W e then pro ceed to calculate the principal di- 10 FIG. 10. Periodic dynamics within b ehavioural states. (a) Periodic oscillations in the 3rd and 4th p ostural eigenmo des during a typical high-frequency running sequence. (b) Average p ow er sp ectral density for all instances of this b ehaviour (the b ottom-most region in (h)). (c)-(d) Phase reconstruction of the data in (a) for mo des three and four (respectively) . (e)-(f ) Probabilit y densities of pro jections along the 3rd and 4th mo des (resp ectively) for all instances of the b ehaviour shown in (a)-(d). The black line is the phase-av eraged curve (via Eqn E1). (g) Comparison b et ween the phase-av eraged curves for 7 differen t lo comotion gaits. Line colours are proportional to the mean gait frequency . (h) Lo comotion gaits from Fig 7(c), colour-co ded by mean frequency . The colour scale here is the same as in (g). (i) 3-D plots of the phase-av eraged tra jectories for 5 different b ehaviours. The first three p ostural modes are plotted here. (j) Regions corresp onding to the orbits shown in (i) (co ded by color). rections of v ariation in these data using PCA, as seen in Fig 3. Lastly , the question remains of ho w many mo des to k eep in our analysis, a task made more am biguous due to the smo othness of the eigen v alue sp ectrum. Our ap- proac h to determining the truncation p oint is to com- pare the PCA eigenv alues with a null model based on the noise prop erties of our data set. Specifically , we assume that the noise is due to finite data collection. Although additional errors in image segmen tation and registration assuredly exist in our data set, this set null mo del pro vides an upper-b ound on the num b er of sta- tistically meaningful eigenmo des. T o calculate this truncation p oin t, we tak e our data matrix, X , and shuffle each of its columns independently from one another, hence eliminating all meaningful cor- relations betw een them. Given finite sampling (ev en if v ery large), how ever, there will still remain some resid- ual correlations, resulting in off-diagonal non-zero terms in the cov ariance matrix. Hence, if w e diagonalise this new cov ariance matrix, the largest eigenv alue pro vides a resolution limit for our ability to distinguish signal from finite sampling noise. Performing this analysis, w e find that only 50 modes hav e eigen v alues larger than this largest sh uffled eigenv alue. These 50 mo des account for sligh tly more than 93% of the observed v ariation in the data. App endix C: W av elet calculations W e use the Morlet contin uous w a v elet transform to pro vide a multiple time scale representation of our p os- tural mode dynamics. More explicitly , w e calculate this transform, W s,τ [ y ( t )], via W s,τ [ y ( t )] = 1 √ s Z ∞ −∞ y ( t ) ψ ∗  t − τ s  dt, (C1) 11 FIG. 11. Comparison b et ween male and female b ehaviours. (a) Measured behavioural space PDF for male (left) and female (righ t) flies. (b) Difference b et ween the tw o PDFs in (a). Here we observe large dimorphisms b etw een the sexes, particularly in the ”Lo comotion Gaits” and ”Idle and Slo w Mo vemen ts” regions. (c) PDFs for b ehaviours in the “Wing Mov ements” p ortion of the b ehavioural space (the lo wer left of the full space). These PDFs (male on the left and female on the right) are normalised so that they each integrate to one. The black lines are the b oundaries found from a watershed transform and are included to guide the eye. (d) Difference b etw een the t wo normalised behavioural spaces in (c). Dashed lines enclose regions in whic h the median male and the median female PDF v alues are statistically different via the Wilcoxan rank sum test ( p < . 01). (e) Zo om-in on the boxed region in (d). Both of these regions corresp ond to left wing grooming, but with b ehaviours within the male-preferred region incorp orating an additional leg kick (Supplementary Movies S12-13). (f ) Average p erio dic orbits for p ostural eigenmodes 1, 2, 6, and 7. The area surrounding the lines represents the standard error of the mean at each p oint along the tra jectory . Average p erio dic orbits for all of the first 25 p ostural mo des are shown in Fig S6. with ψ ( η ) = π − 1 / 4 e iω 0 η e − 1 2 η 2 . (C2) Here, y i ( t ) is a p ostural time series, s is the time scale of in terest, τ is a p oin t in time, and ω 0 is a non-dimensional parameter (set to 5 here). The Morlet w av elet has the additional prop ert y that the time scale, s , is related to the F ourier frequency , f , b y s ( f ) = ω 0 + p 2 + ω 2 0 4 π f . (C3) This can b e derived b y maximizing the resp onse to a pure sine wa v e, A ( s, f ) ≡    W s,τ [ e 2 π if t ]    , with resp ect to s . Ho w ev er, A ( s, ω ) is disprop ortionally large when re- sp onding to pure sine w av es of low er frequencies. T o 12 correct for this, w e find a scalar function C ( s ) such that C ( s ) A ( s, ω ∗ ) = 1 for all s, (C4) where ω ∗ is 2 π times the F ourier frequency found in Eq. C3. F or a Morlet w a v elet, this function is C ( s ) = π − 1 4 √ 2 s e 1 4  ω 0 − √ ω 2 0 +2  2 . (C5) Accordingly , we can define our p ow er sp ectrum, S (k , f ; t ), via S (k , f ; τ ) = 1 C ( s ( f ))    W s ( f ) ,τ [ y k ( t )]    (C6) Last, we use a dyadically-spaced set of frequencies b e- t w een f min = 1 Hz and the Nyquist frequency ( f max = 50 Hz) via f i = f max 2 − ( i − 1) / ( N f − 1) log 2 f max f min (C7) for i = 1 , 2 , . . . , N f (and their corresp onding scales via Eq. C3). This creates a wa velet spectrogram that is resolv ed at multiple time-scales for each of the first 50 p ostural mo des. App endix D: t-SNE implemen tation F or our initial embedding using t-SNE, we largely fol- lo w the metho d introduced in [18], minimizing the cost function C = D K L ( P || Q ) = X ij p ij log p ij q ij , (D1) where p ij = 1 2 ( p j | i + p i | j ), q ij = (1 + ∆ 2 ij ) − 1 P k P ` 6 = k (1 + ∆ 2 k,` ) − 1 , (D2) and ∆ ij is the Euclidean distance b et w een p oints i and j in the embedded space. The cost function is optimised through a gradien t descent pro cedure that is preceded b y an early-exaggeration p erio d, allowing for the system to more readily escape lo cal minima. The memory complexit y of this algorithm prev ents the practical n umber of p oin ts from exceeding ≈ 35 , 000. Al- though improving this num b er is the sub ject of current researc h [43], our solution here is to generate an em b ed- ding using a selection of roughly 600 data p oints from eac h of the 59 individuals observed (out of ≈ 360 , 000 data p oints p er individual). T o ensure that these p oints create a representativ e sample, we p erform t-SNE on 20,000 randomly-selected data p oints from each individ- ual. This em b edding is then used to estimate a probabil- it y density by con v olving each p oint with a 2D gaussian whose whose width is equal to the distance from the p oin t to its N embed = 10 nearest neighbours. This space is segmented by applying a w atershed transform [34] to the in verse of the PDF, creating a set of regions. Finally , p oin ts are group ed by the region to whic h they b elong and the num ber of p oints selected out of each region is prop ortional to the integral ov er the PDF in that re- gion. This is p erformed for all data sets, yielding a total of 35,000 data points in the training set. Giv en the em b edding resulting from applying t-SNE to our training set, we wish to embed additional p oints in to our b ehavioural space b y comparing each to the training set individually . Mathematically , let X b e the set of all feature vectors in the training set, X 0 b e their asso ciated embeddings via t-SNE, z b e a new feature v ector that we would like to embed according to the mapping b etw een X and X 0 , and ζ be the embedding of z that w e would like to determine. As with the t-SNE cost function, we will em b ed z b y enforcing that its transition probabilities in the t wo spaces are as similar as p ossible. Lik e b efore, the tran- sitions in the full space, p j | z , are given by p j | z = exp  − d ( z , j ) 2 / 2 σ 2 z  P x ∈ X exp  − d ( z , k ) 2 / 2 σ 2 z  , (D3) where d ( z , j ) is the Kullbac k-Leibler divergence betw een z and x ∈ X , and σ z is once again found by constraining the en trop y of the condition transition probability dis- tribution, using the same parameters as for the t-SNE em b edding. Similarly , the transition probabilities in the em b edded space are giv en by q j | ζ = (1 + ∆ 2 ζ ,j ) − 1 P x 0 ∈ X 0 (1 + ∆ 2 ζ ,x 0 ) − 1 , (D4) where ∆ ζ ,x 0 is the Euclidean distance b etw een ζ and y ∈ X 0 . F or each z , w e then seek the ζ ∗ that minimises the Kullbac k-Leibler div ergence betw een the transition probabilit y distributions in the tw o spaces: ζ ∗ = arg min ζ D K L ( p x | z || q y | ζ ) (D5) = arg min ζ X x ∈ X p x | z log p x | z q y ( x ) | ζ . (D6) As b efore, this is a non-conv ex function, leading to p o- ten tial complexities in p erforming our desired optimiza- tion. How ever, if w e start a lo cal optimization (us- ing the Nelder-Mead Simplex algorithm [44, 45]) from a weigh ted av erage of points, ζ 0 , where ζ 0 = X x ∈ X p x | z y ( x ) , (D7) this point is almost alw ays within the basin of attraction of the global minimum. T o ensure that this is true in all cases, how ever, we also p erform the same minimisation pro cedure, but starting from the p oint y ( x ∗ ), where x ∗ = arg max x p x | z . (D8) 13 This returned a b etter solution appro ximately 5% of the time. Because this embedding can b e calculated indep en- den tly for each v alue of z , the algorithm scales linearly with the num b er of p oints. W e also make use of the fact that this algorithm is embarrassingly parallelizable. Moreo v er, b ecause w e ha v e set our transition en tropy , H , to b e equal to 5, there are rarely more than 50 p oints to whic h a given z has a non-zero transition probability . Accordingly , we can sp eed up our cost function ev al- uation considerably by only allo wing p x | z > 0 for the nearest 200 p oints to z in the original space. Lastly , w e find the space of behaviours for the female data sets by embedding these data into the space cre- ated with the male training set. W e find that the me- dian re-embedding cost (Eqn. D5) for the female cost is only 1% more than the median re-em b edding cost for the male data (5.08 bits vs. 5.12 bits) indicating that the embedding works well for b oth sexes. App endix E: Phase-av eraged orbits After applying the Phaser algorithm, w e find the phase-a v eraged orbit via a von Mises distribution w eigh ted av erage. More precisely , w e construct the a v- erage orbit for eigenmode k , µ ( k ) ( φ ) via µ ( k ) ( φ ) = X i y ( k ) i exp[ κ cos( φ − φ i )] P j exp[ κ cos( φ − φ j )] , (E1) where y ( k ) i is the pro jection on to the k th eigenmo de at time p oint t i , φ i is the phase associated with the same time p oint, and κ is related to the standard deviation of the v on Mises distribution ( σ 2 v M ( κ ) = 1 − I 1 ( κ ) I 0 ( κ ) , where I ν ( x ) is the mo dified Bessel function of ν th order). Here w e find the v alue of κ ≈ 50 . 3, which is the κ resulting in σ v M = . 1. Because phase reconstruction only is unique up to an additive constant, to compare phase-a v eraged curves of differen t b ehavioural b outs, an additional alignmen t needs to o ccur. This is performed b y first finding the maxim um v alue of cross-correlation b etw een the phase- a v eraged curves for eac h mo de. Then, the phase offset b et w een that pair of 50-dimensional orbits is given b y the median of these found phase shifts. 14 App endix F: Supplementary T ables T ABLE S1. Parameters used in eigen-decomposition P arameter Description V alue N θ Num b er of angles used in Radon transforms 90 M Num b er of p ostural eigenmo des (found, not defined a priori ) 50 T ABLE S2. Parameters used in w a velet analysis P arameter Description V alue N f Num b er of frequency channels 25 ω 0 Non-dimensional Morlet wa v elet parameter 5 f min High-pass frequency cut-off 1 Hz f max Lo w-pass frequency cut-off 50 Hz T ABLE S3. Parameters used in t-SNE implementation P arameter Description V alue H T ransition entrop y 5 N train T raining set size 35,000 N neighbor s Maxim um non-zero re-embedding transitions 200 N embed Nearest neighbors for training set calculation 10 T ABLE S4. Parameters used in behavioural segmentation P arameter Description V alue σ Width of gaussian kernel density estimator for embedded p oints 1.5 τ min Minim um b eha vioural time scale . 05 s σ v Width of gaussian for embedded-space velocity calculations . 02 s σ vM Standard deviation of von Mises smo othing for phase av eraging . 1 radians 15 App endix G: Supplementary Figures 0 10 20 30 40 50 60 0 0.5 1 1.5 2 2.5 3 Time (minutes) Velocity (mm/s) FIG. S1. Flies adapt to the arena in about five minutes. Median mo vemen t v elo cit y for male flies ( N = 59) during the 65 minute filming p erio d. Pink regions represent the 25% and 75% quantiles at eac h p oint in time. Note that almost all adaptation o ccurs within the first five minutes placed in the dish. This region is excluded from our analyses 16 FIG. S2. Radon vs. pixel representation of images. A) Probability density function of pixel standard deviations. B) Probabilit y densit y function of Radon pixel standard deviations. Note the clear minim um that exists in B), allowing for an effectiv e reduction in the num ber of pixels necessary to represent the data. FIG. S3. Comparison b etw een embedding in to 3-D (A) versus 2-D (B) via t-SNE. Other than the embedding dimension, all other parameters remain constan t. Color labels are proportional to the logarithm of the normalizing amplitude ( P k,f S ( k , f ; t )) for each point in the training set. There is a 2% improv ement in the error function (Equation D1) in the 3-D case as compared to the 2-D embedding (2.9 versus 3.3 bits out of a total of 20.6 bits for the transition matrix P ). 17 0 10 20 30 40 50 60 − 0.3 − 0.2 − 0.1 0 0.1 0.2 0.3 0.4 0.5 Individual # Inter − individual JS divergence − Intra − individual JS divergence FIG. S4. Inter- v ersus intra- individual v ariation. Plotted is the inter-individual v ariation minus the intra-individual v ariation for each of the 59 male flies. Here, w e measure the difference b etw een tw o b ehavioural spaces as the Jensen-Shannon (JS) div ergence b etw een their resp ective probability densities. In tra-individual v ariation is measured as the JS divergence b etw een maps generated from the first 20 minutes and the last 20 min utes of an individual data set. In ter-individual v ariation is measured as the median JS divergence b etw een the b ehavioural space b etw een individuals. Error bars are calculated via N-1 b ootstrapping on the inter-individual v ariance. Accordingly , a p ositive v alue on this plot implies that the tw o p ortions of the data set are more similar to each other than they are to other individuals’ b ehavioural spaces. W e find that 48 of the 59 individuals display significantly less intra- than inter-individual v ariations ( µ = . 16 ± . 02). 18  = 0.3  = 0.5  = 1  = 1.5  = 2  = 3.5 − 5 0 5 x 10 − 4 FIG. S5. Effect of v arying the smo othing parameter, σ , on the iden tified male/female wing mov ement distinctions. Plots are of the difference b etw een the median female region-normalized probability density function and the median male region- normalized probability density function for v arious v alues of σ . Lines within the regions encircle regions where the p -V alue of the Wilcoxon rank sum test are less than .01. 19 0 2 4 6 − 100 − 50 0 50 100 150 Mode #1 Projection 0 2 4 6 − 100 − 50 0 50 100 Mode #2 0 2 4 6 − 150 − 100 − 50 0 50 100 Mode #3 0 2 4 6 − 100 − 50 0 50 100 Mode #4 0 2 4 6 − 100 − 50 0 50 100 Mode #5 0 2 4 6 − 15 − 10 − 5 0 5 10 Mode #6 Projection 0 2 4 6 − 10 − 5 0 5 10 Mode #7 0 2 4 6 − 40 − 20 0 20 40 Mode #8 0 2 4 6 − 40 − 20 0 20 40 Mode #9 0 2 4 6 − 40 − 20 0 20 40 Mode #10 0 2 4 6 − 60 − 40 − 20 0 20 40 Mode #11 Projection 0 2 4 6 − 15 − 10 − 5 0 5 10 Mode #12 0 2 4 6 − 10 − 5 0 5 10 Mode #13 0 2 4 6 − 20 − 10 0 10 20 Mode #14 0 2 4 6 − 20 − 10 0 10 20 Mode #15 0 2 4 6 − 15 − 10 − 5 0 5 10 Mode #16 Projection 0 2 4 6 − 40 − 20 0 20 40 Mode #17 0 2 4 6 − 15 − 10 − 5 0 5 10 Mode #18 0 2 4 6 − 10 − 5 0 5 Mode #19 0 2 4 6 − 20 − 10 0 10 20 Mode #20 0 2 4 6 − 20 − 10 0 10 20 Mode #21 Projection Phase 0 2 4 6 − 40 − 20 0 20 40 Mode #22 Phase 0 2 4 6 − 20 − 10 0 10 20 Mode #23 Phase 0 2 4 6 − 20 − 10 0 10 20 Mode #24 Phase 0 2 4 6 − 10 − 5 0 5 10 Mode #25 Phase FIG. S6. Postural space p erio dic orbits for left wing gro oming b eha viours in Fig 11 in the main text. Blue curves are the a verage orbits for region ( i ), and red curves are av erage orbits for region ( ii ). Line thicknesses represent that standard error of the mean at each phase. 20 App endix H: Supplementary Mo vie Legends Supplementary movies c an b e obtaine d by emailing GJB (gb erman@princ eton.e du) Mo vie S1. Ra w video data of a b ehaving fly (left) and the corresp onding segmented and aligned data (righ t). Mo vie S2. Dynamics in b ehavioural space. Ra w video of a behaving D. melano gaster (middle) is displa y ed alongside co ordinates of the fly’s p osition within the filming apparatus (left) and its p osition in the embedded b ehavioural space (right). The red circles represent the p ositions in the appropriate co ordinate system and the trailing lines are the p ositions trav ersed in the previous .5 s. The ligh t blue shading indicates that a particular behaviour is b eing p erformed, and the blue text b elow the video of the fly gives a coarse lab el for the b ehaviour. The first p ortion of the mo vie is 5 s, pla y ed at real time (indicated b y “Real Time” ab ov e the fly video), and the subsequent p ortion of the movie is slow ed do wn by a factor of 5 for clarity (indicated b y “Slow ed 5 × ”). Mo vies S3-11. Eac h mo vie is a mosaic of m ultiple instances of sp ecific regions in b ehavioural space as display ed in Fig 9 and T able S5. Ev ery mo vie con tains m ultiple segmen ts from many different individuals and are slow ed by a factor of 4 for clarity . T ABLE S5. Behavioural Movies Mo vie Lab el Mo vie S3 Idle Mo vie S4 Right wing gro oming Mo vie S5 Left wing gro oming Mo vie S6 Left wing and legs gro oming Mo vie S7 Wing waggle Mo vie S8 Ab domen gro oming Mo vie S9 Running Mo vie S10 F ront leg gro oming Mo vie S11 Head gro oming Mo vie S12. Comp osite movie (slo wed b y a factor of 4) of randomly ch osen instances of flies from the male-preferred b eha vioural region in Fig 11 of the main text. Mo vie S13. Comp osite mo vie (slo w ed b y a factor of 4) of randomly c hosen instances of flies from the female-preferred b eha vioural region in Fig 11 of the main text. REFERENCES [1] Altmann, J., 1974 Observ ational Study of Behavior: Sampling Metho ds. Behaviour 49 , 227–267. [2] Lehner, P . A., 1996 Handb ook of etholo gic al metho ds . Cam bridge, UK: Cambridge Universit y Press, 2nd edition. [3] Gould, J. L., 1982 Etholo gy: The Me chanisms and Evolution of Behavior . New Y ork, NY: W. W. Norton & Company . [4] Stephens, G. J., Osb orne, L. C. & Bialek, W., 2011 Searching for simplicit y in the analysis of neurons and b ehavior. Pr o c. Nat. A c ad. Sci. 108 , 15565–15571. (doi:10.1073/pnas.1010868108). [5] Glimc her, P . W. & Dorris, M., 2004 Neuronal studies of decision making in the visual-saccadic system. In The Co gnitive Neur oscienc es III (ed. M. S. Gazzaniga), pp. 1215–1228. MIT Press. [6] Manoli, D. S., Meissner, G. W. & Baker, B. S., 2006 Blueprints for b ehavior: genetic sp ecification of neural circuitry for innate b ehaviors. T r ends In Neuroscienc es 29 , 444–451. (doi:10.1016/j.tins.2006.06.006). [7] Osborne, L. C., Lisb erger, S. G. & Bialek, W., 2005 A sensory source for motor v ariation. Natur e 437 , 412–416. (doi: 10.1038/nature03961). [8] Stephens, G. J., Johnson-Kerner, B., Bialek, W. & Ryu, W. S., 2008 Dimensionality and dynamics in the b ehavior of C. elegans. PL oS Comp. Bio. 4 , e1000028. (doi:10.1371/journal.p cbi.1000028). [9] Stephens, G. J., de Mesquita, M. B., Ryu, W. S. & Bialek, W., 2011 Emergence of long timescales and stereot yp ed b eha viors in Caenorhab ditis elegans. Pr oc. Nat. A c ad. Sci. 108 , 7286–7289. (doi:10.1073/pnas.1007868108). [10] Desrochers, T. M., Jin, D. Z., Go o dman, N. D. & Graybiel, A. M., 2010 Optimal habits can develop sp ontaneously through sensitivity to lo cal cost. Pr oc. Nat. A c ad. Sci. 107 , 20512–20517. (doi:10.1073/pnas.1013470107). [11] Bro wn, A. E. X., Y emini, E. I., Grundy , L. J., Jucik as, T. & Schafer, W. R., 2013 A dictionary of b ehavioral mo- tifs reveals clusters of genes affecting caenorhab ditis elegans lo comotion. Pr o c. Nat. A cad. Sci. 110 , 791–796. (doi: 10.1073/pnas.1211447110). [12] Dank ert, H., W ang, L., Ho opfer, E. D., Anderson, D. J. & Perona, P ., 2009 Automated monitoring and analysis of so cial b eha vior in Drosophila. Natur e Metho ds 6 , 297–303. (doi:10.1038/nmeth.1310). 21 [13] Branson, K., Robie, A. A., Bender, J., Perona, P . & Dickinson, M. H., 2009 High-throughput ethomics in large groups of Drosophila. Natur e Metho ds 6 , 451–457. (doi:10.1038/nmeth.1328). [14] Kabra, M., Robie, A. A., Rivera-Alba, M., Branson, S. & Branson, K., 2013 JAABA: interactiv e machine learning for automatic annotation of animal b ehavior. Natur e Metho ds 10 , 64–67. (doi:10.1038/nmeth.2281). [15] de Chaumont, F., Coura, R. D.-S., Serreau, P ., Cressant, A., Chabout, J., Granon, S. & Olivo-Marin, J.-C., 2012 Computerized video analysis of so cial interactions in mice. Natur e Metho ds 9 , 410–417. (doi:10.1038/nmeth.1924). [16] Kain, J., Stokes, C., Gaudry , Q., Song, X., F oley , J., Wilson, R. & de Biv ort, B., 2013 Leg-tracking and automated b eha vioural classification in Drosophila. Natur e Communic ations 4 , 1910. (doi:10.1038/ncomms2908). [17] Simon, J. C. & Dickinson, M. H., 2010 A new cham b er for studying the b ehavior of Drosophila. PLoS ONE 5 , e8793. (doi:10.1371/journal.p one.0008793). [18] v an der Maaten, L. & Hinton, G., 2008 Visualizing data using t-SNE. J. Mach. L e arning R ese ar ch 9 , 85. [19] Cann y , J., 1986 A computational approach to edge detection. IEEE T r ans Pattern Analysis and Machine Intel ligenc e 8 , 679–714. (doi:10.1109/TP AMI.1986.4767851). [20] De Castro, E. & Morandi, C., 1987 Registration of translated and rotated images using finite F ourier transforms. IEEE T r ans on Pattern Anal and Mach Int 5 , 700–703. [21] Reddy , B. S. & Chatterji, B. N., 1996 An FFT-based technique for translation, rotation, and scale-inv arian t image registration. IEEE T r ans on Image Pr o c essing 5 , 1266–1271. (doi:10.1109/83.506761). [22] Wilson, C. A. & Theriot, J. A., 2006 A correlation-based approach to calculate rotation and translation of moving cells. IEEE T r ans on Image Pr o c essing 15 , 1939–1951. (doi:10.1109/TIP .2006.873434). [23] Guizar-Sicairos, M. & Th urman, S. T., 2008 Efficien t subpixel image registration algorithms. Opt L ett 33 , 156–158. (doi:10.1364/OL.33.000156). [24] Revzen, S. & Guck enheimer, J. M., 2012 Finding the dimension of slow dynamics in a rh ythmic system. Journal of the R oyal So ciety Interfac e 9 , 957–971. (doi:10.1098/rsif.2011.0431). [25] Ristroph, L., Berman, G. J., Bergou, A. J., W ang, Z. J. & Cohen, I., 2009 Automated hull reconstruction motion tracking (HRMT) applied to sidewa ys maneuvers of free-flying insects. J Exp Bio 212 , 1324–1335. (doi:10.1242/jeb.025502). [26] F ontaine, E. I., Zabala, F., Dickinson, M. H. & Burdick, J. W., 2009 Wing and b o dy motion during flight initiation in Drosophila revealed by automated visual tracking. J Exp Bio 212 , 1307–1323. (doi:10.1242/jeb.025379). [27] Y e, L. & Keogh, E., 2011 Time series shapelets: a nov el tec hnique that allo ws accurate, in terpretable and fast classification. Data Mining and Know le dge Disc overy 22 , 149–182. (doi:10.1007/s10618-010-0179-5). [28] Goupillaud, P ., Grossman, A. & Morlet, J., 1984 Cycle-o ctav e and related transforms in seismic signal analysis. Ge o ex- plor ation 23 , 85–102. [29] Daubechies, I., 1992 T en L e ctur es on Wavelets . Philadelphia, P A: SIAM. [30] Co x, T. F. & Cox, M. A. A., 2000 Multidimensional Sc aling . Bo ca Raton, FL: Chapman and Hall, 2 edition. [31] T enenbaum, J. B., de Silv a, V. & Langford, J. C., 2000 A global geometric framework for nonlinear dimensionality reduction. Scienc e 290 , 2319–2323. (doi:10.1126/science.290.5500.2319). [32] Ro weis, S. T. & Saul, L. K., 2000 Nonlinear dimensionality reduction by lo cally linear em b edding. Scienc e 290 , 2323–2326. (doi:10.1126/science.290.5500.2323). [33] Co ver, T. M. & Thomas, J. A., 2006 Elements of Information The ory . Hoboken, NJ: Wiley-Interscience, 2nd edition. [34] Mey er, F., 1994 T op ographic distance and watershed lines. Signal Pro c essing 38 , 113–125. [35] F ull, R. J. & Ko ditschek, D. E., 1999 T emplates and anchors: neuromechanical hypotheses of legged lo comotion on land. Journal of Exp erimental Biolo gy 202 , 3325–3332. [36] Holmes, P ., F ull, R. J., Ko ditschek, D. & Guck enheimer, J., 2006 The dynamics of legged lo comotion: Mo dels, analyses, and challenges. Siam R eview 48 , 207–304. (doi:10.1137/S0036144504445133). [37] Strauss, R. & Heisenberg, M., 1990 Co ordination of legs during straight w alking and turning in Drosophila melanogaster. J Comp Physiol A 167 , 403–412. [38] W osnitza, A., Bo ck em¨ uhl, T., D ¨ ubb ert, M., Scholz, H. & B¨ uschges, A., 2013 Inter-leg co ordination in the con trol of w alking sp eed in Drosophila. Journal of Exp erimental Biolo gy 216 , 480–491. (doi:10.1242/jeb.078139). [39] Revzen, S. & Guck enheimer, J. M., 2008 Estimating the phase of synchronized oscillators. Phys R ev E 78 , 051907. (doi:10.1103/Ph ysRevE.78.051907). [40] Mendes, C. S., Bartos, I., Ak ay , T., Mark a, S. & Mann, R. S., 2013 Quantification of gait parameters in freely w alking wild type and sensory deprived Drosophila melanogaster. eLife 2 , e00231–e00231. (doi:10.7554/eLife.00231.027). [41] Le Bourg, E., 1987 The rate of living theory . Sp on taneous locomotor activity , aging and longevity in Drosophila melanogaster. Exp erimental ger ontolo gy 22 , 359–369. [42] Wilco xon, F., 1945 Individual comparisons by ranking metho ds. Biometrics Bul letin 1 , 80–83. [43] v an der Maaten, L., 2013 Barnes-Hut-SNE. arXiv p. 1301.3342. [44] Jongen, H. T., Meer, K. & T riesch, E., 2004 Optimization The ory . Boston, MA: Kluw er Academic Publishers. [45] Lagarias, J. C., Reeds, J. A., W right, M. H. & W right, P . E., 1998 Conv ergence prop erties of the nelder-mead simplex metho d in low dimensions. SIAM J. Optim. 9 , 112–147.