Tracking changes in behavioural dynamics using prediction error

Automated analysis of video can now generate extensive time series of pose and motion in freely-moving organisms. This requires new quantitative tools to characterise behavioural dynamics. For the model roundworm Caenorhabditis elegans, body pose can be accurately quantified from video as coordinates in a single low-dimensional space. We focus on this well-established case as an illustrative example and propose a method to reveal subtle variations in behaviour at high time resolution. Our data-driven method, based on empirical dynamic modeling, quantifies behavioural change as prediction error with respect to a time-delay-embedded ‘attractor’ of behavioural dynamics. Because this attractor is constructed from a user-specified reference data set, the approach can be tailored to specific behaviours of interest at the individual or group level. We validate the approach by detecting small changes in the movement dynamics of C. elegans at the initiation and completion of delta turns. We then examine an escape response initiated by an aversive stimulus and find that the method can track return to baseline behaviour in individual worms and reveal variations in the escape response between worms. We suggest that this general approach—defining dynamic behaviours using reference attractors and quantifying dynamic changes using prediction error—may be of broad interest and relevance to behavioural researchers working with video-derived time series.


Introduction
Behaviour mediates individual interaction with the outside world. A systematic description of behaviour is therefore key to linking dynamic internal (e.g. neural) states, with external biotic and abiotic conditions in natural situations. However, even in cases where behaviour can be easily observed, finding a simple quantitative dynamic description is challenging. Increasingly, motion and pose of freely-moving organisms can be automatically tracked and quantified using computer vision [1]. The model roundworm, C. elegans, presents a particularly wellestablished, promising and tractable case: its free motion on a surface can be automatically quantified from video in exceptional detail using only a few coordinates [2][3][4], and its neuronal network is topologically stereotyped and contains only 302 neurons [5]. Here we address the challenge of quantitatively defining and comparing behavioural dynamics in time series data, focusing on C. elegans locomotion as a case study.
On agar, C. elegans locomotion is characterized by forward motion with intermittent brief reversals. 35 years [6] of studies involving genetic perturbations, genetic sensors, electrophysiology, laser ablation, electron microscopy, computational modeling and other tools have mapped the connectome underlying this behaviour, clarified the contributions of allelic differences to behavioural variation, and have quantitatively dissected how locomotory decisions arise from signaling between network components (reviewed in [7]). Here we focus at the organism level, on the locomotion itself, following a series of papers over the last 12 years that have developed a quantitative framework for describing and analyzing worm movement.
This quantitative framework for C. elegans movement on a surface relies on the observation that freely moving worm body poses can be represented sufficiently as coordinates in a 4 to 6 dimensional space of eigenworms (eigenvectors of worm body segment angles) [2][3][4]. This low dimensional representation can be linked back to traditional categorical descriptions such as direction of motion and turn type or can provide a finer-grained dictionary of pose motifs [3]. From such approaches, previously unclassified distinctions in behaviour and novel characteristics of mutants [3,4] demonstrate the power of this near-complete quantitative description of movement.
Categorizing worm behaviour by pose or motion [3,[8][9][10][11][12] can provide essential information regarding worm strain or internal state. This information can be exploited through the use of parametric or non-parametric models. For example, timing of reversals in direction [13] and escape speed [14] have been successfully predicted using parametric models and behavioural data alone (i.e. without using neural activity data). Recently, nonparametric models, which may make fewer assumptions, have shown great success whether including neural input [15] or using body motion alone [16,17]. [16] used a hybrid approach in which eigenworm time series were dynamically partitioned and approximated locally with a sequence of linear models, revealing larger variation within behavioural types than previously thought. Recent papers have exploited delay-coordinate embeddings to reconstruct an approximate state space (in the dynamical systems sense) in which present states maximally determine future states. [17] used this to obtain a comprehensive characterisation of the nonlinear dynamics of C. elegans locomotion. [15], by incorporating neural dynamics, were able to predict behavioural transitions up to 30 seconds in advance. A defining feature underlying the success of these nonparametric, embedding-based approaches is their ability to exploit the nonlinear dynamics evident in natural data of this kind, without requiring that a simple functional form can be found to accurately describe those dynamics [18].
Because even genetically identical individual worms can show distinct yet self-consistent patterns of motion [15,19], it can be difficult to measure subtle changes in behaviour that may reflect important changes in internal state. Categorical definitions of movement or pose behaviour [8], however fine-grained [3], necessarily collapse subtle variation onto a dictionary of stereotypes. Non-categorical distinctions of behaviour must necessarily confront the problem of measuring distances between time series-a highly non-trivial problem. The recent success of dynamic embedding representations of worm motion [15,17], suggests that this problem could be addressed within the framework of empirical dynamic modeling (EDM) [18].
EDM relies on reconstructing attractors that encode system dynamics directly from time series data (for a simple introduction, see Introduction to Empirical Dynamic Modeling). In the context of worm motion, attractors can be reconstructed from data agnostically, i.e. describing the dynamics in a given set of data without requiring prior assumptions about worm, strain, or motion type. These attractors thereby provide a definition of worm dynamics for a given set of observational data to which other worm dynamics can be easily compared.
In this manuscript, we assess the feasibility and robustness of extending eigenworm representations of motion to dynamic attractors based on time-delay-coordinate embeddings. We then exploit this observed robustness to propose a prediction-based scheme for detecting dynamic change. We validate our scheme using the fluctuations of motion dynamics inherent to delta turns [4] and show that our scheme provides a novel, informative and nuanced view of the return to baseline worm motion following an aversive stimulus.

Method overview
To quantitatively assess changes in worm motion dynamics, we require a quantitative representation of those dynamics. We construct this representation based on a time-delay embedding of a reference time series, or library, of observed worm poses (Fig 1). Each point in this embedding corresponds to a short sequence, of length E, of successive worm body poses in the library, at regular time intervals τ ( Fig 1A). Note that because we encode each worm pose here using five eigenworm coefficients, the embedding space has 5E dimensions in total. The dynamic information about the order in which these pose sequences occur within the library defines a trajectory through the embedding space ( Fig 1B). With an appropriate choice of E and τ, the trajectory unfolds into a nearly-repeating but non-intersecting structure that we here call an attractor.
The dynamic information contained in such an attractor can be exploited to make predictions. Nearby points (where each point is a sequence of E worm poses) on the attractor belong to locally similar trajectories ( Fig 1B). This means the trajectory of an out-of-sample point can be estimated by following nearby points on the library attractor as they move forward in time [18,20]. An out-of-sample point originating from dynamics that are similar to those encoded in the library attractor should permit a better prediction than an out-of-sample point originating from different dynamics ( Fig 1C). Thus, out-of-sample prediction accuracy may function as a useful measure of dynamic similarity (c.f. [21]).
We wish to maximise the temporal resolution of our method, in order to precisely localise behavioural changes in time. Therefore we set the embedding-time-delay interval τ and the from a reference library of worm behaviour are embedded in their dynamic context: short sub-sequences of E worm poses. B) These embedded points, together with their temporal sequence, define a trajectory. An appropriately-embedded time series resolves into an attractor. By following nearby points on the attractor, predictions can be made for out-of-sample data. C) Out-of-sample data from dynamics similar to the reference library result in good predictions (red), whereas different dynamic regimes will be predicted poorly (orange).
https://doi.org/10.1371/journal.pone.0251053.g001 prediction horizon T p equal to the sampling interval of the data. For the forward prediction, we use the S-Map method [20] as it allows the sensitivity to local attractor structure to be tuned, thereby optimally exploiting nonlinearity in the data. This leaves us with two free parameters that define our measure: the number of successive worm poses E, and the S-Map nonlinearity parameter θ (see also Methods and materials).

