Statistical Modelling of Cortical Connectivity Using Non-invasive Electroencephalograms

The following experimental protocol is in accordance with all local, national, and international ethics guidelines for human research. The data used to test the protocol have been acquired with authorization of the Ethical Committee of region Tuscany-protocol 2018SMIA112 SI-RE.

NOTE: The scripts used for implementing the analyses described are available at https://github.com/conorkeogh/NetworkAnalysis.

1. Raw Data Collection

1. Prepare subject conditions.
   1. To ensure consistency across recordings, conduct all EEG recordings in a dedicated recording environment. Remove all equipment or stimuli that is not directly relevant to the task to be performed during recording from the environment to avoid distraction.
      NOTE: If resting state recordings are to be performed, remove all sources of distraction from the room and expose subjects to the recording environment prior to the recording session to remove novelty from the environment.
   2. Provide the subject with clear instructions regarding the task to be performed. Once the equipment has been set up, leave the subject alone in the recording environment to get used to the environment prior to beginning recording to minimize movement and distraction.
   3. If the subject has intellectual disabilities, allow him/her the necessary time to accustom to the environment to limit any stress. Sometimes this may require multiple visits and an extended stay in the recording room.
2. Mount the electrodes.
   1. Attach the electrode cap to the patient's head, taking care to ensure correct alignment. Inject conductive gel into each of the electrode ports, beginning at the scalp and slowly withdrawing to the cap surface to establish electrical contact with the scale and improve the signal-to-noise ratio.
   2. Attach electrodes to the electrode cap using a predetermined electrode montage based on the 10−20 system. Attach appropriate ground electrodes (e.g., to the mastoid processes).
3. Set up the EEG.
   1. Connect all electrodes to an electrophysiological recording system. Link the recording system with an appropriate digital recording environment.
   2. Examine all the recording channels to ensure the offset is within an appropriate range and to avoid excessive channel noise. If a channel has an excessive offset or noise, additional conductive gel can be added in order to improve the electrical connection, taking care to avoid causing bridging between electrode sites.
   3. Instruct the subject that recording has started and to avoid all unnecessary movements. Conduct a short test recording to verify appropriate recording quality.
4. Prepare the behavioral task for recording.
   1. Clarify all the task-related instructions with the subject. Reiterate the importance of avoiding all unnecessary movements.
   2. Explain that recording will begin on a clearly agreed signal (e.g., a knock on the recording environment door). Leave the subject in the recording environment. Start recording. Give the agreed signal to the subject.
   3. Following the completion of the task or period of resting state, stop recording, visually examine the data to ensure quality, and save the data.

2. Data Preprocessing

NOTE: The data preparation and feature extraction pipeline is illustrated in Figure 1.

1. Prepare the software.
   1. Load the EEG data to be analyzed into a data analysis environment. Load any additional script libraries necessary, such as EEGLab²².
2. Convert all recordings to the same data format if necessary, with all channels in their corresponding locations.
   1. Discard the beginning and end of each recording (e.g., 5 min) to reduce the contamination of movement artifacts. Split the data into epochs based on task or, if it is a resting state recording, predetermined duration (e.g., 10 min). See NetworkAnalysis_Demonstration.m (section Feature Extraction) and Supplementary Figure 1 for a demonstration of implementation.
      NOTE: Selection of the epoch length can have important effects on the measures of coherence. Epochs of a sufficient length should be used to ensure that true relationships between the signals emerge in the computations to avoid unnoticed artifacts or transient, spurious synchronizations having an excessive weighting. However, in this work there was no statistically significant difference in the overall network structure when ten-minute epochs were compared to an average of ten one-minute epochs following thorough artifact rejection.
