The data were collected at Jagiellonian University in Krakow, Poland. The study included two groups of subjects: patients with early onset relapsing–remitting multiple sclerosis (RRMS) and healthy subjects. In the group of MS patients, there were two subgroups: those treated with IFN-β and those treated with DMF. The total number of participants for this analysis is 175. To further break it down, 39 patients were on IFN-β, 53 were on DMF, and there were 83 healthy controls. The IFN-β group consisted of participants between 22 and 63 years old (M = 39.15, SD = 7.909), and there were 24 females (61.5%) and 15 males (38.5%). The DMF group contained participants between 18 and 54 years old (M = 32.11, SD = 7.250). This group had 33 females (62.3%) and 20 males (37.7%). The participants in the control group were between 21 and 61 years old (M = 36.22, SD = 8.498). There were 53 females (63.9%) and 30 males (36.1%). There were two rounds of data collection (first measurement and second measurement). The data for the second measurement were obtained 1 year after the data for the first measurement were collected. MS patients’ Expanded Disability Status Scale (EDSS) scores (Kurtzke, 1983) ranged from 1 to 4 in the first measurement and from 1 to 4.5 in the second measurement. The number of relapses in the year prior to each measurement ranged from 0 to 2. A Wilcoxon signed-rank test indicated that there was no significant difference between EDSS scores in the first and second measurements, z = −0.958, p = 0.338. The median EDSS score was 1 in both the first and second measurements. Similarly, there was no significant difference in the number of relapses in the year prior to each measurement between the first and second measurements, z = −0.915, p = 0.360. The median number of relapses in the year prior was 0 in both the first and second measurements. The control group did not undergo a second round of data collection because there should not be significant changes in resting state EEG in healthy subjects within 1 year (Kondacs and Szabó, 1999).

For this study, data were collected during a resting state task. The resting state task included a six-minute procedure without any stimuli. In the first 3 minutes, subjects were asked to have their eyes open while focusing on a fixation point, and they had to keep their eyes closed in the last 3 minutes. Commands were pre-recorded and played by speakers. A 256-channel dense array EEG system (HydroCel Geodesic Sensor Net, EGI System 300; Electrical Geodesic Inc., OR, USA) was used to collect the data. The researchers decided to remove channels located on the cheeks (E225, E226, E227, E228, E229, E230, E231, E232, E233, E234, E235, E236, E237, E238, E239, E240, E241, E242, E243, E244, E245, E246, E247, E248, E249, E250, E251, E252, E253, E254, E255, and E256) due to many artifacts of low interest in the signal.

The EEG data underwent pre-processing using MATLAB’s EEGLAB software to ensure data quality and integrity (Delorme and Makeig, 2004). The initial pre-processing stage involved discarding 5 seconds of data that followed sound commands—eliminating these potential artifacts or confounding effects because the experimental instruction allowed for a more precise analysis of the EEG signals. A high pass filter was employed to exclude any signals below the frequency of 0.5 Hz. Adding on, a notch filter to remove power line interference and its harmonics was integrated to reject 50 Hz and its multiplicities from the signal. Independent component analysis (ICA) was conducted. Fifty principal components were used for the analysis to identify and reject artifact components, such as components related to eye movements, muscle activity, or other sources of artifact. Every removed channel was interpolated to estimate the missing values based on surrounding electrodes and provide comprehensive coverage of all channels. Each subject had a sampling rate of 250 Hz for this study.

A commonly used linear analysis with applications in neurophysiological data, lag-1 autocorrelation (AC1), was carried out to validate the use of nonlinear analysis (Meisel et al., 2017; Scheffer et al., 2009). AC1 is a reliable measure of the rate at which the autocorrelation function decays (Huang et al., 2018). The autocorrelation function (ACF) is defined in Equation 1, where x t represents the envelope signals, N is the length, μ is the mean, and v is the variance:

To obtain lag-1 autocorrelation, we set s = 1 (Meisel et al., 2017). Higher AC1 values indicate greater predictability in the signal, whereas lower AC1 values suggest less predictability (Huang et al., 2018).

