A summary of relevant information is presented in Supplementary Table 1, which displays the study category, nonlinear dynamical analysis, number of participants, the number of EEG channels, experiment, and major findings.

A variety of machine learning models to diagnose MS have been developed using MRI data, fMRI data, EEG data, etc. Since the symptoms, severity, and progression of MS vary enormously among individuals, each patient’s prognosis and subsequent treatment decisions should be tailored to their initial conditions. A machine learning algorithm can enable the search and analysis of large datasets about potential biomarkers to help find a clinically useful course of treatment (Hossain et al., 2022).

Seven articles used feature extraction, feature selection, and feature classification in their studies (Ahmadi and Pechenizkiy, 2016; Torabi et al., 2017; Kotan et al., 2019; Raeisi et al., 2020; Karaca et al., 2021; Karacan et al., 2022; Mohseni and Moghaddasi, 2022).

To gain a comprehensive understanding of the methods employed, Table 2 provides a detailed overview of the feature extraction, selection, and classification methods each article employed, along with the nonlinear analyses used to analyze the data. We have also provided a list of definitions related to the concepts of chaos and complex systems in Table 3 as a reference for common terminology. This comprehensive summary highlights the diverse approaches adopted in the studies. The following is a description of the chaotic measures used in studies combining nonlinear dynamical methods with classification:

Wavelets were developed in the 1980s as a substitute for Fourier transform. Currently, wavelets are used for analyzing signals, image processing and recognition, turbulence, etc., in a variety of different fields such as physics, fluid dynamics, and others (Pavlov et al., 2012). Continuous wavelets are useful in explaining nonstationary signals where their spectral composition and statistical characteristics fluctuate over time. According to Pavlov et al. (2012), wavelets:

There is a smoothing effect in MS data processing that poses a challenge. Large signals can exhibit these localized features as peaks or noise. To separate the true features from noise without bias, a few studies (Karaca et al., 2021; Karacan et al., 2022) have used wavelet transforms for MS EEG to capture these localized features in different resolutions. The continuous wavelet transform (CWT) is a mathematical operation that involves multiplying the original signal by a wavelet function ψ that is scaled and shifted across all time. The CWT coefficients are obtained using Equation 1 where f (t) represents the original signal, a is the scaling coefficient, and b is the shift coefficient (Karaca et al., 2021).

The wavelet function, ψ_a,b (t), is the complex conjugate of the wavelet function and is defined as follows:

Nonlinear and nonstationary signals can be analyzed using the empirical mode decomposition (EMD) method. EMD decomposes nonlinear signals into intrinsic mode functions (IMFs). IMFs refer to output signals whose amplitudes and frequencies are slowly shifting. The components of the main signal can be reflected in IMFs. It is important to note that some of these components are of limited use, while others reflect the characteristics of the original signal. Thus, selecting the appropriate IMFs is crucial (Kotan et al., 2019).

Empirical mode decomposition (EMD) breaks down the non-periodic and non-stationary signal, X_DFT (t), into a limited set of IMFs and a residue, r_N (t), as shown in Equation 3.

In Equation 3, N represents the total number of IMFs [IMF_j (t) denotes the j^th IMF], and r_N (t) corresponds to the residue obtained by selecting N IMFs (de Santiago et al., 2018).

The IMFs must satisfy two key conditions: firstly, the number of extremes (maxima and minima) and the number of zero crossings should either be equal or differ by no more than one throughout the dataset. Secondly, the mean value of the envelope formed by the local maxima and the envelope formed by the local minima should be zero at each point. In essence, the IMFs exhibit nearly periodic behavior with a mean of zero (de Santiago et al., 2018).

To decompose the signal, X (t) = X_DFT (t), into its constituent IMFs, the following four-step method is employed:

This process is repeated iteratively until the resulting signal, c(t), satisfies the criteria of an IMF. At this stage, c (t) becomes IMF₁, and the residue, r (t) = X (t) − c (t), replaces the original input signal for the subsequent step (1). Therefore, X (t) = r (t) (de Santiago et al., 2018).

