Sleep hd-EEG recordings from 12 FEP, seven BD patients, and 13 control subjects were analyzed, see Methods for details and Figure 1. Raw data are extremely fine-grained: the sampling frequency of 500 Hz during an average of 8 h of sleep, multiplied by the 64 electrodes that comprise the EEG headset yields approximately, 1,000,000,000 measurements per patient. Clearly, these measurements are not all independent of each other, but they encode information that spans several sleep phases and brain regions. Therefore, we aim at finding a satisfactory compromise between compression and information.

To do so, we process the raw sleep EEG records through our pipeline as described in Methods in detail and illustrated in Figure 1. The first step is to remove artifacts from the data. Eye-movement artifacts are well known to influence raw sleep EEG data. To mitigate their impact on our results, we use a fast linear regression model to correct for eye movements, see Methods for details and Figure 1B: in this illustrative figure, the horizontal electro-oculogram potential (HEOG) well correlated with channels AF8, F7, FP2 and FP1 in the top plot. After the correction step (bottom plot in pannel B), this dependence was almost completely eliminated. After splitting the signal into different frequency bands (see Methods for details), we compute time- and band-specific electrode-to-electrode correlations of the form C i j b ( t ) , represented as a heatmap in Figure 1C. Finally, we construct time-varying multilayer networks using an innovative strategy that takes into account the whole dataset (and not each time snapshot individually), maximising the total amount of information contained in the time-varying dataset. Figure 1D offers a visual representation of the final output we obtain after processing the raw EEG data: a set of time-varying multilayer networs, where different layers correspond to different frequency bands, network nodes represent electrodes and edges represent high EEG correlations.

Networks offer a simplified and effective representation of interactions between nodes, but deciding the correlation threshold beyond which edges are added to the network is a nontrivial subject. In order to make an informed choice, here we introduce the Integrated Jensen-Shannon Divergence (IJSD), I ( θ ) = ∑ t = 1 θ D ( ρ t − 1 , ρ t ) , (1) a measure of the total information change over time (Grosse et al., 2002), computed as the sum of the Jensen-Shannon divergence of each epoch with respect to the previous one. Here ρ _t are the density matrices associated to each network in the framework of spectral entropies (Domenico and Biamonte, 2016), see Methods for details. The value of I depends on θ in non-trivial ways, but the limit cases are clear: if θ is too low (high), all edges are present (absent) at all time steps, so there is no information change over time and thus I = 0 for both θ = 0 and θ = 1. It is only for intermediate values of the correlation threshold θ that the sequence of multilayer networks can display richer temporal variations, yielding a higher information change. This can be clearly seen in Figure 2 panels (A, B, C), which show the value of I as a function of θ for one BD, one FEP and one control example. As anticipated, I ( θ ) = 0 for both θ = 0 and θ = 1, with a clear maximum at around θ ∼ 0.7 for most frequency bands.

We implement two strategies to choose the optimal correlation threshold θ* from the analysis of the information content quantified by IJSD. In the first approach, we set a global absolute value for the correlation threshold, while in the second approach that value is relative to each network and chosen to maintain a constant edge density, keeping only the interactions with highest absolute correlation. In both cases, the adjacency matrices can be build as A i j b ( t ) = 1 if  | C i j b ( t ) | ≥ θ * 0 else (2) that is, we place edges for both large positive and large negative correlations.

The optimal correlation threshold θ* for fixed-threshold networks is computed as the average of the band- and patient-specific optimal values that result from optimizing each case separately, θ * = θ b , p * b , p (3) θ b , p * = a r g m a x θ ∈ [ 0,1 ] I b , p ( θ ) (4) where I θ ( b , p ) denotes the IJSD of patient p at frequency band b. In other words, for each patient p we compute a band-specific optimal threshold θ b , p * . The group averages and variability of these are shown in Figure 2. Taking the average of all θ b , p * , we reach an overall value of θ* = 0.72, shown as a black solid line in Figure 2. Overall, the figure shows that a single global threshold can reasonably accommodate for the band- and patient-specific optimal values.

The second approach consists in keeping the same fraction of edges in all networks, yielding what we call fixed-density networks. The optimal density value in this case is set so that it coincides with the average density of the fixed-threshold networks. This second approach takes into account that different patients, time point or bands might have different intrinsic correlation levels, and presents additional advantages from the network analysis point of view.

