This section outlines the methodologies employed to compute the causal theoretical quantifiers for the LFPs acquired from iEEG. Accurately quantifying the information content within observed neural activity is vital for the analysis of neural systems. Integrating “permutation patterns” with various metrics offers a more comprehensive understanding of the characteristics of a time series (Olivares et al., 2020).

The initial step in quantifying information using causal entropy measures for a time series involves associating it with a probability distribution function (PDF) extracted using the Bandt and Pompe (BP) methodology. Each time series, χ(t) = {x_t; t = 1, ⋯ , M}, consists of M measurements of the observable χ. To analyze the time series, it is divided into n = M−(D−1)τ overlapping segments, where the chosen embedding dimension is D and the embedding delay is τ. This method is based on the construction of a histogram of ordinal patterns. A brief description is given here, as it has been studied extensively, but the reader is referred to the following references for a formal explanation, simple examples, and similar applications (Bandt and Pompe, 2002; Jauregui et al., 2018; Pessa and Ribeiro, 2021; Zanin and Olivares, 2021; Amigó and Rosso, 2023; Guisande et al., 2023).

The partitions are represented by a D-dimensional vector, which is used to determine the permutation of index numbers. All possible permutations of order D are considered, and the relative frequency of each permutation is calculated to obtain the distribution of ordinal patterns. Thus, the frequency of BP can be calculated using the following equation:

where p_j(Π_j) represents the relative frequency of the j−th ordinal pattern, Π_j, and π_i is the sequence of all ordinal patterns for all partitions. Note that the estimated PDF is discrete since it is computed using a histogram. It is also important to emphasize that for reliable statistics, the analyzed time series, χ(t), needs to be much longer than the total number of possible ordinal patterns (M >> D!. where D! represents the total number of possible ordinal patterns of size D).

This study employs the method of ordinal patterns to identify and quantify the existence of ordinal structures in time series data. The objective is to compare computationally simulated data, based on the power law (k-noise) and experimental data by analyzing the dynamics of physiological signals in normal brain regions across biological sexes, while taking into account the causality of the signals. Using the resulting probability distribution, the causal entropy (also known as permutation entropy) and causal complexity (also known as permutation complexity) of Rényi are computed. In this study, any mention of entropy or complexity refers specifically to these two quantifiers.

The Rényi entropy serves as a measure of uncertainty within a probability distribution. It offers a unique advantage over Shannon entropy, especially in the context of evaluating brain activity, due to its exceptional versatility in capturing diversity within complex systems. While Shannon entropy provides a single information index, Rényi entropy employs a parametric family of indices, allowing for a more comprehensive sensitivity to both rare and abundant elements (Jost, 2007).

The Rényi entropy is defined as (Rényi, 1961; Pessa and Ribeiro, 2021):

In the context of BP, N represents the total number of possible states in the probability distribution, corresponding to the number of potential ordinal patterns (N = D!).

The q factor in the Rényi entropy formula functions as a weighting coefficient for the probability distribution. It indicates the order of entropy and determines the specific type of Rényi entropy being calculated. By adjusting the value of q, the Rényi entropy formula can reveal various aspects of the distribution, such as inclination toward rare events or focus on common events. Essentially, q facilitates the analysis of different characteristics of the probability distribution based on its value.

If the probabilities are uniform, all Rényi entropies of the distribution are equal, given by S_q = ln N. For non-uniform distributions, the entropies weakly decrease with respect to q. As q approaches 0, it converges to the max-entropy (S_{q → 0} = ln N), assigning equal weight to all possible events regardless of their probabilities. When q = 0, it represents the logarithm of the size of the support of χ. As q approaches ∞, it converges to the min-entropy (Sq→∞=-ln maxj(pj)), considering only events with the highest probability. The intermediate case at q = 1 corresponds to Shannon entropy (Sq=1=-∑j=1Npjln(pj)); (Zmeskal et al., 2013; Zhou and Zheng, 2022).

The adaptability of Rényi entropy makes it a useful measure for quantifying various levels of information and capturing the dynamics of non-stationary processes. By utilizing ordinal patterns, it enables the tracking of changes in entropy distribution over time (A-iyeh and Peters, 2016; Shalymov and Fradkov, 2016; Jauregui et al., 2018). Furthermore, Rényi entropy proves highly effective in describing multifractal systems (Jizba and Arimitsu, 2001) and maintains a close relationship with the scale-free exponent, making it a valuable tool for investigating the scale-free dynamics of the brain.

