Caenorhabditis elegans is a free-living soil nematode that has been a workhorse genetic model system. The nematode’s transparent tissue, simple anatomy, and fast reproduction contribute to both ease in culture and a literal window into the internal workings of a living organism. Its completely sequenced genome contains many genes that are homologous to human genes, and importantly the ability to manipulate genes with relative ease makes it an extremely attractive model system. For neuroscience in particular, C. elegans presents a unique opportunity with its simple nervous system (just 302 neurons) that is complex enough to exhibit many sensory modalities, including mechanosensation, chemosensation, and response to heat, osmolarity, and smell.

Behavior characterization in C. elegans was historically qualitative, mainly relying on experimentalists specifying end-point assessment (e.g. whether the worm chemotaxes to a particular source of odor within a certain amount of time), or experimentalists using heuristics to assess behavior (e.g. naming worms genes “unc” for uncoordinated). In the last decade, machine vision tools first replaced human identifications of worms from images and videos, which allows much larger dynamic datasets to be annotated and analyzed. In recent years, further development in quantitative behavior characterization tools such as tracking (Stirman et al., 2011; Swierczek et al., 2011; Husson et al., 2012; Yemini et al., 2013; Porto et al., 2019), eigenworms (Stephens et al., 2008), behavior “dictionaries” (Brown et al., 2013), and t-SNE (Berman et al., 2016; Liu et al., 2018) have moved the field away from merely describing the outcome to understanding the types of behavior the brain of this simple system can generate. While many of these techniques do well in quantitatively describing behavior and distinguishing differences in behavior, behavioral dynamics are rich and opportunities abound in exploring behavioral dynamics using other mathematical tools.

Persistent homology has been used to analyze time series data in many different settings. Some earlier work was theoretical and studied the interaction between persistence and sliding window embeddings—which we used in this research—as well as proposed possible applications (Firas and Elizabeth, 2015; Perea and Harer, 2015; Perea, 2016). Research into gene expression has used persistent homology to detect patterns or classify whether a signal is periodic (Dequéant et al., 2008; Perea et al., 2015). Frequently, persistence has been used to study neural data (Petri et al., 2013; Backholm et al., 2015; Stolz et al., 2017; Sizemore et al., 2019; Helm et al., 2020; Lütgehetmann et al., 2020), and in many cases neural data from C. elegans, but the analysis tends to rely on clique complexes as the topological space of interest instead of sliding window embeddings.

C. elegans (N2 strain) were cultured at 20°C under standard conditions on agar plates seeded with OP50 E. Coli. Animals were age-synchronized via hatch-off and cultured on plate until they reached day 1 of adulthood. For behavior experiments on agar, animals were prepared, imaged, and tracked as previously described (Porto et al., 2019). For behavior experiments in methylcellulose media, synchronized populations were then washed off of culture plates with M9 buffer. Unless otherwise noted, video data was collected on a dissecting microscope (Leica MZ16) using a CMOS camera (Thorlabs DCC3240M), with a frame rate of 30 frames per second and a magnification of ×1.2.

Behavior data was collected with animals confined with microfluidic devices. In these devices the cavities in which worms are loaded have only slightly greater depth than the width of an adult worm, which restricts worms to the focal plane of the microscope and to almost entirely 2-dimensional behavior. Microfluidic devices were fabricated as described previously (Chung et al., 2011). Methylcellulose solutions were prepared at concentrations of 0.5 % ,   1 % ,   2 % ,   and   3 % weight in volume of M9 buffer. To ensure that single animals could be isolated in single chambers of the unbonded microchamber microfluidic device, we first picked animals onto a room-temperature, unseeded plate. To ensure that animals were fully immersed in methylcellulose mixture, we used a glass pipet to aspirate a small amount of methylcellulose solution, and then aspirated animals from the unseeded plate one at a time into the methylcellulose solution. Then, single animals surrounded by methylcellulose mixture were pipetted into individual chambers of an unbonded PDMS chamber device. The device could then be flipped over onto a sterile 10 cm Petri dish and gently pressed down until the individual chamber walls came into contact with the Petri dish, preventing animals from leaving their chambers. Animals were then imaged in devices for about 5 min at 30 frames per second, resulting in time series data with 10,665 points.