Sample entropy (SampEn), initially developed by Richman and Moorman (2000) to measure regularity, was used to analyze the EEG signals across all groups (Duran et al., 2013; Richman and Moorman, 2000). Greater entropy values indicate that the system is complex, irregular, and unpredictable, often associated with a healthy system. Conversely, low entropy values indicate a more deterministic and predictable system, meaning the EEG signals show more regular patterns and less complexity (Duran et al., 2013; Pincus, 2006; Delgado-Bonal and Marshak, 2019). Two notable parameters are used in calculating SampEn: m and r. The parameter m represents the length of the subseries, and r represents the similarity criterion (Ramdani et al., 2009). Following the guidance of Costa et al. (2005) and Duran et al. (2013) selected m = 2 and r = 0.15 as the parameters, and it was noted that the selection of the parameters does not negatively impact the overall pattern of the results (Costa et al., 2005; Duran et al., 2013). Thus, others typically default to the parameters Duran et al. (2013) used, as they are considered standard and, therefore, were deemed appropriate for this study. Following the guidance outlined by Ramdani et al. (2009), the equation for sample entropy is as follows (Richman and Moorman, 2000; Ramdani et al., 2009):

With time series x₁, x₂, … x_N, subsequences of length m are first defined in Equation 2:

After, the quantity is calculated by the following:

The Heaviside function is defined by Θ , and | | · | | ∞ represents the maximum norm, which is y j m − y i m ∞ = max 0 ≤ k ≤ m − 1 | X j + k − X i + k | . To explain, Equation 3 calculates the sum of the quantity of vectors, y j m , that are within the radius, r, from y i m that exist in the reconstructed phase space. Identical matches are excluded and are represented by j ≠ i . Also, N – m represents the total amount of vectors in the (m + 1) dimensional state space.

Equation 4 calculates the density:

Calculations in the (m + 1) space to extend the template matching process are as follows:

In Equation 5, the number of sequences y j m + 1 within radius r of y i m + 1 is calculated, with the term y j m + 1 − y i m + 1 representing the maximum difference between the two subsequences. After calculating the individual template matches A i m r , they are all averaged across all vectors to give A m r , as shown in Equation 6. Then, the total amount of template matches in a m-dimensional and m + 1 dimensional phase space with r is represented by Equations 7 and 8:

The sample entropy can then be calculated as follows in Equation 9:

The sample entropy MATLAB script provided by Richman and Moorman (2000) was used in conjunction with an unpublished modified script from Amon (2021) to conduct the analysis (Richman and Moorman, 2000; Amon, 2021).

Higuchi’s fractal dimension (HFD) was also employed to analyze the EEG signals. It is another method frequently used in nonlinear analysis, and it details the time series’ complexity and self-similarity (Accardo et al., 1997). Following the outline of the computation summarized in Hernandez et al. (2023), the calculation of HFD involves analyzing a time series data sequence, denoted as X (1), X (2), …, X (N), where N represents the total number of samples (Hernandez et al., 2023). The selection of a scale factor, m, begins the process. This scale factor, m, defines the length of the subseries under investigation. The selection of k is also necessary to commence the process, as this is the index of the subseries. The cumulative length, L(m, k), is calculated by comparing the absolute differences between adjacent data points within the subseries, as shown in Equation 10 (Porcaro et al., 2020):

N is the length of the original time series X and N − 1 i n t N − m k normalizes the function. The average cumulative length across all subseries is calculated to acquire L k , the average length for the given scale factor, as represented in Equation 11:

The Higuchi’s fractal dimension is then computed by taking the logarithm of L k , as defined in Equation 12:

The resulting fractal dimension value represents the fractal dimension of the time series, providing insight into its complexity (Porcaro et al., 2020). The method for calculating Higuchi’s fractal dimension was adopted from Jesús Monge-Álvarez¹.