In one of the selected articles, Kotan et al. (2019) used the EMD approach to distinguish between MS and healthy controls. Considering that one of the major causes of disability in MS patients is visual deficits, Raeisi et al. (2020) proposed a method based on visual stimulation and EMD for classification.

Torabi et al. (2007) employed the following feature extraction methods to find the most appropriate combination of nonlinear dynamical methods and classifier methods to differentiate MS patients from healthy subjects (Torabi et al., 2007).

Lempel-Ziv (LZ) complexity is a measure used to analyze the complexity of discrete-time physiologic signals, such as frequency, number of harmonics, frequency variability of signal harmonics, and signal bandwidth (Aboy et al., 2006). The calculation of the LZ complexity steps involves graining the time series X (x₁, x₂, …, x_n) into subsequences using a threshold value, which is always the average quantity. The signal is divided into two parts based on the threshold, assigning 1 to data larger than the threshold and 0 to data smaller than the threshold (Zhao et al., 2023).

This graining process creates a sequence P ={s₁, s₂, …, s_n}. Subsequences S and Q are then connected to form SQ. For example, if S = {s₁, s₂, …, s_r} and Q = {s_j, …, s_{j + m}}, then SQ is SQ = {s₁, s₂, …s_r, s_j, …, s_{j + m−1}}. If Q is not a subset of SQ, it becomes a new subsequence, and if Q is a subset of SQ, it is reconstructed as Q = {s_j, …, s_{j + m + 1}}. The number of times Q becomes a new subsequence is indicated as c (n). Subsequently, the complexity is normalized by the following equation:

However, since coarse graining ignores details of the original time series, a refined graining method is developed using median and quartiles as boundary points to divide the signal into four sections:

M denotes the median. Lastly, the refined Lempel-Ziv complexity (Equation 7) can be defined as:

Lempel-Ziv (LZ) has been applied to analyze the reactions of neurons in the primary visual cortex when exposed to various stimuli (Szczepaski et al., 2003). Also, it has been employed to evaluate the entropy of neural discharges, commonly referred to as spike trains (Amigó et al., 2004). LZ complexity has several advantages, and they are as follows:

In EEG signals, the reliance between amplitude and frequency plays an important role in the calculation since the binarization is led by slow rhythms. Thus, slow rhythms, or high amplitudes, impact the mean or median more than faster rhythms or low amplitudes (Ibáñez-Molina et al., 2015). LZ complexity is a suitable measure to aid in further understanding the complex and random nature of the EEG signals obtained from subjects with MS (Torabi et al., 2017).

Entropy is a measure used to examine the ambiguity of an information source and the probability distribution of the source’s samples. Entropy can be an indicator of the complexity of a system. There are several types of entropy analyses, such as approximate entropy (ApEn) and sample entropy (SampEn) (Ma et al., 2018). Quantifying irregularity in time series is done using ApEn (Pincus, 1991). However, ApEn requires a large dataset, and if the data length is short, the method is not robust enough (Richman and Moorman, 2000). SampEn is defined as the negative logarithm of the conditional probability that two sequences of m points remain similar when m+1 is reached, with each vector being counted over all others except itself. Thus, unlike ApEn, SampEn is relatively consistent irrespective of dataset size (Richman and Moorman, 2000). Mohseni and Moghaddasi (2022) used sample entropy in the feature-extraction step of their proposed method to create an MS diagnostic tool (Mohseni and Moghaddasi, 2022). See reference Li et al. (2018) for more information on entropy, its measures, and variants.

This measure appraises long-term memory processes in a time series. To analyze a time sequence X (x₁, x₂, …, x_n) with continuous values, the first step is to take the logarithm of the sequence and then perform a single differentiation of M_i.

This results in a logarithm sequence. The logarithm sequence is then divided into “A” adjacent subsets, with each subset having a length of h=(n- 1)A. Within each subset, the mean value is denoted as e_a and the standard deviation as S_a, where a = 1, 2,…, A. Within each subset, the accumulated intercept of the mean for each previous k point is calculated as follows:

The fluctuating range of each subset can be obtained by R_a:

Next is rescaling the range:

By increasing the value of h, we can obtain the rescaled range of subsets with varying lengths. The Hurst exponent describes the proportional relationship between (R/S)_h and h as (R/S)_h = c × h^HE, where c is a constant and HE is the Hurst exponent. Thus, the Hurst exponent can be estimated by plotting the logarithm of (R/S)_h against the logarithm of h. The slope of the fitted line in the plot corresponds to the Hurst exponent (Zhao et al., 2023).

The value of the exponent ranges between 0 and 1. There are three classification categories:

Kotan et al. (2019) extracted the Hurst exponent from the collected EEG signals in MS patients to examine the self-similarity of the data.

The Lyapunov exponent is a measure used to evaluate the average exponential divergence or convergence of nearby trajectories in phase space. A positive Lyapunov exponent indicates chaos in the system (Yakovleva et al., 2020).

In formula 12, the equation k(t) represents the measure of divergent distance, where K represents the initial distance and λ₁ denotes the largest Lyapunov exponent (Zhao et al., 2023).

In this approach, a proper feature selection is necessary for improving modeling efficiency and laying the groundwork for subsequent steps (Wang et al., 2013). As a final step, different classification methods will be used to identify MS patients from controls, including k-nearest neighbors (KNNs), support vector machines (SVMs), artificial neural networks (ANNs), and random forests (RFs), among others (Ahmadi and Pechenizkiy, 2016; Torabi et al., 2017; Raeisi et al., 2020).

Fractal analysis is a measure of the complexity and self-similarity of a signal that is investigated in the phase space via the attractor dimension or other correlated parameters, and it is able to analyze information at shorter intervals when compared to other linear and nonlinear analyses (Accardo et al., 1997). More specifically, fractal analysis quantitatively assesses the roughness of neural structures, approximation of time series, and representation of patterns. It is able to differentiate different brain states in the physiopathological spectrum (Di Ieva et al., 2015). According to Accardo et al. (1997), there are several ways to calculate the fractal dimension (i.e., box counting, Walker’s Ruler, and a modified correlation dimension D2 analysis) (Mandelbrot and Mandelbrot, 1982; Pickover and Khorasani, 1986; Tirsch et al., 1991; Gonzato et al., 1998). Higuchi’s is the most widely used method for determining fractal dimension.

The calculation of Higuchi’s fractal dimension involves analyzing a time series data sequence, denoted as X(1), X(2), …, X (N), where N represents the total number of samples. The process begins by selecting a scale factor, m, which determines the length of the subseries for analysis and k represents the index of the subseries being analyzed. For each subseries, the cumulative length, L(m, k), is computed using a specific formula that compares the absolute differences between adjacent data points within the subseries (Porcaro et al., 2020):

N represents the original time series X’s length, and N-1i⁢n⁢t⁢(N-mk) is a factor used to normalize the function. The cumulative lengths are then averaged across all subseries to obtain L (k), the average length for the given scale factor:

Finally, the Higuchi fractal dimension is calculated by taking the logarithm of L (k):

The resulting fractal dimension value represents the fractal dimension of the time series, providing insight into its complexity.

The fractal dimension is always expected to be between 1 and 2 since the dimension of a plane is 2, and the dimension of a line is 1. The fractal dimension increases when an EEG line fluctuates as more of the plane is covered (Accardo et al., 1997). A fractal dimension analysis of an MRI scan clearly shows the changes in white matter and gray matter in the early stages of MS. Clinical decision-making could be supported by the use of this approach as an early diagnostic biomarker (Di Ieva et al., 2015). Torabi et al. (2017), Kotan et al. (2019), and Mohseni and Moghaddasi (2022) extracted nonlinear features such as fractal dimensions to explore the nonlinear nature of EEG signals and sub-bands in MS patients (Torabi et al., 2017; Kotan et al., 2019; Mohseni and Moghaddasi, 2022). Porcaro et al. (2019) used fractal dimensions to analyze task-based EEG to develop a better treatment method (Porcaro et al., 2019).