To extract midline data from videos, we first found masks for each frame to isolate the worm from the background using a combination of Otsu thresholding (Otsu, 1979), image smoothing using a Gaussian kernel, and size filtration. Otsu thresholding is a thresholding algorithm based on the gray-level histogram of an image. The threshold is identified by the grayscale pixel value that minimizes the intra-class variance of background and foreground pixels. We then broadly followed the method used in Stephens et al. (Stephens et al., 2008) to represent the worm’s posture in “worm-centric” coordinates. Briefly, we found the midline of the worm in each frame by thinning the mask to a single line and interpolating between pixels of this line such that the midline was represented by 101 evenly spaced points. We calculated the tangent angle between each pair of adjacent points along the midline so that the animal’s posture could be represented as a vector of angles, and then transformed those vectors with PCA so users could balance accuracy requirements and resource limitations via truncation of the data. We replaced frames in which animals were self-occluded with the data from the most recent non-self-occluded frame. We used untruncated PCA data for most computations because it has the same persistence output as the raw angle data. We used truncated PCA data (the first five principal components) for the computations in Section 3.1.

The videos for this study were selected from a much larger set of data based on how well they could be segmented and skeletonized. Some videos have subsequences that are difficult to automatically skeletonize because the animals self-occlude, i.e. bend in such a way as to cross over themselves. Thus, this dataset is likely biased toward less complex behaviors like thrashing and there are some cases where there are multiple videos of the same animal. The resulting data set has 40 samples of 10,665 points each with 10 samples for each viscosity condition.

Sliding window embeddings turn time series data into point cloud data in a way that does not forget the temporal information of the time series. There are some additional benefits to sliding window embeddings, including that they “separate” points that intersect each other in a time series, such as in Examples 2.6 and 2.7.

Definition 2.1:

A time series is a sequence of vectors ( x t ) t ∈ T   =   ( x t , x t + 1 , x t + 2 , … ) where each x t is in the same finite-dimensional vector space V and T is a totally ordered set.

Remark 2.2:

The totally ordered set T, which indexes the time series, can be ℤ , ℕ , or a finite set like [ N ]   =   { 1,2 , … , N } . For many applications including the ones in this paper, the indexing set is finite and will be omitted in notation for brevity, as in ( x t ) t .

Given any time series we can construct a new time series called a sliding window embedding, which is also known as a time delay embedding with a lag or delay time of 1.

Definition 2.3:

Given a time series τ   =   ( x t ) t with vectors x t   ∈   V , a sliding window embedding of window length l of τ is a new time series, τ l   =   ( x t ˜ ) t , with x t ˜   =   [ x t   x t + 1   …   x t + l − 1 ]   ∈   V l , where [ ⋅ ] is concatenation of vectors.

That is, the t t h vector in the new time series is the concatenation of l consecutive vectors in the original time series and has dimension equal to l ⋅ dim   ( V ) .

Remark 2.4:

If the original time series has N points, then the sliding window embedding of window length l has N − l + 1 points, as one can see in Example 2.5 .

Example 2.5:

Consider the time series τ   =   ( [ 1,2 ] , [ 3,4 ] , [ 5,6 ] , [ 7,8 ] , [ 9,10 ] ) in ℝ 2 . The sliding window embedding of τ of window length l   =   3 is τ 3   =   ( [ 1,2,3,4,5,6 ] ,   [ 3,4,5,6,7,8 ] , [ 5,6,7,8,9,10 ] )   ⊆   ℝ 6 , which has 5 − 3 + 1 = 3 points.

We applied persistent homology (Section 2.3) to sliding window embeddings of C. elegans video data in order to quantify behavior. Degree one persistent homology detected cycles in these sliding window embeddings which we show correspond to particular behaviors.