3. Perform artifact rejection by visually inspecting the epoch data and rejecting visually unsuitable data.
   NOTE: As the modelling technique described relies on the relationships between signals, it is essential to ensure thorough rejection of artifacts. These can corrupt the channel data, leading to artificial increases (if the artifact is represented on multiple channels) or decreases (if the artifact is represented only on some channels) of measures of coherence.
   1. Identify bad channels in the recordings.
      1. High pass filter data at 0.5 Hz to remove the baseline drift due to the floating ground of the acquisition system.
      2. Select all channels meeting the appropriate statistical criteria (e.g., those with a standard deviation greater than three times or less than one third of the average channel standard deviation).
         NOTE: Removing the channels with data that are unlikely to have originated from neural sources avoids spurious relationships being introduced into the network models.
      3. Examine these channels to determine whether they are suitable.
      4. Reject the epochs with unsuitable channels if possible. Alternatively, exclude the bad channels and interpolate the data at these channels (e.g., using EEGLab's spline interpolation algorithm).
         NOTE: Interpolation across a large number of channels or with only a small number of recording channels may generate unsuitable data for analysis. Further, this introduces no new information into the dataset and can result in artificially high measures of coherence between interpolated signals and the signals from which they are derived.
   2. Perform the independent component analysis on the remaining epochs (e.g., using EEGLab's ICA function). Visually inspect the derived components and reject visually unsuitable data.
   3. Apply the appropriate statistical thresholds to identify the potential artifacts not immediately evident on visual inspection (e.g., based on extreme values or abnormal spectra). Examine these and determine whether rejection is appropriate.
   4. Repeat the independent component analysis and artifact identification on the surviving epochs.
   5. Identify the data epochs to be saved for further analysis. Discard all the rejected data epochs. Identify all the epochs to be taken forward for further analysis.
      NOTE: Where only one epoch per subject is required, select the first suitable epoch for further analysis.
4. To prepare the data, correct the baseline of the recordings by subtracting the mean of all channels from the recordings to avoid the impact of baseline wandering during prolonged recordings. Re-reference all the channels to an appropriate reference (e.g., the ground electrode or the average of all channels). See NetworkAnalysis_Demonstration.m, NetworkAnalysis_Preprocess.m, and Supplementary Figure 2 for examples of implementation.
   NOTE: Reference selection can have important effects on the network measures. As the reference data are "subtracted out" of all analyzed channels, any neural data that are represented on the reference channel will be subtracted out and thus not contribute to the model generation. It is common practice to use reference signals recorded over bony prominences without immediately underlying neural structures, such as the mastoid process. However, these can be corrupted by neural data due to volume conduction effects through the scalp and therefore distort network measures differentially based on location relative to the reference. As a result, for resting state data it is best to use an average of all scalp channels as reference. This means that all data are not referenced relative to a specific spatial location, distorting measures, because all channels contribute to the reference. This can have effects such as dampening the apparent overall activity and can distort measures by subtracting out signals that are very strongly represented on some channels and thus contribute heavily to the average. This is a greater issue for activity- and event-related signals, but is typically not the case with resting state data.
   1. Digitally filter all the channels to isolate frequencies of interest (e.g., 1 Hz-50 Hz). See NetworkAnalysis_Demonstration.m, NetworkAnalysis_Preprocess.m, and Supplementary Figure 3 for examples of implementation.
      NOTE: Ensure the use of appropriate frequency limits and filter parameters for the intended analysis to avoid the distortion of frequencies at the extremes of the examined range and aliasing effects. Zero phase-shift 4^th-order Butterworth filters perform appropriately. Appropriate filtering ensures that the activity of interest is isolated for modelling. Even with a broad range (e.g., 1 Hz-50 Hz), this ensures that high frequency artifacts and low frequency baseline wandering are not interpreted as coherent between channels, distorting measures.

3. Feature Extraction

1. Assess spectral power.
   1. Calculate overall power spectra by performing a Fourier transform of each channel being analyzed across the whole frequency range to be assessed (e.g., 1 Hz-50 Hz).
   2. Assess the activity in individual frequency bands: isolate the theta band at 4 Hz-8 Hz. Isolate the alpha band at 8 Hz-12 Hz. Isolate the beta band at 12 Hz-30 Hz. Isolate the delta band at 0.5 Hz-4Hz. Isolate the gamma band at >30 Hz (e.g., 30-50 Hz). See NetworkAnalysis_Demonstration.m, NetworkAnalysis_FeatureExtraction.m, and Supplementary Figure 4 for examples of the implementation of spectra derivation and isolation of frequency bands.
      NOTE: EEG data are traditionally divided into frequency "bands" for investigation. These are primarily named based on the order in which they were discovered, and the specific bandwidths vary somewhat. The functional significance of oscillations at specific frequencies remains an area of active investigation. It is thought that oscillations within specific bands may be related to specific neural activities, such as the emergence of a high-amplitude alpha wave in the occipital region with the eyes closed, although the exact relationship between neural functions and oscillatory activity in EEG recordings remains unclear.
   3. Evaluate overall power across the whole scalp by calculating the mean of individual channel spectra. Normalize the power in individual bands with respect to the overall power to give a measure of relative power and allow more accurate comparisons between conditions.
2. Perform network mapping.
   1. Evaluate the interactions between the first electrode pair by deriving a measure of interelectrode coherence:
      
      See NetworkAnalysis_Demonstration.m, NetworkAnalysis_FeatureExtraction.m, and Supplementary Figure 5 for examples of implementation.
      1. Calculate the cross-spectrum of the two channels:
         
         1. Compute the Fourier transform of each signal, X and Y
            
         2. Calculate the cross-spectrum:
            
            Where: t is the sampling interval, T is the length of recording, X is the Fourier transform of x, and Y^* is the complex conjugate of Y.
         3. Ignore the negative frequencies and correct measures. The second half of the computer frequency axis can be disregarded in the case of real-valued signals, and the power measures multiplied by two to correct for this.
            NOTE: This is equivalent to the Fourier transform of the cross-correlation of x and y.
      2. Normalize the cross-spectrum by the power spectra of both channels: .
         1. Compute the Fourier transform of each signal:
            
         2. Calculate the power spectrum:
            
            Where: t is the sampling interval, T is the length of recording, X is the Fourier transform of x, and X^* is the complex conjugate of X.
         3. Ignore the negative frequencies and correct measures: the second half of the computer frequency axis can be disregarded in the case of real-valued signals, and the power measures multiplied by two to correct for this.
         4. Use the calculated power spectra to normalize the cross-spectrum and derive a measure of coherence:
            
            NOTE: This generates C, a measure of the coherence between signals x and y at the frequencies . This is a measure of the relationship between these signals at the frequencies examined, measured on a scale from 0 to 1. Where there is a constant phase relationship between the two signals examined at all timepoints, the coherence will have a value of 1, indicating a strong relationship between signals at those frequencies, implying that activity in one signal is functionally related to activity in the other (i.e., that there is communication between the two). Where there is no phase relationship between the two signals, the coherence will have a value of 0, indicating that the signals are not related.
   2. Repeat this procedure for each unique pair of electrodes to develop a measure of phase stability between the signals at each electrode pair, building up a model of functional connectivity across all electrodes.
      NOTE: For a montage of n electrodes, this will produce coherence measures. This represents mapping the measured time series data onto a high-dimensional plane based on the relationships between recorded signals, allowing the nature of these interactions to be investigated.

4. Data Visualization

1. Perform spectral power analysis.
   1. Examine the power matrices.
      1. Map the measurements of the spectral power to be visualized onto a two-dimensional data structure where each column is an electrode location, each row is a frequency band, and each cell is the spectral power at that location, within that band.
      2. Identify the maximum and minimum power levels across all conditions to be compared. Set these at the maximum and minimum for all conditions. Map the spectral power values between the identified maximum and minimum to colors. Export a color map visualizing the spectral power at each frequency band at each electrode location (Figure 2).
   2. Perform topographic mapping.
      1. Create a data structure containing the labels of each of the 10-20 system electrode locations used, in order corresponding to that of the data structure to be mapped. Using EEGLab's topoplot() function, the spectral power data, the identified maximum and minimum, and the channel list, generate a plot mapping the distribution of spectral power across the scalp.
2. Assess coherence.
   1. Examine the coherence matrices.
      1. Map the measurements of the interelectrode coherence to be visualized onto a two-dimensional data structure where each column is an electrode location, each row is an electrode location, and each cell is the coherence between the corresponding electrode pair.
      2. Map the coherence values between 0 and 1 to colors. Export a color map visualizing the interelectrode coherence between each electrode pair within the frequency limits used (Figure 3). Repeat this procedure for each frequency band to be investigated. See Supplementary Figure 6 and produce_plots.r for examples of implementation. See Figure 3 for example output.
   2. Perform network visualization.
      1. To visualize higher-order interactions between cortical areas and map out network dynamics, calculate how each electrode pair's coherence measure covaries with those of every other unique electrode pair across the overall spectrum and within specific bands.
      2. Map these covariance measures to colors. Export a color map visualizing the network dynamics within and across-frequency bands. See produce_plots.r for examples of implementation. See Figure 4 for example output.

5. Analyzing Network Models

NOTE: The application of modern statistical methods to the models derived allows for taking advantage of the relationships modelled in the high-dimensional network feature space to investigate cortical function. A number of approaches that offer advantages over traditional comparisons of the individual measures or averages of the coherence measures can be taken. Some of the potential approaches these network models facilitate are outlined below. These are discussed only superficially as indicative of the potential applications of network modelling, because a thorough discussion of each technique is beyond the scope of the present work.

1. Perform dimensionality reduction.
   NOTE: Comparisons at the individual variable level fail to take advantage of the relationships represented by the models created, while performing comparisons on all of the measures in the dimensional constructs created is problematic due to the huge number of comparisons required and the failure to integrate the high-level information contained in the statistical models. Mapping the high-dimensional data onto a lower-dimensional space while maintaining the information generated by the model generation process allows for the performance of meaningful comparisons while taking full advantage of the models' data-rich structure.
   1. Derive measures for comparison between the groups that represent the overall network dynamics within the statistical models generated using the principal component analysis. See NetworkAnalysis_Demonstration.m and Supplementary Figure 7 for an example of implementation.
      1. As above, construct a covariance matrix for the pairwise coherence measures. This will generate a dimensional covariance construct where . This model is therefore extremely high-dimensional and allows visualization of high-level network relationships as outlined above.
      2. Decompose the covariance matrix into eigenvectors and corresponding eigenvalues. This allows identification of the axes within the model feature space that contain the greatest variance, without being bounded by the existing measures.
      3. Rank the eigenvectors by corresponding eigenvalue to identify those accounting for the greatest proportion of variance within the model.
   2. Compare the first principal components derived from the network models. See NetworkAnalysis_Demonstration.m and Supplementary Figure 7 for an example of implementation.
      NOTE: The first principal component accounts for the greatest degree of variance within the model. Therefore, comparison of this measure allows for comparison of the overall network dynamics throughout the model as a whole between groups with a single statistical test, allowing for simultaneous analysis of the complex relationships being modelled and avoiding the issues associated with many comparisons.
2. Perform a region of interest analysis. The models derived represent network connectivity across the whole cortex, across all frequency bands. If there is interest in specific anatomical areas or in functions within specific bands, these regions of the model can be isolated and analyzed separately.
   1. Choose an anatomical region of interest.
      NOTE: Restricting analysis to specific anatomical areas allows for evaluation of the network activity within or between specific cortical areas in order to identify the relationships that may not be apparent on analysis of the model as a whole.
      1. Identify the coherence data within the model relating to the anatomical areas of interest.
      2. Derive a covariance matrix and perform principal component analysis as described above to compute measures of overall network architecture within the regions of interest.
      3. Compare the measures of network dynamics within the anatomical regions of interest between the groups as outlined above.
   2. Choose a functional region of interest.
      NOTE: Restricting analysis to specific frequency bands allows for assessment of the network activity within specific oscillatory frequencies (Figure 4).
      1. As with anatomical analyses, isolate the coherence data within the frequency bands of interest. See NetworkAnalysis_Demonstration.m and Supplementary Figure 8 for examples of implementation, using interactions within the overall spectrum only as an example.
      2. Perform principal component analysis to derive measures of overall network activity within the bands of interest.
      3. Compare the measures between groups to evaluate the network differences at specific oscillatory frequencies.
3. Use machine learning.
   NOTE: Modern statistical learning approaches can be applied to the models generated in order to further interrogate the high-level relationships represented within them.
   1. Use supervised learning.
      NOTE: Using data with predefined classes, the models of cortical networks can be used to derive classifiers that can be used to identify signatures within the complex relationships represented by the models to classify new data, opening the possibility for investigating novel diagnostic and prognostic biomarkers, etc. Furthermore, which features within the models drive these classifications in order to gain insights into the underlying mechanisms can be investigated.
      1. Derive the classifiers. Using pre-labelled data, a classifier can be derived to predict the class of a set of data based on the network models.
         1. Divide the data into a set of subject data for training and a set for testing the classifier.
         2. Train a classification algorithm such as a support vector machine or a random forest on the labelled training data.
         3. Assess the performance of the model-trained classifier on the test data.
            NOTE: These approaches allow use of the statistical models as inputs to derive novel biomarkers.
      2. Perform sequential elimination.
         NOTE: Using the model to train a classifier, data can be removed iteratively and the training process can be repeated to identify which components of the model are driving its predictive ability, allowing for investigation of the underlying mechanisms.
         1. Train a classifier on the model as described above.
         2. Remove the model feature with the lowest variability between groups.
         3. Repeat the training process and assess performance.
         4. Repeat the iterative feature removal until the features that contribute most to the performance are identified. These are the model components responsible for the ability to differentiate between classes.
   2. Perform unsupervised learning.
      NOTE: Using the models alone, insight can be gained into the groups being investigated. By modelling the data as high-dimensional constructs based on the relationships between recordings, relationships between groups that were not seen at the level of individual recordings may become apparent. Unsupervised techniques such as clustering algorithms allow for the investigation of relationships within the models without being restricted by the predefined classes.
      1. Using a distance metric such as Euclidean distance, compute the measures of distance between subjects within the space defined by the network model. See NetworkAnalysis_Demonstration.m and Supplementary Figure 9 for an example of implementation.
      2. Using a clustering algorithm such as k-nearest neighbors, identify the groups within the data based on the model parameters (Figure 5).
      3. Repeat this procedure using a sequential elimination procedure as described above to investigate how individual features contribute to the groupings within the model.
         NOTE: This allows use of the derived models to identify the groups within the data that were not apparent otherwise. This may allow for the derivation of disease subtypes, pathological groupings, etc., that are only evident at the network level.