Recurrence is defined as stretches of long or short repeated patterns that exist in living and non-living systems. As signals become more complex, recurrence becomes more uncommon (Webber and Zbilut, 2005). Analyzing recurrence patterns requires a mathematical or quantification analysis, which paves the way for using RQA (Webber and Marwan, 2015). Five variables are considered in RQA: percent recurrence (REC) or recurrence rate (RR), percent determinism (DET), maximal length in the diagonal direction (Dmax), Shannon entropy of the frequency distribution of the diagonal line lengths (ENT), and trend (TND) (Webber and Marwan, 2015). These measures are mainly concerned with the lengths, numbers, and distributions of the diagonal lines in recurrence plots. Three additional variables were added to examine intermittency and chaos-order transitions: laminarity (LAM), the average length of vertical structures (trapping time or TT), and the maximal length of the vertical structures (Vmax). These RQA parameters can be obtained from the recurrence plots (RP) of the EEG signals using the following formulas (Khodabakhshi and Saba, 2020):

Where N is the total count of states and x_i is the considered state.

Here, l_min, refers to the length of the smallest diagonal line, while P (l) represents the distribution of frequencies for different lengths (l) of diagonal lines. The D_max, ENT, and LAM, respectively, are given by:

In this context, P (v) refers to the histogram that captures the distribution of lengths, (v), of vertical lines. In addition, TT and V_max are calculated by:

Several studies have compared EEGs from patients with those from healthy controls using RQA identifiers, seeking to identify a relationship between baseline EEG and MS status (Carrubba et al., 2012, 2019). RQA is useful in various nonlinear datasets and can also be coupled with principle component analysis (PCA) (Webber and Marwan, 2015).

Systems in real life that have many components interacting nonlinearly usually require correlation analysis. Similarly, correlation analysis has evolved into complex network analysis (Xie et al., 2019). Finding synchronization between the signals of a neural system, such as the brain, has a critical role in characterizing the system activities and the integration of information within and across a disorder. The visibility graphs (VGs) algorithm enables the mapping of time series to complex networks. In the case of an EEG signal {x⁢(t)}t=1N with N data samples, each sample is represented as a node in a histogram-based graph. The height of each histogram bar corresponds to the value of the respective data node. Connections between nodes exist if the tops of the bars are visible. For two nodes (t_i, x_i) and (t_j, x_j), an edge is formed between them if there exists any data node (t_k, x_k) between (t_i, x_i) and (t_j, x_j) that satisfies the convexity criterion (Lacasa et al., 2008):

The horizontal visibility graphs (HVGs) algorithm is a modified version of the VG algorithm. In HVG, two data nodes (t_i, x_i) and (t_j, x_j) are considered to have horizontal visibility if they meet the condition stated in Equation 24 (below):

Therefore, (t_k, x_k)represents a data node located between (t_i, x_i) and (t_j, x_j). The complex network is represented by an adjacency matrix A = (a_ij)_{N × N}, where a_ij = 1 if nodes t_i and t_j are connected, and a_ij = 0 if they are not (Ahmadi and Pechenizkiy, 2016).

The application of synchronization measures, such as VGs and HVGs, to capture any intrinsic interactions between two-time series that are not stationary can thus prove beneficial (Dong et al., 2019). Thus, synchronization measures that capture the key features of the system can significantly contribute to our understanding of it. Based on EEG signals from healthy controls and MS patients, the authors suggest the successful application of the HVG method to construct a synchronization matrix of the brain network (Ahmadi and Pechenizkiy, 2016).

Mutual information is a measure of nonlinear dependence between two random variables. It is reliant on the random variable that diverges from random chance. If both random variables are independent, the joint entropy of their variables equals the sum of the marginal entropies. Dependence is observed when the joint entropy is less than the sum of the marginal entropies. If the joint distribution of two random variables is given by p (X, Y) and their factored marginal distributions are p (X) and p (Y), the formula for mutual information, I (X; Y) is as follows (Moermans et al., 2018):