The human brain, with its intricate network of interconnected neurons, exhibits multifractal characteristics across various spatial and temporal scales (França et al., 2018). Hence, one of the primary reasons for selecting Rényi's entropic quantifiers as a central component of this approach is their effectiveness in describing multifractal systems. This approach can provide valuable insights into brain systems across different spatial and temporal scales, as well as varying levels of complexity (Tozzi et al., 2018).

Following the methodology of Martin et al. (2006), the Rényi statistical complexity is defined as:

where H_q is defined as

It is important to note that D_q (Jensen-Rényi divergence) is determined by

and Dq* (a constant normalization representing the maximum value of D_q) is given by

To model and analyze patterns of brain activity exhibiting scale-free structures, correlated noise (k-noise) time series following a power law distribution (f^−k) were computationally simulated. These time series simulate the intrinsic dynamics of scale-free neuronal activity, allowing modeling of the diverse and evolving nature of the brain across different time and frequency scales.

The computational generation process involved three essential steps: First, a set of k values representing the exponent of the power law distribution was defined using a logarithmic scale. Random noise data points were then generated and transformed from the time domain to the frequency domain using the Fast Fourier transform (FFT). The FFT decomposed the signal into its frequency components, which were then multiplied by a function that followed a specific power law distribution. This operation adjusted the amplitudes of different frequencies to match the desired power law distribution. Finally, the inverse Fourier transform was applied to obtain the resulting time series with power law characteristics.

After generating the time series, the BP methodology was used to associate these data points with a PDF. Subsequently, quantitative measures such as Rényi entropy and Rényi complexity were calculated for analysis purposes.

The time series corresponding to the LFPs utilized in this study were sourced from the MNI Open iEEG Atlas database. This database comprises recordings of intracranial activity in typical brain regions during various states, including quiet wakefulness with eyes closed (W), non-REM sleep stage N2 (N2), non-REM sleep stage N3 (N3), and REM sleep (R; Frauscher et al., 2018a,b; von Ellenrieder et al., 2020).

A total of 106 patients with focal epilepsy were included in the atlas, and recordings were made from 1,772 channels. Only channels located in gray matter and considered “normal” (far from epileptic regions) were utilized. Various types of intracerebral electrodes were employed, including Dixi, homemade MNI, and AdTech electrodes, as well as AdTech subdural strips and grids. Importantly, all recorded signals from these electrodes were incorporated into the analysis in this work without any differentiation between them.

The dataset contains patient information, including sex, channel type, hemisphere, channel name, channel position, and channel region. Additionally, all signals were resampled to 200 samples per second to ensure consistency. Power-line interference was minimized using an adaptive filter. Furthermore, all channels were zero-padded to a length of 68 s (13,600 samples) to maintain uniformity across segments, regardless of their number.

To facilitate the comparison of patient activity and the accumulation of results from multiple subjects, the electrodes were placed in a common stereotactic space. For more detailed information about the database and acquisition methods, readers can refer to the following references (Frauscher et al., 2018a,b; von Ellenrieder et al., 2020).

The signals were classified based on hemisphere, region, and biological sex, allowing for the analysis of each region in both hemispheres separately for female and male patients. To ensure a sufficient sample size for comparing behavior between males and females within each region, a criterion was established, requiring a minimum of five patients of each sex. This criterion was consistently applied, regardless of the type of electrodes used in the signals. Cases that did not meet this minimum requirement were excluded from the analysis.

To mitigate the effect of age, five males and five females were selected from each region of interest, specifically those that minimized differences in mean age and standard deviation between groups. The analysis was limited to the left hemisphere as it met more criteria in this region.

The LFPs from each channel for selected patients were analyzed, associating them with a PDF using the BP method. Subsequently, H_q and C_q were computed for q in the range [0.1, 7] with a step size of 0.01. The results for males and females were separately averaged within each region, and the standard deviation was calculated.

The Table 1 provides a detailed overview of mean ages and standard deviation values for specific brain regions in both females and males within the left hemisphere. Figure 1 schematically illustrates the identified Regions of Interest (ROIs) in the left hemisphere. The red shading indicates the specific region studied for each case. All highlighted areas belong to the left hemisphere. This letter nomenclature for the ROIs is consistent throughout the study and is used in subsequent figures where the results for each region are analyzed. These ROIs were reconstructed using the nodes provided by the MNI Open iEEG Atlas, which follows the Desikan-Killiany parcellation.

