Introduction {#Sec1}
============

Quantitative comparison of the behavior of animals 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^[@CR1]--[@CR4]^. 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 speed 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^[@CR2]^". 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^[@CR5]^. Trait variability between individuals such as locomotor handedness is also influenced by alellic differences^[@CR6]^. Consequently, the above mentioned gross approximations are at best imperfect representations of the behavior that neither capture nor do justice to the intra and inter subject variability^[@CR7]^. 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 useful, but also somewhat unsatisfactory as a measure since it omits and obscures a large portion of the behavioral data. In an ideal scenario, an entire behavioral episode is represented in as complete a form as possible during quantification.

Studies over the past few years have addressed these (or related) issues by applying tools from computational geometry, computer vision, and machine learning to describe behaviors^[@CR8]^. These methods have proved effective in describing animal behavior more comprehensively. For example, in examining the entire behavioral repertoire of the hydra^[@CR9]^, for automated annotation of behavior in fruit flies^[@CR10]^, for discovering generative rules governing *Drosophila* locomotion^[@CR11]^, for identifying the temporal features that explain spontaneous swimming behavior in worms^[@CR12]^, for identifying the dynamics of shoaling in fish^[@CR13]^, and for detecting sub-second modules in mice behavior^[@CR14]^. Such studies are providing new insights like the continuity between behavioral states^[@CR15]^, between the generative rules of locomotion of vertebrates and invertebrates^[@CR11]^, and objective means for comparison of new data-sets to perform spatiotemporal mapping of posture of non-stereotyped actions^[@CR7]^. Here, we report on two new applications of computational geometric methods that can be applied broadly for the representation and quantitative analysis of dissimilarity in the locomotory behavior 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 applies a technique for automated curve alignment developed by Sebastian *et al*.^[@CR16]^ that can be efficiently calculated with a method called *dynamic time warping*^[@CR17]^. Dynamic time warping has been applied widely to time series analysis^[@CR18],[@CR19]^ as well as the study of biological phenomena, such as speech patterns^[@CR20]^ and cardiac filling phases^[@CR21]^. 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^[@CR22]^, whales^[@CR23]^, and songbirds^[@CR24]^. The same technique has also been adapted to compare locomotion and infer postures in humans^[@CR25],[@CR26]^. Motivated by a classic result in differential geometry called the Frenet-Serret Formula^[@CR27]^, we introduce the notion of a *behavior curve* associated to an animal's trajectory and derive a measure that we call the Behavioral Distortion Distance (BDD). This can be efficiently calculated via dynamic time warping to provide meaningful quantitative comparisons for entire episodes of locomotion.

The second application we considered 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^[@CR28]--[@CR32]^ used previously in cortical cartography^[@CR33],[@CR34]^, 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^[@CR35]^. We compared medaka fish from a preexisting data-set that had been exposed either to an alarm substance or to a control^[@CR36]^. Though the species where such a phenomenon has been described have expanded considerably^[@CR37]--[@CR40]^ since its discovery in minnows approximately 80 years ago^[@CR41]^, alarm behavior to conspecific injury has been mainly studied in the Ostariophysi species of fish^[@CR42],[@CR43]^. 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.^[@CR44]--[@CR48]^. As such, it is difficult to predict what an alarmed behavior in response to the chemosensory alarm substance will look like in a species where this has not been described before, or when they have multiple defensive behaviors^[@CR49]^. We therefore chose as our test case a single published study (at the time we started this project) from one of the authors describing an alarm response in medaka^[@CR36]^. We asked if the individuals in the two conditions (alarmed and normal) were distinguishable without us defining what action sequence or locomotory pattern constitutes an alarmed state.

We found that BDD and CSD can differentiate between alarmed and normal states. The BDD between most pairs of alarmed individuals were also high, which suggests that swimming trajectories can become so unpredictable when medaka are alarmed that they may not match a second individual responding to the same stimulus. On the other hand, the CSD between most pairs of alarmed individuals were low, highlighting the existence of a discernible spatiotemporal pattern in spite of low predictability in their swimming trajectory. 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 comparison of locomotion of any kind, in any animal.

Results {#Sec2}
=======

The preexisting data-set we used for this report recorded the swimming motion of 36 medaka over a 10-minute interval^[@CR36]^. 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 (minute 2 in a 10 minute observation), while the other half were exposed to a control substance (tank water or NSS). Mathuru^[@CR36]^ discovered a statistically significant difference in the average speed of the experimental subjects (SS) compared to the control subjects (NSS) during the post-stimulus (2--10 minute) time interval (Mann-Whitney U = 49.0, p = 0.0003) and a statistically significant difference in intervals of immobility (defined as displacement less than 3 mm in 1 second by the experimenter; Mann-Whitney U = 13.0, p = 0.000002).

Here, we explored the use of two additional measures that rely on notions of "distance" between the swimming patterns of pairs of subjects using the tracked data generated from the automated tracker in the previous study^[@CR36]^. The BDD or the Behavioral Distortion Distance, defined in the methods section, measures the differences in local locomotory behavior (independent of location and relative time) whereas the CSD or the Conformal Spatiotemporal Distance, also defined in the methods section, measures the spatiotemporal location pattern of a subject over a specified duration. For all of the results listed in this article, BDD refers to the behavioral distortion distance indexed by speed and curvature, or BDD~Θ~ for Θ = {*s*, *κ*} in the notation of the Methods section.

Analysis of pairwise differences in behavior {#Sec3}
--------------------------------------------

The BDD and CSD between all pairs of the 36 test subjects over five 2-minute time intervals and the entire post-stimulus (minutes 2 to 10) time interval are illustrated in Fig. [1](#Fig1){ref-type="fig"}. The pairs of subjects are split into three categories: pairs of control subjects (NSS-NSS), mixed pairs with one control and one experimental subject (NSS-SS), and pairs of experimental subjects (SS-SS). The BDD and CSD for each pair of subjects and specified duration are illustrated in Figs [2](#Fig2){ref-type="fig"} and [3](#Fig3){ref-type="fig"}, respectively, and the summary statistics are given in Tables [1](#Tab1){ref-type="table"} and [2](#Tab2){ref-type="table"}.Figure 1Panels (A--C) show the distributions, means, and standard deviations of BDD for pairs of fish over five 2 minute time intervals broken down into three categories: pairs of control fish (NSS-NSS; Green), mixed pairs with one control and one alarmed fish (NSS-SS; Orange), and pairs of alarmed fish (SS-SS; Purple). Panel (D) shows the distribution, mean, and standard deviation of BDD for pairs in each category over the entire post-stimulus time interval. Panels (E--G) show the distributions, means, and standard deviations of CSD for all pairs of fish over five 2-minute time intervals broken down, again, into three categories: pairs of control fish (NSS-NSS; Green), mixed pairs with one control and one alarmed fish (NSS-SS; Orange), and pairs of alarmed fish (SS-SS; Purple). Panel (H) shows the distribution, mean, and standard deviation of CSD for pairs in each category over the entire post-stimulus time interval. The distributions over the pre-stimulus interval (minutes 0 to 2) are all plotted in blue.Figure 2The Behavioral Distortion Distances (BDD) between all pairs of the thirty-six fish in the dataset over five 2-minute time intervals and the entire post-stimulus time interval. The alarmed (SS) individuals are labeled magenta and the control (NSS) individuals are labeled cyan.Figure 3The Conformal Spatiotemporal Distances (CSD) between all pairs of the thirty-six fish in the dataset over five 2-minute time intervals and the entire post-stimulus time interval. The alarmed (SS) individuals are labeled magenta and the control (NSS) individuals are labeled cyan.Table 1The means and standard deviations for the BDD between pairs of subjects over five 2-minute time intervals and the post-stimulus time interval, broken down into three categories: pairs of control subjects (NSS-NSS), mixed pairs with one control and one experimental subject (NSS-SS), and pairs of experimental subjects (SS-SS).Time Interval\
(minutes)MeanStandard DeviationAllNSS-NSSNSS-SSSS-SSAllNSS-NSSNSS-SSSS-SS0--20.0610.0620.0610.0610.00760.00740.00760.00762--40.0730.0650.0760.0740.01450.01050.01460.01454--60.1090.0660.1150.1400.10570.01150.10550.13946--80.1080.0690.1110.1400.10710.01320.10680.14228--100.1030.0660.1050.1360.10800.01140.10770.14452--100.0750.0610.0800.0780.02070.00850.02230.0200Table 2The means and standard deviations for the CSD between pairs of subjects over five 2-minute time intervals and the post-stimulus time interval, broken down into three categories: pairs of control subjects (NSS-NSS), mixed pairs with one control and one experimental subject (NSS-SS), and pairs of experimental subjects (SS-SS).Time Interval\
(minutes)MeanStandard DeviationAllNSS-NSSNSS-SSSS-SSAllNSS-NSSNSS-SSSS-SS0--20.01390.01330.01400.01420.003880.003830.003940.003752--40.01250.01350.01250.01140.003730.003480.003710.003744--60.01160.01370.01170.009570.003940.003820.003610.003656--80.01180.01350.01190.00990.004100.004130.003960.003528--100.01190.01370.01190.01020.003920.004090.003640.003522--100.01270.01350.01290.01140.004080.004090.003930.00406

Figure [1](#Fig1){ref-type="fig"} reveals that while the pre-stimulus interval BDD are similar between the three categories, the mean BDD between NSS-SS pairs and between SS-SS pairs are both larger than the mean BDD in the pre-stimulus period, as well as the mean BDD between NSS-NSS pairs for each of the post-stimulus time intervals. In a comparable manner, the mean CSD between NSS-SS pairs and between SS-SS pairs are both smaller than the mean CSD in each of the post-stimulus time intervals, while it remains the same between NSS-NSS pairs. Figure [2](#Fig2){ref-type="fig"} provides a colored visual aide to compare BDD between all pairs, while Table [1](#Tab1){ref-type="table"} summarizes this information. It shows that the mean BDD ranged between 0.061 and 0.069 for NSS-NSS pairs for both the pre- and post-stimulus intervals. However, the ranges increased in the post-stimulus period to between 0.074 and 0.141 for NSS-SS pairs and SS-SS pairs. In a similar manner, Fig. [3](#Fig3){ref-type="fig"} provides a colored visual aide to compare CSD between all pairs and Table [2](#Tab2){ref-type="table"} and shows that the mean CSD between NSS-NSS pairs over the pre- and post-stimulus time intervals ranged between 0.0133 and 0.0137, while, for NSS-SS pairs and SS-SS pairs, the range decreased from approximately 0.014 to between 0.0095 and 0.0125. Kruskal-Wallis H statistical tests confirmed that the differences over the five 2-minute time intervals for NSS-NSS pairs for both the BDD (Fig. [1A](#Fig1){ref-type="fig"}, H-statistic = 6.759, *p* = 0.149) and the CSD (Fig. [1E](#Fig1){ref-type="fig"}, H-statistic = 0.897, *p* = 0.925) were statistically insignificant. However, they did reveal statistically significant differences for the BDD between NSS-SS pairs (Fig. [1B](#Fig1){ref-type="fig"}, H-statistic = 387.579, *p* = 1.34 × 10^−82^), between SS-SS pairs (Fig. [1C](#Fig1){ref-type="fig"}, H-statistic = 259.470, *p* = 5.92 × 10^−55^); and the CSD between NSS-SS pairs (Fig. [1F](#Fig1){ref-type="fig"}, H-statistic = 76.781, *p* = 8.37 × 10^−16^), between SS-SS pairs (Fig. [1G](#Fig1){ref-type="fig"}, H-statistic = 124.940, *p* = 4.70 × 10^−26^). This suggests that the BDD and CSD over each of the specified 2-minute time intervals throughout the duration of the assay (10-minute) is unchanging when control substance or tank water is introduced, but statistically significant differences emerge when one or both members of the pair are exposed to *Schreckstoff*.

To evaluate if the exposure to *Shcreckstoff* is the reason for the differences in the BDD and CSD measures above, we performed pairwise comparisons between the pre-stimulus time interval (0--2 minutes) and each of the four 2-minute post-stimulus intervals via Wilcoxon signed-rank T tests. This revealed statistically significant increases in the BDD between NSS-SS pairs (Fig. [1B](#Fig1){ref-type="fig"}, all T-statistics \<4600, all *p*-values \< 10^−37^) and the BDD between SS-SS pairs (Fig. [1C](#Fig1){ref-type="fig"}, all T-statistics \<1400, all *p*-values \< 10^−15^). Statistically significant decreases in the CSD between NSS-SS pairs (Fig. [1F](#Fig1){ref-type="fig"}, all T-statistics \<18000, all *p*-values \< 10^−6^) and between SS-SS pairs (Fig. [1G](#Fig1){ref-type="fig"}, all T-statistics \<2600, all *p*-values \< 10^−8^) were also detected. Therefore, the BDD increases while the CSD decreases between pairs of medaka after stimulus delivery when at least one individual is exposed to *Schreckstoff*. Next, we also conducted pairwise comparisons of the distributions of BDD and CSD between the three categories of pairs over the entire post-stimulus time interval. In Fig. [1D](#Fig1){ref-type="fig"}, the mean BDD between NSS-SS pairs and SS-SS pairs are both larger than the mean BDD between NSS-NSS pairs. To verify that these differences are statistically significant, we conducted Mann-Whitney U tests between the distributions of BDD for NSS-NSS pairs and NSS-SS pairs (Fig. [1D](#Fig1){ref-type="fig"}, U-statistic = 10730.0, *p* = 1.49 × 10^−23^) and NSS-NSS pairs and SS-SS pairs (Fig. [1D](#Fig1){ref-type="fig"}, U-statistic = 5282.0, *p* = 1.05 × 10^−16^). In Fig. [1H](#Fig1){ref-type="fig"}, the mean CSD between SS-SS pairs is less than that between NSS-SS pairs, which is less than that between NSS-NSS pairs. Again, to verify that these differences are statistically significant, we conducted Mann-Whitney U tests between the distributions of CSD for NSS-NSS pairs and SS-SS pairs (Fig. [1H](#Fig1){ref-type="fig"}, U-statistic = 15126.0, *p* = 9.84 × 10^−6^) and for NSS-SS pairs and SS-SS pairs (Fig. [1H](#Fig1){ref-type="fig"}, U-statistic = 30470.0, *p* = 5.24 × 10^−5^). The Mann-Whitney U tests did not reveal statistically significant differences between the distributions of BDD for NSS-SS pairs and SS-SS pairs (Fig. [1D](#Fig1){ref-type="fig"}, U-statistic = 25437.0, *p* = 0.64), nor did they reveal statistically significant differences between the distributions of CSD for NSS-NSS pairs and NSS-SS pairs (Fig. [1H](#Fig1){ref-type="fig"}, U-statistic = 26729.5, *p* = 0.17). Therefore, this suggests that the BDD detected different levels of variation in behavior between pairs of control subjects and pairs with at least one experimental subject. On the other hand, CSD detected different levels of variation in behavior between pairs of experimental subjects and pairs with at least one control subject over the post-stimulus time interval. Moreover, the variations in behavior measured by BDD between pairs of experimental subjects were as large as those between mixed pairs with one control and one experimental subject.

We also compared the distributions of BDD and CSD between the three categories of pairs and the set of all pairs (again over the entire post-stimulus time interval). Mann-Whitney U tests revealed that the mean BDD between NSS-NSS pairs is significantly lower than that of all pairs (U-statistic = 68674.5, *p* = 3.33 × 10^−16^) and the mean BDD between NSS-SS pairs is significantly higher than that of all pairs (U-statistic = 87351.5, *p* = 2.63 × 10^−4^). Additionally, Mann-Whitney U tests revealed the mean CSD between SS-SS pairs is significantly lower than that of all pairs (U-statistic = 57300.5, *p* = 2.85 × 10^−4^). The Mann-Whitney U tests did not reveal statistically significant differences between the distributions of BDD for SS-SS pairs and all pairs, nor did they reveal significant differences between the distributions of CSD for NSS-NSS pairs and all pairs, nor NSS-SS pairs and all pairs. Therefore, the variations in behavior detected by BDD between NSS-NSS pairs and NSS-SS pairs are significant compared to the set of all pairs, not just to SS-SS pairs and each other. Similarly, the variations in behavior detected by CSD between SS-SS pairs are significant compared to the set of all pairs, not just to the other two categories of pairs (NSS-NSS and NSS-SS).

Finally, to further illustrate that the differences reported in the previous paragraphs are unlikely to be the results of random fluctuations, we conducted an additional analysis based on randomized (also known as permutation) tests^[@CR50]^. For this, we compared the mean BDD between NSS-SS pairs to a distribution of mean BDD between randomly relabeled pairs. Specifically, we uniformly sampled 100000 of the 9075135300 ways to divide the 36 test subjects into two 18-element subsets and calculated the mean BDD between the subsets in each of those partitions. Then, we compared the mean BDD between the NSS-SS partition to the distribution of mean BDD for randomized partitions. The corresponding distribution statistics, standard scores, and *p*-values are listed in Table [3](#Tab3){ref-type="table"}.Table 3The mean BDD between NSS-SS pairs compared to the distribution of mean BDD between pairs from randomly generated partitions of the 36 assay subjects into two 18-element subsets, with corresponding standard scores and *p*-values, over each of the specified time intervals.Time Interval\
(minutes)Randomized PartitionsNSS-SS PartitionDistribution MeanDistribution Std. Dev.Mean BDDZ-scoreP-value0--20.0610.000230.061−0.480.6282--40.0730.000460.0766.934.21 × 10^−12^4--60.1090.000690.1158.423.80 × 10^−17^6--80.1080.000470.1116.672.64 × 10^−11^8--100.1030.000420.1054.280.0000192--100.0750.000540.0799.242.36 × 10^−20^

Together, these results and statistical analyses suggest that the variability in dynamics and trajectory increase, while the variability in spatiotemporal preference decreases, when the pair of medaka being compared contains at least one alarmed individual.

Individual variability and classification {#Sec4}
-----------------------------------------

Following our analyses of pairwise comparisons, we asked if BDD and CSD could yield meaningful information about individual behavioral episodes, without direct comparison to other episodes in a larger data-set. One approach for this is to measure the distance between a subject and itself over different (non-overlapping) time intervals. This is not practical for CSD, since that method only captures meaningful information over extended periods of time. We calculated the mean BDD from each subject to itself over 1000 non-overlapping sub-intervals of the post-stimulus time interval, sampled uniformly at random. The resulting values, which we call the Intra-Individual BDD (IIBDD) for each subject, are listed in Supplementary Table [1](#MOESM1){ref-type="media"}. A Mann-Whitney U test showed that the mean IIBDD of experimental subjects is significantly higher than that of control subjects (U-statistic = 23.5, *p*-value = 1.25 × 10^−5^). Therefore, the variation in behavior measured by BDD from an experimental subject to itself is, on average, significantly larger than that from a control subject to itself.

This led us to ask if IIBDD could be used to classify normal and alarmed behaviour by itself. We conducted three logistic regressions on the thirty-six test subjects using average speed, and immobility, the parameters used in the previous report, and IIBDD. The logistic regression for average speed had an intercept of 4.449 and a coefficient of −0.168 (*χ*^2^ = 15.562, *p*-value = 7.98 × 10^−5^); the logistic regression for immobility had an intercept of −3.005 and a coefficient of 0.302 (*χ*^2^ = 30.659, *p*-value = 3.08 × 10^−8^); and the logistic regression for IIBDD had an intercept of −20.680 and a coefficient of 244.250 (*χ*^2^ = 25.494, *p*-value = 4.44 × 10^−7^). The data used for these regressions are in Supplementary Table [1](#MOESM1){ref-type="media"} and the corresponding classifications are given in Table [4](#Tab4){ref-type="table"}. The logistic regressions show that all three parameters can be used to classify normal and alarmed behavior, with duration of immobility and IIBDD exhibiting better goodness of fit than average speed. Therefore, IIBDD allows for classification of a new animal's behavior as alarmed or normal if the new experiment replicates the conditions of the original experiment.Table 4Logistic regression results for classifying normal (NSS) versus alarmed (SS) behavior given the average speed, duration of immobility, and Intra-Individual BDD (IIBDD), respectively. Regression outputs less than 0.5 are classified as NSS and those greater than 0.5 are classified as SS. Incorrect classifications are indicated with bold font.SubjectSpeed ClassificationImmobility ClassificationIIBDD ClassificationNSS 01NSS (0.189)NSS (0.069)NSS (0.065)NSS 02NSS (0.222)NSS (0.074)**SS** (**0**.**754**)NSS 03**SS** (**0**.**764**)**SS** (**0**.**658**)NSS (0.437)NSS 04NSS (0.218)NSS (0.137)NSS (0.152)NSS 05**SS** (**0**.**760**)**SS** (**0**.**804**)NSS (0.411)NSS 06NSS (0.209)NSS (0.045)NSS (0.063)NSS 07NSS (0.309)NSS (0.080)NSS (0.038)NSS 08NSS (0.029)NSS (0.056)NSS (0.009)NSS 09NSS (0.141)NSS (0.052)NSS (0.032)NSS 10NSS (0.251)NSS (0.074)NSS (0.009)NSS 11NSS (0.147)NSS (0.045)NSS (0.006)NSS 12**SS** (**0**.**593**)NSS (0.120)NSS (0.447)NSS 13**SS** (**0**.**591**)NSS (0.347)NSS (0.257)NSS 14**SS** (**0**.**523**)NSS (0.267)**SS** (**0**.**806**)NSS 15NSS (0.119)NSS (0.049)NSS (0.031)NSS 16NSS (0.108)NSS (0.074)NSS (0.193)NSS 17NSS (0.332)NSS (0.056)NSS (0.085)NSS 18NSS (0.118)NSS (0.056)NSS (0.297)SS 01SS (0.507)SS (0.998)SS (0.988)SS 03SS (0.866)SS (1.000)SS (0.889)SS 04**NSS** (**0**.**380**)SS (0.791)SS (0.966)SS 05SS (0.864)SS (0.998)SS (1.000)SS 06SS (0.960)SS (1.000)SS (0.613)SS 07SS (0.636)SS (0.874)SS (0.847)SS 08**NSS** (**0**.**405**)SS (0.952)**NSS** (**0**.**323**)SS 09SS (0.816)SS (0.927)SS (0.518)SS 10**NSS** (**0**.**140**)**NSS** (**0**.**188**)SS (0.534)SS 14SS (0.757)**NSS** (**0**.**493**)SS (0.635)SS 15SS (0.782)SS (0.998)SS (0.993)SS 16SS (0.872)SS (1.000)SS (1.000)SS 17SS (0.779)SS (0.813)SS (0.986)SS 20SS (0.907)SS (1.000)**NSS** (**0**.**266**)SS 21**NSS** (**0**.**132**)**NSS** (**0**.**200**)SS (1.000)SS 22SS (0.962)SS (1.000)SS (1.000)SS 23SS (0.970)SS (1.000)SS (1.000)SS 24SS (0.607)SS (0.512)**NSS** (**0**.**349**)Accuracy75%86%86%

Discussion {#Sec5}
==========

Complex behaviors with common action patterns or modules, albeit in exaggerated or in compressed forms, can be precluded from identification as different. To address this issue we investigated two new applications of known techniques in computational geometry for curve alignment, which we call the Behavioral Distortion Distance or BDD, and for surface alignment that we call the Conformal Spatiotemporal Distance or CSD. These applications were motivated by the aim of devising generic methods that can be used to detect differences in a locomotory pattern without the need for experimenters to define the action sequence that constitutes a particular behavioral state. We applied them here to reexamine a report of an alarm response in medaka^[@CR36]^. As the pattern of change in locomotion when a fish is alarmed is not universal and can vary in a species or context dependent manner^[@CR38],[@CR39]^, it serves as a good test case for the two applications. We find that both BDD and CSD are suitable for quantitative comparison of such behaviors without the need to rely on domain specific expert knowledge of the underlying ethology. They confirm that indeed medaka behavior changes after exposure to *Schreckstoff* derived from conspecifics. The methods section describes the mathematical framework and definitions applied to develop these applications in detail.

The section analyzing pairwise differences shows that the fish in the two groups have similar locomotory behavior in the first 2-minutes, or in the pre-stimulus period. Therefore, there is no intrinsic behavioral difference between the groups at the onset. A significant change in BDD and CSD between the pre- and post-stimulus time intervals is noticeable only when the pairs had at least one fish exposed to *Schreckstoff*. This can be visualized more clearly in Figs [2](#Fig2){ref-type="fig"} and [3](#Fig3){ref-type="fig"}. In such cases, the changes in BDD and CSD were complementary in the sense that the BDD was larger and the CSD was smaller. Notably, we did not observe significant differences in BDD between pairs of experimental fish and mixed pairs with one control fish and one experimental fish. This suggests that the behaviors of pairs of alarmed fish are as different from one another as they are from fish exhibiting normal behavior. In particular, BDD did not identify a stereotypical alarm response. While the CSD between mixed pairs was not significantly different from that of control pairs and experimental pairs, the mean CSD between experimental pairs was significantly lower than control (in fact, all) pairs. This suggests that alarmed fish exhibited a lower level of variability in their action pattern. One possible interpretation for this, following a visual inspection of the heatmaps of these individuals, is that this reduced variability could be due to periods of immobility represented as hot spots in their heatmaps.

Our analysis of the medaka alarm response using BDD and CSD brings together two seemingly contrasting ideas of stereotypic innate behavior and individual variability. Combining our observations, we conclude that the alarm response to *Schreckstoff* in medaka comprises increased variation in the local geometry and kinematics of a trajectory with a common spatiotemporal pattern. This is consistent with the ethogram of periods of erratic swimming interspersed with periods of immobility. Increased variation in the local geometry suggests swimming trajectories become less predictable over time in alarmed fish that likely serves as an advantageous anti-predator strategy. Yet, the ability of CSD to capture a commonality mirrors a description in Schleidt^[@CR5]^ of "observer's eye and mind perform a complex task of averaging and assigning weights to many different characteristics", essentially integrating information over the experimental time period.

Following our analysis of pairwise BDD and CSD, we investigated the ability of Intra-Individual BDD (IIBDD, the average BDD from a subject to itself over many randomly sampled non-overlapping time intervals) to classify normal and alarmed behaviors. We found that IIBDD performed comparably in this task to the duration of immobility via logistic regressions.

For our implementation of BDD, we took into account kinematic factors of medaka swimming such as the instantaneous speed (rate of displacement), as well as intrinsic geometric properties of the trajectory such as curvature (the rate of change of turning while traveling at unit speed; see Fig. [4](#Fig4){ref-type="fig"} and Supplementary Video [S1](#MOESM2){ref-type="media"}). Our choice of these factors was dependent on the information available in the original data-set, however, additional factors such as posture, if available, can be easily incorporated to calculate the distance between two generic behavioral curves. When comparing two behavioral curves, multiple optimal solutions are possible. In our context, such cases are considered equivalent. For instance, if we compare an episode in which an animal turned exactly once to an episode where a second animal turned exactly twice with the same curvature for each turn and speed, there will be two minima for BDD - aligning the single spike in curvature from the first episode to either of the two spikes in curvature from the second episode. Both these possible alternative alignments will give the same BDD score. We consider the two alternative alignments to be equivalent as they do not reveal any additional information about the behavior of the two animals other than the fact that the two episodes exhibited different behaviors.Figure 4The trajectory, speed, and curvature plots of the first 15 seconds of the fish NSS 01. The same color scale is used in each plot so corresponding points on each curve are colored with the same color.

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 when comparing two conditions. As demonstrated here, it is possible to overcome this limitation to an extent by (1) statistical analysis of all pairs and (2) intra-individual comparisons. 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 and CSD. In medaka, a higher BDD value and a lower CSD value represented an alarmed state. This may not be a characteristic feature of alarmed states in all species. However, if there is any substantial change in the pattern of locomotion compared to normal state in a species, it will likely be revealed in a change in BDD and CSD values in spite of individual differences in the execution of the motor action.

More generally, we have implemented the curve alignment process so that it can be used to compare any two curves that change over a common length of time. Our implementation can generate parametric curves from tracking data normally obtained for analysis of behavior (that is, *x*- and *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. They will be particularly useful for characterizing the degree of similarity when comparing two or more conditions - such as when different treatments of a drug are applied to different groups, or in studying differences in exploratory pattern in different groups of animals. As demonstrated above, another useful aspect of these measures is that once a data-set is generated, the IIBDD from a smaller subset in a subsequent experiment is also informative provided the conditions of the original experiment are replicated when generating new data. These applications can also be applied to animals moving in 3-dimensional arena. We expect these methods will be more widely used after this demonstration of their utility. Accordingly, our software package is publicly available at the GitHub repository: <https://github.com/mechunderlyingbehavior/locomotion.git>.

Methods {#Sec6}
=======

Behavioral data {#Sec7}
---------------

The data used for testing the applications comes from Mathuru^[@CR36]^. 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. Briefly, fish behavior was recorded on an external Agent v5 HD web camera placed approximately 50 cm in front of the tank. Recorded videos were digitized at 20 frames/second. The frames were converted to 8-bit grayscale thresholded images by subtracting the background on ImageJ (from <https://imagej.nih.gov/ij/>). The spatial resolution of the videos was 2.4 pixels/mm. Fish position was tracked automatically on MetaMorph 6.3 (from <https://www.moleculardevices.com>) using the "track objects" algorithm. The original tracking data used in our model is uploaded to public data repository and is provided with this paper. Experiments in the original dataset were performed in accordance to the guidelines recommended by the Institutional Animal Care and Use Committee (IACUC) of the Biological Resource Center at A\*STAR. Approved experimental protocols (IACUC 161110) were followed.

Curve alignment and behavioral distortion distance (BDD) {#Sec8}
--------------------------------------------------------

Our first method compares the behavior of two animals by comparing the intrinsic geometric properties and dynamics of their trajectories. It is an application of the *edit distance* between curves introduced by Sebastian *et al*.^[@CR16]^, which can be calculated via a technique known as *dynamic time warping* introduced by Sakoe and Chiba^[@CR20]^.

### Theoretical framework and definitions {#Sec9}

We assume that the trajectory of each fish is a continuous, smooth parametric curve given by a vector-valued function $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\bf{r}}:[{t}_{0},{t}_{f}]\to {{\mathbb{R}}}^{3}$$\end{document}$ 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 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\bf{r}}:[{t}_{0},{t}_{f}]\to {{\mathbb{R}}}^{3}$$\end{document}$ it is well known in differential geometry^[@CR27]^ that the unit tangent, normal, binormal vectors, defined by$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$T(t)=\frac{{\bf{r}}^{\prime} (t)}{||{\bf{r}}^{\prime} (t)||},\,N(t)=\frac{T^{\prime} (t)}{||T^{\prime} (t)||},\,{\rm{and}}\,B(t)=T(t)\times N(t),$$\end{document}$$respectively, at each point of the trajectory form an orthonormal basis of $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${{\mathbb{R}}}^{3}$$\end{document}$. (Here, *g*′(*t*) denotes the derivative of a function *g* with respect to *t*, \|\|**v**\|\| denotes the magnitude of a vector $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\bf{v}}\in {{\mathbb{R}}}^{3}$$\end{document}$, and × denotes the cross product.) In other words, we can recoordinatize $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${{\mathbb{R}}}^{3}$$\end{document}$ 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, given by the famous Frenet-Serret formula^[@CR51],[@CR52]^, which states that if $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\bf{r}}:[{t}_{0},{t}_{f}]\to {{\mathbb{R}}}^{3}$$\end{document}$ 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$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\frac{d}{dt}\,[\begin{array}{c}T(t)\\ N(t)\\ B(t)\end{array}]=s(t)\,[\begin{array}{ccc}0 & \kappa (t) & 0\\ -\kappa (t) & 0 & \tau (t)\\ 0 & -\tau (t) & 0\end{array}]\,[\begin{array}{c}T(t)\\ N(t)\\ B(t)\end{array}]$$\end{document}$$where$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$s(t)=\Vert {{\bf{r}}}^{{\rm{^{\prime} }}}(t)\Vert ,\,\kappa (t)=\frac{||{{\bf{r}}}^{{\rm{^{\prime} }}}(t)\times {{\bf{r}}}^{{\rm{^{\prime} }}{\rm{^{\prime} }}}(t)||}{||{{\bf{r}}}^{{\rm{^{\prime} }}}(t{)||}^{3}},\,{\rm{a}}{\rm{n}}{\rm{d}}\,\tau (t)=\frac{det({{\bf{r}}}^{{\rm{^{\prime} }}}(t),{{\bf{r}}}^{{\rm{^{\prime} }}{\rm{^{\prime} }}}(t),{{\bf{r}}}^{{\rm{^{\prime} }}{\rm{^{\prime} }}{\rm{^{\prime} }}}(t))}{||{{\bf{r}}}^{{\rm{^{\prime} }}}(t)\times {{\bf{r}}}^{{\rm{^{\prime} }}{\rm{^{\prime} }}}(t{)||}^{2}}.$$\end{document}$$

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 *s*(*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, speed, and curvature of the Medaka NSS 01 fish (Supplementary Video [1](#MOESM2){ref-type="media"}) over the first fifteen seconds of its experiment are shown in Fig. [4](#Fig4){ref-type="fig"} and in Supplementary Video [2](#MOESM3){ref-type="media"}.

Depending on the application, there may be additional factors that are relevant to behavior, such as position, posture, or distance to an object of interest. To accomodate such factors, we introduce the following notion:

#### **Definition 1.** {#FPar1}

Let $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\bf{r}}\,:[{t}_{0},{t}_{f}]\to {{\mathbb{R}}}^{3}$$\end{document}$ be the trajectory of an animal and let $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\Theta =\{{\theta }_{1},\ldots ,{\theta }_{m}\}$$\end{document}$ be a set of *m* behavioral factors, such as *s*, *κ*, and *τ*. The behavior curve of the trajectory indexed by Θ is the parametric curve $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${C}_{\Theta }\,:[{t}_{0},{t}_{f}]\to {{\mathbb{R}}}^{m}$$\end{document}$ given by $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${C}_{\Theta }(t)=({\theta }_{1}(t),\ldots ,{\theta }_{m}(t))$$\end{document}$.

Our objective is to find an acceptable notion of distance between the behavior curves of two episodes. Because it is virtually impossible to deliver an external chemical stimulus simultaneously to two freely swimming animals (as the stimulus takes time to disperse through the tank water), 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. To solve this curve alignment problem, we apply the *edit distance* developed by Sebastian *et al*.^[@CR16]^, which measures the least amount of energy required to bend, stretch, and/or compress one parametric curve onto another by finding an alignment that minimizes a specified energy. In our context, an *alignment* between two parametric curves $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${C}^{1}\,:[{t}_{0},{t}_{f}]\to {{\mathbb{R}}}^{{d}_{1}}$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${C}^{2}\,:[{t}_{0}^{\ast },{t}_{f}^{\ast }]\to {{\mathbb{R}}}^{{d}_{2}}$$\end{document}$ is a function$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\phi \,:[{u}_{0},{u}_{f}]\to [{t}_{0},{t}_{f}]\times [{t}_{0}^{\ast },{t}_{f}^{\ast }]\,{\rm{with}}\,\phi (u)=({\phi }_{1}(u),{\phi }_{2}(u))$$\end{document}$$satisfying the conditions that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\phi ({u}_{0})=({t}_{0},{t}_{0}^{\ast })$$\end{document}$, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\phi ({u}_{f})=({t}_{f},{t}_{f}^{\ast })$$\end{document}$, and *φ*~*i*~ is monotone increasing and differentiable for each *i* ∈ {1, 2}, and the *energy* we seek to minimize is given by$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$E(\phi )=\frac{1}{{u}_{f}-{u}_{0}}{\int }_{{u}_{0}}^{{u}_{f}}\sqrt{\mathop{\sum }\limits_{i=1}^{d}\,{({C}_{i}^{1}({\phi }_{1}(u))-{C}_{i}^{2}({\phi }_{2}(u)))}^{2}}\,du.$$\end{document}$$

We can now define the behavioral distortion distance:

#### **Definition 2.** {#FPar2}

Let $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${{\bf{r}}}_{1}\,:[{t}_{0},{t}_{f}]\to {{\mathbb{R}}}^{d}$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${{\bf{r}}}_{2}\,:[{t}_{0}^{\ast },{t}_{f}^{\ast }]\to {{\mathbb{R}}}^{d}$$\end{document}$ be the trajectories of two animals and let $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\Theta =\{{\theta }_{1},\ldots ,{\theta }_{m}\}$$\end{document}$ be a set of behavioral factors. The behavioral distortion distance between **r**~1~ and **r**~2~ with respect to Θ is$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${{\rm{BDD}}}_{\Theta }({{\bf{r}}}_{1},{{\bf{r}}}_{2})=E({\phi }_{{\rm{\min }}})$$\end{document}$$where *φ*~min~ is an alignment between the behavior curves $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${C}_{\Theta }^{1}$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${C}_{\Theta }^{2}$$\end{document}$ corresponding to **r**~1~ and **r**~2~, respectively, that minimizes *E*.

An example optimal alignment between the behavior curves with Θ = {*s*} for the first 15 seconds of observation for the Medaka NSS 01 fish and Medaka SS 01 fish is illustrated in Fig. [5](#Fig5){ref-type="fig"}. Supplementary Video [3](#MOESM4){ref-type="media"} is an animation example that illustrates this alignment over a 30 second observation interval of the same pair of animals.Figure 5An alignment between the speed curves of the NSS 01 and SS 01 medaka over the first 15 seconds of the experiment. (**A**) The optimal alignment *φ*~min~; (**B**) the speed curve for NSS 01; (**C**) the speed curve for SS 01; (**D**) the two curves plotted simultaneously with pairs of matched points indicated by gray line segments; (**E**) the two curves plotted with parameterized or "warped" time.

### Implementation {#Sec10}

We obtain the trajectories from the "track objects" method in the MetaMorph software as sequences of *x*- and *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 Mathuru^[@CR36]^, 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.

We use *dynamic time warping* (DTW) to calculate BDD, the same approach Sebastian *et al*.^[@CR16]^ use to calculate edit distance. Our implementation of the DTW algorithm is an extension of the standard DTW method in the Machine Learning module, mlpy^[@CR53]^, in Python (see Muller^[@CR54]^ for a description of how DTW works). Specifically, we modified the existing DTW method in mlpy to compare multivariate time series (instead of only univariate time series) by changing the distance function from 1-dimensional Euclidean distance to *m*-dimensional Euclidean distance for a specified *m*. While our method is able to handle input data of an arbitrarily high dimension (any number of behavioral factors), the analyses in this study have only used speed and curvature. The speed and curvature values are generated from the *x*- and *y*-coordinates using the formulas$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$s(t)=\sqrt{x^{\prime} {(t)}^{2}+y^{\prime} {(t)}^{2}}\,{\rm{and}}\,\kappa (t)=\frac{|x^{\prime} (t)y^{\prime\prime} (t)-y^{\prime} (t)x^{\prime\prime} (t)|}{v{(t)}^{3}},$$\end{document}$$where the derivatives are computed by applying the gradient function in the NumPy^[@CR55]^ Python package to the *x* and *y* time series. We smoothed the trajectories using the Savitzky-Golay function in the SciPy^[@CR56]^ Python package (see Savitzky and Golay^[@CR57]^ for a description of the Savitzky-Golay smoothing filter). That function takes two integer parameters, namely the order of the polynomials used for smoothing and length (number of frames) of the smoothing window centered at each frame. Since the choice of parameters is application dependent, we searched for the parameters that simultaneously minimized mean square error (deviation from raw trajectory) and jerk (rate of change of acceleration) for every order between 1 and 9 and every window length up to 600 frames (30 seconds). We encountered numerical instability for orders above 7. The products of mean square error and mean jerk with respect to window length are plotted for each order between 1 and 7 in Fig. [6](#Fig6){ref-type="fig"}. We used the order 5 Savitzy-Golay smoothing with a window length of 53 frames.Figure 6Products of mean square error and mean jerk of Savitzky-Golay smoothing with respect to window length for each odd order between 1 and 7. The plots for the pairs of orders 2 and 3, 4 and 5, and 6 and 7 were indistinguishable in the graph, so the even orders were omitted.

Furthermore, to account for variations in mean values of behavioral factors between animals and differences in magnitudes between behavioral factors themselves, we normalize all of the data using the sigmoid function$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rm{sig}}\,(x)=\frac{1}{1+{e}^{-(x-\mu )/\sigma }},$$\end{document}$$where *μ* and *σ* are the mean and the standard deviation of the respective input sequence. We chose the sigmoid function for our particular application because the distribution of the curvatures of the trajectories is positively skewed, as can be seen from Fig. [4](#Fig4){ref-type="fig"}. Since the range of the sigmoid function is \[0, 1\], the value of BDD~Θ~ is bounded above by the square root of the number of factors in Θ.

Surface alignment and conformal spatiotemporal distance (CSD) {#Sec11}
-------------------------------------------------------------

The second application we considered 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. To do this, we apply the *symmetric distortion energy* between surfaces introduced by Koehl and Hass^[@CR30],[@CR58]^.

### Theoretical framework and definitions {#Sec12}

To examine the overall spatiotemporal patterns of an animal within an enclosure, we study its spatiotemporal distribution function, that is, the continuous function $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\rho :\Omega \to {\mathbb{R}}$$\end{document}$ 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 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${{\mathbb{R}}}^{3}$$\end{document}$ with points {(*x*, *y*, *ρ*(*x*, *y*)\|(*x*, *y*) ∈ Ω)}, as illustrated in Fig. [7](#Fig7){ref-type="fig"}. We assume that this surface is smooth and Riemannian, that is, it is equipped with a well-defined notion of distance and angles.Figure 7The approximately uniform Delaunay triangulations corresponding to the heat maps of the NSS 01 and SS 01 medaka, their conformal flattenings to the unit disk, the symmetric distortion minimizing Möbius transformation from the unit disk to itself, and corresponding alignments.

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*^[@CR59]^, which states that 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 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\lambda }_{f}:{S}_{1}\to {\mathbb{R}}$$\end{document}$ 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. Koehl and Hass^[@CR58]^ introduced the *symmetric distortion energy* of a conformal diffeomorphism *f* : *S*~1~ → *S*~2~ with dilation functions $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\lambda }_{f}\,:{S}_{1}\to {\mathbb{R}}$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\lambda }_{{f}^{-1}}\,:{S}_{2}\to {\mathbb{R}}$$\end{document}$ as$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${E}_{sd}(f)=\sqrt{{\int }_{{S}_{1}}\,{({\lambda }_{f}(z)-1)}^{2}\,dA}+\sqrt{{\int }_{{S}_{2}}\,{({\lambda }_{{f}^{-1}}(z)-1)}^{2}\,dA}.$$\end{document}$$

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 3.** {#FPar3}

Let *S*~1~ and *S*~2~ be the surfaces of spatiotemporal distribution functions $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rho }_{1}\,:\Omega \to {\mathbb{R}}$$\end{document}$ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\rho }_{2}\,:\Omega \to {\mathbb{R}}$$\end{document}$, respectively, and for every diffeomorphism *f* : *S*~1~ → *S*~2~ and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$z\in \Omega $$\end{document}$ let *z*~*f*~ denote the unique point such that $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$f((z,{\rho }_{1}(z)))=({z}_{f},{\rho }_{2}({z}_{f}))$$\end{document}$. The *conformal spatiotemporal distance* between *S*~1~ and *S*~2~ is given by$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\text{CSD}({S}_{1},{S}_{2})=\frac{1}{2||\Omega ||}(\sqrt{{\int }_{\Omega }{({\rho }_{1}(z)-{\rho }_{2}({z}_{{f}_{\ast }}))}^{2}\,dA}+\sqrt{{\int }_{\Omega }{({\rho }_{2}(z)-{\rho }_{1}({z}_{{f}_{\ast }^{-1}}))}^{2}\,dA}\,)$$\end{document}$$where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$||\Omega ||$$\end{document}$ denotes the area of Ω and $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${f}_{\ast }\,=\,{{\rm{argmin}}}_{{\rm{conformal}}f}{E}_{sd}(f)$$\end{document}$.

### Implementation {#Sec13}

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 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${{\mathbb{R}}}^{3}$$\end{document}$. 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^[@CR60],[@CR61]^ 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^[@CR62]^. In particular, we use an algorithm developed by Collins and Stephenson^[@CR63]^ to find a discrete conformal mapping $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${f}_{i}:{S}_{i}\to {\mathbb{D}}$$\end{document}$ 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* : *S*~1~ → *S*~2~ as $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$f={f}_{2}^{-1}\circ \,m\,\circ \,{f}_{1}$$\end{document}$ where $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$m\,:{\mathbb{D}}\to {\mathbb{D}}$$\end{document}$ 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 $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$m\,:{\mathbb{D}}\to {\mathbb{D}}$$\end{document}$ whose composition $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${f}_{2}^{-1}\circ m\circ {f}_{1}$$\end{document}$ minimizes the symmetric distortion energy. This is illustrated in Fig. [7](#Fig7){ref-type="fig"}.

A well known application of Schwarz's Lemma^[@CR64]^ in complex analysis is that every Möbius transformation from $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${\mathbb{D}}$$\end{document}$ 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^[@CR58]^. 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^[@CR29]^. We use the following discrete analogue of the symmetric distortion energy presented in Koehl and Hass^[@CR58]^: 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$$\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$${E}_{sd}(f)\approx \sqrt{\sum _{{e}_{ij}\in {T}_{1}}\,{(\frac{\ell (f({e}_{ij}))}{\ell ({e}_{ij})}-1)}^{2}\frac{{A}_{ij}}{3}}+\sqrt{\sum _{{e}_{ij}\in {T}_{2}}\,{(\frac{\ell ({f}^{-1}({e}_{ij}))}{\ell ({e}_{ij})}-1)}^{2}\frac{{A}_{ij}}{3}}$$\end{document}$$where *e*~*ij*~ denotes an edge in *T*~1~ or *T*~2~, $\documentclass[12pt]{minimal}
                \usepackage{amsmath}
                \usepackage{wasysym} 
                \usepackage{amsfonts} 
                \usepackage{amssymb} 
                \usepackage{amsbsy}
                \usepackage{mathrsfs}
                \usepackage{upgreek}
                \setlength{\oddsidemargin}{-69pt}
                \begin{document}$$\ell $$\end{document}$ 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.

Statistical analyses {#Sec14}
--------------------

All of the Kruskal-Wallis H, Mann-Whitney U, and Wilcoxon T tests presented in the Results section were executed in SciPy^[@CR56]^. A p value of 0.01 was considered significant. When multiple comparisons were made Bonferroni correction for multiple comparisons was applied to determine statistial significance. When reporting, the P values up to the third decimal place are presented instead of an asterisk sign.

Supplementary information
=========================

 {#Sec15}

Supplementary Table S1 Supplementary Video S1 Supplementary Video S2 Supplementary Video S3

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information
=========================

is available for this paper at 10.1038/s41598-019-52300-8.

This work was supported by Yale-NUS College through grants R-607-265-225-121 and IG16-LR003. The authors acknowledge the support from Haroun Chahed, Goh Rui Zhe and Maharshee Karia in the initial stages of the project.

A.S.M. and M.T.S. conceived and designed the project, and wrote the manuscript with assistance from S.G. A.S.M. provided the data-set to test the methods and interpreted the results. M.T.S. formulated the models and conducted the computational experiments. S.G. developed and implemented the software package with assistance from M.T.S.

The authors declare no competing interests.