Typically, the fractal dimension ranges between 1 and 2, where higher HFD values indicate greater complexity and lower values suggest reduced complexity (Accardo et al., 1997; Scarpa et al., 2017).

Currently, no standard method is used to select the most appropriate value for the k_max parameter (Kesić and Spasić, 2016). The method selected in this paper is a common method used by Doyle et al. (2004) and Wajnsztejn et al. (2016). They considered the most appropriate k_max parameter to be where HFD approaches a local maximum or asymptote (saturation point) (Wanliss et al., 2021; Doyle et al., 2004; Wajnsztejn et al., 2016). According to Figure 1, the data reaches a local maximum at k_max = 70. Therefore, k_max = 70 was the parameter chosen for this study.

For the analysis, each participant’s EEG signal was divided into short 15-s time windows with 50% overlap. This was decided after following the advice of several articles that have opted to divide EEG signals into short time windows for computational efficiency (Mohseni and Moghaddasi, 2022; Ramanand et al., 2004; Er et al., 2021; Kesić and Spasić, 2016). The 50% overlap was chosen to prevent any discontinuity at the frame’s beginning or end (Er et al., 2021).

Several statistical analysis techniques were used to understand the data and answer the research questions comprehensively. Descriptive statistics provided a summary of the data. Levene’s and Mauchly’s tests were conducted to test for homogeneity and sphericity. Although homogeneity was violated in most cases, it was not violated in the second measurement of AC1. There was no indication of a violation of sphericity. Given the sample size (n > 30) and following guidance from Hair et al. (2010) and Byrne (2010), parametric tests were utilized, as skewness (between −2 and + 2) and kurtosis (between −7 and + 7) were within acceptable ranges (Hair et al., 2010; Byrne, 2010). A paired samples t-test was used to compare the means within subjects, and mixed analysis of variance (ANOVA) was used to investigate the main effects of time and group. Welch’s ANOVA was employed to analyze the means between subjects to address the violation of homogeneity, and standard ANOVA was used to evaluate the means between subjects in the second measurement of AC1, where homogeneity was not violated. Games-Howell post hoc test was completed to identify which groups demonstrated significant differences. An alpha level of 0.05 was used as the threshold for determining the effect’s significance.

Lag-1 autocorrelation (AC1) was carried out to assess the linearity of the dataset. The mean AC1 value of the IFN-β group was 0.800 (SD = 0.044) in the first measurement and 0.815 (SD = 0.042) in the second measurement. For the DMF group, the mean AC1 value was 0.812 (SD = 0.052) in the first measurement and 0.805 (SD = 0.050) in the second measurement. The mean AC1 of the control group was 0.806 (SD = 0.034). Paired samples t-test revealed no significant differences in the means within the IFN-β group (t(38) = −1.676, p = 0.102) and DMF group (t(52) = 0.901, p = 0.372). According to the mixed factorial ANOVA, time did not have a significant effect, F(1, 172) = 0.727, p = 0.395. However, a significant interaction effect of time and group was reported F(2, 172) = 3.396, p = 0.036, highlighting a significant change in the pattern over time across groups. Due to the violation of homogeneity in the first measurement, F(2, 172) = 3.344, p = 0.038, Welch’s ANOVA was conducted for between-subjects comparison. No significant differences were reported in the first measurement, F(2, 82.498) = 0.651, p = 0.524. Since the data in the second measurement, F(2, 172) = 1.636, p = 0.198, did not violate homogeneity, standard ANOVA was carried out. Like in the first measurement, no significant differences were reported, F(2, 172) = 0.728, p = 0.484.

To assess the complexity of the EEG data, box plots with 95% confidence intervals were created to understand the distribution and central tendency of the SampEn and HFD values across different groups and measurements (Figure 2). Referring to the point plots in Figure 3, both treatment groups at each measurement had recorded relatively high mean SampEn values and HFD values compared to the control group. Summary statistics are shown in Table 2. A paired samples t-test was employed to evaluate the significance of the difference within each treatment group.