Parameter robustness
A natural concern is whether such a measure will be robust enough to parameter variation to actually be useful. We assess this on a suite of publicly available worm motion data from four strains: wild type (N2) and three diverse mutants (octr-1(ok317), unc-80(e1069) and dpy-20 (e1282)IV), using library and prediction intervals of 1000 samples each (Fig 2, see Methods and materials).
Average prediction error remains stable as a function of both E and θ, over a wide range of parameter settings, for almost all worms of each type (Fig 2). Within each type of worm, the shape of the error as a function of the parameters shows strong consistency between individuals. Generally, the optimal parameter settings occur at E > 1, and the wildtype shows evidence of nonlinearity, i.e., the optimal value of θ > 0, which is consistent with previous observations [16,17]. The mutant worms on the other hand generally do not show clear evidence of nonlinearity on average, though the sensitivity to θ is very weak over the parameter range explored. This more linear characteristic of the mutant data may be attributed to at least two causes: 1)

Fig 2. Prediction error is robust to parameter variation.
Top row: robustness to E (at θ = 2) for data from a wild type and three mutants (headings; see text). Bottom row: robustness to θ (at E = 5) for the same four worm types as top row. Thin lines are individual worms, thick lines are means. Note that there is a coordinate system variation between the wildtype and mutant worm data, which hinders direct comparison of error magnitudes.
https://doi.org/10.1371/journal.pone.0251053.g002 the mutant data do not contain self-intersecting worm poses, whereas these were reconstructed in a second processing step for the N2 wildtype data [4]; 2) the sampling rate of the mutant worm data is almost twice that of the N2 wildtype data (30 Hz vs. 16 Hz), which will make the data appear locally far more linear at the single timestep scale. Indeed, looking at an additional wildtype from the same data source as the mutant time series (i.e., sampled at 30 Hz without self-intersecting turns) reveals a very similar E, θ robustness profile as the other mutants (S1 Fig). Despite all of these differences between the acquisition and pre-processing of the time series however, the prediction error remains stable over a single wide parameter range across all of these data. Furthermore, the prediction error is robust in detail through time, with the error time series showing strong correlation across parameter settings (S2 Fig). For all subsequent analyses, we therefore choose mid-range parameters E = 5 and θ = 2, to allow for any transient highly nonlinear behaviour that may occur, without substantially sacrificing average prediction error. In all subsequent analyses also, we work only with the N2 wildtype time series from [4]. We therefore convert the observed and predicted eigenworm coefficients to worm body section angles [2] prior to calculating error, to allow the use of intuitive error units: RMS error of body section angles in radians.