An example of temporal iEEG series fragments used to analyze the different brain states studied is shown in Figure 2.

An embedding dimension D = 3, 4, 5, and 6 for BP was chosen to investigate how the choice of embedding dimension affects the results. This selection satisfies the condition M≫D!, where the length of the atlas iEEG data was M = 13600.

In this study, the delay time of τ = 1, as recommended by Bandt and Pompe in their original article (Bandt and Pompe, 2002), was used to analyze causal relationships on a small time scale for theoretical scenarios.

For the analysis of iEEG, two time delays were implemented for each mode in the construction of the probability functions. On one hand, τ = 1 was implemented, aligning with the traditional BP methodology. On the other hand, previous studies have demonstrated a strong relationship between the parameter τ and the intrinsic temporal scales of the analyzed system (Soriano et al., 2011; Zunino et al., 2012, 2022). Hence, the embedding time (τ) was set equal to the shorter characteristic time (τ_s) to account for the short intrinsic temporal scales of each mode. Although investigating different values of τ might have provided additional insights, the primary goal of this study was to match this parameter with the shorter characteristic time (τ_s) observed in the analyzed iEEG and to maintain the traditional delay of τ = 1.

A reliable determination of a system's time delay can be achieved through permutation entropy and statistical complexity. Importantly, these metrics reach extremes when the embedding delay τ aligns with the characteristic delay τ_s of the system. The detection of this parameter is more sensitive using statistical complexity (Zunino et al., 2010a,b). The statistical complexity (MPR; Lamberti et al., 2004; Martin et al., 2006) was computed for the parameter τ, ranging from 1 to 30, with dimensions D set to 3, 4, 5, and 6, across all channels of the atlas in all states. The analysis showed that the highest value is achieved when τ = 1. Statistical complexity may be computed using:

The normalized Shannon entropy H is defined as:

where S_max = S[P_e] = log₂N, and

where N is the number of possible states of the physical system under consideration, and p_j is the probability of each state. The disequilibrium, QJ=Q0J, is determined by the Jensen-Shannon divergence, given by:

with a normalization constant Q₀. Here, P_e represents the uniform distribution (Pe=1N). To obtain the statistical complexity, QJ and H are multiplied together.

Moreover, a widely utilized conventional tool for measuring the characteristic time of a signal is autocorrelation (Lepri et al., 1994; Schäfer and Kratky, 2008; Zunino et al., 2010b). Autocorrelation involves the cross-correlation of a signal with itself, facilitating the identification of repetitive patterns and the detection of signal periodicity, particularly in the presence of noise. Let χ(t) denote the signal, and the continuous autocorrelation R_χχ(τ) at a delay τ is defined as (Zoughi and Boostani, 2010):

Where χ¯ represents the complex conjugate and * denotes convolution. If the function is real, χ¯(t)=χ(t). For a wide-sense stationary process, it is defined as (Zoughi and Boostani, 2010):

iEEG signals are not inherently stationary; however, by considering epochs, stationarity can be achieved (Zoughi and Boostani, 2010). In this study, it is assumed that 68-s samples can be considered stationary epochs.

To test the stationarity of the dataset, an Augmented Dickey Fuller (ADF) test (Mushtaq, 2011; Avramidis et al., 2021) was conducted on all continuous signals within the dataset, using a significance level of 0.05. All signals showed evidence of stationarity. In the W mode, no signals contained NaN values, so all signals were considered stationary when subjected to the test.

Furthermore, to determine τ, the autocorrelation of each signal was calculated, and the distances between the peaks were identified. Subsequently, these distances were averaged to derive a value of τ for each signal. To ascertain τ_s for each mode, the τ values of all channels were averaged. Figure 3A schematically illustrates this process for an arbitrary signal from the set. In this figure, the autocorrelation is visually represented in blue, and the peaks of local maxima are highlighted with red dots.

When finding the characteristic time (τ_s) of the system by analyzing the peaks of the average statistical complexity of all signals in each state, only the first maximum at τ = 1 is found. Although this method allowed the calculation of τ_s for single signals and sets of several hundred channels, averaging the complexities of all signals attenuates the peaks, making detection impossible.

