Recurrence eigenvalues of movements from brain signals

The ability to characterize muscle activities or skilled movements controlled by signals from neurons in the motor cortex of the brain has many useful implications, ranging from biomedical perspectives to brain–computer interfaces. This paper presents the method of recurrence eigenvalues for differentiating moving patterns in non-mammalian and human models. The non-mammalian models of Caenorhabditis elegans have been studied for gaining insights into behavioral genetics and discovery of human disease genes. Systematic probing of the movement of these worms is known to be useful for these purposes. Study of dynamics of normal and mutant worms is important in behavioral genetic and neuroscience. However, methods for quantifying complexity of worm movement using time series are still not well explored. Neurodegenerative diseases adversely affect gait and mobility. There is a need to accurately quantify gait dynamics of these diseases and differentiate them from the healthy control to better understand their pathophysiology that may lead to more effective therapeutic interventions. This paper attempts to explore the potential application of the method for determining the largest eigenvalues of convolutional fuzzy recurrence plots of time series for measuring the complexity of moving patterns of Caenorhabditis elegans and neurodegenerative disease subjects. Results obtained from analyses demonstrate that the largest recurrence eigenvalues can differentiate phenotypes of behavioral dynamics between wild type and mutant strains of Caenorhabditis elegans; and walking patterns among healthy control subjects and patients with Parkinson’s disease, Huntington’s disease, or amyotrophic lateral sclerosis.