Validation and dynamic sensitivity
In addition to being stable, we also require our measure to be meaningful: to detect subtle variations in worm motion dynamics that may potentially have an exogenous origin, such as a change in (latent) internal state or external conditions. We propose here that delta turns [4] may provide an ideally constrained example of such exogenous interference. During delta turns, where the worm must intersect and then cross over (or under) its own body, the initiation and completion of the turn may result in both worm body shape deformation (thereby temporarily modifying the coordinate system of the dynamics) and stick-slip-type dynamics that are likely to be highly variable. Furthermore, because the traditional algorithms used to extract worm pose from image data fail when the worm body intersects itself, a supplementary inference method is needed [4]. The uncertainty inherent in making such an inference results in a small increase in worm pose estimation error [4]. This error is exogenous to the natural dynamics of worm motion. Taken together, we thus we expect that a method for detecting variation in worm dynamics should show a distinct spike during the transition from non-intersecting pose to crossing pose (at delta turn initiation), and from crossing to non-intersecting (at delta turn completion).
In Fig 3, the time series of a foraging worm from the [4] data (see Methods and materials) do not show obvious differences between the reference library and prediction intervals ( Fig  3A). Accordingly the prediction based on our attractor performs well at almost all time points, with the exception of four pronounced pairs of sharp peaks and one isolated sharp peak ( Fig  3B). The location of these peaks does not correspond to any clear signal change in the original time series (arrows, Fig 3A). The time points where the magnitude of the third eigenworm coefficient |a 3 |>15 correspond to four tight turns [4] that upon visual inspection are confirmed to be delta turns (red points, Fig 3B). The four delta turns are each book-ended by pronounced spikes in error, that correspond to the transitions during initiation and completion of the delta turn described above (Fig 3C). Closer inspection of the localisation of prediction error along the worm body confirms that the head-end of the worm is responsible for the error peak during delta turn initiation, and the tail-end during delta turn completion in all four cases, exactly as would be expected (Fig 4). The origin of the isolated spike in error around 8 seconds is thought to arise from a dynamic sequence during a 'w'-shaped reverse motion that is sufficiently different from motions sampled in the library (see also Discussion). Note that values in the prediction set are only predicted one time step ahead based on the previous embedded value from the prediction set; i.e., we are not "forecasting" 60 seconds into the future (see Method Overview and Methods and materials). Error peaks in (B) indicated by arrows do not correspond to obvious features in the time series. B) Prediction error corresponding to prediction set in (A). Turns (red dots), characterised here by third eigenworm coefficient |a 3 |>15, are book-ended by periods of anomalous dynamics. The mean error for a constant predictor (i.e. the prediction that the worm pose is the same as the previous time step) is higher than the vertical axis range (0.126). C) A closer look at the worm pose time series (black) and predictions (orange) reveals that the anomalous dynamics correspond to delta turn initiation and completion, when the worm head and tail respectively transition between self-intersection and no self-intersection. https://doi.org/10.1371/journal.pone.0251053.g003

PLOS ONE
Tracking changes in behavioural dynamics using prediction error