To overcome this limitation, autocorrelation is used. This is a technique that measures the similarity between a signal and a delayed version of itself. Analysis of the autocorrelation function helps identify repeating patterns or cycles within the signal. Autocorrelation is applied individually to each signal to find peaks in the autocorrelation function that correspond to the periodicity or characteristic times of the signal. Calculating the average distance between all of these peaks helps to estimate the characteristic short period for each channel (MathWorks, 2024). Finally, the characteristic times are averaged across all channels for each of the four modes. When using these short periods with the Bandt and Pompe (BP) method, smaller-scale causal relationships are taken into account. These short times were preferred in part to preserve this small-scale causality when sampling the signal to create the probability function of ordinal patterns, given the limited length of the signals that ensure stationary epochs to work with.

An example of a REM mode signal is shown in Figure 3B. This is a zoom-in of a section where the signal exhibits a smooth and regular periodicity, like a sinusoidal wave. Its characteristic time surpasses that considered between the sawtooth teeth seen when all local maxima are considered. Nevertheless, this order is on the scale of τ = 800, which, when multiplied by the sampling frequency, yields ~4 s. This constitutes a considerably large separation for constructing the distribution of ordinal patterns when analyzing 60-s signals. Therefore, shorter times were chosen. Analyzing correlations between widely separated regions relative to signal length is beyond the scope of this study. Using τ significantly larger than D could result in the loss of significant causal information on small time scales.

This yielded values of τsW=10.8, τsN2=12.6, τsN3=16.6, and τsR=10.25, for Wakefulness with Closed Eyes, non-REM Sleep Stage N2, non-REM Sleep Stage N3, and REM Sleep modes, respectively. Figure 4 presents the MPR complexity curves as functions of τ across dimensions D = 3, 4, 5, and 6 for each mode. The corresponding characteristic times τ_s, computed via autocorrelation, are also indicated.

Regarding the Rényi parameter, in this study, quantifiers H_q and C_q were calculated for q in the range [0.1, 7] with a step size of 0.01. As a specific case to visually illustrate the differences in mean values of these quantifiers across regions, the value of q that maximizes Rényi complexity, calculated analogously to τ, was chosen. For D = 6, this resulted in qmaxW=1.4, qmaxN2=0.93, qmaxN3=0.81, and qmaxR=1.06.

The need to introduce the concept of Rényi entropy-complexity causality space (H_q × C_q × q) arises from the necessity to distinguish between different curves of Rényi complexity-entropy. These curves are formed by points derived from different parameters q, within the same distribution. Sometimes these curves overlap when projected onto a single plane. This is true even though their values differ when compared at certain q values.

When computing generalized quantifiers such as the pair H_q and C_q, it is necessary to generalize the concept of the entropy-complexity causality plane (H × C) proposed by Rosso et al. (2007) in order not to lose information in the graphical representation. By introducing the dimension corresponding to the parameter q, it is possible to construct a space, which is analogously called the Rényi Entropy-Complexity Causality Space (H_q × C_q × q), formed by the superposition of multiple causality planes for successive values of q.

The bounds of maximum and minimum complexity can be calculated using the methodology of Martin et al. (2006). These bounds indicate the theoretical upper and lower bounds of complexity and serve as a benchmark for placing different dynamical systems in each q plane within the space. They depend on the choice of parameters of the BP D and the Rényi parameter q.

In this work, the expressions corresponding to Rényi entropy and complexity were modified for the calculation of maximum and minimum complexity bounds in the Ordpy library (Pessa and Ribeiro, 2021)). This library enables the implementation of the BP methodology for Rényi entropy and complexity. Additionally, readers can find the modified bounds functions for Rényi complexity as supplementary material to accompany this Python library.¹

This space is constructed by stacking multiple causal planes (H_q × C_q) within the interval q∈{0.5, 0.7, 1, 1.4, 2, 3, 4, 5, 6, 7}, following the approach by Martin et al. (2006). The Logistic Map (r = 4), the Hénon Map (a = 1.4 and b = 0.3), Schuster Map (parameter z∈{3/2, 2, 5/2}), white noise, and correlated noise for k were represented. Below are the equations of the maps used to generate the temporal series (Pessa and Ribeiro, 2021).

Logistic Map:

Hénon Map:

Schuster Map:

The f^−k behavior of each region under study was explored in the frequency domain by calculating the Power Spectral Density (PSD).