The values for mutual information can range from 0 to infinity, where 0 represents total independence between the random variables, and infinity represents two correlated and continuous random variables (Smith, 2015). The effect of treatment on fatigue in MS was evaluated by calculating mutual information between the somatosensory and motor cortex in each hemisphere (Porcaro et al., 2019). Other studies have also used mutual information to evaluate and compare two EEG signals (Lenne et al., 2013; Tramonti et al., 2018).

Coherence has been extensively employed across a range of research disciplines, including but not limited to time-series-based studies, physics, and image processing, for the purpose of measuring the linear synchronization between two-time series. In simple terms, coherence denotes the interrelationship of two-time series at a specific frequency. One important application of coherence is in linear filtering, particularly in the context of EEG analysis that involves the presence of noise (Lenne et al., 2013).

In this context, for EEG data, coherence calculations can determine neural population synchronization levels among different brain regions (Nunez et al., 1999). The classic expression for coherence is:

In this equation, s_xyω indicates the cross-spectrum between x and y signals in the ω frequency (Jouzizadeh et al., 2021).

Coherence seeks to understand the information transfer between a known variable and an unknown variable. The measure increases with dependency between the two variables (Lenne et al., 2013). In progressive MS patients, EEG coherence shows a significant decrease between the anteroposterior and interhemispheric areas in alpha and theta bands, which correlates with cognitive dysfunction and subcortical lesion severity found on MRI results (Leocani et al., 2000; Lenne et al., 2013).

Figure 6 lists the name and frequency of each method used in the studies included in this review. The most common method of chaos quantification was fractal dimension analysis.

Kotan et al. (2019) compared three methods for selecting IMFs for EEG signals; power-based, correlation-based, and power-spectral density-based combined with Hurst exponent and Higuchi fractal dimensions. Then, three classifiers: KNN, multilayer perceptron neural networks, and RF, were employed to detect cognitively impaired and cognitively intact MS patients using the nonlinear features of Hurst exponent and fractal dimension. The authors concluded that IMF selection methods have a variable effect on accuracy based on classifier selection. The findings of their study indicate that different methods for selecting IMFs have varying effects on the classification accuracy of EEG signals collected from multiple sclerosis (MS) patients. Therefore, it is recommended to employ various IMF selection methods along with different classifiers and different types of signals to better understand their effects on classification accuracy.

Using RQA, Carrubba et al. (2019) developed a method for detecting MS from EEG signals. They assessed RQA application by statistically comparing the values of its quantifiers of EEGs from patients having MS with values from the EEGs of healthy subjects. They found that optimal embedding dimension and a time delay of 5 points maximize RQA’s ability to detect deterministic activity in EEG time series (Carrubba et al., 2019). In the MS patients, the RQA quantifier values were significantly higher than those in the healthy controls, indicating that the disease is associated with detectable changes in the EEG. Accordingly, a decrease in complexity associated with MS was observed in EEG recurrence plots (Carrubba et al., 2010). In another study by Carrubba et al. (2012), the same RQA method was employed. They proposed that EEG signals represent an instantaneous sum of contributions from several neural networks that elicit intra- and internetwork interactions. Therefore, the authors used RQA quantifiers to detect any changes from the “unknown, but certain laws” that brains without MS usually follow. The authors examined the RQA quantifiers, percent recurrence (%R), and percent determinism (%D). They observed an increase in %R in MS patients and no difference in %D in MS patients compared with healthy subjects, suggesting that %R is an indicator of nonlinearity (Carrubba et al., 2012). Therefore, due to their ability to reveal evoked potentials, %R and %D are possible indicators for MS diagnosis.

Lenne et al. (2013) also analyzed functional information about interhemispheric and intra-hemispheric cortical communication. As part of their methodology, they compared the coherence between different EEG bands and the mutual information between bipolar EEG signals between patients and healthy controls. The main outcome of cortical communication impairment in MS was a significant decrease in mutual information between brain areas. This study indicated that averaged interhemispheric mutual information in the resting state may indicate neurological dysfunction in MS patients. In addition, these findings indicate that this nonlinear measure may serve as an indicator of MS patients’ information processing deficits. Using mutual information, one can thus characterize the connectivity of brain information transmission (Lenne et al., 2013).

