Computational Geometric Tools for Modeling Inherent Variability in Animal Behavior

A fundamental challenge for behavioral neuroscientists is to represent inherent variability among animals accurately without compromising the ability to quantify differences between conditions. We developed two new methods that apply curve and shape alignment techniques to address this issue. As a proof-of-concept we applied these methods to compare normal or alarmed behavior in pairs of medaka (Oryzias latipes). The curve alignment method we call Behavioral Distortion Distance (BDD) revealed that alarmed fish display less predictable swimming over time, even if individuals incorporate the same action patterns like immobility, sudden changes in swimming trajectory, or changing their position in the water column. The Conformal Spatiotemporal Distance (CSD) technique on the other hand revealed that, in spite of the unpredictability, alarmed individuals share an overall swim pattern, possibly accounting for the widely held notion of “stereotypy” in alarm responses. More generally, we propose that these new applications of known computational geometric techniques are useful in combination to represent, compare, and quantify complex behaviors consisting of common action patterns that differ in duration, sequence, or frequency.


Introduction
Quantitative description of the behavior of an animal is a key element of neuroscience. As important as this is, even the most comprehensive ethograms fail to capture the inherent variability in the behaviors of individuals, whether they consist of simple or complex sets of actions [1,2,3,4]. This is further compounded by the error-prone nature of subjective assessment during human observations. A pertinent issue for behavioral studies, therefore, is the development of appropriate methods for the quantitative comparison of behaviors between two animals or a single individual under two or more conditions. Researchers typically quantify parameters (such as the average velocity or average height climbed) based on expert knowledge of the ethology of that species or the kinematics of a behavior. These measures are useful, yet they can be inadequate as descriptors of behavior for a variety of reasons. One notable reason is that they are gross approximations of the behavioral action being measured. Averaging of observational data may not be representative of the kinematics and can be confounded by what has been described elsewhere as "the failure of averaging" [2]. Furthermore, there can be considerable variability in the execution of even those behaviors that are considered "stereotypical" or instinctive.
As elaborated at great length by Schleidt in 1974, individual variability can exist in the form of "incompleteness of the elements" that make up an innate behavior or in the form of changes in the characteristics of those "elements" such as their duration, sequencing, or frequency within a given episode of executing the behavior [5]. Consequently, the above mentioned gross approximations are at best an imperfect representation of the behavior that neither capture nor do justice to the intra and inter subject variability [6]. Another reason measures fail to be adequate descriptors of behavior is that even if specific parameters are defined as relevant in a behavioral readout, they may not be available when examining a new species or when studying previously uncategorized behaviors. Quantification of infrequent events (such as episodes of immobility or jumps) is also somewhat unsatisfactory as a measure since it omits and obscures a large portion of the behavioral data. In the ideal scenario, an entire behavioral episode would be represented in as complete a form as possible during quantification.
Studies over the past few years have started to address these (or related) issues by applying tools from computational geometry, computer vision, and machine learning to describe behaviors [7]. These methods have proved useful, for example in examining the entire behavioral repertoire of the hydra [8], for automated annotation of behavior in fruit flies [9], identifying the temporal features that explain spontaneous swimming behavior in worms [10], identifying the dynamics of shoaling in fish [11], and detecting sub-second modules in mice behavior [12]. Such studies are providing new insights like the continuity between behavioral states [13] and objective means for comparison of new datasets to perform spatiotemporal mapping of posture of non-stereotyped actions [6].
Here, we report the development of two new methods for the representation and quantitative analysis of locomotion of pairs of animals moving freely in a two dimensional space. We applied these methods to compare swimming behavior of pairs of fish. The first among these methods is based on an approach to automated curve alignment developed by [14] that can be efficiently calculated with a computational technique called dynamic time warping [15]. Dynamic time warping has been applied in the past to analyze biological phenomena, such as speech patterns [16] and cardiac filling phases [17]. In recent years, others have employed techniques based on dynamic time warping to examine and compare temporal patterns and acoustic features in bottle-nosed dolphins [18], whales [19], and songbirds [20]. The same technique has also been adapted to compare locomotion and infer postures in humans [21,22]. Motivated by a classic result in differential geometry called the Frenet-Serret Formula [23], we introduced the notion a behavior curve associated to an animal's trajectory and derived a measure that we call the Behavioral Distortion Distance (BDD), which can be efficiently calculated via dynamic time warping to provide meaningful quantitative comparisons for entire episodes of locomotion.
The second method we developed examines the spatiotemporal kinetics encoded in heat maps. When comparing heat maps between animals that represent the time spent at each location in an arena, it is not useful to look at a straightforward average since hot spots from different animals will average out and disappear. Instead, one must first align the heat maps so that corresponding hot spots coincide in order to reveal if there are common patterns and similarities. We implemented a well-established algorithm [24,25,26,27,28] used previously in cortical cartography [29,30], for finding conformal (angle preserving) alignments between surfaces -a heat map can be interpreted as a surface in 3-dimensional space where different intensities correspond to different heights in the z-direction -that uses discrete geometric objects called circle packings. We call the resulting measure between heat maps the Conformal Spatiotemporal Distance (CSD).
As a proof-of-concept, we asked if these methods can quantitatively distinguish normal swimming behavior from a behavior described as "alarmed" in fish. Though the species where such a phenomenon has been described have expanded considerably [31,32,33] since its discovery in minnows approximately 80 years ago [34], alarm behavior to conspecific injury has been mainly studied in the Ostariophysi species of fish [35,36]. Even among Ostariophysi, species differ in the expression of alarmed behavior. Many different parameters, in captivity and in the wild, have been described for different fish including a change in inter-individual distance in shoals, freezing, darting, change in the vertical position in a water body, etc. [37,38,39,40,41]. As such, it is difficult to predict what an alarmed behavior will look like in a species where this has not been described before. We therefore chose as our test case the single published study from one of the authors (at the time we started this project) describing an alarm response in medaka [42]. We asked if the individuals in the two conditions (alarmed and normal) were distinguishable without an experimenter's subjective criteria in defining what action sequence or pattern constitutes an alarmed state.
We found that individuals in each condition cluster together. Unexpectedly, the BDD between most alarmed individuals were also high. This suggests that swimming trajectories can become so unpredictable when medaka are alarmed that they may not match a second episode or a second individual expressing the same behavior. On the other hand, the CSD between most alarmed individuals were low, highlighting the existence of a discernible pattern in spite of low predictability. Overall, our results suggest that BDD and CSD together are useful aides to represent long episodes of swimming behavior. Our results also indicate that these computational geometric methods can be applied more universally for quantitative analysis of locomotion behavior of any kind, in any animal.

