Deep learning-assisted comparative analysis of animal trajectories with DeepHL
DeepHL system architecture
The DeepHL system consists of three server computers. The first one is a web server that receives a trajectory data file from a user and provides analysis results to the user (Intel Xeon E5-2620 v4, 16 cores, 32 GB RAM, Ubuntu 14.04). The second one is a storage server that stores data files and analysis results. The third one is a GPU server that analyzes data provided by the user (Intel Xeon E5-2620 v4, 32 cores, 512 GB RAM, four NVIDIA Quadro P6000, Ubuntu 14.04). Supplementary Information, Algorithm, provides a complete description of the DeepHL method. DeepHL is accessible on the Internet through http://www-mmde.ist.osaka-u.ac.jp/maekawa/deephl/. Supplementary Information, User guide to DeepHL, provides a user guide to DeepHL. In addition, Supplementary Information, Usage of Python-based Software, and Supplementary Software 1 present the Python code of DeepHL.
Preprocessing
An input trajectory is a series of timestamps and X/Y coordinates associated with a class label. To perform position- and rotation-independent analysis, we convert the series into time series of speed and relative angular speed and then standardize them (Supplementary Information, Algorithm). Note that the absolute coordinates of wild animals, which can relate to the distance from a nest or feeding location, for example, are important in understanding behavior of the animals. Hence, DeepHL allows the original coordinates to be input to DeepHL-Net along with the speed and relative angular speed. In addition, other biological time-series sensor data measured by the user can be fed into DeepHL-Net when these time-series data are included in a data file uploaded by the user. For example, a time series of the heading direction of animals obtained from digital compasses can be useful for behavior understanding. Moreover, primitive features usually used in trajectory analysis can be easily fed into DeepHL-Net. DeepHL automatically computes the travel distance from the initial position, the straight-line distance from the initial position, and the angle from the initial position (Supplementary Table 1) as primitive features. Using the web interface of DeepHL, the user can easily select primitive features and other sensor data to be fed into DeepHL-Net (Supplementary Information, User guide to DeepHL). See Supplementary Information, Effect of input features, for effects of input features on classification accuracy. Normally, the inputs of DeepHL-Net are two-dimensional time series, that is, speed and relative angular speed. When we input an additional time series (such as the original coordinates) into DeepHL-Net, the additional time series are added as additional dimensions of the inputs.
Multi-scale layer-wise attention model (DeepHL-Net)
Here, we explain DeepHL-Net shown in Fig. 2f in detail. The input of the model is a time series of primitive features, that is, an lMAX × Nf matrix, where lMAX is the maximum length of the input trajectories and Nf is the dimensionality of the time series, that is, the number of the primitive features. Because the lengths of observed trajectories are not identical to each other in many cases, we fill in missing elements in the matrix with −1.0 and mask them when we train DeepHL-Net. In each 1D convolutional layer of the convolutional stacks, we extract features by convolving input features through the time dimension using a filter with a width (kernel size) of Ft. We use different filter widths in the four convolutional stacks (3%, 6%, 9%, and 12% of lMAX) to extract features at different levels of scale. We use a stride (step size) of one sample in terms of the time axis. We also use padding to allow the outputs of a layer to have the same length as the layer inputs. In addition, to reduce an overfitting, we employ a dropout, which is a simple regularization technique in which randomly selected neurons are dropped during training44. The dropout rate used in this study is 0.5.
In each LSTM layer of the LSTM stacks, we extract features considering the long-term dependencies of the input features. LSTM is a recurrent neural network architecture with memory cells, and it permits us to learn temporal relationships over a long time scale. LSTM learns long-term dependencies by employing memory cells that hold past information, updating the cell state using write, read, and reset operations with input, output, and forget gates (see Supplementary Information, Algorithm). In addition, we employ dropout to reduce overfitting. The attention information of each layer is computed by using Eq. (1), and then it is multiplied by the layer output. Here, the softmax and tanh functions in Eq. (1) are defined as follows:
$$,{text{softmax}},({x}_{j})=frac{exp ({x}_{j})}{{sum }_{i}exp ({x}_{i})},$$
(2)
$$tanh ({x}_{j})=frac{exp ({x}_{j})-exp (-{x}_{j})}{exp ({x}_{j})+exp (-{x}_{j})}.$$
(3)
Note that parameters in Eq. (1) for each layer, that is, Wa and ba, as well as parameters in the convolutional and LSTM layers are estimated during the network training phase. Here, we introduced the tanh activation function into Eq. (1) to smooth out the output attention values. When an outlying large value is included in WaZT + ba at time t, attention values other than time t become extremely small without using the tanh function. When we visualize a trajectory using such attention values, only a single data point is colored in red, making it difficult for a user to identify important segments.
Training and testing of DeepHL-Net
The DeepHL user can select the parameters of DeepHL-Net used in the analysis, that is, the number of convolutional/LSTM layers and the number of neurons in each layer (default: four layers with 16 neurons). Then, DeepHL-Net is trained on 80% of randomly selected trajectories to minimize the binary classification error of the training data, employing backpropagation based on Adam45 (Supplementary Information, Algorithm). (Note that each trajectory has a class label for binary classification.) Then, the trained DeepHL-Net is tested using the remaining 20% of trajectories to compute the classification accuracy, providing an indication of the degree of difference between the two classes.
Computing the score of each layer
To screen the layers in DeepHL-Net, we compute a score for each layer according to Eq. (4)
$$s({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})={s}_{mathrm{fc}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})+{s}_{mathrm{it}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}}).$$
(4)
Here, ({A}_{i,{C}_{mathrm{A}}}) is a set of attention vectors calculated from trajectories belonging to class A using the ith layer. In addition, ({A}_{i,{C}_{mathrm{B}}}) is a set of attention vectors calculated from trajectories belonging to class B using the ith layer. As mentioned in the main text, an attention vector from a discriminator layer should have large values within limited segments. Therefore, ({s}_{mathrm{fc}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})) in Eq. (4) calculates the averaged variance of the attention values normalized by the average length of the trajectories, as described in Eq. (5). When the layer focuses on a part of a trajectory, the variance increases
$${s}_{mathrm{fc}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})=sqrt{frac{1}{| {A}_{i,{C}_{mathrm{A}}}cup {A}_{i,{C}_{mathrm{B}}}| cdot l({A}_{i,{C}_{mathrm{A}}}cup {A}_{i,{C}_{mathrm{B}}})}sum _{{bf{a}}in {A}_{i,{C}_{mathrm{A}}}cup {A}_{i,{C}_{mathrm{B}}}}V({bf{a}})}.$$
(5)
Note that V(⋅) calculates the variance and l(⋅) calculates the average length of the trajectories. We take the square root of the average variance to derive the average standard deviation. Using (l({A}_{i,{C}_{mathrm{A}}}cup {A}_{i,{C}_{mathrm{B}}})), which calculates the average length of ({A}_{i,{C}_{mathrm{A}}}cup {A}_{i,{C}_{mathrm{B}}}), we normalize the computed variance. Because the softmax function in Eq. (1) ensures that all values sum to 1, resulting in a larger variance for longer trajectories, we normalize the average variance using the average length.
In addition, as mentioned in the main text, the distribution of attention values by the layer for one class should be different from that for another class. Therefore, ({s}_{mathrm{it}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})) calculates the difference between the distributions of the attention values of classes A and B as follows:
$${s}_{mathrm{it}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})=(1-,{mathrm{Intersect}},(h(A_{{i,{C}}_{mathrm{A}}}),h({{A}}_{{i,{C}}_{mathrm{B}}}))).$$
(6)
Here, h(⋅) calculates a normalized histogram of attention with 200 bins, and Intersect(⋅ , ⋅) calculates the area overlap between two histograms, and is described as follows:
$${mathrm{Intersect}},(H_{1},H_{2})=mathop{sum}limits_{i}min (H_{1}(i),H_{2}(i)),$$
(7)
where H1(i) shows the normalized frequency of the ith bin of histogram H1. As described in Eq. (4), the final score is calculated as the sum of the two scores of ({s}_{mathrm{fc}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})) and ({s}_{mathrm{it}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})).
Here, ({s}_{mathrm{fc}}({A}_{i,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}})) in Eq. (4) is used to find a layer that focuses only on a portion of a trajectory. Owing to the term, only a small important portion of trajectories is highlighted in many cases, as shown in Figs. 3, 5, and 6, especially for the trajectories of beetles. However, substantial portions of several trajectories of the normal mice are highlighted, as shown in Fig. 4d. Because the characteristics of the normal mouse trajectories are the distance from the initial position, the segments in the trajectories far from the initial position are highlighted.
Computing the correlation between attention values and handcrafted features
To help the user understand the meaning of the highlights, DeepHL automatically computes the Pearson correlation coefficients between the attention values of each layer and handcrafted features computed by DeepHL, as shown in Supplementary Table 1. In addition, the correlation coefficients with sensor data and handcrafted features included in a trajectory data file are automatically computed. Computing the correlation with environmental sensor data can reveal the relationship between a behavior and environmental conditions. If a specific behavior is exhibited only when the temperature is high, for example, we can infer that the behavior relates to the high temperature condition. Furthermore, DeepHL automatically computes the moving average, moving variance, and derivative of each of the above features/sensor data, and then computes the correlation coefficients with the attention values, which are presented to the user (Supplementary Fig. 1).
Computing the difference between distributions of each handcrafted feature for the two classes within highlighted segments
To help the user understand the meaning of the highlights, DeepHL automatically computes the difference between distributions of each handcrafted feature for two classes within highlighted segments. The difference is computed as follows:
$${mathrm{diff}}({A}_{i,{C}_{mathrm{A}}},{F}_{j,{C}_{mathrm{A}}},{A}_{i,{C}_{mathrm{B}}},{F}_{j,{C}_{mathrm{B}}})=1-,{mathrm{Intersect}},(h(m({{A}}_{{i,{C}}_{mathrm{A}}},{{F}}_{{j,{C}}_{mathrm{A}}})),h(m({{A}}_{{i,{C}}_{mathrm{B}}},{{F}}_{{j,{C}}_{mathrm{B}}}))).$$
(8)
Here, ({F}_{j,{C}_{mathrm{A}}}) is a set of time series of the jth handcrafted feature calculated from trajectories belonging to class A. In addition, m(⋅ , ⋅) is a masking function that extracts feature values within highlighted segments. Because the softmax function in each attention layer ensures that all attention values in a sum of 1, we consider an attention value larger than c/(# time slices) as a potential attended value (c = 1.2 in our implementation).
Data acquisition of worms
Data acquisition was performed according to Yamazoe-Umemoto et al.22. In brief, several worms were placed in the center of an agar plate in a 9-cm Petri dish, 30% 2-nonanone (v/v, EtOH) was spotted on the left side of the plate, which was covered by a lid and placed on the bench upside down. Then, the images of the plate were captured with a high-resolution USB camera for 12 min at 1 Hz. Because the worms do not exhibit odor avoidance behavior during the first 2 min because of the rapid increase in odor concentration46, the data for the following 10 min (i.e., 600 s) was used. From the images, individual worms were identified and the position of the centroid was recorded by an image processing software Move-tr/2D (v. 8.31; Library Inc., Japan). The number of recorded trajectories is 325 (Supplementary Table 2). The comparison was between the naive worms (control class) and the worms after preexposure to the odor (preexposed class).
DeepHL analysis of worms
A multivariate time series of movement speed, relative angular speed, distances from the initial position, and angle from the initial position extracted from the time series of trajectories was fed into DeepHL-Net, yielding a binary classification accuracy of 93.9%, where 20% of the data are used as test data. The discriminator layer used in this investigation has the highest score of all layers. As shown in Fig. 3d, which was calculated from the moving variance of the speed within highlighted segments, we can state that the changes in the speed of preexposed worms is larger than those of control worms. Figure 3e shows spectrograms of the speed calculated from entire trajectories (Fig. 3c) with a 128-s wide sliding window shifted in 1-sample intervals. In addition, Fig. 3f shows histograms of the dominant frequency of speed calculated from entire trajectories using the 128-s wide sliding window shifted in 1-sample intervals. These results also indicate the difference in the frequency of speed between the preexposed and control worms. Our investigation revealed that the dominant frequency of speed significantly differs between the preexposed and control worms using GLMM with Gaussian distributions (t = −6.60; d.f. = 322.8; p = 1.68 × 10−10, effect size(r2) = 0.232). The p value is two sided. Individual factors were treated as random effects. The number of data points for the control class is n = 76, 784 and that for the preexposed class is n = 75, 750. We used GLMM with Gaussian distributions because the objective variable has a continuous value and we used the lmerTest package (v. 2.0–36) of R (v. 3.4.3) for the analysis.
Data acquisition of mice
We collected 52 trajectories of normal mice and unilateral 6-hydroxydopamine (OHDA) lesion mouse models of PD while they freely moved for 10 min in an open field (60 × 55 cm2, wall height = 20 cm; normal: 22, PD: 30). The trajectories were detected by the animal’s head position, which was captured by an overhead digital video camera (60 fps). Two sets of small red and green light-emitting diodes were mounted above the animal’s head so that it could be located in each frame. Custom softwares based on Matlab (R2018b, Mathworks, MA, USA) and LabVIEW (Labview 2018, National Instruments, TX, USA) were used for tracking. We then created 30-s segments by splitting each trajectory because training a DNN requires a number of trajectories. We used 966 segments in total (normal: 374, PD: 592) collected from nine C57BL/6J mice (normal: 5, PD: 4). Note that we excluded 30-s segments that contain no movements of a mouse.
DeepHL analysis of mice
Movement speed, relative angular speed, travel distances, straight-line and travel distances from the initial position, and angle from the initial position were fed into our model. The accuracy for the binary classification of normal and 6-OHDA model mice was 74.7%, where 20% of the data are used as test data. The score of the discriminator layer was the highest of all LSTM layers and the sixth highest of all layers. Our investigation revealed that the behavior of visiting locations far away from the initial position can be characteristic of normal mice.
To evaluate PD symptoms from animal behaviors, previous studies have exclusively focused on the movement speed of animals in the open-field tests (frequency and bout duration of ambulation as well as immobility or fine movement) because typical symptoms in the animal model of PD are thought to be slowness of movement and a paucity of spontaneous movements. As shown in Fig. 4e–g, we found significant differences in average movement speed during ambulation periods, average movement speed during fine movement periods, and average maximum distance within a ±60-s window in a session. These differences were derived from the findings of DeepHL using the two-sided Wilcoxon rank-sum test (W = 544, p = 3.486 × 10−5, effect size (Cliff’s delta) = −0.648; W = 511, p = 5.869 × 10−4, effect size (Cliff’s delta) = −0.548; W = 521, p = 2.666 × 10−4, effect size (Cliff’s delta) = −0.579). The 95% confidence intervals are [1.222, 3.481], [0.139, 0.468], and [13.726, 43.175], respectively. We used the exactRankTests package (v. 0.8–29) of R (v. 3.2.3). Note that these behavioral features are extracted from original 10-min trajectories.
The maximum distance, which was derived from a finding of DeepHL, is more useful for evaluating the PD symptoms than conventional measures based on the movement speed. Note that the new feature is designed based on an insight drawn from an analysis by deep learning. These results suggest that DeepHL helps find a novel measure not directly linked to the movement speed, that is, a straight-line distance within a certain time window. When the aim of an animal is to visit all locations in an area, the travel distance over a short duration commonly becomes longer. Besides, it is well known that rodents, including mice and rats, spontaneously prefer to explore an environment, particularly in novel places. Thus, DeepHL may capture the fact that the abnormal behavior of the 6-OHDA lesion model of PD hinders such spontaneous behavioral traits of normal mice. Indeed, the 6-OHDA lesion mouse model appears to remain in the same place. Although this hypothesis should be verified based on the causality between behavioral traits and neural activity patterns underlying PD symptoms using neuronal recording together with its optogenetic manipulation in the basal ganglia and motor cortex23, it is beyond the scope of this study.
Behavioral features of mice
According to Kravitz et al.23, ambulation was defined as periods when the velocity of the animal’s center point averaged >2 cm/s for at least 0.5 s. Immobility was defined as continuous periods of time during which the average change of the trajectory was More