Application example
We use the measure of worm dynamics variation to examine the return to baseline behaviour following an aversive stimulus. In this example, described in [2], worms were subjected to a brief infra-red laser pulse to the head after 10 seconds of recording. The worms then exhibited a characteristic escape response: rapid reverse crawling motion, followed by a tight (omega) turn, and then forward motion in the new direction away from the location of the laser pulse [2].
We use the first 80% of time points before the stimulus as a reference library of baseline behaviour, to allow us to compare out-of-sample prediction error before and after the stimulus. We do not perform any additional optimisation of the parameters of our algorithm, taking instead the same parameters used earlier: E = 5 and θ = 2 (see Methods and materials). We do this despite the fact that these escape response datasets have a different sampling rate from the foraging worm datasets (20 Hz vs. 16 Hz) in order to highlight the practical robustness of our method. Because the escape response contains behaviours that are almost surely not contained in the pre-shock reference library of worm dynamics, we perform an approximate classification of worm behaviour across the entire time series to aid interpretation. We identify forwards and backward crawling motion based on the sign of the apparent phase velocity in the (a 1 , a 2 ) plane [2] (provided the amplitude in this plane is sufficiently large, see Methods and materials), and we conservatively identify as putative tight turns time points where |a 3 |>10. Fig 5 shows the typical escape-response sequence of reverse-turn-forward motion following the stimulus. This sequence is also reflected in the prediction error time series: at stimulus onset, there is a sharp increase in the prediction error that persists throughout the rapid reverse crawling and subsequent turn, and then eventually returns close to baseline during