Individual Variability in Swimming Patterns
We used a data-set that recorded the swimming motion of thirty-six medaka over a 10 minute interval [42]. The observational tanks were featureless small rectangular tanks with dimensions 20 cm × 12 cm × 5 cm (L × H × W ). While the tanks restricted movement compared to the animals' natural habitat, they allowed multiple degrees of freedom in swimming motion and can be accepted as reasonably "free" spaces for motion. In particular, the fish could swim in any manner across the tank approximately six to seven body lengths long, four body lengths deep, and one and a half body lengths wide. Half of the fish were exposed to an alarm substance (Schreckstoff or SS) at the same time point (the 2 minute mark) in the assay as depicted in (Figure 1A; SS) while the other half were exposed to a control substance (tank water; NSS). We used an automated tracker to record each animal's position in the tank (Supplementary Video S1) and developed a measure of difference in behavior, or the BDD between two trajectories over a common length of time. Our model takes into account kinematic factors, such as the velocity (rate of displacement) of an animal at a given moment in time, as well as intrinsic geometric properties of the trajectory, such as curvature (the rate of change of turning while traveling at unit speed; see Figure 5 and the Methods section for details). Additionally, we included the vertical height in the tank as a parameter in our model because it is known to change for alarmed individuals in other species [43,44,45,32,37]. For all of the results listed in this article, BDD refers to the behavioral distortion distance indexed by height, velocity, and curvature, or BDD Θ for Θ = {y, v, κ} in the notation of the Methods section.
To evaluate the differences in behavior between distinct individuals, we first set out to understand the amount of variability to expect within a single individual's behavior over different periods of time. In particular, we sought to identify (1) the average BDD from an individual to itself at (non-overlapping) time intervals and (2) the length of time at which we see the least amount of variation in BDD from an individual to itself. An answer to (1) would provide a reasonable benchmark for the BDD between different individuals and an answer to (2) would suggest the time length over which we could obtain the most consistent/reliable BDD between different individuals.
As a preliminary step, we calculated the BDD from each fish to itself over a thousand randomly generated pairs of non-overlapping intervals of various lengths between 0 and 5 minutes. Figure 1B shows an example of the BDD values calculated for a representative alarmed (SS) or non-alarmed (NSS) fish. Each point in the plot corresponds to a pair of intervals with the length of the intervals plotted along the x-axis and their corresponding BDD along the y-axis. The calculations revealed that the BDD values of both the alarmed and non-alarmed fish exhibited significant variation for shorter periods of time, but quickly stabilized as the length of the time intervals increased. This is expected since the set of pairs of non-overlapping intervals with longer duration is smaller than the set of such intervals with shorter duration. For instance, there is only one non-overlapping pair of intervals of length 5 minutes, so only one BDD per fish of this duration can be generated. The stability we observed, however, made us optimistic that the questions raised by (1) and (2) could be answered precisely.
We ran the following, more structured, computation: For each fish in the experiment and interval length d between 0.25 and 4.75 minutes at 0.25 minute increments, we randomly generated a hundred pairs of non-overlapping time intervals of length d and calculated the BDD for the fish over each pair of intervals. We then calculated the mean and standard deviation of the values for each fish and choice of d. Since the size of the set of pairs of non-overlapping intervals of length d varies as a function of d 2 , we normalized the standard deviations of the BDD for each d by the size of the corresponding domain. The mean and normalized standard deviations (standard deviation of BDD per unit area of input) are shown in Table 3 with corresponding box plots in Figure 1C.
The results of this computation show that the mean BDD for the control group (NSS) is fairly level for d ≥ 0.5 after an initial substantial decline from d = 0.25. This indicates that the BDD from a control fish to itself over durations of 30 seconds or more is consistent. On the other hand, the normalized standard deviations of BDD for the control group (NSS) has a slight "J" shape with an initial decline from d = 0.25 followed by a rapidly increasing incline for d ≥ 1 with minimum between d = 0.5 and d = 0.75. This indicates that the least variation in BDD for individuals occurs over 30 to 45 second intervals, which means we can expect to obtain the most consistent comparisons of BDD over such durations.
This analysis also revealed that the mean BDD values for alarmed fish were significantly higher than that of non-alarmed fish ( Figure 1C) for the majority of lengths of comparison. As seen in Figure 1C, the first quartile of BDD values for alarmed fish is greater than the third quartile of the BDD values for non-alarmed fish for every duration except d = 4.75.
In other words, more than 75% of the BDD values for fish exposed to the alarm substance were strictly greater than at least 75% of the BDD values for fish exposed to control substance. This suggests that it is possible to predict whether or not a fish can be considered alarmed with a specified level of confidence by looking only at its BDD values with itself, that is, without the need to compare it to any other fish [46].