The results reported by Leocani et al. (2000) align closely with those reported by Lenne et al. (2013). The authors observed a reduction in coherence in EEG activity between the anteroposterior and interhemispheric areas in the cognitively impaired MS patients. Some findings revealed that there was unusual synchronization and hyperconnectivity in certain frequency bands, which could be due to the compensatory mechanisms or pathological abnormalities associated with MS. These findings emphasized that coherence analysis could be used as an indicator for cognitive dysfunction connected to multiple sclerosis (MS) while understanding changes that take place in the brain. Thus, cognitive impairment in MS is contingent on the corticocortical connection tied to the demyelination and/or the axonal loss that lies beneath the cortex in the white matter (Leocani et al., 2000).

Tramonti et al. (2018) used weighted symbolic mutual information (wSMI) to evaluate non-random joint fluctuations between two EEG signals after participants underwent gait rehabilitation. Their analysis provided strong results for understanding behavioral changes and showed that phase synchronization is correlated with functional recovery (Tramonti et al., 2018).

Carrubba et al. (2010) successfully applied chaotic measures for the diagnosis of MS. By introducing a subliminal stimulus in patients with MS and those without and through RQA quantifiers, they tested whether cognitive processing was altered in patients with MS, finding that MS patients had 27% onset responses in EEG signals and subjects without MS had 85% onset responses. Using this study’s methodology, one can assess the degree of synchronization between brain networks, which is exactly the high-level brain function they believed to be impaired (Carrubba et al., 2010).

Mohseni and Moghaddasi (2022) also employed nonlinear feature extraction and classification. Through segmenting random signals in EEG signals, they proposed a model for diagnosing MS. While the alpha, beta, and gamma sub-bands have an impact on signal analysis, the robustness of the technique with random selections of 15-to-30-second segments of the EEG signal and the detecting power of the algorithm makes it appropriate for generalization (Mohseni and Moghaddasi, 2022).

Jouzizadeh et al. (2021) applied coherence analysis to understand the functional connectivity of MS. They discovered differences in the functional connectivity between patients with MS and healthy subjects. The information transfer and increase in communication distances were observed due to the deficit in global efficiency and the increase in path size. To add, it was determined that MS has a large impact on the brain since the connectivity patterns in the frontal, parietal, and occipital regions were all affected, which indicates there would be issues with the cognitive, sensory, and visual processing areas. The authors noted that the results produced similar results to those generated by fMRI, and it was also determined that betweenness centrality and small-world propensity are effective indicators in differentiating an individual with MS from an individual without MS (Jouzizadeh et al., 2021).

Buyukturkoglu et al. (2017) aimed to investigate the functional connectivity patterns in the brain related to fatigue in multiple sclerosis (MS) patients. They focused on the concept of functional connectivity, which refers to the degree of coherence or synchronization between different brain regions. The main finding of the study was that a potential EEG-Neurofeedback system for MS fatigue should train patients to voluntarily decrease coherence in the beta frequency band between homologous temporoparietal cortices. By targeting and modulating this specific brain organization feature, which becomes more altered as fatigue symptoms worsen, the researchers believe that the symptom of fatigue can be ameliorated. The choice of the beta frequency band was based on its involvement in motor processing and its potential relevance to the experience of fatigue. Additionally, the researchers considered that training in the alpha and/or beta bands might yield better results compared to slow delta or high gamma bands, as the latter can be affected by eye artifacts, which are common during visual tasks involved in EEG-Neurofeedback training (Buyukturkoglu et al., 2017).

With that, the studies in this section call attention to the variety of methods and techniques used for analyzing EEG signals collected at a resting state. Methods and techniques, such as IMF selection methods, recurrence quantification analysis (RQA), coherence analysis, mutual information, and mutual information, were all used in their respective experiments to extract nonlinear features and evaluate the functional connectivity in MS patients. These studies show how these measures can be used for detecting changes caused by MS in EEG signals and recognizing signs of cognitive dysfunction, neurological impairment, and information processing deficits accompanied by MS. The variety in the findings highlights the importance of identifying the appropriate analysis to understand the EEG activity in subjects with MS.