PLOS ONE
Tracking changes in behavioural dynamics using prediction error  stimulus at 10s, despite a very small library size (first 8 seconds, grey). Different rows correspond to different worms (see Materials and methods). Background colours indicate approximate behavioural classification into tight turns (red, |a 3 |>10), forward and backward motion (blue and yellow, positive and negative phase velocity in the a 1 , a 2 plane, respectively), and unclassified (white, when amplitude in the a 1 , a 2 plane or the estimated phase velocity are too small), intended as an approximate guide only (see Methods and materials). Grey lines subsequent forward motion (e.g. Fig 5A and 5B). In some cases (e.g. Fig 5D and 5E), the prediction error remains high, suggesting that the movement dynamics of the worm towards the end of the time series remains substantially different from the movement dynamics prior to the shock. However, care is required when interpreting this result. By definition, baseline is whatever dynamics occupy the reference library. If the reference library does not adequately sample the behavioural regime of interest, then it can not adequately represent it. For these escape response data, there are very few time points available before the shock from which to construct a reference library, and so the reference library is short: 160 time points, or 8 seconds (by contrast, for the foraging worms we used 1000 time points, or 62.5 seconds). Nonetheless, even using a limited library, in Fig 5A, 5B, 5F, 5G and 5H, the prediction error clearly indicates that the forward crawling motions at the beginning and the end of the sequence are similar (despite differences in phase velocity) and that they are substantially different from the bulk of the escape response.
Comparing between escaping worms reveals fine-grained variations in escape response. If we use, for example, the worm in Fig 5A as the reference library for the worm in Fig 5B, and vice versa (Fig 6), we see a dramatic difference in the escape response of the worm in Fig 5A  around 15 seconds (Fig 6B). This behavioural sequence, that is different from any sequence exhibited by the reference worm, results in a peak error several times higher than other major peaks in either error time series. This anomalous behaviour is also visible in the approximate classification scheme: the apparent (a 1 , a 2 ) plane phase velocity changes sign between the escape reversal and the turn, resulting in a transient 'forward motion' classification. Importantly, however, detecting this difference using our method does not require any special foresight about which data features (e.g. phase velocity) should be used for detection. In addition to this most obvious feature, we see clearly that the forward motions of each worm before and after the escape response sequence are broadly similar (low error, Fig 6A and 6B). The error during the initial reversal shows similar behaviour in each case, with an initial peak of similar magnitude, that then declines throughout the reversal sequence. This highlights a challenge in the measurement of worm motion behaviour: in most dynamic representations [3,17] rate and shape are coupled. That is, if a worm were to run through exactly the same sequence of poses at a different rate, this would map into a different position in the dynamic representation (embedding) space. This is also true of our method. Hence, during the rapid reversal, which is initially fast and gradually decays in speed and phase velocity (Fig 5 and [14]), we see a concurrent decrease in prediction error as motion slows (except at the end of the reversal in Fig 6A, where the reference worm's behaviour differed). Note however, that the prediction error during reversal when predicting from one worm to the other is still substantially lower than the prediction error during reversal in Fig 5, where no reversal time sequence was included in the library.

Discussion
The approach presented here exhibits several advantages in comparison with other methods that can be used to detect change in worm motion dynamics. Figs 3-5 show that the output signal of the algorithm lets us localize changes in worm behaviour to within a few time samples. This property is a direct consequence of our use of prediction. Template-matching approaches to measure motion change [3] are similar on a local scale to time-delay indicate the apparent relative phase velocity in the a 1 , a 2 plane used in the classification scheme (see Methods and materials). The dashed horizontal lines indicate the mean error for the constant predictor (i.e. predict that the worm pose is the same as at the previous time step) provided that error lies within the vertical axis range.
https://doi.org/10.1371/journal.pone.0251053.g005 embeddings. However, because these approaches do not use the dynamic trajectory information in an embedded space, they do not fully exploit the time resolution of the data. Comparing predictive output at a single time point focuses this information to single-sample resolution. In order for such fine scale prediction to be meaningful, however, it has to be able to exploit relevant information globally in time-achieved here using the attractor. This data-driven focusing of relevant information across time scales is a natural advantage of the EDM approach.
A second advantage of our approach is its flexibility and sensitivity. These arise from the use of a user-specified library for the construction of the attractor, allowing the relevant dynamic reference information to be tailored with great precision. Although here we have focused on dynamic change within and between individual worms, the approach can be easily and naturally extended to composite libraries encompassing broader classes of behaviour.
This flexibility and sensitivity comes with a cost. An inappropriate choice of reference library can lead to output that carries little useful information (see Fig 5D and 5E). We suggest, therefore, the use of an independent method of behavioural categorisation (such as the one used here) in conjunction with our method, to aid the user in distinguishing dynamic differences due to behaviours that were absent or under-sampled in the reference library from dynamic differences due to more subtle changes in state (see S3-S14 Figs). In some applications, however, it may not be necessary to make this distinction. For example, methods that look at average behavioural differences over an extended time window to give a global measure of behavioural difference will usually conflate these two types of difference but may still yield useful information [3].
Because our method uses a fixed sampling interval to define the embedding, pose sequences performed at different rates will be represented at different locations in the embedding space.
Although the method appears robust (see Fig 2 and the adoption of parameters from Fig 3 in Figs 5 and 6 despite a change in sampling rate) it may be desirable in some instances to disentangle the rate of progression through a pose sequence from the pose sequence itself. This might allow, for example, the differences in worm escape response reversal to be determined more precisely (see Fig 6). To achieve such a separation would require the time delay embedding construction to be adjusted according to, e.g. phase velocity or worm speed: a process that would likely benefit from an additional dimensionality reduction step in time series embedding such as that used by [17]. Here, we have focused on the simplest embedding scheme to aid interpretability.
In cases where the sampling rate of the data is very high (e.g. 30 Hz and above), it can be expected that the constant predictor will become quite accurate, as the pose change between successive frames becomes very small. (Note the substantially lower constant predictor error in Figs 5 and 6 where the sampling rate is 20 Hz, than in Figs 3 and 4 where the sampling rate is 16 Hz.) In some such cases, it may be desirable to sacrifice temporal resolution for increased sensitivity, by setting τ and T p to twice the sampling interval of the data. In this initial study we have focused on maximum temporal resolution as a base case.
More generally, the simple, flexible, sensitive and robust characteristics of our behavioural change measure when analyzing C. elegans pose dynamics result from the ability of EDM tools to exploit the exceptional richness and resolution of time series data, without imposing strong model assumptions. The recent success of other empirical nonlinear dynamics approaches [15,17] suggests that such techniques may prove increasingly useful in future C. elegans behavioural analyses.
Although we have focused on C. elegans here, the embedding and prediction tools that underlie our approach are applied widely to time series data in general. The key observation that prediction error can be used to measure meaningful behavioural change can therefore also be explored in the increasing number of other high-resolution behavioural time series extracted from video (e.g. [1]).
In our robustness analysis we include time series of the first five eigenworms from three additional mutant strain datasets from OpenWorm [22], collected as per [19]. The mutations selected are octr-1(ok317) (strain VC224), unc-80(e1069)V (strain CB1069), and dpy-20(e1282) IV (strain CB1282). Twelve individual worms were analyzed for octr-1(ok317) and dpy-20 (e1282)IV, while nine individuals were analyzed for unc-80(e1069)V. Embeddings were constructed using the discontinuous method (see below), starting from sample 1001, and generating embeddings of 1000 points each for the library and prediction. See S1 File for lists of the specific files included.
Embedding construction. We use two slightly different methods to construct embeddings: a continuous method and a discontinuous method.
The continuous method is suitable for time series data that do not contain substantial breaks (i.e. NaN values). It generates embeddings from fixed data intervals (e.g. from sample 12001 to 13000 inclusive). This allows precise delineation of behaviours of interest, and simple identification of prediction error with time. Note that because this method only has access to exactly the data interval specified, the embedding coordinates for the first few time points will be undefined, as time lags from before the interval are not available. The continuous method is used for all figures where error time-series are plotted for the N2 wildtype worms: S2-S14 Figs; Figs 3-6.
The discontinuous method can account for breaks (NaN values) in the time series. It is used here for all figures that include the mutant worm data (which contain NaN values during self-intersecting worm poses): Fig 2 and S1 Fig The discontinuous method takes as input only a starting sample (e.g. sample 1001), and proceeds to construct an embedding from there, accepting only embedding points for which all lag coordinates and target coordinates are available. To avoid embedding mismatches due to changing E during robustness testing, we construct a master embedding for the highest value of E, and use subsets of this embedding for lower values of E.
Prediction and error metric. We use the rEDM package in R. Using the S-map [20] feature in rEDM for each eigenworm coefficient and the E � 5 dimensional embedding, we make out-of-sample predictions of the first 5 eigenworm coefficients. In all cases we used τ = T p = Δt, where Δt is the sampling interval of the data. For quantification of prediction error for any given worm at a single time point, we convert predicted and observed eigenworm coefficients to 100 worm body section angles [2], then calculate RMS error in radians, unless otherwise specified (cf. Robustness calculations, below).
Robustness calculations (Fig 2). We use the first 1000 embedded points (see Data and Embedding construction, above) to as the reference library attractor and predict from the remaining 1000 embedded points. The vertical axis of Fig 2 shows the mean over time of the single time point RMS error between the observed and predicted eigenworm coefficients, for each worm. The thick lines described as "means" are simply constructed by taking the arithmetic mean of the values in the thin (individual worm) lines.
Time dependence of error (Figs 3 to 6). Fig 3 uses S-map analysis as above to make predictions on the eigenworm coefficients of foraging worm 1 from [4], with parameters E = 5, θ = 2. The reference library is constructed from frames 10001-11000, while the eigenworm coefficient predictions are made using frames 11001-12000. Predicted and observed worm poses are displayed at points of interest by converting the body angles calculated using the eigenworms to coordinates in the plane assuming fixed length of each worm segment, and aligning predicted and observed worms at the head. Fig 4 incorporates the same data and method as Fig 3, but with RMS error in predictions calculated over each of five equally-sized non-overlapping regions of the worm as indicated in the figure. These equally-sized regions are each composed of twenty body angles extracted from the eigenworm coefficients.
In Fig 5, we use the first eight seconds of pre-stimulus behaviour to predict post-stimulus escape response while tracking prediction error. For approximate behavioural classification, see below.
In Fig 6, we proceed as per Fig 5, but instead of using the first 8 seconds of the predicted worm time series as the library, we use the entire time series of a different escaping worm, as specified in the text.
Approximate classification of worm behaviours. Given time series a 1 (t), a 2 (t), a 3 (t) corresponding to coefficients of the first, second and third eigenworms respectively, sampled at regular time intervals Δt. Let c(t) be the class label assigned to the worm behaviour at time t. We define our approximate classification scheme as follows [2,4]. Note that in the main text we have reversed the sign of the phase velocity relative to what is described here, such that forward motion corresponds to positive phase velocity, to aid intuition. where o a 1 ;a 2 ðtÞ is a local estimate of phase velocity in the (a 1 , a 2 ) plane given by: α(t) ≔ arctan(a 2 (t + Δt)/a 1 (t + Δt)) − arctan(a 2 (t − Δt)/a 1 (t − Δt)) jo a 1 ;a 2 ðtÞj ¼ min½aðtÞ; p À aðtÞ�=2Dt signðo a 1 ;a 2 ðtÞÞ ¼ signðaðtÞÞ if α(t)<π − α(t), and signðo a 1 ;a 2 ðtÞÞ ¼ À signðaðtÞÞ otherwise.
The minimum phase velocity threshold � = 0.1 rad/s, and minimum phase plane amplitude δ = 3 (no units) were chosen to improve stability.