For the PSD estimation, Welch's method was employed on iEEG signals. This involved computing the magnitude of the discrete time Fourier transform for 59 overlapping blocks, each lasting 2s with a 1s step, and weighted by a Hamming window. Averaging these transformed blocks yielded the spectral density for each channel. To ensure independence from the amplitude of the iEEG signal, the spectral density in each channel was subsequently normalized to total power. These parameter choices were guided by prior work from the developers of the MNI Open iEEG Atlas, as referenced in Frauscher et al. (2018a,b); von Ellenrieder et al. (2020).

The power law coefficient k was determined by fitting an exponential function to the data's PSD within the 10–40 Hz range. The lower limit of 10 Hz was selected instead of 1 Hz due to a shoulder between 1 and 10 Hz observed in most channels, as indicated by the atlas authors (Frauscher et al., 2018a). The upper limit of 40 Hz is set by the data's sampling rate. The highest frequency unaffected by anti-aliasing downsampling filters is 80 Hz, but accurate computation of the scale-free spectrum is limited to half that frequency (Frauscher et al., 2018a). Each data channel's PSD was adjusted using the least squares method with the formula A·f^−k, where A represents the amplitude and k is the exponent. Subsequently, the parameters from each channel in all regions were averaged individually for males and females. The mean fits and their corresponding coefficient of determination (R-square, denoted as R below) were plotted on the mean PSD. To compare, a boxplot was generated for males and females, showing the median, interquartile range (IQR), outliers, and a notch representing the 95% confidence interval of the median power law coefficients (k) from the fits.

In Figure 5, four panels show embedding dimensions D = 3, 4, 5, and 6, with a fixed time delay of τ = 1. At the top section of each panel, classical dynamical systems are shown to illustrate their positions in the (H_q × C_q) causality planes corresponding to q∈0.7, 1, 1.4, 4, presented in reading order. These systems are bounded by the maximum (Cqmax) and minimum complexity (Cqmin) lines, shown as solid black lines. Logistic Map (r = 4), the Hénon Map (a = 1.4 and b = 0.3), the Schuster Map (parameter z∈3/2, 2, 5/2), white noise, and correlated noise for k∈[1, 5], with increasing values of k from right to left, were represented.

Below each panel, the (H_q × C_q × q) space is shown schematically. The maximum complexity (Cqmax) is shown in red, and the minimum complexity (Cqmin) is shown in blue. At q = 1, the classical Jensen-Shannon complexity-entropy causal plane (Rosso et al., 2007) is evident. For q > 1, Cqmax increases and Cqmin decreases, expanding the range of possible values. Conversely, for q < 1, these bounds approach each other, reducing the area between the curves. The space (H_q×C_q×q) is defined by the volume formed between these curves as the number of q values considered increases.

The position of the systems within the space defined between the bounds depends on the embedding dimension D used. To distinguish system features, one should choose an appropriate τ, optimized according to some criterion (e.g., maximizing complexity), and implement the D that maximizes the differences between the dynamical systems under study. Chaotic systems such as the Hénon map and the logistic map separate at low dimensions and converge as D increases. The same trend is observed for the three implemented cases of the Schuster map. In the case of correlated noise, this behavior depends on the values of k studied.

In Figure 6, computationally simulated correlated noise (f^−k with k ranging from 1 to 5 spaced logarithmically) is analyzed.

Rényi entropy-complexity causality (H_q and C_q) calculations have been performed using embedding dimensions of D = 3, 4, 5, and 6, with a fixed time delay of τ = 1. This is presented in four panels, one for each dimension. Time series were generated with a length of M = 10000, and the Rényi parameter varies within the range q∈[0.1, 7].

For each dimension, three different views of the Rényi Entropy-Complexity Causality Space (H_q × C_q × q) are presented, with the projection onto the H_q × C_q plane on the right. A graduated color scale is used to represent different values of the exponent k.

In the three-dimensional representations, the entropy-complexity curves of the different series do not intersect, allowing the dynamics of the time series with different noise values of k to be distinguished in space. In the case of projections onto the plane (H_q × C_q), it is more difficult to separate the curves because they touch at points where the values of (H_q, C_q) are the same, even though they come from different q parameters.

Care must be taken when choosing the embedding dimension D, as the behavior in space appears to be opposite to that of the plane, which is obtained by projecting three-dimensional curves onto it. This can be seen especially for the intermediate values analyzed, where the curves open up in space as D decreases, but when projected onto the plane, the overlap is greater.

