Constituting the most employed metrics, linear modelbased approaches are sufficient for identifying the wide range of oscillatory interactions that take place under the hypothesis of oscillatory phase coupling
. Linearmodelbased approaches allow the frequency domain representations of multiple interactions in terms of transfer functions, partial coherence, and partial power spectrum decomposition
. This feature is extremely helpful in the study of brain signals that usually exhibit oscillatory components in wellknown frequency bands, resulting from the activity of neural circuits operating as a network
It is important to distinguish between time and frequencydomain techniques, as the latter reveal connectivity mechanisms related to the brain rhythms that operate within specific frequency bands
^{[21]}^{[39]}. While approaches such as correlation, mutual information, Granger causality, and transfer entropy are linked to a timedomain representation of the data, some others, such as coherence, directed transfer function, directed coherence, and partial directed coherence, assume that the acquired data are rich in individual rhythmic components and exploit frequencydomain representations of the investigated signals. Although this can be achieved through the application of nonparametric techniques (Fourier decomposition, wavelet analysis, Hilbert transformation after bandpass filtering
^{[72]}), the utilization of parametric AR models has collected great popularity, allowing one to evaluate brain interactions within specific spectral bands with physiological meanings
^{[21]}. Furthermore, timefrequency analysis approaches, which simultaneously extract spectral and temporal information
^{[73]}, have been extensively used to study changes in EEG connectivity in the timefrequency domain
^{[74]}^{[75]}^{[76]}, and in combination with deep learning approaches for the automatic detection of schizophrenia
^{[77]} and Knearest neighbor classifiers for monitoring the depth of anesthesia during surgery
^{[78]}.
4. Functional Connectivity Estimation Approaches
4.1. TimeDomain Approaches
Several timedomain approaches devoted to the study of FC have been developed throughout the years. Despite phasesynchronization measures, such as the phase locking value ^{[79]} and other modelfree approaches ^{[80]} being still abundantly used in brainconnectivity analysis, linear methods are easier to use and sufficient to capture brain interactions taking place under the hypothesis that neuronal interactions are governed by oscillatory phase coupling ^{[6]}.
In a linear framework, ergodicity, Gaussianity, and widesense stationarity (WSS) conditions are typically assumed for the acquired data, meaning that the analyzed signals are stochastic processes with Gaussian properties and preserve their statistical properties as a function of time. These assumptions are made, often implicitly, as prerequisites for the analysis, in order to assure that the linear description is exhaustive and the measures can be safely computed from a single realization of the analyzed process. Under these assumptions, the dynamic interactions between a realization of M Gaussian stochastic processes (e.g., M EEG signals recorded at different electrodes) can be studied in terms of timelagged correlations. In the time domain, the analysis is performed via a linear parametric approach grounded on the classical vector AR (VAR) model description of a discretetime, zeromean, stationary multivariate stochastic Markov process, $\mathbf{S}={\left[{X}_{1}\dots {X}_{M}\right]}^{\u22ba}$ . Considering the time step n as the current time, the dynamics of $\mathbf{S}$ can be completely described by the VAR model ^{[21]}^{[24]}:
$${\mathbf{S}}_{n}=\sum _{k=1}^{p}{\mathbf{A}}_{k}{\mathbf{S}}_{nk}+{\mathbf{U}}_{n},$$
where ${\mathbf{S}}_{n}={[{X}_{1,n}\dots {X}_{M,n}]}^{\u22ba}$ is the vector describing the present state of S, and Spn=[S⊺n−1…S⊺n−p]⊺ describes its past states until lag p, which is the model order defining the maximum lag used to quantify interactions; Ak is the M×M coefficient matrix quantifying the timelagged interactions within and between the M processes at lag k; and U is a M×1 vector of uncorrelated white noise with an M×M covariance matrix Σ=diag(σ211,…σ2MM). Multivariate methods based on VAR models as in (1) depend on the reliability of the fitted model, and especially the model order. While lower model orders can provide inadequate representations of the signal, orders higher than are strictly needed tend to provide overrepresentation of the oscillatory content of the process and drastically increase noise ^{[81]}. One should pay attention to the procedure for selecting the optimum model order, which can be set according to different criteria, such as the Akaike information criterion (AIC) ^{[82]} or the Bayesian information criterion (BIC) ^{[83]}.
It should be noted that, in multichannel recordings such as with EEG data, the analysis can be multivariate, which means taking all the channels into account and fitting a full VAR model, as in (1), or it can be done by considering each channel pair separately, which means fitting a bivariate AR model (2AR) in the form of (1) with M=2:
$${\mathbf{Z}}_{n}=\sum _{k=1}^{p}{\mathbf{B}}_{k}{\mathbf{Z}}_{nk}+{\mathbf{W}}_{n},$$ where Z=[XiXj]⊺, i,j=1,…,M (i≠j), is the bivariate process containing the investigated channel pair, with Zn=[Xi,nXj,n]⊺ and Zpn=[Z⊺n−1…Z⊺n−p]⊺ describing, respectively, the present and p past states of Z; Bk is the 2×2 coefficient matrix quantifying the timelagged interactions within and between the two processes at lag k, and W is a 2×1 vector of uncorrelated white noises with 2×2 covariance matrix Λ. The pairwise (bivariate) approach typically provides more stable results, since it involves the fitting of fewer parameters but leads to loss of information due to the fact that only a pair of time series is taken into account ^{[84]}. Indeed, since there are various situations that provide significant estimates of connectivity in the absence of true interactions (e.g., cascade interactions or common inputs) ^{[21]}^{[84]}, the core issue becomes whether the estimate of pairwise connectivity reflects a true direct connection between the two investigated signals or is the result of spurious dynamics between multiple time series. To answer this question, it is recommended to take into account the information from all channels when estimating the interaction terms between any pair of time series. Even if at the expense of increased model complexity resulting in a more difficult model identification process, moving from a pairwise to a multivariate approach can significantly increase the accuracy of the reconstructed connectivity patterns. This would allow distinguishing direct from indirect interactions through the use of extended formulations obtained through partialization or conditioning procedures ^{[16]}^{[21]}^{[85]}.
4.2. FrequencyDomain Approaches
To examine oscillatory neuronal interactions and identify the individual rhythmic components in the measured data, representations of connectivity in the frequency domain are often desirable. The transformation from the time domain to the frequency domain can be carried out by exploiting parametric (modelbased) or nonparametric (modelfree) approaches. Nonparametric signalrepresentation techniques are mostly based on the definition of the power spectral density (PSD) matrix of the process as the Fourier transform (FT) of ${\mathbf{R}}_{k}$ , on the wavelet transformation (WT) of data, or on Hilbert transformation following bandpass filtering ^{[72]}. In general, they bypass the issues of the ability of linear AR models to correctly interpret neurophysiological data and the selection of the optimum model order. The latter choice can be problematic, especially with brain data, because it strongly depends on the experimental task, the quality of the data and the model estimation technique ^{[86]}. However, the nonparametric spectral approach is somewhat less robust compared to parametric estimates, since it can be characterized by lower spectral resolution; e.g., it has been shown to be less efficient in discriminating epileptic occurrences in EEG data ^{[87]}.
On the other hand, parametric approaches exploit the frequencydomain representation of the VAR model, in the multivariate (1) or in the bivariate (2) case, which means computing the model coefficients in the Zdomain and then evaluating the model transfer function
$\mathbf{H}\left(\omega \right)$ on the unit circle of the complex plane, where
$\omega =2\pi \frac{f}{{f}_{s}}$ is the normalized angular frequency and
fs is the sampling frequency
^{[21]}. The
M×M PSD matrix can then be computed using spectral factorization as
$$\mathbf{P}\left(\omega \right)=\mathbf{H}\left(\omega \right)\Sigma {\mathbf{H}}^{*}\left(\omega \right),$$
where * stands for the Hermitian transpose ^{[21]}; note that $\Sigma $ is replaced by Λ in the case of (2), i.e., when M=2. It is worth noting that, while the frequencydomain descriptions ubiquitously used are based on the VAR model representation, their key element is actually the spectral factorization theorem reported above and that approaches other than VAR models can be used to derive frequencydomain connectivity measures ^{[88]}.
4.3. InformationDomain Approaches
The statistical dependencies among electrophysiological signals can be evaluated using information theory. Concepts of mutual information, mutual information rate, and information transfer are widely used to assess the information exchanged between two interdependent systems ^{[60]}^{[69]}, the dynamic interdependence between two systems per unit of time ^{[22]}^{[23]}, and the dynamic information transferred to the target from the other connected systems ^{[13]}^{[70]}, respectively. The main advantage of these approaches lies in the fact that they are probabilistic and can thus be stated in a fully modelfree formulation. On the other hand, their practical assessment in the information domain is not straightforward because it comprises the estimation of highdimensional probability distributions ^{[89]}, which becomes more and more difficult as the number of network nodes increases in multichannel EEG recordings. Nevertheless, informationbased metrics can also be expressed in terms of predictability improvement, such that their computation can rely on linear parametric AR models, where concepts of prediction error and conditional entropy, GC and information transfer, or TD and mutual information rate, have been linked to each other ^{[19]}^{[20]}^{[90]}^{[91]}. Indeed, it has been demonstrated that, under the hypothesis of Gaussianity, predictability improvement and informationbased indexes are equivalent ^{[19]}. Based on the knowledge that stationary Gaussian processes are fully described in terms of linear regression models, a central result is that for Gaussian data, all the information measures can be computed straightforwardly from the variances of the innovation processes of full and restricted AR models ^{[19]}.
4.4. Other Connectivity Estimators
Brain connectivity can be estimated through a large number of analyses applied to EEG data. Multivariate timeseries analysis has traditionally relied on the use of linear methods in the time and frequency domains. Nevertheless, these methods are insufficient for capturing nonlinear features in signals, especially in neurophysiological data where nonlinearity is a characteristic of neuronal activity
^{[80]}. This has driven the exploration of alternative techniques that are not limited by this constraint
^{[80]}^{[92]}. Moreover, the utilization of AR models with constant parameters, and the underlying hypotheses of Gaussianity and WSS of the data, can be key limitations when stationarity is not verified
^{[93]}. A number of approaches have been developed to overcome this issue, providing timevarying extensions of linear modelbased connectivity estimators using adaptive AR models with timeresolved parameters, in which the AR parameters are functions of time
^{[94]}^{[95]}^{[96]}.
5. EEG Acquisition and PreProcessing
The acquisition and conditioning of the EEG signal represent two important aspects with effects on the entire subsequent processing chain. The main steps of acquisition and preprocessing are indicated in Figure 1. An example of application of the preprocessing pipeline to experimental EEG is shown in Figure 2.
Figure 1. Main steps in EEG acquisition and preprocessing. In general, source localization is not mandatory, as represented by the dashed round brackets.
Figure 2. Schematic representation of the preprocessing pipeline applied to EEG signals acquired on the scalp (s253 recording of the subject 2 that could be found in
https://eeglab.org/tutorials/10_Group_analysis/study_creation.html#descriptionofthe5subjectexperimenttutorialdata, accessed on 15 February 2023). (
A) Unipolar EEG signals are acquired using a mastoid reference (Ref, in red). For clarity, only a limited number of the recorded signals, among the original 30 channels, is plotted. The average rereferencing process and the preprocessed signals are illustrated below. Notably, red arrows indicate blinking artifacts that are clearly visible. (
B) The rereferenced signals are filtered using a 1–45 Hz zerophase passband filter, followed by independent component analysis (ICA) to extract eight independent components (ICs), shown on the right. (
C) The first IC, suspected to be an artifact, is analyzed, with a scalpshaped heatmap assessing its localization in the frontal area and its temporary coincidence with the artifacts shown in panel (
A). After removing the first IC, the cleaned signal is plotted at the bottom of the panel.
Sampling frequency, the number of electrodes, and their positioning, each have an important role in assessing connectivity; too low of a sampling frequency cannot be used to analyze highfrequency bands for the Nyquist–Shannon sampling theorem ^{[97]}, and the electrode density defines with which accuracy and reliability further processing will be performed ^{[40]}.
An EEG signal is a temporal sequence of electric potential values. The potential is measured with respect to a reference, which should be ideally the potential at a point at infinity (thus, with a stable zero value). In practice, an ideal reference cannot be found, and any specific choice will affect the evaluated connectivity
^{[6]}^{[98]}. Unfortunately, as for most of the FC analysis pipeline, there is no gold standard for referencing, and this is clearly a problem for crossstudy comparability
^{[99]}.
5.1. Resampling
Nowadays, EEG signal is usually acquired with a sampling rate (SR) equal or superior to 128 Hz, since this value is the lower base2power that permits one to capture most of the information content in EEG signal and a big part of the γ band. It is noticeable that highfrequency oscillations (HFOs), such as ripples (80–200 Hz) and fast ripples (200–500 Hz), will not be captured with these SRs ^{[100]}^{[101]}^{[102]}^{[103]}. Moreover, there are studies suggesting that low SRs affect correlation dimensions, and in general, nonlinear metrics ^{[104]}^{[105]}.
5.2. Filtering and Artifact Rejection
Filtering EEG is a necessary step for the FC analysis pipeline, not only to extract the principal EEG frequency waves, but particularly to reduce the amounts of noise and artifacts present in the signal and to enhance its SNR. Research suggests the use of finite impulse response (FIR) causal filters that could be used also for realtime (RT) applications, or infinite impulse response (IIR) filters, which are less demanding in terms of filter orders but distort phase, unless they are applied with reverse filtering, thereby making the process noncausal and not applicable for RT applications. In general, if sharp cutoffs are not needed for the analysis, FIR filters are recommended, since they are always stable and easier to control ^{[106]}.
If one is interested in investigating FC in the
γ
band, electrical line noise at 50 or 60 Hz could be a problem, since it is not fully removable with a lowpass filter. Notch filters are basically bandstop filters with a very narrow transition phase in the frequency domain, which in turn leads to an inevitable distortion of the signal in the time domain, such as smearing artifacts ^{[106]}. To avoid this problem, some alternatives have been developed. A discrete Fourier transform (DFT) filter is obtained by subtracting from the signal an estimation of the interference obtained by fitting the signal with a combination of sine and cosine with the same frequency as the interference. It avoids potential distortions of components external to the powerline frequency. It assumes that the powerline noise has constant power in the analyzed signal segment. As this hypothesis is not strictly verified in practice, it is recommended to apply the DFT technique to short data segments (1 s or less).
Another proposed technique is CleanLine, a regressionbased method that makes use of a sliding window and multitapers to transform data from the time domain to the frequency domain, thereby estimating the characteristics of the power line noise signal with a regression model and subtracting it from the data ^{[107]}. This method eliminates only the deterministic line components, which are optimal, since EEG signal is a stochastic process, but in the presence of strong nonstationary artifacts, it may fail ^{[108]}.
It is pretty normal that EEG signals can be corrupted by many types of artifacts, defined as any undesired signal which affects the EEG signal whose origin cannot be identified in neural activity. Generators of these undesirable signals could be physiological, such as ocular artifacts (OAs such as eye blink, saccade movement, rapid eye movements), head or facial movements, muscle contraction, or cardiac activity ^{[109]}^{[110]}. Powerline noise, electrode movements (due to nonproperly connected electrodes), and interference with other electrical instrumentation are nonphysiological artifacts ^{[111]}. Artifact management is crucial in the analysis of connectivity. In fact, their presence in multiple electrodes can result in overestimation of brain connectivity, skewing the results ^{[112]}^{[113]}.
The most effective way for dealing with artifacts is to apply prevention in the phase of acquisition to avoid as much as possible their presence, for example, by making acquisitions in a controlled environment, doublechecking the correct positioning of electrodes, or instructing the patients to avoid blinking the eyes in certain moments of the acquisition. In fact, eye blinks are by far the most common source of artifacts in the EEG signal, especially in the frontal electrodes ^{[114]}^{[115]}. They are shown as rapid, highamplitude waveforms present in many channels, whose magnitudes exponentially decrease from frontal to occipital electrodes. However, saccade movements produce artifacts, generating an electrooculographic signal (EOG). Acquisition of this signal, concurrently with the EEG signal, is known to be a great advantage in identifying and removing ocular artifacts, since vertical (VEOG), horizontal (HEOG), and radial (REOG) signals diffuse in different ways through the scalp ^{[116]}.
Once the parts of the signals corrupted by the artifacts have been identified, it is still common practice to eliminate these portions, avoiding alterations on the signal that could lead to spurious connectivity. As for channel rejection, however, it is preferable to retain as much information as possible ^{[117]}.
5.3. Bad Channel Identification, Rejection, and Interpolation
It could happen that some EEG channels present a high number of artifacts (eye blink, muscular noise, etc.) or noise, due to bad electrodescalp contact. In these cases, the rejection of these channels could be an option. However, it is necessary to check whether the deleted channels are not fundamental and the remaining channels are sufficient to carry on the analysis, considering also that deleting channels will result in an important loss of information that is likely not recovered anymore. Some authors discussed the criteria for detection of bad channels and suggested considering the proportion of bad channels with respect to the total to assess the quality of the dataset (for example, imposing a maximum of 5%) ^{[118]}.
Identification of bad channels could be performed visually or automatically. Visual inspection requires certain experience with EEG signals to decide if a channel is actually not recoverable and needs to be rejected
^{[119]}^{[120]}. Indeed, this process is highly subjective. Automatic detection of bad channels could be performed in various ways. A channel correlation method identifies the bad channels by comparing their divergence from the Pearson correlation distribution among each pair of channels with the other couples. This method assumes high correlations between channels due to the volumeconduction effect proper of the EEG signal, which is discussed in detail in the source localization section. High standard deviation of an EEG signal can be an indicator of the presence of a great amount of noise. By setting a proper threshold, the standard deviation can be a useful index for identifying bad channels.
5.4. ReReferencing
As described in the previous section, many decisions are made during the acquisition of the signals; however, some of them are revertible. Rereferencing and downsampling is one of them. It consists of changing the (common) reference of each EEG channel into another one by performing for each channel the addition of a fixed value (notice that it is a fixed value in space, i.e., common for all channels, but not in general constant in time).
Rereference with respect to a new reference channel can be performed by simply subtracting each channel with the new reference. In most cases, average reference (CAR) is considered to be a good choice, especially if electrodes cover a large portion of the scalp with the assumption that the algebraic sum of currents must be zero in the case of uniform density of sources
^{[121]}. It consists of referencing all the electrode signals with respect to a virtual reference that is the arithmetic mean of all the channel signals. Due to charge conservation, the integral of the potential over a closed surface surrounding the neural sources is zero. Obviously, the EEG channels give only a discrete sampling of a portion of that surface, so that the virtual reference can be assumed to be only approximately zero. However, it could be much better than using a single location as reference, such as an EEG channel or the mastoid (both LM and LR), which have been shown to generate larger distortion than the average reference
^{[112]}.
6. Source Connectivity Analysis
Measuring information dynamics from EEG signals on the scalp is not straightforward due to the impact of volume conduction, which can modify or obscure the patterns of information flow across EEG electrodes ^{[122]}. This effect is due to the electrical conductance of the skull, which serves as a support to diffuse neural activity in all the directions and for this reason is also known as field spread problem. The neural current sources are related to the Poisson equation to the electrical potential, which diffuses across the scalp and can be measured by many electrodes, also pretty far from the original source. This is the reason why interpreting scalplevel connectivity requires caution. In fact, the estimated FC between two electrodes could be reflecting the activation of a single brain region, rather than two functionally connected regions. Even though this effect can be compensated for when working with scalp EEG signals ^{[122]}^{[123]}, it is often recommended to use source signal reconstruction to obtain a more accurate representation of the underlying neural network. This is because the sourcebased network representation is considered a more accurate approximation of the real neural network structure ^{[124]}.
The general pipeline that source localization algorithms follow is based on this twostep loop ^{[125]}:

Forward problem—definition of a set of sources and their characteristics and simulation of the signal that would be measured (i.e., the potential on the scalp) knowing the physical characteristics of the medium that makes it diffuse;

Inverse problem—comparison of the signal generated by the head model with the actual measured EEG and adjustment of the parameters of the source model to make them as similar as possible.
The first part, also known as estimation problem, can be carried out by defining a proper head model using boundary element models (BEMs) or finite element models (FEM). The second part of the pipeline is also called the appraisal problem, and it is not trivial at all, being one of the fundamental challenges in EEG processing analysis. This process is in fact formally an inverse illposed problem, where there exists an infinite number of combinations of sources that can explain the acquired signal ^{[126]}.