Behavioral Distortion Distances (BDD) between Alarmed and Non-Alarmed Medaka
With the least variation in BDD for individuals occurring over 30 second durations, we calculated the BDD between each pair of the thirty-six medaka over 30 second intervals throughout the 10 minute experiment. The results are shown in Figure 2 and detailed in Table 1.
As seen in the Table 1, the mean BDD within the control fish (NSS-NSS) was fairly constant throughout the experiment ranging between 0.173 and 0.187, only slightly higher than the mean BDD of 0.168 we observed in the intra-individual comparisons. The mean BDD between alarmed fish (SS-SS) and between the control and alarmed fish (NSS-SS) were also comparable to the control fish, but only for the first two minutes. They exhibit a significant increase after the 2 minute mark. We interpret this to be the consequence of exposure to the alarm substance delivered at the 2 minute mark as the conditions for each group were identical up to the that point. A similar pattern appears in the standard deviations of BDD, except that the values do not exhibit a significant jump until after the 4 minute mark.
In addition, we calculated the average BDD from each fish to all the other fish over each of the time intervals and sorted the fish from lowest to highest for each interval ( Figure 2). For the first 2 minutes, the two groups (NSS and SS) were evenly mixed with respect to this ordering. As time progressed, the two groups became more and more distinguished with most of the control fish at the lower end of the spectrum and most alarmed at the higher end. This pattern was most     Table 1 shows the distribution split, that is, the percentage of SS fish among the lower half in the ordering versus the percentage in the upper half of the order, and the corresponding likelihood of such a split if the fish had been ordered at random. Between 3 and 3.5 minutes, the likelihood of obtaining the ordering we found is on the order of 1 in 200000. This suggests that fish exposed to alarm substance have a distinct swimming behaviour a short while after exposure to Schreckstoff, but not immediately.