This section presents the results of the H_q × C_q × q curves for both sexes, considering the standard deviation and the exponents of the fits made to the PSD considering the notches. Emphasis is placed on the differences found in the H_q × C_q × q curves, relating them to the differences found in the k exponent of the power law (f^−k).

Figure 7 provides an illustrative comparison of curves within the causal space (H_q × C_q × q) for both sexes, showcasing the ROI of region B for the R mode (REM sleep). The curves within the interval q ∈ (0, 7) do not show complete separation; however, for values greater than 1, separation occurs when considering one standard deviation, while they begin to converge as values approach 0. In the upper panel, Figure 7A displays the curves for males (in cyan) and females (in violet) across the full range, with a zoom to three q intervals: 0.1–1, 1–2, and 3–7 for D = 6. A schematic of the theoretical bounds is included in the center, illustrating that for values < 1, the space rapidly contracts, bringing all systems closer together. The upper corner presents the plane (H_q × C_q), where the curves almost entirely overlap across the plane when considering deviations. Below, in Panels B–D, the same curves are shown in the space (H_q × C_q × q) for D = 5, D = 4, and D = 3, respectively. In all cases, a time delay of τ = 1 was used. It is observed that similar behaviors are found in other dimensions.

Figure 8 present, from top to bottom, the mean Power Spectral Density (PSD) within the frequency range of 10 − 40Hz, along with the mean fitting and R values for males and females. Additionally, boxplots depicting the power law fitting exponents are shown, as well as the curves within the Rényi Entropy Complexity Space for D = 6 with delay times of τ = 1 and τ = τ_s.

Results for females are depicted in violet and for males in light blue. The top panel corresponds to the Wakefulness with Closed Eyes state, and the bottom panel to the non-REM Sleep Stage N2. Each column represents an analyzed Region of Interest (ROI) following the convention of Figure 1.

In the Appendix Figure A1 display, for interested readers, the medians of PSD for males and females are shown, along with the Interquartile Range (IQR) within the frequency range of 0.5 − 30Hz, corresponding to the δ, θ, α, and β bands.

Figure 9 are analogous to the aforementioned, but the panels correspond to non-REM Sleep Stage N3, and the bottom one to REM Sleep.

Table 2 provides a summary of the sex differences found across all the analyzed regions. The rows show the k exponent from the f^−k fit over the PSD, analyzing the notches of the boxplots and the cases of τ = 1 and τ = τ_s, marking the standard deviations of H_q and C_q as shaded ellipses. A “1” indicates that differences were found, while a “0” indicates that no differences were found. Strong sky blue shading highlights regions where differences in exponent were found for both time delays used; lighter shading indicates regions where only one of the τ showed curve separation.

The cases where differences were found in the exponents and curves simultaneously are detailed below. In the W state, differences were found in τ = 1 in B. In N2, they were found simultaneously in τ = 1 and τ = τ_s in B, in τ = 1 in I and J, and τ = τ_s in E. In N3, they were found simultaneously in τ = 1 and τ = τ_s also in F and I and in τ = 1 in B and E. In R, the greatest number of differences were found, simultaneously in τ = 1 and τ = τ_s in A, J, and L and in τ = 1 in B, and τ = τ_s in E.

The ROI that shows the most differences in the curves at the same time as in the exponent is B. On the other hand, in the τ curves, differences were found simultaneously with differences in the k exponent with a higher frequency with respect to τ_s, the frequency of occurrence of this event being 12 and 8, respectively.

In Figure 10, four panels present a visual scheme illustrating the regions studied using D = 6 and q = q_max across the four states investigated: Wakefulness with Closed Eyes, non-REM Sleep Stage N2, non-REM Sleep Stage N3, and REM Sleep from top to bottom. Each panel displays four different views of a brain diagram, color-coded to represent the findings in the regions of interest (ROI) as delineated in Figure 1. The colors indicate the absolute difference in average Rényi entropy and complexity between men and women, ranging from the minimum difference to the maximum (indicated by the scale to the right of each panel). In each instance, the regions are labeled in the order of decreasing differences in the reading direction. For the Wakefulness state, the regions with the most significant differences between men and women in Rényi entropy are B, L, and I, while in complexity, they are J, A, and D. For the N2 and N3 stages of non-REM sleep, the regions with the greatest disparities in entropy and complexity for q_max in descending order are H, E, and G for both metrics. In the case of REM Sleep, the regions with the most notable differences in entropy are G, B, and E, while for complexity, they are B, E, and G.