The cycles that persistent homology detects may consist of collections of points that trace out a closed curve. A cycle is “large” or highly persistent if it encloses an area that could fit a large ball; a cycle that is tall and skinny has small persistence.

Below we see two examples where a time series exhibits a single periodic behavior but persistent homology will detect either two or zero non-trivial cycles. In contrast, the persistent homology of a sliding window embedding detects exactly one non-trivial cycle in both examples.

Example 2.6:

Figure 1A displays one period of a periodic time series in ℝ 2 with the property that if successive points are connected by line segments then the path of the time series self-intersects. To discover this figure-eight-shaped loop, one might try to use persistent homology (Section 2.3). However, persistent homology would detect two distinct loops, each comprising half of the period. See Figure 2D for an illustration of these loops.

Figures 1B,C show two-dimensional PCA projections of sliding window embeddings of the figure-eight for l = 10 and l = 20 , respectively, using the first and third principal components. In these point clouds the time series draw out simple closed curves, and in fact in each of these cases persistence detects a single loop.

Notice that as the window length increases, the “size” of the loop increases. This increase in the loop’s persistence makes it easier for persistent homology to robustly detect it.

Example 2.7:

Figure 3A shows a 1-dimensional time series that is a discretization of a sine wave. This periodic behavior creates no loops — in fact, because the points take values in ℝ , the time series cannot produce degree 1 homology. However, a sliding window embedding, in this case of window length 4, creates a loop that is detected by persistent homology. That loop in ℝ 4 is projected down to two dimensions in Figure 3B.

In this section we provide an overview of persistent homology and how it may be used to produce quantitative summaries of the shape of a collection of points such as the sliding window embedding discussed above.

Definition 2.9:

A simplicial complex on a set of vertices V is a collection K of non-empty subsets of V such that if τ   ∈   K and τ ′   ⊂   τ , then τ ′   ∈   K . An element τ   ∈   K is called a simplex. An n-simplex or simplex of dimension n is a simplex τ   ∈   K with size | τ |   =   n + 1 . The 1-skeleton of a simplicial complex K is the set of simplices with dimension at most one. A filtered simplicial complex or filtration is a collection { K r } r ∈ S of simplicial complexes K r where S   ⊆   ℝ such that K r   ⊆   K s for all r , s   ∈   S with r   ≤   s .

Definition 2.10:

Let X   ⊂   ℝ d be a finite set and let r ≥ 0 . The Vietoris-Rips complex of X at scale r, denoted ℛ r ( X ) , is the simplicial complex with vertex set X and whose simplices are given as follows. A subset { x 0 , … , x n }   ⊂   X is an n-simplex in ℛ r ( X ) if and only if | x i − x j |   ≤   r for all i , j   ∈   { 0 , … , n } .

Definition 2.11:

The Vietoris-Rips filtration of a finite set X   ⊂   ℝ d is the collection ℛ ( X )   : =   { ℛ r ( X ) } r ≥ 0 .

Note that while the Vietoris-Rips complex of X is parameterized by the non-negative reals, the finiteness of X guarantees that ℛ ( X ) consists of only finitely many distinct simplicial complexes.

Example 2.12:

Figure 2A shows the 1-skeleton of the Vietoris-Rips complex of a pointcloud in ℝ 2 at four scales. Notice that each simplicial complex includes into the next.

The persistent homology of a Vietoris-Rips filtration can be represented by a multiset in ℝ 2 called a persistence diagram in which each point gives the scale of the appearance and disappearance of a topological feature (such as a loop) in the filtration.

Example 2.13:

Figures 2B,C show the persistent homology in degree 1 of Example 2.12. Notice that both Figure 2B—the persistence diagram—and Figure 2C—the persistence landscape (see Definition 2.14)—show two cycles, but because they are born and die at exactly the same radius parameters they are plotted in the same place. The two cycles are shown in Figure 2D.

It is difficult to apply standard tools of statistics and machine learning directly to persistence diagrams, which, for example, need not have unique averages (Mileyko et al., 2011). A solution is to map persistence diagrams into a vector space or Hilbert space. One such mapping is the persistence landscape. See (Bubenik, 2015) for the following definitions and results.