Conformal Spatiotemporal Distances (CSD) between Alarmed and Non-Alarmed Medaka
We conducted a similar analysis for the Conformal Spatiotemporal Distance (CSD) method we developed to compare heat maps representing the amount of time a fish spends in each location of the observational tank throughout the experiment or a specified portion of it (see Figure 7 and Methods). We noted that in order for discernible patterns to emerge in these heatmaps, longer intervals than those used for BDD were required. We used 2 minutes intervals over the course of the 10 minute experiment along with the entire 10 minute interval( Figure 3) and detailed in Table 2.
The mean CSD within the control fish (NSS-NSS) was essentially constant over the five 2 minute intervals ranging between 0.022 and 0.023. Once again, like the BDD,the mean difference between alarmed fish (SS-SS) and between the control and alarmed fish (NSS-SS) were comparable to control fish for the first 2 minutes, but exhibited a noticeable change after this point. However, unlike the BDD, there was drop in CSD for alarmed fish.The most substantial drop occurred during the 4-6 minute interval. The standard deviations did not exhibit any clear changes for any of the three comparisons.
As with BDD, we calculated the average CSD from each fish to all the other fish over each of the time intervals and sorted the fish from lowest to highest for each interval Figure 3. Initially, there were twice as many control fish than alarmed in the first half of this ordering, but as time progressed, that trend reversed with most of the alarmed at the lower end of the spectrum. This pattern was most apparent during the 6-8 minute interval and was mirrored by the entire experiment (0-10 minute). However, as Table 2 shows the likelihood of obtaining the ordering we found for the entire interval duration was less than 1 in 50, we consider CSD less capable of identifying a "typical" alarmed behavior. While the CSD did not perform as well as BDD when it came to classifying the state, it did capture a level of similarity in the behavior between alarmed fish that was not captured by BDD.

