Abstract
In recent years, functional connectivity in the developmental science received increasing attention. Although it has been reported that the anatomical connectivity in the preterm brain develops dramatically during the last months of pregnancy, little is known about how functional and effective connectivity change with maturation. The present study investigated how effective connectivity in premature infants evolves. To assess it, we use EEG measurements and graph-theory methodologies. We recorded data from 25 preterm babies, who underwent long-EEG monitoring at least twice during their stay in the NICU. The recordings took place from 27 weeks postmenstrual age (PMA) until 42 weeks PMA. Results showed that the EEG-connectivity, assessed using graph-theory indices, moved from a small-world network to a random one, since the clustering coefficient increases and the path length decreases. This shift can be due to the development of the thalamocortical connections and long-range cortical connections. Based on the network indices, we developed different age-prediction models. The best result showed that it is possible to predict the age of the infant with a root mean-squared error () equal to 2.11 weeks. These results are similar to the ones reported in the literature for age prediction in preterm babies.
1. Introduction
The brain can be seen as a complex network of interacting regions and hierarchical communications, which are constrained by the anatomy, but not limited to it [1]. The neuronal clusters can actually work together and communicate to perform a joint task beyond their structural locations. The clinical literature [2] distinguishes this type of connectivity from the anatomical one, which is often called structural. The consequence of this functional infrastructure is the generation of complex electrophysiological patterns, which are temporally correlated, by distant cerebral areas [3]. Those spatiotemporal patterns are dynamic; they change according to the individual development trajectory [4]. In particular, the last trimester of gestation is a period of brain development, which includes both anatomical rewiring [5] and electrophysiological modifications [6]. Different authors illustrated that the cortical regions undergo differentiation, folding, and gyrification, while the subcortical areas experience synaptogenesis and myelination as well as neural pruning to establish thalamocortical connections or long distance cortical connections [7, 8]. Based on MRI scans of preterm babies, Dubois et al. [8] showed that the white matter volume and the inner cortical surface increase with gestational age. Furthermore, the same author [8, 9] demonstrated that fractional anisotropy (FA) of the brain fiber bundles increases with postbirth maturation, although the different connection pathways seem to develop in an asynchronous way. According to Batalle et al. [5], FA is a measure of anatomical directivity that describes the connectivity strength among brain regions together with the density and the percentage of connecting streamlines. Hüppi and Dubois [9] eventually argued that diffusion tensor imaging parameters, such as anisotropy, can be structural markers of network functional organization [5]. A possible implication is that the temporal correlation among brain regions is also expected to change since the functional connectivity is related to the structural topology [1]. However, less is known about functional connectivity and maturation [10]. Therefore, functional connectivity in developmental science received increasing attention in recent years [1]. Moreover, different tools to describe functional connectivity became available in the last two decades. According to Bullmore and Sporns [11], brain networks can be represented as a graph, where the nodes are brain regions and the edges are the connection strengths. The nodes are defined by the neural activity, while the edges define the presence of a coupling between time series and the associated weight represents the coupling degree. Friston [12] recognizes two main types of functional connectivity: the actual functional connectivity, which investigates the temporal correlation among signals, and the effective connectivity, which measures the causality of the coupling. While the former does not imply any directionality, the second one highlights which area or node causes the activity of other ones. Functional connectivity methods generate undirected graphs, while the effective connectivity methods are associated with directed graphs [11]. In different studies of the preterm infant brain [5, 13], the neural activity has been measured using fMRI. Although this neuroimaging technique can investigate the subcortical structures, it is an expensive method and it is not suitable to measure effective connectivity due to fMRI low temporal resolution [14]. In contrast, EEG is a suitable measurement for effective connectivity and has been employed in recent papers to study the brain connectivity maturation [14, 15]. The main drawbacks of those studies are the limited investigated maturation period or the absence of graph-theory metrics to investigate the connectivity changes. Studies about the long term development of effective connectivity based on EEG have not been performed yet. In particular, there is a lack of research that investigates the effective connectivity from birth until full-term age in premature infants. Furthermore, in the literature, there are a few works that study the connectivity maturation in the preterm brain from a graph point of view [5, 13]. The main aim of the present study was a thorough investigation of the effective connectivity evolution during the preterm maturation after birth, exploiting different network metrics (as in [5]). In particular, we measured the effective connectivity by means of transfer entropy [16] and Granger causality [17]. On the obtained directed graph, different graph metrics were computed to track their evolution from 27 to 42 postmenstrual age (PMA) weeks. The network analysis was complemented by a regression study to predict the age of the patient and assess the predictive power of the connectivity analysis (as performed in one of our previous studies [15]). The paper is organized as follows: in Section 2, we introduce the simulated and the real dataset alongside the methods used to estimate effective connectivity. In Section 2, the applied graph-theory metrics are also discussed. Section 3 reports the results, while Section 4 includes our discussion on the reported results.
2. Methods
2.1. Data
Different datasets were used. A simulated dataset was also employed in order to show how connectivity analysis performs in a controlled case. In particular, the main objective of the simulation was to illustrate the meaning of the most common graph indices used in the literature. A dataset comprised of EEG signals was the main part of the maturation study.
2.1.1. Simulated Data
The first experiment of the study consisted of a simulation based on a linear Gaussian regression model (derived from [18]), expressed by the following equations:
The parameter is a scalar value that varies between 0 and 1 and it influences the strength of the coupling among the variables. The variables represent white noise with unit variance. The associated graph is reported in Figure 1. In order to give the reader a clear insight of the graph-theory measures, the objective of the simulation was twofold: firstly, we investigated the influence of the strength of the coupling, provided by the parameter , and show the effect of the weakening of causality among the variables. Secondly, we investigated the influence of the noise using different levels of signal-to-noise ratio (SNR), which were obtained varying the variance of white noise in (1). In the latter objective, time series were simulated using the AR model (1) and the coupling was estimated with the effective connectivity methods discussed in the following paragraphs. We also studied the impact of filtering on the graph based metrics by computing the scores from filtered signals, using a sampling frequency which was set at 200 Hz and a band-pass filter between 1 and 80 Hz was applied on the simulated data. The reason to investigate the impact of filtering is explained in Section 2.3. The band-pass frequencies were obtained from [19].
2.1.2. EEG Data
The second dataset comprised 25 preterm infants, with gestational age (GA) ≤ 32 weeks, who were recruited for a larger EEG study to assess brain development and automatically detect quiet sleep epochs [20, 21]. Each patient was enrolled with informed parental consent at the Neonatal Intensive Care Unit (NICU) of the University Hospital of Leuven, Belgium. All the included subjects presented good outcome at 2 years. EEG was recorded with 9 electrodes (, , , , , , , , and reference ) according to the 10–20 international system (BRAIN RT, OSG Equipment, Mechelen, Belgium). Throughout connectivity analysis, the channel was disregarded. The measurements were performed twice during their stay in the unit, resulting in 88 recordings ranging from postmenstrual age of 27 weeks to 42 weeks. Mean time of recording EEG was 4 h 55 min. Monopolar signals were recorded at a sampling frequency of 250 Hz. Two independent EEG readers annotated the data for the different sleep stages, namely, quiet sleep (QS) and nonquiet sleep (NQS) epochs. The reason for this binary manual classification is due to the difficulty to discriminate active sleep from the awake state in the very young child.
2.2. Connectivity Analysis
The literature about effective connectivity presents different methods to assess coupling among time series, such as directed transfer function (DTF), partial directed coherence, or Granger causality test (GCT) [22]. With these approaches, the coupling is frequently defined as G-causality [23], in the sense that a time series Granger-causes another one when the past values of a time series significantly explain the evolution of the other one. Even though Granger himself pointed out the limited applicability to the bivariate case, the recently published conditional multivariate Granger causality (cMVGC) [24] overcomes the problem. According to [23, 25], these methods are based on the same multivariate or vector autoregressive (VAR) modeling. However, they can show different aspects in the causality analysis. While the DTF estimates the reachability from one channel to another (defined also as G-influentiability by [23]), both PDC and cMVGC look into the direct active link between two channels, which is defined as G-connectivity [23]. Those two methods investigate the same connectivity model in two different domains, which are the frequency (PDC) and time (cMVGC) domain. Since in this paper we are mainly interested in the time-domain coupling, we opted for the cMVGC in its autoregressive (AR) and information dynamics implementations, respectively, described in [18, 24].
2.2.1. Granger Causality
The first method employed to assess the G-connectivity among EEG channels was Granger causality (GC) with a VAR modeling framework. In the present study, we followed the formulation by [24]. A th order multivariate autoregressive model VAR() takes the form where is a collection of process , while and are, respectively, the regression coefficient matrices and the stochastic process residuals. In the conditional scenario [24], where we want to know the pairwise influence of process over considering the presence of the other variables, can be rewritten aswhere and are the two time series of interest, while represents the third set of variables involved in the analysis. In our study, and can be two variables in (1), for example, and , or two EEG channels, for example, and . Consequently, will be a vector time series which contain the remaining in (1) or the remaining EEG channels. Based on the VAR() model (2) and the split defined in (3), we may actually explain the future of based on the following full and reduced regressions: In both cases we consider the conditioning variable , although only the first model considers explicitly the influence of . The difference is also highlighted by the two regression coefficient matrices and . In order to test whether the coefficient matrices are significantly different from zero, the Granger causality is defined as the log-likelihood ratio where and are the covariance matrices of the residuals and . Based on this multivariate framework, the G-causality basically quantifies the reduction of the prediction error if the variable is added to explanatory variables of , one of which is the variable [24]. Unlike the classical Granger causality test limited to the bivariate case, we defined conditional multivariate Granger causality. A special case of the conditional scenario is the pairwise conditional G-causality, where the pairwise effective couplings among all pair of variables and contained in are measured. Since we take in account the spurious effect due to the presence of other variables (i.e., the coupling between and conditioned by the presence of the other variables), the pairwise conditional causalities are defined as which is an matrix, where is the number of processes, and contains all the pairwise coupling estimations. Those quantities may be considered a weighted directed graph, also known as G-causal graph, and the matrix as the associated adjacency matrix . The G-causality was computed with the multivariate Granger causality toolbox [24] implemented in MATLAB (Mathworks, Natick, MA, USA).
2.2.2. Transfer Entropy
The second method applied to assess the G-connectivity among EEG channels was transfer entropy. According to the information dynamics framework by [16], the expected Kullback-Leibler divergence defines the transfer entropy (TE) from process to process as follows:where is the conditional probability that at times is explained by past values of and with respective memory order and . is the conditional probability that at times is only explained by past values of . is the joint probability distribution among the three variables , , . Transfer entropy inherently implies directionality (i.e., effective connectivity), since it is asymmetric and contains transition probabilities [16, 26]. In order to estimate the information dynamics coupling, Montalto et al. [18] pointed out that TE is equivalent to the difference of two conditional entropies (CE)
With respect to (7), the variable is here introduced. In a multivariate analysis, we cannot discriminate the exclusive relationship between two processes [17, 24]. However, it is possible to investigate the contribution of time series in the evolution of time series with respect to all the other agents involved in the analysis. Consequently, we can actually define as a vector variable that does not contain neither nor as past values and illustrates the transfer of information from to taking into account other time series involved. If we assume that , , have a joint Gaussian distribution, the two CE in (8) can be expressed as linear regression of past values of the vector variables involved in the multivariate system as follows:
The first equation explains with a regression on the vector , where , , approximate respectively , , with a vector of size as follows:The second equation explains with a regression on the vector , which only contains , . Equations (9) are usually referred to as full and restricted regressions. The reader can easily recognize the same structure of Granger causality equations (4), which is asymptotically equivalent to transfer entropy [16]. At the light of the previous results, can be rewritten as where and are the variances of the white noise residuals in (9). Except for a scalar factor, (11) is the same expression of (5). It is important to notice that also TE estimates the G-connectivity in a conditioned framework. In the present study, the TE entropy has been computed using the MUTE toolbox [18] implemented in MATLAB (Mathworks, Natick, MA, USA).
2.2.3. Graph-Theory Indices
The effective connectivity methods generate an adjacency matrix that describes a weighted directed graph [24], called G-causal graph. Unlike a binary network, those graphs tend to be complete or full without a specific thresholding and the adjacency matrix is asymmetric. Since the maximal number of couplings is (where is the number of signals), graph-theory measures can be used to summarize the causal connectivity [27]. Although there is no minimal theoretical number of graph nodes, 8 EEG channels can be considered quite limited for a graph analysis. However, there are studies in the neural processing literature where graph theory was applied on a limited number of time series in order to give a concise view of a high-density (or complete) network [27, 28], in particular if the sources-level is concerned. Bullmore and Sporns [11] illustrated a list of measures to describe the graph topology, such as the average characteristic path length, the clustering coefficient, and the diameter. The path length represents the minimum number of edges to reach each other node in the graph starting from one specific node. The average path length is the mean of the nodes’ shortest paths and can be considered a measure of network integration capacity [5]. The latter can be also investigated by means of the clustering coefficient. In case of a weighted network, a clustering coefficient is defined as average intensity of all triangles around each node [29]. A common overall measure is the average of nodes clustering coefficients, defined in a similar way to the average path length. The last index to evaluate the integration capacity of the graph is the diameter. Let the node eccentricity be the maximum graph distance from one node to any other node in the network [30]. The diameter is then defined as the maximum eccentricity in the graph. All those indices have been computed thanks to the brain connectivity toolbox [31] implemented in MATLAB (Mathworks, Natick, MA, USA). In addition, we computed also the causal density as the sum of all significant couplings of the adjacency matrix, as defined by [32]. Alongside those network metrics, the spectral indices such as the spectral radius, the spectral gap, and the algebraic connectivity were considered. The spectral radius and the spectral indices were the maximal absolute eigenvalue and the difference between the first and the second absolute eigenvalues of the adjacency matrix, respectively [33, 34]. The spectral radius, commonly known as “Page Rank,” represents the dominance degree of a node in the network: the higher the value, the higher the centrality of the dominant node in the network. In other words, the dominant node behaves as the center of a hub [33]. Another way to look into the change of dominance is the spectral gap: if the difference between the first and second eigenvalues in the graph matrix decreases, the node with highest eigenvalue (the spectral radius) is less dominant with respect to the other nodes in the network. However, Estrada and Hatano [34] argued that the spectral gap behaves as a measure of clustering in the undirected graph: in case of lower values of spectral gap, the graph can present small clusters in the network. The last spectral measure is the algebraic connectivity, which is the second smallest eigenvalue of the Laplacian matrix and it illustrates how easily a graph can be divided into clusters [30]. In case of a directed graph, the Laplacian matrix can be obtained as follows:where is the adjacency matrix obtained by the effective connectivity tools described above and is the degree matrix of the associated graph [35]. An overview of all graph indices is reported in Table 1.
2.3. Algorithmic Pipeline and Statistical Analysis
According to different authors [19, 36], filtering could add spurious connectivity in the effective connectivity analysis or make the estimation of the underlying Granger causality VAR model unstable. Although a theoretical invariance of causality estimation has been demonstrated under filtering, the G-causality works in practice if, and only if, the data are stationary. The use of filtering as mean to reach stationarity with filtering has already been investigated by [36]. However, three main reasons can undermine this approach. The first reason is the increase of the estimated VAR model order due to the fitting of filtered process of theoretical infinite order with a numerical finite order. This could lead to a poor, and not robust, parameter estimation or even unstable VAR model, which is the second reason to avoid band-pass filtering. The last reason is the ill-conditioning of the Toeplitz matrix of the autocorrelation sequence , which is necessary for VAR model estimation. All the related theoretical details are explained in [36]. The practical downsides of the band-pass filtering are the dramatic increase of the VAR model order compared to the unprocessed data or the increase of false positive detection of connectivity links with both the GC [36] and the PDC [19]. In particular, the amount of false detection increases with narrowing of the filtering frequency band or the increase of the filter order. In general there is no distinction between FIR and IIR filtering for both the authors. However, in both studies [19, 36], notch filtering and differentiation (or high pass filtering at 1 Hz) were pointed out as methods to keep the VAR model order low and reduce the false detection, even compared to unprocessed data. Those approaches help to achieve stationarity or keep the VAR model order low for nearly nonstationary processes like EEG. In addition, the presence of trends or seasonality can add unit-roots to the time series (poles on the unit circle or outside it in the complex plane), which violates the hypothesis of covariance stationarity. In the presence of unit-roots, the impulse response of the VAR model would be oscillatory or diverge to infinity. Consequently, it is suggested to eliminate trends and seasonality by first-order differentiation or differentiation at various lags (notch filtering). Theoretical and numerical details of the different type of filtering are reported in [24, 36]. Given the stated literature results, on the one hand, we decided to investigate the impact of filtering on the simulated data and, on the other hand, we applied only notch filtering at 50 Hz and 100 Hz and differentiation on the EEG data in order to reduce nonstationarities in the time series. However, the EEG can be affected by muscle artifacts, which can spread out among the different electrodes and bias the connectivity analysis. To mitigate this effect, we applied canonical correlation analysis (CCA) to remove the artifacts caused by the change in the EMG signals [37]. To investigate the effect of the CCA on the analysis, we compared the output of the connectivity analysis in two scenarios: the first one did not consider the usage of CCA; the second one applies CCA and reconstructed the EEG removing the 3 sources with lowest autocorrelation. In the first case, the authors segmented the EEG in 5 s windows and computed the effective coupling with both listed methods. In the second scenario, CCA was applied on 5 s EEG segments. However, before performing any connectivity analysis, we recombined the segments in 30 s intervals in order to increase the rank of the reconstructed EEG matrix. This step was required since both GC and TE request data that do not present collinearity (i.e., the time series matrix cannot be rank-deficient [24]). The window length was suggested as a further step to keep time series stationarity [24]. For each recording, we computed the adjacency connectivity matrix for EEG segments and we averaged over the QS epochs and NQS epochs. This averaging step did not consider coupling values that were not statistically significant. According to [18, 24], the GC and TE follow an asymptotic -distribution, like statistics. Therefore, an -test was implemented to test the significance of coupling among EEG channels and all the couplings with were set to zero. Consequently, coupling graphs were obtained, where 88 is the total number of recordings, 2 is the employed causality methods, 2 is the considered number of sleep states. On the average matrices we computed the network indices described above, which results in a tensor , where 7 is the number of graph features and 4 is the number of combinations considering the number of connectivity methods and sleep states involved (TE in QS, TE in NQS, GC in QS, and GC in NQS). For each feature, we evaluated the maturation trend in three distinctive age groups (≤31, , ≥37 PMA weeks) as median (IQR), where IQR is InterQuartile Range. Besides, we computed the Pearson correlation coefficient between the variable age and each single feature and its statistical significance were computed. Those results were meant to give a general overview of the feature maturation trend for each connectivity method, for each sleep state, with or without CCA preprocessing. In addition, we computed an ordinary least squares (OLS) regression for each single feature in case of GC during QS and the associated confidence interval at , in order to give a visual representation of any network index prediction power. In particular, we split randomly the single graph index dataset 100 times in 70% training set and 30% testing set and we computed the prediction error on the test set as root mean-squared error, , as shown by [38]. Furthermore, each slice of the tensor (a matrix ) was used to predict the patient’s age at the moment of the recording with a multivariate linear regression. Similar to the single feature approach, we randomly split the dataset 100 times, as in [38] and we assessed in the test set. At each iteration, index and the -test statistics were computed.
3. Results
In the following paragraphs, the results obtained in both the simulated and the real dataset are discussed. In particular, we will show how the network indices behave in both examples, the simulated data and the EEG maturation dataset. The last part summarizes the predictive power of those graph metrics to infer the postmenstrual age of the subject.
3.1. Simulation Study
Figures 2(a), 2(b), and 2(c) display the integration graph indices for the AR model defined in (1), whose network is depicted in Figure 1. In the original model (), the graph presents two distinct clusters, which are loosely connected by the edges between node 1 and nodes 4 and 5. This two-hubs structure is reflected by Figures 2(a) and 2(b). When the intracluster connectivity is high (), the clustering coefficient reaches its highest level, while the path length is at its lowest level. Those results are in line with a rich-club or small-world network [11]. However, when the coefficient starts decreasing, the clustering coefficient proportionally decreases, while the path length increases. The spectral radius decreases with vanishing values of , as the clustering coefficient does. Figures 2(d), 2(e), and 2(f) display the ratio between the network indices estimated from the simulated time series via transfer entropy and the original measures (in particular, in (d), , in (e), , in (f), ). Different noise levels have been used. The results are quite straightforward for the clustering coefficient and the spectral radius: the higher the signal-to-noise ratio, the higher the values of two indices are and the closer they are to the original values (Figures 2(d) and 2(f)). In the case of the path length, the absolute value is decreasing with higher SNR. However, the estimated value becomes similar to the original one for very low noise variance (Figure 2(e)). Figures 2(d), 2(e), and 2(f) show also the effect on band-pass filtering on the graph indices estimation. The most remarkable results are related to the path length and the clustering coefficient. The estimated clustering coefficient tends to underestimate the original one, while the estimated path length tends to be persistently higher than the original one. On the contrary, the ratio for the spectral radius in case of filtering behaves similar to the one obtained using the raw data.
(a)
(b)
(c)
(d)
(e)
(f)
3.2. EEG Data: Graph Indices
Figure 3 shows graph metrics for the adjacency matrices obtained using the EEG measurements and Granger causality in QS (GC in QS). Figures 3(a), 3(b), 3(c), and 3(d) show the scatter plots with the fitted OLS regression model, while Figures 3(e) and 3(f) display the clustering coefficient and path length dynamics in three distinct age groups. As already mentioned in Section 2.3, Figure 3 gives a visual representation of the trends for the graph features, while Tables 2 and 3 provide a complete overview for all coupling methods and all sleep states. Each single feature has a significant trend with age, although the Pearson correlation coefficient increases when CCA is used as a preprocessing step. Specifically, the trend for the clustering coefficient, the spectral radius, and the spectral gap is negative, while the path length is increasing with age. This result is persistent in each method and each sleep state. The connectivity weakening for GC in QS is also reported in Figure 4, which shows the average connectivity graph for three distinct age groups. The three panels show how the coupling among time series decreases by the reduction in arrows width and the color shift from red to blue. Table 4 reports the results for age prediction with a multivariate linear model, combining all the network features. All the models can predict the age of the infant recording with a between 2 and 3 weeks and the CCA models always outperform the model without CCA as a preprocessing step. Furthermore, the explained variance () is higher with the models that include CCA. It is also interesting to notice that best prediction results are obtained with GC during QS ( PMA weeks, PMA weeks). Table 4 does not report the results for one single model estimation, but the median and InterQuartile Range (IQR) of the evaluation parameters for 100 bootstrap iterations. In each single iteration, the model proved to be significant as reported by the value column in Table 4 ().
(a)
(b)
(c)
(d)
(e)
(f)
4. Discussion
In the present study, we quantified the effective brain connectivity in preterm infants to track their maturation. Although there are some studies that investigate connectivity in the neonates [4, 5] and its change in the first days of life (short maturation period) [14], this is the first study to track maturity using effective EEG-based connectivity in preterm patients on a wide maturation period, from birth to full-term age. To the best of our knowledge, MRI [8] and fMRI [10] have been the leading method to track maturity. However, fMRI is only suitable for functional connectivity, as discussed above. In this article, we compared two well-known methods to estimate coupling between processes, like GC and TE. Based on the obtained connectivity matrices, we estimated integration and spectral network indices for directed weighted graphs. Those features were used to predict the age of the patient during the recording. The same process was applied on a simulated dataset to investigate how the network measures behaved in a controlled case.
4.1. Simulated Dataset
According to [11], graphs with high clustering coefficient and low path length behave like a rich-club network, while graphs with low clustering coefficient and high path length denote a random network, where the number of edges for each node is normally distributed. Figure 1 portraits two club networks, where the nodes are connected to each other with a short distance. This leads to a high clustering coefficient and low path length when the coupling coefficient is equal to 1. However, when decreases, the intracluster connectivity weakens and the graph becomes more similar to a random network. This type of graph is characterized by low clustering coefficient and high path length: indeed, the nodes are less connected among each other in a more homogeneous network. This result is also supported by the direct proportionality between spectral radius and coupling coefficient. Another interesting point is related to the filtering. Figure 2 illustrates clearly how the clustering coefficient and the path length estimation can be highly affected by the filtering. Those results are in line with the analysis by [36], which demonstrated that careless filtering can add spurious connectivity in the time courses. In our simulation, the effect of filtering weakens the intracluster connectivity (adding intercluster connectivity). The net effect is a decrease in clustering coefficient and an increase in path length. Therefore, we decided only to apply a notch filter and differentiation as preprocessing steps on the EEG.
4.2. EEG Data
In the literature, a number of studies can be found to have assessed the brain maturation in children and adolescents by graph theory [1, 4, 10] and a few papers focused on preterm brain maturation by network metrics [5, 13]. The first result obtained in our study is the change in effective connectivity with age. Although Schumacher et al. [14] used a different method, they also concluded that there is a change in effective connectivity mainly driven by postnatal maturation. In addition to that, we have been able to demonstrate the existence of a relationship between connectivity and PMA. In this study, we also observed a change in graph parameters that suggest that the EEG-scalp network moved from a rich-club network to a more random network. The integration and spectral indices decreased with age, except for the path length, which reflects a segregation of nodes due to a higher graph distance as well as less intense triangle patterns around the nodes themselves. Hub-networks have high clustering coefficients since they have central club nodes, which are surrounded by triangle patterns. On the contrary, a random network presents nodes, which are connected to any other node in the network with a weak coupling. The net effect is a low clustering coefficient and high path length. This emergence of a normal-distributed network is also confirmed by the decrease of the spectral gap, spectral radius, and the algebraic connectivity. The latter two indices emphasize how easily the graph can be clustered and a negative trend would suggest the absence of modules or groups in the graph. On the contrary, a negative trend of the former index should suggest an increase of modularity, as shown by [34]. However, since the spectral gap is also inversely proportional to the path length, its reduction just shows that the spectral radius is less dominant with respect to the other eigenvalues [33] and it becomes another measure of modularity like the other two spectral indices. The results of the multivariate model further highlight the shift from a rich-club network at younger age to a more random one at full-term age. In particular, the lowest on the test set is around 2.11 weeks compared to other studies [15, 39]. Finally, it is important to notice that those negative trends for graph features are consistent for each effective connectivity method and each sleep state. At first sight, the results we obtained seem to be in contradiction with maturation trends that can be found in children or adolescents, where a shift from random to a rich-club network has been discovered. However, it should be taken in account that, on one side, only 9 electrodes on the scalp were used, due to the small size of the preterm brain and, on the other hand, there is a fast development of the brain during this monitoring period with different trajectories for the different cerebral regions [40]. This composite evolution is mainly driven by the different disappearance timing of cortical subplate in the various brain areas [40]. In particular, two main changes took place. The first one is the faster development of the thalamocortical connections compared to the corticocortical ones [5]. The shift from one type of connections to the other happens only at later age [41]; therefore, the maturation process might lead to “separation” of the nodes, as reflected by the results in this study. In simple terms, the EEG-scalp connections became weaker in favour of a connection strengthening between the cortical and subcortical areas. The second important change is the negative correlation between age and short-range corticocortical connections, as shown by [5] via fMRI. In [5], the study results showed that long distance connections develop faster than the short ones. Furthermore, the development is characterized by a strengthening of the former connections and weakening of the short-range couplings [1]. Consequently, the EEG electrodes/nodes (which measure short-range connectivity) tend to separate each other, with an increase of the path length and a reduction of the clustering coefficient. This hypothesis is also supported by the decrease in fronto-frontal and occipito-occipital functional coupling measured by fMRI [5]. This segregation can be emphasized by the fact that there are a few electrodes on infant scalp. However, a study with high-density EEG on preterm infants [13] found an increased modularity on the scalp EEG network and a reduced clustering coefficient in the postcentral network, while the clustering coefficient increases in the precentral network. This result could confirm the segregation or the more local integration of the brain network due to the pruning of short-range connections, as also shown by [42] in the comparison between adults and children. It is important to notice that both models with and without CCA found the same trends, but source filtering increased the prediction power of the model. It is possible that the EMG artifacts disturbed the connectivity analysis and biased the prediction model in the first considered scenario.
5. Conclusions
In the present study, we investigated effective EEG-based brain connectivity in premature infants, whose PMA ranged from 27 to 42 weeks. Results showed that the EEG-graphs changed with age in terms of topology. In particular, the clustering coefficient and the spectral radius decreased with maturation, while the path length increased. This perspective suggests that the EEG graph shifted from a small-world network to a random network. This apparent nodes’ segregation can be a consequence of the thalamocortical connections development and the strengthening of the long-range cortical connections. The lowest age-prediction error was 2.11 PMA weeks (obtained with GC in QS), which is in line with literature results. Application of source filtering methods, like CCA, can improve the performance of the connectivity analysis.
Disclosure
This paper reflects only the authors’ views and the Union is not liable for any use that may be made of the contained information.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper. The received funding stated in the Acknowledgments does not lead to any conflicts of interest.
Acknowledgments
This research is supported by Bijzonder Onderzoeksfonds KU Leuven (BOF): The Effect of Perinatal Stress on the Later Outcome in Preterm Babies (no. C24/15/036); iMinds Medical Information Technologies (SBO, 2016); Belgian Federal Science Policy Office, IUAP no. P7/19/(DYSCO, “Dynamical Systems, Control and Optimization,” 2012–2017); Belgian Foreign Affairs-Development Cooperation (VLIR UOS Programs (2013–2019)); and EU: the research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007–2013)/ERC Advanced Grant: BIOTENSORS (no. 339804). A. Caicedo is a post doc fellow at Fonds voor Wetenschappelijk Onderzoek-Vlaanderen (FWO), supported by Flemish government. M. Lavanga is a SB Ph.D. fellow at Fonds voor Wetenschappelijk Onderzoek-Vlaanderen (FWO), supported by Flemish government.