Introduction
Caenorhabditis elegans is a nematode or roundworm about 1 mm in length [1] and commonly used as a model organism in the study of genetics because of its powerful genetics and fully characterized simple nervous system [2,3]. The nervous system of the C. elegans hermaphrodite consists of 302 neurons that form 118 morphologically distinct neuron classes [4]. These neurons activate many distinct stimulus modalities and then combine them to produce distinct patterns of behavior [5,6].
The study of the movement of C. elegans is important because the knowledge gained from understanding the mechanism underlying the moving of these worms can be useful for discovering new characteristics in behavioral genetics. It is known that behavior is a visual display of sensitive and integrative information of nervous system function and plays as an effective measure for evaluating the effects of mutation or efficacy of drug treatment for animals [7]. Ultimately, such knowledge is expected to provide potential alternatives for better diagnosis and therapeutics of human disease [8].
An interesting study reported in [9] carried out an experiment on the movements of C. elegans. An explanatory summary of the study can be found online (https:// ilove sympo sia. com/ 2008/ 10/ 02/ eigen worms/). In this study of the dimensionality and dynamics of C. elegans, the shape of a worm was modeled as a vector of 100 angles measured between 101 adjacent segments. Consequently, a 100 × 100 correlation matrix was then constructed and its 100 eigenvalues were computed. The study discovered that only 4 of the 100 eigenvalues accounted for more than 90% of the variability of worm shapes. These 4 eigenvectors define 4 fundamental worm shapes, which are called eigenworms. As a finding, the worm motion could be categorized into four characteristic shapes whose combination can describe any shape that the worm can form. The first two eigenworms, which correlate with a wave propagating along the body of a worm, contribute to the forward motion. The third eigenworm correlates with the curvature whose variations navigate the movement of a worm along its trajectory. Finally, the fourth eigenworm correlates with the movement of the worm's head as it searches for food and traverses.
As an important step toward the elucidation of behavior in C. elegans, an ability to measure the difference between patterns of motion of these worm types is necessary to allow for pattern prediction that can provide insights into the control of the dynamics of movement [10]. It is therefore of interest to explore methods for quantifying and differentiating complexity in time and vision series of eigenworm motions between wild type and mutant strains of C. elegans [10,11].
Neurodegenerative diseases affect many physical activities of the patients, particularly balance and movement. Many of these diseases are thought to be gene inheritance disorders, but the cause is largely unknown [12]. Diseases of degenerative nerves can be life-threatening and currently have no cure. Therefore, effective treatments, including medications and surgery, may help the patients improve mobility and relieve pain. Gait impairment is a common symptom in neurodegenerative disorders. Specifically, gait variability, which is the stride-to-stride fluctuations measured with time, has been known to be associated with neurodegeneration [13].
Several studies on gait dynamics in human neurodegeneration have been reported in literature [14][15][16][17][18][19]. These studies attempted to discover new features of gait time series mainly used for pattern classification of healthy control (HC), Parkinson's disease (PD), Huntington's disease (HD), and amyotrophic lateral sclerosis (ALS) subjects. However, quantitative characterization of gait and its impairments has not been well investigated. New computational methods that can provide scalar values to represent some attributes of the gait dynamics of the HC and diseases can be useful because these quantitative descriptors can used for early disease detection or physiological makers.
In this paper, the method for computing the largest eigenvalue of a convolutional fuzzy recurrence plot of time series [20] is investigated for measuring the quantity of the complexity of behavioral and moving patterns of eigenworm behavioral phenotypes of different C. elegans strains, and patients with neurodegenerative disorders (PD, HD, and ALS) and HC subjects. Here, based on the concepts of nonlinear dynamics and chaos theory, time series of the time series are reconstructed into phase space sets and their matrices of spatial recurrence features extracted. A convolutional kernel is then iteratively applied on these recurrence features to encode feature invariance. Finally, the largest eigenvalues of the deepest convolutional recurrence matrices are determined and used as indicators of complexity of motion or walking patterns in different worm types or human subjects, respectively.

Eigenworm data
The time-series dataset used in this study was described in [21]. The data relate to 258 traces of worms converted into eigenworm time series. Worm motions on an agar plate were recorded and a range of human-defined features measured [22]. The worm outline was extracted, where each frame of the worm motion was captured and represented with amplitudes when the shape was projected onto the eigenworms. There are five classes of C. elegans: N2, goa-1, unc-1, unc-38, and un63. The N2 is wild type (normal), and the other four are mutant strains. This dataset, which is publicly available online [23], consists of the time series of the first dimension or first eigenworm. Table 1 shows the number and length of the time series of each of the five worm types.

Gait in neurodegenerative disease data
This third-party public database [24] includes gait signals recorded from patients with PD, HD, ALS, and HC subjects. The database also provides clinical information for each subject, including age, gender, height, weight, walking speed, and a measure of disease severity or duration. The signals were obtained using force-sensitive resistors, with output proportional to the force under the foot. Stride-to-stride measures of footfall contact times were derived from these signals. Detailed information about these data was provided in [25,26]. Both time series of the left swing interval (LSI) and right swing interval (RSI) are used in this study. Both LSI and RSI were measured in seconds. Because these original time series have not been filtered, they are filtered using 1-D median filtering that applies a third-order one-dimensional median filter to the input time series and considers the signal to be zero beyond the endpoints. To ensure the time series of all subjects are of the same length, all the time series have a truncated length of the first 120 time points, which is equivalent to the original shortest timeseries length. Table 2 shows the numbers of HC, PD, HD, and ALS subjects and truncated length of the LSI and RSI time series.

Eigenvalues of recurrence dynamics
A sequence of values in time can be transformed into an object in space, which is called the phase space of a dynamical system, and the object in the phase space is referred to as the phase space set. Such a transformation is particularly useful because some properties underlying complex data can be more easily discovered from the phase space set than from the original time series [27]. According to the embedding theorem [28,29], it is always possible to reconstruct the phase space from a time series. A popular method for the phase-space reconstruction is the technique of time delay [28,29], which is applied for reconstructing the phase space of the eigenworm time series as follows. Let t = (t 1 , t 2 , . . . , t N ) be a time series of scalar values generated by the eigenworm motion. Given m as an embedding dimension and τ as a time delay, the time series t can be reconstructed as . Now it can be seen that the time series t can be transformed into X that is expressed as a matrix, where each row is a phase-space vector Based on the phase space set reconstructed from a time series, a method for studying behaviors of systems of nonlinear dynamics is the recurrence plotting that examines the revisit or recurrence of the trajectory in the phase space [30]. A recurrence plot (RP) is a binary symmetrical matrix showing either pairs of time points at which the trajectory is at the same place (representing with a black dot) or not (representing with a white dot). As an extended algorithm of an RP, the construction of a fuzzy recurrence plot (FRP) [31] was introduced to overcome the difficulty for determining the threshold for similarity and limited binary expression of recurrence imposed by the RP method. An FRP represents the recurrence of the trajectory in the phase space as a grayscale image taking values in [0, 1], where 0 is a black pixel, 1 a white pixel, and other values gray pixels. The construction of an FRP works by partitioning X into a number of clusters, denoted as c, using the fuzzy c-means (FCM) algorithm [32]. The FCM assigns a fuzzy membership grade, denoted as µ ij taking values in [0, 1], to each x i , i = 1, 2, . . . , M , with respect to each cluster center v j , j = 1, 2, . . . , c . The fuzzy membership variable µ ij expresses the degree that x i possibly belongs to v j .
The FRP method applies the two following properties: • Reflexivity: • Symmetry: The fuzzy relation between the pairs of the phase-space vectors x i and x k ; i, k = 1, 2, . . . , M can be inferred by transitivity using the max-min composition operator [33] as As a result, an FRP is a square grayscale image or matrix defined as By using on the concept of convolutional filtering in deep learning for invariant feature extraction of image objects, where convolutional layers are building blocks for the design of convolutional neural networks [34], a deep convolutional FRP, denoted as cFRP, can be performed by a series of convolutions of the FRP and a convolution kernel as [20] where w is a convolution kernel and ⊛ denotes the convolutional operator.
In this study, the kernel is used as a sharpening-effect 3 × 3 matrix, whose elements are given as After each step of the convolution, the rectified linear unit (ReLU) is then applied to eliminate negative values of the cFRP. The ReLU for a cFRP returns a zero or positive value to each element cFRP(i, k) as The next step is to reduce the size of the convolved FRP by using a pooling operator, which aims to reduce the dimensionality of a convolutional feature map but still keeps the useful information. Different types of pooling include the max, average, or sum operator. A popular choice for the pooling of a convolutional matrix in deep learning is the max pooling operator. Let a set of pooling regions be � = (� 1 , � 2 , . . . , � Q ) , where � q = (ω q,1 , ω q,2 , . . . , ω q,m×m ) , q = 1, 2, . . . , Q . The number of pooling regions Q within a convolutional matrix is determined by the pool size m × m and stride that is the step size for traversing the convolutional FRP. The max pooling that operates on a pooling region of size m × m , denoted as P max , is defined as Finally, the eigenvalues of a cFRP having a certain small square matrix size can be determined, where the largest eigenvalue, denoted as max , is used as the characteristic value of recurrence dynamics.
Procedure for computing the largest eigenvalue of a cFRP: 1 Input: an N × N FRP 2 Given w, ReLU, pool size, and stride 3 Given n as the desired final n × n convolved FRP 4 While N > n 5 Perform convolution of the FRP and w. 6 Apply ReLU on the convolved FRP. 7 Perform max pooling on the convolved FRP. 8 End While loop. 9 If N = n , compute the largest eigenvalue of the n × n cFRP. Figure 1 shows the iterative process for computing the final convolved FRP whose largest eigenvalue is determined.

Eigenworm behavioral phenotypes
To determine the largest eigenvalues from the time series of the C. elegans behavioral phenotypes presented in the foregoing section, the embedding dimension m = 4 was selected based on the identification of the four principal dimensions of the eigenworms, time delay τ = 1, and number of clusters c = 3, 5, and 7. Figure 2 shows some time series and FRPs of the five C. elegans classes, where the FRPs were constructed with m = 4, τ = 1, and c = 3. It can be seen that, the texture of the FRP of N2 (wild type or normal) is distinctive from the four mutant classes (goa-1, unc-1, unc-38, and unc-63), while the texture patterns of (goa-1 and unc-1) and (unc-38 and unc-63) are similar to each other. To compute the largest eigenvalues, the size of the final cFRPs for the five C. elegans classes was set as n = 2, that is a 2 × 2 matrix, to allow the deepest feature extraction based on the deep-learning approach. Table 3 shows the average largest eigenvalues and standard deviations (SDs) of the C. elegans behavioral dynamics with embedding dimension = 4 and different numbers of clusters. Table 3 also shows the p-values, 95% To compare the results obtained from the largest eigenvalues with one of the most popular methods for quantifying the fluctuation or predictability of time series, the method of sample entropy (SampEn) [35] was used to compute the SampEn values of the time series of the five C. elegans types. Input parameters for computing SampEn values of the time series were selected as: embedding dimension m = 4, time delay τ = 1, and threshold δ = 0.1σ , 0.2σ , and 0.3σ , where σ is the standard deviation of the time series (a practical and welladopted selection of δ was suggested to be between 0.1σ and 0.25σ [36,37]). Table 4 shows the average SampEn values and SDs of the C. elegans dynamics with different values for δ as well as the p-values, 95% CIs, and 99% CIs.

Gait in neurodegenerative diseases
To determine the largest eigenvalues from the time series of gait patterns of the HC, PD, HD, and ALS cohorts, the embedding dimension m = 1 was selected based on the one dimension of the signals, time delay τ = 1, and number of clusters c = 3. Figures 3 and 4 show some time series and FRPs of the HC, PD, HD, and ALS subjects. Being similar to the case of the C. elegans, to compute the largest eigenvalues, the size of the final cFRPs for the HC, PD, HD, and ALS time series was set as n = 2, that is a 2 × 2 matrix, to allow the deepest feature extraction based on the deep-learning approach. Table 5 shows the average largest eigenvalues with standard deviations (SDs), p-values, 95% confidence intervals (95% CI), and 99% confidence intervals (99% CI) of the LSI and RSI of the HC, PD, HD, and ALS subjects. Table 6 shows the average SampEn values with standard deviations (SDs), p-values, 95% confidence intervals (95% CI), and 99% confidence intervals (99% CI) of the LSI and RSI of the HC, PD, HD, and ALS subjects. Input parameters for computing SampEn values of the time series were selected as: embedding dimension m = 2 (m = 1 resulting in some SampEn values = ∞ ), time delay τ = 1, and threshold δ = 0.3σ.

Discussion
The average largest eigenvalues of the eigenworm recurrence dynamics that were computed with the embedding dimension of 4 show consistently with evidence of statistical significance that the average largest eigenvalues of N2 (wild type/normal) are smaller than those of the four mutant worms with c = 3 and 5. The mutant type goa-63 has the average largest eigenvalues for c = 3 and 5. The goa-1 has the largest average largest eigenvalue for c = 7, but the average largest eigenvalue of N2 (5.2862) is slightly larger than unc-1 (5.2543). The mutant type unc-38 has the smallest average largest eigenvalues among the other three mutant types for both c = 3 and 5. However, the upper bounds of 95% and 99% CIs of the four mutants are higher than those of the wild type for all three different numbers of clusters. In general, these results suggest the use of c being either 3 or 5 for producing consistent In contrast to the maximum eigenvalues, the average SampEn values of the wild type for all three thresholds are smaller than those of mutant types goa-1 and unc-1 but larger than those of mutant types unc-38 and unc-63.
Such results suggest the SampEn method fails to differentiate between the wild type and mutant strains of C. elegans.
For the analysis of gait dynamics, the results shown in Table 5 consistently indicate that the average eigenvalue of the HC is lower than those of the patients with neurodegenerative disorders using either LSI or RSI data.  While the eigenvalues of PD, HD, and ALS are closer together, the eigenvalues of HD and ALS are closest for the case of LSI; but PD and HD are closest for RSI (see phylogenetic trees in Figure 5a, b, which were constructed using the unweighted pair group method with arithmetic mean (UPGMA) [38]). In both cases of LSI and RSI, the HC is clearly separated from the neurodegenerative−disorder cohort (Figure 5a, b). On the contrary, values of the SampEn computed using either LSI and RSI time series as shown in Table 6 do not differentiate gait patterns between the HC and neurodegenerative-disorder subjects. These results can be visualized with the phylogenetic trees in Fig. 5c, d constructed using the UPGMA, where HC and HD are located in the same group (Fig. 5c) for the case of the LSI, and HC and ALS in the same group (Fig. 5d) for the RSI. Furthermore, for SampEn analysis of gait patterns of HC, PD, HD, and ALS subjects, it was found that for m = 1, SampEn = ∞ for some LSI and RSI time series. This was because no similarity between sub-segments of the time series had been detected, giving a conditional probability of zero and resulted in an infinite value of SampEn. Therefore, m=2 was then chosen.
Findings from the eigenvalues of gait suggest that LSI and RWI are potential markers of neurodegenerative disorders, which can help identify early stage of the diseases for optimal treatment and intervention. In deep learning, which is the state-of-the-art of artificial intelligence, convolutional layers, which perform convolution operations, are the critical basic units implemented in convolutional neural networks. In functional analysis, a convolution is a mathematical operation on two functions or in this case it applies a kernel to an input, resulting in a transformed feature that expresses how the property of the original input is modified by the kernel. The iterative operation of the same kernel to the subsequent inputs results in a powerful feature map that can be used for classifying images of different objects. This application of the same kernel across an input image is a powerful idea [39].
In this study, the filter kernel was designed to sharpen edges or fuzzy information of the input FRP. The kernel, which was repeatedly applied across the FRP after the operations of the ReLU and max pooling, allowed the creation of deep features in the down-sampled FRPs. This ability of the AI method is commonly referred to as translational invariance that has been proven to be robust for addressing prediction problems [39]. This deep-learning mechanism offers an insight into the capability of the cFRP-based eigenvalues for identifying the dynamics of movements controlled by brain signals.

Conclusion
The proposed algorithm for computing the largest eigenvalues of convolutional fuzzy recurrence plots of time series can differentiate normal from mutant types of Caenorhabditis elegans, and gait patterns between healthy control subjects and patients with neurodegenerative disorders. Methods for prediction of worm types by using time series of spatial movements recorded from the eigenworms are needed to assist life-science and neuroscience researchers to gain deeper understanding of the genetic basis of behavior and facilitate studies to probe how molecular, cellular, and systems-level approaches can use sensory inputs to infer neural circuits and behaviors [6] as well as other complex neural actions [40]. The eigenvalues of left and right swing intervals of gait dynamics can be useful for gaining insight into the dynamics of the sub-phases of the stride as well as effect of the diseases on gait asymmetry, and potentially used as markers of disease progresses.
The method for computing eigenvalues of a convolutional FRP appears to be a useful computational tool for nonlinear time-series analysis to quantify other types of behavioral dynamics and movements produced from brain signals. There are several other different filter kernels that can be explored and implemented either separately or parallelly to extract deep feature maps of FRPs Table 5 Average largest eigenvalues and confidence intervals of left and right swing intervals to enhance the discriminatory power of the convolutional eigenvalue method for intraclass identification. Such a future study, if being successful, will be very useful because it can overcome the difficulty for using machine learning when training data are limited and the cost for data acquisition can be expensive.