Definition 2.14:

For a   <   b let f a , b   :   ℝ   →   ℝ be the piecewise-linear function given by f a , b ( t )   =   { t − a ,   if   a ≤ t ≤ a + b 2 b − t ,   if   a + b 2 ≤ t ≤ b 0 ,   otherwise. Given a persistence diagram Dgm p ( K ) , the corresponding kth persistence landscape is the function λ k   :   ℝ   →   ℝ given by defining λ k ( t ) to be the kth largest value of f a , b ( t ) over all points ( a , b ) ∈ Dgm p ( K ) . The persistence landscape is the sequence ( λ k ) k . The parameter k is called the depth of the persistence landscape. For a point cloud X, we will denote by PL   ( X ) the persistence landscape obtained by applying degree 1 persistent homology to the Vietoris-Rips filtration of X.

Persistence landscapes have unique averages, satisfy the law of large numbers and central limit theorems, and can be discretized for computations. Because the sequence of functions that make up a landscape are nested, they can all be graphed on the same plot as in the right column of Figure 4.

While the persistence landscape is defined to be an object in a space of continuous functions, it can be discretized and turned into a finite-dimensional vector. Through discretization, each depth of the landscape transforms from a continuous function on ℝ to a vector where the i^th entry in the vector corresponds to the function value at the i^th discrete parameter value. The vectors for each depth of the landscape are concatenated together to produce a single high-dimensional vector. These discrete landscapes can be computed directly, which we did for the computations outlined in Section 2.4.

This vector space (in fact, Hilbert space) setting lets us use linear algebra-based statistical and machine learning techniques such as principal component analysis (PCA). The principal components from PCA on discretized landscapes can be converted into a format much like a persistence landscape—a sequence of continuous functions on ℝ —but the principal components are not themselves persistence landscapes because the functions fail to be nonnegative.

In this section we give details for our analysis of C. elegans data. The input consist of piecewise linear midlines of C. elegans from video recordings as described in Section 2.1. These midlines were parameterized by the 100 angles between adjacent segments and then were transformed using PCA, so each sample input to our system was a time series τ of 100-dimensional vectors measured in radians. See Figure 1 in (Stephens et al., 2008) and the accompanied description for more details on this parameterization of the C. elegans midlines or see our short summary of the procedure in Section 2.1.

The time domain of this time series was divided into overlapping patches of a given size called the patch length, resulting in a collection { τ i } i of smaller time series. For our experiments, a patch length of 300 was chosen, with adjacent patches overlapping by half of the patch length. The sliding window embeddings of the τ i were then computed with window length parameter l   =   20 , resulting in a new collection { ( τ i ) l } i of time series of length 300   −   l   +   1   =   281 . This analysis is not particularly sensitive to the choices of the hyperparameters patch length and window length; only extreme changes in either parameter lead to significant changes in results. The hyperparameter choices were motivated by the timescales at which C. elegans complete meaningful behaviors: for patch length, 150 frames of 30 fps video is 5 s of behavior; for window lengths, 20 frames is 0.67 s and corresponds to roughly one period of forward crawling in adult C. elegans submerged in the 0.5 % methylcellulose environment. Strategies for choosing appropriate window lengths are described in (Perea and Harer, 2015). The method for cross validation of the choice of window length is described at the end of this section.

Persistence diagrams were computed for each of the patches ( τ i ) l . This step accounts for the vast majority of the computational resources of the pipeline, and the computational costs are made worse by the concatenation of vectors in a sliding window embedding. This is where we greatly benefit from the preprocessing that turns video data, which is extremely high-dimensional (see (Tralie and Perea, 2018) Section 3.1), into a 100-dimensional time series. On a 2017 15-inch MacBook Pro with a 2.8 GHz Intel Core i7 processor and 16 GB of RAM, this step took 22588.632 s, or about 6 h and 15 min.

For each ( τ i ) l , a (discretized) persistence landscape PL   ( ( τ i ) l ) was computed from the persistence diagram of the Vietoris-Rips complex ℛ ( ( τ i ) l ) . The grid of parameter values on which the persistence landscape was evaluated to produce its discretization was chosen to include all of the bars of the barcode and to be sufficiently fine to produce nice visualizations. Since the persistence landscapes are piecewise linear with slope bounded by ± 1 , the step size of this discretization bounds the error and there is eventually little to be gained from a finer discretization. The step size of the discretization we used was 0.1. The maximum depth was chosen to include all nonzero depths of the persistence landscapes.

The collection { PL   ( ( τ i ) l ) } i of persistence landscapes was then assembled into a single summary for a given video by averaging the persistence landscapes across patches to result in a single average persistence landscape associated to each video. These steps are summarized in Figure 5. For each environment viscosity, the average persistence landscapes for each video were averaged to give an average persistence landscape of the class.

The persistence landscapes for each class were used for analysis of the sliding window embedding as follows. Distances between each class’ average persistence landscapes were computed using the usual Euclidean distance. The pairwise distances were visualized using multidimensional scaling to give a 2-dimensional visualization of the similarities between the classes; see Figure 6A.

Principal component analysis (PCA) was applied to the set of average persistence landscapes for each video and the first two principal components were plotted as sequences of functions. We plotted the PCA projection of these video average persistence landscapes together with the average of each class; see Figure 7. These plots visualize some of the similarity between classes and variation within classes.

Next, we further studied the variation within classes in two ways. First, the standard deviations of each coordinate were computed for the average persistence landscapes of the videos in each class. These were visualized as sequences of functions in Figure 8B to show the variation in different parts of the average persistence landscapes. Second, we applied PCA to the average persistence landscapes of the videos in each class and plotted the cumulative variances explained by the first n principal components for n = 1 , … , 10 (10 is the number of samples in each class). The first three principal components were also computed. See Figure 9.

We conducted a permutation test on the pairwise Euclidean distances between the average persistence landscapes of each videos. We used 10,000 permutations for each permutation test. The approximate p-value equals the percentage of cases in which the distance is at least as large as the observed distance. Results can be found in Table 2A.

We applied multiclass support vector machines (SVM) to classify samples according to viscosity of their environments. We use the ksvm function from the kernlab package in R on the average persistence landscapes of the videos. Accuracy was estimated using 10-fold cross validation with cost set to 10. Cross validation was repeated 20 times and the results were averaged. A confusion matrix for one instance of SVM with 10-fold cross validation was also computed and is shown in Table 2B.

We used support vector regression (SVR) to approximate viscosity of the worm environments given the worm’s behavior data. The goal was to assess how predictive our techniques are. Accuracy was estimated by averaging 10 repetitions of 10-fold cross validation. Results are plotted in Figure 9.

As a final step we applied cross validation to the hyperparameter window length. We cross-validated the choice of window length by running the above pipeline for window lengths of 1, 10, 20, and 30 on a subset of the data. To reduce total computation time, we restricted to the first minute of each video, which corresponds to 1800 frames. We then compared the permutation test and multiclass SVM results to see which hyperparameter choice gave the best results. The results from cross validation of window length can be found in Section 3.2.1.

To help validate our pipeline, we compared our results to those obtained by applying the same computational procedure to a null model given by randomly permuting the frames in each video.

We also studied the effects of using a preprocessing step different from sliding window embeddings: moving average filters. A moving average filter of window length l of time series data creates a new time series. This time series has the same length as the corresponding sliding window embedding and each point in the time series is constructed using the same window of the original time series. The moving average filter, however, takes the average of the vectors in its window, instead of concatenating them.

Furthermore, we compared our results to the ones obtained from two simpler techniques. For the first technique we attempted to characterize C. elegans behavior using the speed of the worm. To do this, we computed 2-norms of the differences between two consecutive frames of angle data and then averaged all of those discrete derivatives, which resulted in a single value per sample. For the second technique, which could be described as measuring the variance of the worm’s pose over a video, we computed standard deviations for the angle data coordinate-wize, which resulted in a vector of length 100 per sample. In each case we performed a permutation test and multiclass SVM on the resulting feature vectors. The results of these experiments appear in Section 3.2.2.

The following results were obtained by carefully analyzing a sample of C. elegans behavior data from a video of a worm crawling on agar. The sample consists of 400 frames of a 30 frames per second video, so roughly 13 s of behavior. Having a solid surface to provide friction forces slower but more complicated behavior than we would see in an aqueous environment, and we take advantage of the resulting clarity of the data.

In the data we analyzed the subject exhibits the following behaviors in chronological order:

1. crawl forward,

Below we analyze the time series τ from this data, the corresponding sliding window embedding τ 20 , the corresponding moving average filter, and the sliding window embedding of the corresponding null model. These comparisons illustrate various strengths of the sliding window embedding: it smooths the noise from the original time series; it retains more geometric data than the moving average filter; and it captures temporal data from the original time series that is destroyed in the null model. Then, using representative cycles we construct synthetic C. elegans midline data that produce a forward crawl and explain how this process gives concrete interpretations of persistence features in terms of synthetic behavior data. See Section 3.1.1.

To visualize the four point clouds on which we will compute persistence, we apply PCA and project onto the first few principal components. Some of these projections are shown in the two left-most columns of Figure 4. In contrast to the three other time series, the null model time series in Figure 4D has no discernible geometric structure. It appears that the corruption of temporal information has destroyed the interpretability of the visualizations of sliding window embeddings. Meanwhile, the original time series τ in Figure 4A has a similar shape to its sliding window embedding τ 20 in Figure 4B, with the caveat that the sliding window embedding has the effect of smoothing the data and making it more robust to noise. The moving average filter in Figure 4C also has this smoothed property. Though the three time series in Figures 4A–C have similar shapes, persistence diagrams, and persistence landscapes, they vary in one important feature: the pause.

In Figure 10A the points in the sliding window embedding that correspond to frames where the worm is performing a specific behaviors are highlighted. The points corresponding to the pause behavior deviate from the path of points corresponding to crawling backwards. This deviation is small compared to the noisiness of the original time series, so the pause deviation does not create a large topological feature in the graph of the original time series. This is reflected in the persistence diagrams and landscapes of Figure 4 as well; the original time series diagram and landscape Figure 4A show only three significant topological features, while the diagram and landscape of the sliding window embedding Figure 4B and moving average filter Figure 4C show four.

The sliding window embedding and moving average filter both smooth the input data and detect the pause behavior, but they are qualitatively different. One piece of geometric information that the sliding window embedding retains that the moving average filter does not is the direction of the time series. Consider plotting a time series and the corresponding reverse-chronological-order time series in the same ambient space. The points of the original time series would align exactly with its reverse, so the two corresponding point clouds are the same. This is also true of a moving average filter on that time series. The sliding window embedding, however, can have distinct point clouds.

Retaining the direction of time in a time series is particularly important for data that has certain types of symmetry. A natural occurrence of such data is the sine wave data in Example 2.7. Because the data in this time series follows a path and then backtracks along that path, it never produces a loop with any significant persistence. There is no loop in the moving average filter of the data, either.

We present our analysis of an experiment where C. elegans are submerged in solutions with varying viscosities. The viscosity of the solution is correlated with how much methylcellulose is added and experimental conditions are labeled with their methylcellulose content, usually in order from low to high methylcellulose and viscosity.

Average persistence landscapes for each class are shown in Figure 8A. The lower viscosity conditions allow for larger depth one landscapes ( λ 1 ) but have relatively few non-zero higher-depth landscapes, which means there are a smaller number of larger loops detected in the sliding window embeddings. This indicates that in lower-viscosity environments, C. elegans exhibit behaviors of higher amplitude but either demonstrate fewer distinct behaviors or have much less variation between repetitions of behaviors. Conversely, the high-viscosity classes show many more cycles in the sliding window embeddings, with each cycle being small compared to the cycles found in the low-viscosity environments. These observations suggest that at high-viscosity, behaviors do not involve large changes in posture and are more varied. From observing the raw video data, it is apparent that in higher-viscosity environments C. elegans can make smaller, tighter body bends, which is consistent with these results.

We also observed that as viscosity increases, the support of the persistence landscapes stretches further to the right and cycles are born at higher radius values. The worms seemed to exhibit less varied behaviors in lower-viscosity environments, so perhaps in such environments they continued “retracing their steps” through the sliding window embedding space which resulted in more densely sampled curves and thus homology classes formed at lower radii.

The pairwise distances between landscapes for each sample are visualized in Figure 11. The normalized pairwise distances between the average persistence landscapes for each class are shown in Table 1. We include the origin—the zero persistence landscape, i.e. the 0 vector—in these distance computations to complete the normalization. Normalization is such that the average distances between each class and the origin is 1. Multidimensional scaling on these distances visualizes the similarities between samples and classes, respectively, and are shown in Figure 6. From the raw distances and the multidimensional scaling of the distances, we can see that the high-viscosity classes ( 2 % and 3 % methylcellulose) are closest together and that this pair, the 0.5 % class, and the 1 % class are roughly equidistant from one another.

Principal component analysis on the average persistence landscapes for each sample gives the graphs in Figure 7. In Figure 7A, the projections of the average persistence landscapes of the samples onto the first two principal components are given by hollow symbols and projections of the average persistence landscapes of the classes are given by solid symbols. Here we see results similar to those from the multidimensional scaling in Figure 6: the low-viscosity class landscapes are far from each other and the high-viscosity classes, while the high-viscosity class landscapes are quite close. We can also see some of the variance within classes. The low-viscosity classes have much more variability than the high-viscosity classes, with the highest-viscosity class, 3 % methylcellulose, having very little variation in these first two principal components.

These conclusions about variation in each of the classes are supported by the standard deviations of each coordinate in the average (discrete) persistence landscapes of each class. In Figure 8B the standard deviations of each coordinate are graphed as sequences of functions so that the standard deviations can be easily matched up with their corresponding locations on the average persistence landscapes. We conclude that there is little variation in the 3 % class, slightly more in the 2 % class, and much more in the 0.5 % and 1 % classes. The 1 % class showed more variation in higher-depth landscapes than the 0.5 % class, suggesting that C. elegans can produce slightly more complex behaviors in a slightly higher viscosity environments. The variances of the 0.5 % and 1 % classes also exhibit a distinct pattern; the 0.5 % samples varied more toward the lower radius parameters (the left side of the graph), whereas the 1 % samples varied more toward higher (more to the right) radius parameters.

In the following analysis we study the complexity of behavior expressed in each class. Figure 9 shows results from using PCA on the videos in each class. The viscosity of the environment is negatively correlated with the percent of variance explained by the first principal component, which suggests that behaviors in low-viscosity environments are simpler than those in high-viscosity environments. Viscosity appears to be correlated with the number of nonzero landscapes, which also suggests that high-viscosity environments allow for more varied behaviors.

We conducted permutation tests between pairs of classes to determine how well average persistence landscapes can distinguish between samples from different classes. The p-values for these computations are shown in Table 2A. The permutation test gives strong evidence of statistical significance, i.e. that the topological summaries of samples from each class are significantly different.

We then used multiclass support vector machines (SVM) to build a classifier for the samples. The estimated accuracy of the classifier, computed by averaging accuracies across 20 instantiations of the multiclass SVM classifier, was 95.125%. This indicates that persistent homology is able to produce meaningful, distinguishing features from the C. elegans videos. A sample confusion matrix for one instance of SVM is shown in Table 2B. Finally, we used support vector regression to estimate the methylcellulose content in the environment for each sample. The results are plotted in Figure 12. There are two outliers on this graph which are estimated as having negative methylcellulose content. The two animals in these samples moved much more quickly than their peers, so we believe that the SVR is picking up on the strong negative correlation between viscosity of the environment and speed, and based on these animals’ fast speed, assigning a methylcellulose content that is so low that it is negative.