Raeisi et al. (2020) used an extension of bivariate empirical mode decomposition (BEMD) to analyze a task-based experiment (three visual stimuli were used to distinguish MS from healthy groups in the study). Each possible pair of EEG recordings was pre-processed to extract five pairs of IMFs by the bivariate empirical mode decomposition (BEMD) method. Phase synchrony between each pair of IMFs was calculated utilizing the mean phase coherence (MFC) method and Hilbert transform. To classify healthy and MS subjects, ReliefF and KNN classification methods were applied. Based on the results, the accuracy, sensitivity, and specificity of the red-green task were 93.09, 91.07, and 95.24%, respectively, while those of the black-white task were 90.44, 88.39, and 92.62%, and those of the blue-yellow task were 87.44, 87.05, and 87.86%. The experimental results indicated that the method could be generalized to automate MS diagnosis systems (Raeisi et al., 2020).

Karacan et al. (2022) focused on classifying the environment according to EEG signals. A cognitive task was performed on a computer by healthy volunteers and the MS patient group, and then the same task was performed in a virtual reality environment. After extracting chaotic entropies and fractal dimensions, the accuracies of different classification methods were compared. A KNN classifier performed best for volunteers with MS with an accuracy of 95.45%. Additionally, after inducing photic stimulation in different EEG subbands [delta (1–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), and beta (13–30 Hz)] and feature extraction, the proposed ensemble subspace KNN classifier algorithm performed the automatic classification of MS with an accuracy of 80%, a sensitivity of 72.7%, a specificity of 88.9%, and a positive predictive value of 88.9%. The authors proposed using photonic stimulation EEGs for the pre-diagnosis of MS (Karacan et al., 2022).

The nonlinear model proposed by Torabi et al. (2017) classified healthy individuals and individuals with MS into two groups. The authors used the nonlinear features of EEG signals to distinguish healthy volunteers and MS patients performing a color-and-luminance-changing task and a direction-changing task. The Katz fractal dimension was the most “informative nonlinear feature” when combined with an SVM. Their model achieved high classification performances of 93.08 and 79.79% for the direction-based and the color-luminance-based tasks, respectively (Torabi et al., 2017).

Porcaro et al. (2019) sought an alternative treatment method to address fatigue brought upon by the side effects MS patients experience with pharmaceutical treatments. They carried out a five-day transcranial direct current stimulation (tDCS) focusing on the somatosensory representation of the body (S1), known as FaReMuS treatment. Utilizing fractal dimension and mutual information, 48% of the variance of fatigue was explained, paving the way for personalized neuromodulation techniques that can be used to address fatigue in patients with MS. The authors observed that the left side of the whole-body somatosensory area (S1) showed a significant change after neuromodulation. The nonlinear analyses helped demonstrate that left S1 was impaired in MS patients prior to treatment, and the difference between S1 in MS patients vs. healthy patients disappeared after treatment, which is a significant step toward identifying an improved treatment (Porcaro et al., 2019).

Tomasevic et al. (2013) The findings of this study, in terms of coherence analysis, indicate that indices related to movement execution are significantly associated with fatigue rather than morpho-structural measures related to the primary sensorimotor network. Specifically, fatigued patients exhibited cortico-muscular coupling at faster frequencies and corrected the pressure exerted during handgrip at higher frequencies. The disruption of primary somatosensory network patterning in MS indicates that intra-cortical synchronization phenomena affecting cortico-muscular coupling also play a significant role in motor control. The study suggests that the quality of communication between the cortex and muscles is impaired in fatigue, while the spectral features of the motor cortex and muscular oscillatory activities remain unaltered (Tomasevic et al., 2013).

The studies in this section underscore the diverse array of methods and techniques used for analyzing EEG signals collected while participants completed a task. The studies provide valuable insight into the growing body of knowledge regarding using EEG in the diagnosis, treatment, and understanding of MS. EEG-based techniques provide essential tools for academics, physicians, and patients in the pursuit of better management and care for people with MS by using modern signal processing, nonlinear analysis, and creative intervention strategies.