We investigate the group differences between BD and control patients, as well as between FEP and control patients. To do so, we need to condense the information contained in our multilayer and time-varying networks into simpler summary statistics. A simple yet useful measure in this case is what we coin as edge presence, which is the fraction of time an edge is present (that is, A i j b ( t ) = 1 ) during one full EEG sleep session. Formally, P i j b = A i j b ( t ) t (5)

Figure 3 shows the group differences of P i j b for each edge (i, j) and each band when comparing BD patients with controls (panels A, C), as well as FEP patients with controls (panels B, D). This analysis is shown both for fixed-density networks (A, B) and for fixed-threshold networks (C, D). In both cases we see differences in the parieto-occipital area, but the signal is stronger for fixed-density networks. If we focus on Figure 3A, for instance, we see that BD patients tend to have a lower edge presence in the parieto-occipital area (strong blue edges). Notice that we employ a colorbar that goes from red to transparent to blue, so that edges that do not have strong differences are effectively not drawn. Overall, the figure shows important differences in the parieto-occipital area, with a similar but stronger signal for fixed-density networks.

Motivated by the results shown visually in Figure 3, we construct a parieto-occipital (PO) specific measure. Selecting the 18 electrodes of that region (see methods for details), we compute the difference of PO presence between the PO area and the rest. P PO = P i j b ( i j ) ∈ P O − P i j b ( i j ) ∉ P O (6)

Additionally, we also consider the average clustering coefficient, the average shortest path length and the betweenness centrality as measures related to clustering and information navigability as candidates to better quantify the differences that we see in Figure 3.

Figure 4 shows boxplots of these four measures comparing, BD and FEP patients with control subjects. Statistically significant differences are marked with a star, see Methods for details. Panels (A, B) show that for bands 1 to 4 (that is, between 1 and 16 Hz), FEP patients have a higher clustering coefficient when compared to controls, while panels (C, D) shows some significant results on the same range of frequencies for the average shortest path length, both for BD and for FEP patients. Turning to betweenness centrality, panel E shows that when using fixed-threshold networks, FEP patients significantly differ from controls in bands 2, 3, and 4 (2–16 Hz) Interestingly, when looking at the parieto-occipital relative presence (panels G, H), we observe a different pattern of marked differences between BD and control patients for lower frequency bands, 0.5–4 Hz. This is consistent with the fact that the more standard network measures used in panels A to F treat all nodes under the same footing, independently of the brain region they correspond to, while PO presence is a tailor-made measure, specifically designed to capture the visual results of Figure 3 taking into account the location of parieto-occipital electrodes.

It is interesting to ask if the network measures shown in Figure 4 are correlated with network edge density, for the case of fixed-threshold networks. Figure 5 shows how indeed edge density is a strong driver of average clustering coefficient, average shortest path length and betweenness centrality for all patient groups, but not of parieto-occipital presence. This is consistent with the fact that, by construction, P _PO is a relative difference of two averages taken on the same network.

Hd-EEG recordings where obtained from San Paolo Hospital in Milano. In particular, the dataset consists of sleep EEG recordings from 12 FEP patients (Eight males and four females, mean age 21.0 ± 3.77), seven BD patients (Three males and four females, mean age 34.57 ± 7.09), and 13 control subjects (Six males and seven females, mean age 25.61 ± 10.64). All participants underwent an in-laboratory sleep hd-EEG recording with a 64-electrode Easycap net designed to enhance electrode contact with the scalp (BrainAmp, Brain Products GmbH, Gilching, Germany). The night of the recording, all subjects were accommodated in a sleep suite and allowed to sleep within 1 h of their usual bedtime. All subjects were recorded throughout the night and until they woke up naturally the next morning. Table 1 shows the average length of recording sessions and total sleep time for each participant group. The headset has 64 unipolar electrodes positioned following the standard 10–20 system, and include two channels that record eye movements (one for vertical movements and one for horizontal movements). All recordings had a sampling frequency of 500 Hz. Data was provided in anonymized form as pairs of .set and .fdt files.

Our preprocessing pipeline transforms the raw EEG recordings into correlation tensors of the form C i j b ( t ) , with (i, j) denoting and edge between electrodes i and j, b a specific frequency band, and t a 30-s epoch. The steps we carry are as follows:

1. Epochs division: divide the raw signal into epochs of approximately 30 s, see below for details, obtaining a signal S _i (t) for electrode i and epoch t.