Discussion
Complex behaviors with common action patterns or modules, albeit in exaggerated or in a compressed forms, can be precluded from identification as different. We designed two new applications of known techniques in computational geometry for curve and surface alignment that are particularly relevant when examining complex behavior and behavioral states. We applied these techniques to reexamine the report of the existence of an alarm response in medaka [42]. Gross parameters such as the total duration of immobility among fish exposed to the alarm substance or Schreckstoff were adequate to differentiate two groups of medaka in the original report, even though individuals showed considerable variability. The analyses presented in this paper provide methods to quantify behaviors that are intrinsically variable and confirm that medaka behavior changes after exposure to Schreckstoff derived from conspecifics. This can be inferred without expert knowledge of the underlying ethology.
This report also provides some new insights into the alarm behavior of medaka that may be generalizable to other fish. First, the BDD analysis revealed that fish exposed to Schreckstoff differed from each other almost as much as they differed from non-alarmed fish. This may appear counter-intuitive at first glance if we expect the alarmed fish to behave in a comparable manner or more "uniformly" when expressing fear, however, these differences in BDD are consistent with the idea that normal kinematics are interrupted when medaka are alarmed. As a consequence, the repetitive aspects in their swimming are reduced and the swimming trajectories become less predictable overall. This is likely to be advantageous as an anti-predator strategy. Second, this analysis also revealed that the BDD from a fish to itself (Figure 1) can be a good indicator of its behavioral state. However, the caveat here is that one would need a length of recording of normal behavior (either from the same animal or comparable individuals) to establish the magnitudes of small and large BDD values.
One question that often arises in such analyses is not just whether two episodes are different, but how exactly do they differ. Imortantly, we can use our model to examine which parameters of the swimming kinematics contribute most to the difference between two behavoiral episodes. We cannot only determine whethor or not the behaviors of two animals were different for a given episode, but we can identify and measure how they were different.  there is a noticeable drop in velocity domination for most alarmed fish (SS-SS pairs; (Figure 4) and a corresponding increase in the domination of parameter vertical height in the tank. When comparing control fish against alarmed (most SS-NSS pairs) however, velocity domination increases with a corresponding decrease in the relevance of vertical height (not shown). This analysis can be interpreted most congruously in the following manner -that most alarmed medaka are fairly similar to each other when it comes to velocity changes, in particular periods of immobility, but they differ from each other primarily because of the location in the water column where they become immobile. In contrast, they differ from control fish in the frequency and the duration of immobility rather than the position occupied in the water column.
The CSD analysis was less effective than BDD for classification; however, it revealed a possible reason for human observers describing alarm behavior as being stereotypic in spite of considerable variability between individuals expressing alarm. As described in [5], "observer's eye and mind perform a complex task of averaging and assigning weights to many different characteristics", essentially integrating information over a large experimental time period. The CSD shows that most alarmed fish ( 60%; Figure 3) share an overall pattern of similarity which might be picked up by the observers as features more frequent in alarmed fish. Our analyses also reveal that control fish can exhibit "an alarmed state" and that alarmed state is an extension of the normal behavior "phase space".
A general caveat of using the BDD and CSD for analysis is that they can only be applied to a pair of curves (resp. heat maps) at a time. It therefore accompanies a degree to difficulty in visualizing a characteristic feature, if any, common to all individuals in a condition when comparing two conditions. As demonstrated here, it is possible to overcome this limitation to an extent by statistical analysis of all pairs. We expect these methods would be able to verify a proposed model for "typical" behavior, however, by showing that the BDD or CSD between the candidate model and a sufficiently large sample of control animals were less than or equal to the intra-individual BDD (resp. CSD) over different time intervals. Another caveat of these measures is the assumption that the behavior being analyzed can be described by the intrinsic geometric and kinematic properties of locomotion. Indeed, behavioural posture is unaccounted in these analyses. However, here this limitation exists because of the recording technique employed. It can be overcome if postural data is identified in videos with higher resolution recordings that could then be incorporated to compute the BDD.
More generally, we have developed the curve alignment process so that it can be used to compare any two curves that change over a common length of time. We have also devised algorithms that generate parametric curves from tracking data normally obtained for analysis of behavior (that is, xand y-coordinates, time-stamp, frame rate, etc.) These new applications can therefore be used to compare locomotion dependent behavior of any animal tracked in a two dimensional arena. It can also be applied to animals moving in three dimensional arena with simple modifications. We expect these methods will be more widely used after this demonstration of their utility. Accordingly, our software will be publicly available on a GitHub repository once this manuscript has been accepted for publication. In the meantime, please contact one of the authors for an up-to-date version of our python package.

Behavioral data
The data used for the development of methods comes from [42]. Details of the experimental setup, including the method for preparation of the alarm substance (Schreckstoff ), can be found in the methods section of the original publication. The original tracking data used in our model is uploaded to public data repository and is provided with this paper.

Curve Alignment & Behavioral Distortion Distance (BDD)
Our first method compares the behavior of two animals by comparing the intrinsic geometric properties and dynamics of their trajectories.

Model
We assume that the trajectory of each fish is a continuous, smooth parametric curve given by a vector-valued function r : [t 0 , t f ] → R 3 where t 0 and t f denote the start and end times of a specified time interval, respectively. To identify the geometric and dynamic features of the trajectories that are most relevant to behavior, we first calculate the rate at which an animal's motion is changing relative to its current frame of reference. We model an animal's frame of reference at a point along its trajectory with its Frenet-Serret frame at that point: For a non-degenerate, continuously differentiable, parametric curve C with the vector-valued function r : [t 0 , t f ] → R 3 it is well known in differential geometry [23] that the unit tangent, normal, binormal vectors, defined by respectively, at each point of the trajectory form an orthonormal basis of R 3 . Here, g (t) denotes the derivative of a function g with respect to t, ||v|| denotes the magnitude of a vector v ∈ R 3 , and × denotes the cross product. In otherwords, we can recoordinatize R 3 according to a fish's position and orientation at each moment in time. Measuring the rate of change of an animal's frame of reference at a given point is equivalent to calculating the rate of change of its Frenet-Serret frame, which is given by the famous Frenet-Serret formula, which we state below as a theorem. Theorem 1 (Frenet-Serret Formula, [47,48]). If r : [t 0 , t f ] → R 3 is a non-degenerate, continuously differentiable, parametric curve and T (t), N (t), and B(t) are its unit tangent, normal, and binormal vectors at time t, then d dt The upshot of the Frenet-Serret formula in our context is that the rate of change of the Frenet-Serret frame is completely determined by v(t), κ(t), and τ (t), which are the speed, curvature, and torsion, respectively, of r at time t. While speed will be familiar to many readers as the rate of displacement, the curvature indicates the rate at which the curve is turning while traveling at unit speed, and the torsion indicates the rate which the curve is twisting while traveling at unit speed. Portions of the trajectory, velocity, and curvature of the Medaka NSS 01 fish over the first thirty seconds of its experiment are shown in Figure 5.
We also considered the height of the fish in the tank as a behavioral factor, since it has been observed in other species that alarmed fish tend to dive more deeply than unalarmed fish. We did not, however, take into account the horizontal position in the tank, since left/right preferences and orientations have not been previously associated meaningfully in fish in the context of alarm behavior. For example, we consider two fish swimming at constant rates in circles of the same radius at constant heights in the tank to exhibit the same behavior even if one is swimming clockwise and the other counter-clockwise when viewed from above. Now that we have identified several appropriate factors for comparing behavior, we must determine how to measure the difference between those factors. To do this, we introduce the following notion: Because it is virtually impossible to deliver an external stimulus simultaneously to two freely swimming animals (the stimulus takes time to disperse through the tank), it is not meaningful to compare two behavior curves at common moments in time. For instance, a slight time lag due to different reaction times could yield arbitrarily large quantitative differences between two animals that exhibit identical sequences of behaviors with that approach. Instead, we must identify an appropriate correspondence between pairs of points from each behavior curve in order to make a meaningful comparison. This is a curve alignment problem, which we address using a method developed by [14]. Definition 3. An alignment between two parametric curves C 1 : with ϕ(u) = (ϕ 1 (u), ϕ 2 (u)) satisfying the conditions that ϕ(u 0 ) = (t 0 , t * 0 ), ϕ(u f ) = (t f , t * f ), and ϕ i is monotone increasing and differentiable for each i ∈ {1, 2}. Such a function associates a unique point on C 2 to each point on C 1 by matching An example alignment between the behavior curves with Θ = {v} for the first 30 seconds of observation for the Medaka NSS 01 fish and Medaka SS 01 fish is illustrated in Figure 6.
For every alignment ϕ between a fixed pair of parametric curves C 1 and C 2 , we can calculate the amount of energy required to bend, stretch, and/or compress C 1 onto C 2 according to ϕ.
To compare the behaviors between two trajectories, we search for an alignment between their corresponding behavior curves that minimizes the distortion energy. Definition 5. Let r 1 : [t 0 , t f ] → R d and r 2 : [t * 0 , t * f ] → R d be the trajectories of two animals and let Θ = {θ 1 , . . . , θ m } be a set of behavioral factors. The behavioral distortion distance between r 1 and r 2 with respect to Θ is where ϕ min is an alignment between the behavior curves C 1 Θ and C 2 Θ corresponding to r 1 and r 2 , respectively, with minimum distortion energy.
The BDD is a modification of the edit distance, d edit , between parametric curves presented in [14], which they show to be a pseudometric, that is, (1) d edit (C 1 , C 1 ) = 0, (2) d edit (C 1 , C 2 ) = d edit (C 2 , C 1 ), and (3) d edit (C 1 , C 3 ) ≤ d edit (C 1 , C 2 ) + d edit (C 2 , C 3 ) for all parametric curves C 1 , C 2 , and C 3 in a common Euclidean space. The fact that BDD is a pseudometric follows from very similar arguments. The edit distance and BDD fail to be metrics, however, since it is possible for distinct trajectories to have (edit or behavioral distortion) distance 0. This is perfectly fine in our context because it is reasonable for animals with distinct trajectories to exhibit identical behaviors.

Implementation
We obtain the trajectories from the "track objects" algorithm in the MetaMorph software as sequences of xand y-coordinates in the real plane. While the observational tank is 3-dimensional, the camera only captures the projection of a 2-dimensional, 20 cm × 12 cm, vertical cross-section of an animals movement. As described in [42], since videos were recorded from the front of the tank, displacement in the orthogonal width could not be recorded, effectively reducing a three dimensional space to a two dimensional observation. Ideally, we would be able to capture 3-dimensional trajectories; however, the BDD has the same derivation for 2-dimensional trajectories. The Frenet-Serret frames of a plane curve are given by the tangent and normal vectors and the Frenet-Serret formula is identical with torsion constantly equal to 0.
To approximate BDD, we use dynamic time warping (DTW) as proposed by [14]. Our implementation of the DTW algorithm is based on the version described in [49] with the distance function modified to account for velocity, curvature, and height. In particuar, we extended the existing DTW method within the Machine Learning module (mlpy) in Python to accommodate our use of multivariate data and different choices of Θ since the original DTW function only accepts two univariate sequences as data. While our method is able to handle input data of an arbitrarily high dimension (3-dimensional data and any number of behavioral factors), in our analyses thus far, we have only used x-coordinates, y-coordinates, velocity, and curvature. The velocity and curvature values are generated from the xand y-coordinates using the formulas where the derivatives are computed by applying the gradient function in the NumPy package in Python to the x and y time series. In addition, we apply a simple moving average smoothing function to the first derivatives x (t) and y (t) before computing the velocity and curvature values in order to reduce the data noise. In our tests, the smoothing window is 10 frames, spanning 5 frames prior up to 4 frames following, for instance, Furthermore, to account for variations in average values of behavioral factors between animals and differences in magnitudes between behavioral factors themselves, we normalize all of the data using the sigmoid function where µ and σ are the mean and the standard deviation of the respective input sequence. Since the range of the sigmoid function is [0, 1], the value of BDD Θ is bounded above by the number of factors in Θ.

Surface Alignment & Conformal Spatiotemporal Distance (CSD)
The second method we developed compares the spatiotemporal distributions between two animals. When comparing heat maps of time spent at each location in an arena between animals, it is not useful to look at a straightforward average since hot spots from different animals will average out and disappear. Instead, one must first align the heat maps so corresponding hot spots coincide in order to reveal if there are common patterns and similarities.

Model
To examine the overall spatiotemporal patterns of an animal within an enclosure, we study its spatiotemporal distribution function, that is, the continuous function ρ : Ω → R where Ω is the enclosure of the animal and ρ(z) is the number of seconds the animal spends at location z within Ω over a specified time interval I. When Ω is a region in the plane, this function can be visualized/interpreted as a heatmap where small values of ρ correspond to "cold" points and large values of ρ correspond to "hot" points. Alternatively, the spatiotemporal distribution function can be visualized/interpreted as the surface in R 3 with points {(x, y, ρ(x, y) | (x, y) ∈ Ω}, as illustrated in Figure 7. We assume that this surface is smooth and Riemannian, that is, it is equipped with a well-defined notion of distance and angles.
We can compare the spatiotemporal distribution functions between two animals by taking an appropriate norm between them, but that will be sensitive to a variety of factors that cannot be controlled, such as the position of an animal when the substance is administered. Mirroring our formulation of BDD, we must build some flexibility into our model by first aligning the corresponding surfaces before calculating a norm between them.
To align a pair of Riemannian surfaces S 1 and S 2 , we look at the space of diffeomorphisms between them. A differentiable map f : S 1 → S 2 is a diffeomorphism if it is a bijection and its inverse f −1 is also differentiable. A diffeomorphism gives a point-to-point correspondence between surfaces in a way that does not introduce any tears or folds, but it can distort the shape of a surface in unusual ways, for instance by stretching it. Ideally, we would like to find a diffeomorphism that does not distort S 1 while mapping it onto S 2 ; such a mapping is called an isometry. Unfortunately, an isometry between S 1 and S 2 may not exist. Any diffeomorphism between surfaces with different areas, for instance, will have to stretch one of the surfaces somewhere. While we cannot hope to find an isometry between S 1 and S 2 , we can always find a map that preserves angles, such a map is called conformal. This is the content of a deep result in geometry called the Uniformization Theorem. Theorem 6 (Uniformization Theorem, [50]). Every simply connected Riemann surface is conformally equivalent to one of three Riemann surfaces: the open unit disk, the complex plane, or the Riemann sphere.
While a conformal mapping f : S 1 → S 2 may distort distances, it will stretch S 1 uniformly in all directions at each point. This defines a function λ f : S 1 → R that maps each point z in S 1 to the factor by which vectors at p are stretched by f . The map f is an isometry if and only if λ f (z) = 1 for all z ∈ S 1 . Accordingly, one can measure how far a conformal diffeomorphism is from being an isometry be calculating how far the stretching factor is from 1 at each point. [51] introduced the symmetric distortion energy of a conformal diffeomorphism f : In our context, an optimal alignment between S 1 and S 2 will be a conformal diffeomorphism that minimizes the symmetric distortion energy. We can then measure the difference between two aligned surfaces by calculating the L 2 -norm between their corresponding functions. Definition 7. Let S 1 and S 2 be the surfaces of spatiotemporal distribution functions ρ 1 : Ω → R and ρ 2 : Ω → R, respectively, and for every diffeomorphism f : S 1 → S 2 and z ∈ Ω let z f denote the unique point such that f ((z, ρ 1 (z))) = (z f , ρ 2 (z f )). The conformal spatiotemporal distance between S 1 and S 2 is given by where ||Ω|| denotes the area of Ω and f * = argmin conformal f E sd (f ).

Implementation
To approximate the spatiotemporal distribution function for each fish, we subdivide the rectangular cross-section of the experimental enclosure into an m × n rectangular grid and assign to each subrectangle the percentage of time spent in that region throughout a specified duration of time. We model the surface corresponding to this function by triangulating each of the rectangular regions and taking the induced triangulation on the corresponding set of points in R 3 . To approximate a conformal mapping from this surface to the unit disk, it will be helpful to have a regular or uniform (all of the edge lengths are approximately equal) triangulation of the surface. Accordingly, we subdivide the edges of our preliminary triangulated surface so that all of the edge lengths are within a specified tolerance of equality and apply the Bowyer-Watson algorithm [52,53] to obtain a Delaunay triangulation of the resulting vertex set.
We apporixmate the set of conformal mappings between two such triangulations with the Discrete Uniformization theorem [54]. In particular, we use an algorithm developed by [55] to find a discrete conformal mapping f i : S i → D from each triangulated surface / heat map to a triangulation of the unit disk. Once we have these conformal flattening maps, f 1 and f 2 , we can express every conformal mapping f : is a Möbius transformation of the unit disk. Finding an optimal alignment between S 1 and S 2 is then reduced to finding a Möbius transformation m : D → D whose composition f −1 2 → m → f 1 minimizes the symmetric distortion energy. This is illustrated in Figure 7.
A well known application of Schwarz's Lemma [56] in complex analysis is that every Möbius transformation from D to itself is determined by which point is sent to the origin and a rotation. For stability reasons, we map the centroid of each surface to the origin [51]. This also reduces the minization problem to a single rotation, which we approximate by brute force, although there are algorithms in the literature for this [25]. We use the following discrete analogue of the  symmetric distortion energy presented in [51]: If T 1 and T 2 are triangulations of smooth Riemannian surfaces S 1 and S 2 , respectively, and f : S 1 → S 2 is a conformal mapping, then where e ij denotes an edge in T 1 or T 2 , denotes the length of an edge (or its image), and A ij is the sum of the areas of the two triangles that contain the edge e ij . Here the sums run over all of the interior edges in T 1 and T 2 respectively.