A mixed factorial ANOVA was conducted for SampEn and HFD to observe the main effects of time and group (control, IFN-β, or DMF). An interaction plot was created to visualize the effects.

Due to the violation of homogeneity, Welch’s ANOVA was performed for the between-subjects effect at the first and second measurements for both SampEn and HFD. A Games-Howell post hoc test was conducted to identify significant differences between groups.

As mentioned, higher entropy values indicate that a system is complex, irregular, and unpredictable, often linked to a healthy system. On the other hand, lower entropy values indicate a more predictable and deterministic system (Duran et al., 2013; Pincus, 2006; Delgado-Bonal and Marshak, 2019). As for HFD, greater values indicate more complexity in the signal (Scarpa et al., 2017). Treatment was expected to have some level of impact on the complexity of the signal (Shalbaf et al., 2012; Thomasson et al., 2000).

In the study, the control, Interferon-β, and dimethyl fumarate groups displayed high SampEn and HFD values at each time measurement, supporting the hypothesis that an increase of complexity was observed. It is shown that both treatment groups displayed higher SampEn and HFD values when compared to the control group, suggesting that the MS patients were found to have a greater number of nonlinear segments. These findings were similar to those of Pezard et al. (2001), who reported higher entropy values compared to the control group when investigating Parkinson’s disease (Pezard et al., 2001). This further reveals MS patients treated with IFN-β and DMF have less predictable and more complex electrical activity compared to the controls (Pezard et al., 2001). The high nonlinearity can also be tied to the dimensionality of the electrical activity. Lachaux et al. (1997) described how dimensionality decreases if nonlinearity increases (Lachaux et al., 1997). This indicates that the MS patients treated with both treatments may have brain dynamics of a lower dimension (Pezard et al., 2001; Stam et al., 1994). Additionally, it has been noted that the increase in the complexity of EEG signals for MS patients is linked to the brain’s compensatory mechanisms and is indicative of the brain’s structural complexity (Wątorek et al., 2024). We can hypothesize that the higher complexity reported in the treatment groups could also be due to the brain’s adaptive response to the effects of the treatments, as they are responsible for the regulation of the immune system and reduction in inflammation (Jakimovski et al., 2018; Linker and Haghikia, 2016; Mills et al., 2018).

There were no significant differences reported in the complexity characteristics of EEG signals between MS patients undergoing treatment with IFN-β and DMF at the first and second measurements, which rejects the hypothesis that patients treated with DMF will exhibit significant differences in complexity characteristics compared to patients treated with IFN-β. However, the second hypothesis was partially supported because the complexity characteristics (SampEn and HFD) of each treatment group differed significantly compared to the control group at each time measurement, as confirmed by Welch’s ANOVA and the Games-Howell post hoc test. These findings are backed by other studies that have concluded that nonlinear EEG measures can be sensitive to treatments (Pezard et al., 1998; Pezard et al., 2001; Wackermann et al., 1993).

In particular, as seen in Tables 6, 8, the mean differences in SampEn between each treatment group and the control group at the first and second measurements were higher than the mean differences observed in the same scenario for HFD. This indicates that SampEn demonstrated the highest sensitivity and the greatest predicted value in evaluating the effects of each treatment group compared to the control group, supporting our third hypothesis. These results suggest that treatments, such as IFN-β and DMF, impact the overall brain dynamics, as reflected by the higher sample entropy and Higuchi fractal dimension values.

Several studies (Hossain et al., 2022; Di Ieva et al., 2015) have investigated using nonlinear analysis in recognizing biomarkers in individuals with MS and healthy controls (Hernandez et al., 2023). Both entropy and fractal dimension have been used to either distinguish between conditions or differentiate between healthy and pathological brains in previous research (Marino et al., 2019; Smits et al., 2016; Zappasodi et al., 2014; Zappasodi et al., 2015; Bauer et al., 2011; Pezard et al., 2001). In this study, we aimed to explore whether sample entropy and HFD are reliable indicators for the progression of MS. We initially hypothesized that MS patients treated with IFN-β and DMF treatments would reveal significant and consistent changes over time relative to the control group.