Epochs division: We divide EEG recordings into epochs of around 30 s following Aboalayon et al. (2016). To be precise, each epoch has a length of 2¹⁴ raw time points which, at a sampling frequency of 500 Hz, corresponds to 32.768 s. This choice is particularly convenient because pure powers of two allow for faster discrete Fourier transform calculations.

Artifact correction: Following Gratton et al. (1983), we correct for eye-movements using a linear regression equation of the form Y = X B (7) where Y corresponds to the EEG data (62 channels in our case), X corresponds to the eye-movement data (Two channels in our case), and B is the regression coefficient matrix to be determined. Solving for B via least squares, the corrected signal X* is computed as X * = ( X − Y B ) T (8)

Bands division: We use seven frequency bands, numbered from 0 to 6 throughout the manuscript, which logarithmically interpolate the 0.5–64 Hz range typical of brain waves.

• Band 0: (0.5, 1) Hz.

We use the Pearson correlation coefficient to measure the strength and direction of dependence between the signals x _i , x _j recorded by two electrodes i, j, C i j = c o v ( x i , x j ) σ x i σ x j . (9)

Repeating this measurement for each band b and timepoint t, we get a full correlation tensor C i j b ( t ) .

We use Jensen-Shannon Divergence (JSD) as a distance measure between networks, in the framework of spectral entropies (Domenico and Biamonte, 2016). For a pair of networks with density matrices ρ and σ, the JSD is defined as J ρ ‖ σ = S ρ + σ 2 − 1 2 S ( ρ ) + S ( σ ) , (10) where S(ρ) is the spectral entropy of the network, S ( ρ ) = log 2 Z + τ ln ⁡ 2 Tr [ L ρ ] , (11) with L denoting the Laplacian, τ diffusion time and the density matrix ρ defined as ρ = e − τ L Z , Z = Tr ( e − τ L ) (12)

Parieto-occipital edge presence: The parieto-occipital area is mapped to the following electrodes: P7, P5, P3, P1, PZ, P2, P4, P6, P8, PO7, PO3, PO4, PO8, O1, OZ, O2, and IZ. From this list, the parieto-occipital presence is computed as explained in the main text, mainly the difference of average presence between nodes in the parieto-occipital area and the rest.

Clustering Coefficient: We use the standard definition of clustering coefficient, c i = 2 ⋅ t i k i ⋅ ( k i − 1 ) (13) as implemented in the networkx python library (Hagberg et al., 2008), where t _i is the number of triangles in which node i is involved and k _i is the node degree. Averaging over all nodes, we define the clustering coefficient of the network as c = 1 N ∑ i = 1 N c i (14)

Betweenness Centrality: We use the convention of Brandes (2008), which defines a node-dependent quantity as follows: c B ( i ) = 2 ( N − 1 ) ( N − 2 ) ∑ j , k ∈ V σ ( j , k | i ) σ ( j , k ) (15) where σ( j, k|i) is the number of shortest path that connect nodes j and k that passes through i and σ( j, k|i) = 0 if i = j, k. σ( j, k) is the total number of shortest path connecting j and k and σ( j, k) = 1 if j = k. By convention the fraction σ ( j , k | i ) σ ( j , k ) is considered zero if both elements are zero. We then average over all nodes to get a single measure for each network: B C = 1 N ∑ i = 1 N c B ( i ) . (16)

Average Shortest Path Length: We start from the standard definition of average shortest path length (ASPL) for a connected graph G a G = 1 N ⋅ ( N − 1 ) ∑ i ≠ j d ( i , j ) . (17) where d(i, j) is defined as the length of the shortest path connecting two nodes, j. If i and j belong to two different connected components d(i, j) is said to be infinite, while d(i, j) = 0 if i = j.

In our setting, networks can have more than one connected component, and we do not want to limit ourselves to the largest connected component as important information could be missed. Hence we employ a weighted version of the ASPL, w a G = ∑ c = 1 n c a c ⋅ w c ∑ c = 1 n c w c (18) where n _c is the number of connected components with more than two nodes and w _c = N _c ⋅ (N _c − 1), N _c is the number of nodes of component c. This formulation takes into account the ASPL of all nodes but effectively gives more weight to the larger components.

Group differences are assessed with a two-sided T-test without assuming equal variances between groups, as implemented in the ttest_ ind function from the scipy Python scientific library. Cases marked as significant (∗) in Figure 4 correspond to a p-value below 0.05.

Data from the SPINDLE-1 study, approved by the Milan Area A Interhospital Ethics Committee (Approval n. 22864). All participants signed an informed consent for participation in the SPINDLE-1 study.