Referencing Figure 2, it was observed that the initial measurements of SampEn and HFD demonstrated more dispersion compared to the second set of measurements. This observation could indicate the progression of MS over time, leading to more consistency in the results. Nevertheless, we determined that the hypothesis could only be partially supported because time and the interaction between time and treatment group significantly impacted only HFD and not SampEn. However, a significant increase from the first measurement to the second measurement was only observed in HFD values of the DMF group. Hence, an increase in signal complexity and positive neurophysiological changes can be attributed to DMF, which is reflected only in HFD. This finding is supported by Viglietta et al. (2015) and Vermersch et al. (2022). Viglietta et al. (2015) concluded that DMF reduces new and enlarging T2 lesions, gadolinium-enhancing lesions activity, and the number of new non-enhancing T2 lesions (Viglietta et al., 2015). Similarly, Vermersch et al. (2022) reported that more pediatric patients treated with DMF did not develop new or newly enlarging T2 lesions compared to those treated with IFN-β (Vermersch et al., 2022). These findings demonstrate the effectiveness of DMF in reducing disease activity and may explain the increase in EEG complexity over time compared to IFN-β. Although SampEn demonstrated the highest sensitivity and greatest predicted value, its responsiveness was limited when time was factored in. This finding signifies how HFD may be more responsive to temporal changes in EEG dynamics than SampEn.

There are a few limitations and opportunities for future research to note in this study. The first limitation is centered on the selection of the kmax parameter. Different methods of kmax parameter selection have been employed previously, but researchers have yet to agree on a universal method (Kesić and Spasić, 2016). Different parameter selection methods could alter the results. However, one of the most common methods was chosen in this study. This method was carried out by selecting the parameter where HFD reached a maximum or asymptote (Wanliss et al., 2021; Doyle et al., 2004; Wajnsztejn et al., 2016). Another possible limitation is the sample size of each treatment group. Increasing the sample size could have enhanced the results reported in this experiment. More specifically, the IFN-β treatment group had the lowest number of participants, and an increase in the number of MS patients on IFN-β could have highlighted clinically significant differences between the treatment groups.

There are several opportunities for future research. First, future studies could expand and balance the sample sizes for each treatment and collect longitudinal EEG data from the control group to strengthen the analysis and validate these findings. The next step in the study could be to analyze the EEG time series using multifractal methodology. This method helps quantify the data’s correlation structure through the set of scaling exponents, providing a deeper understanding of the data’s complexity (Wątorek et al., 2024). Furthermore, there are several methods to characterize complexity. One method is detrended fluctuation analysis (DFA), which is used to evaluate the Hurst exponent and can then be recalculated to determine the fractal dimension (Márton et al., 2014). Another method is the Lyapunov exponent, which is employed to identify chaotic behavior in the data and can be used to quantify data complexity (Yakovleva et al., 2020). The presented study investigates the effects of two immunomodulatory treatments; however, they aren’t the only treatments for multiple sclerosis. MS treatments include immunosuppressants (i.e., fingolimod), immunomodulatory therapies (i.e., IFN-β and DMF), and immune reconstitution therapies (i.e., alemtuzumab and cladribine) (Dobson and Giovannoni, 2019). Future studies could investigate the effects of immunosuppressants and immune reconstitution therapies on the brain’s dynamics via nonlinear analysis. These studies could use nonlinear analysis to investigate how these different treatment groups compare.

As reported by Hernandez et al. (2023), several articles have used machine learning algorithms in studying MS (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). Methods include feature extraction, feature selection, and feature classification, and these methods could allow researchers to swiftly search and analyze large datasets for potential biomarkers (Hernandez et al., 2023; Hossain et al., 2022). In future studies, researchers could build on this study’s approach by developing machine-learning methods that integrate MRI and functional magnetic resonance imaging (fMRI) to compare the efficacy of different MS treatments. This could further enhance the analysis by identifying trends and possible biomarkers more efficiently.