The spontaneous brain activity exhibits long-range spatial correlations detected using functional magnetic resonance imaging (fMRI) signals in newborns when (1) long neuronal pathways are still developing, and (2) the electrical brain activity consists of developmentally unique, intermittent events believed to guide activity-dependent brain wiring. We studied this spontaneous electrical brain activity using multichannel electroencephalography (EEG) of premature and fullterm babies during sleep to assess the development of spatial integration during last months of gestation. Correlations of frequency-specific amplitudes were found to follow a robust bimodality: During low amplitudes (low mode), brain activity exhibited very weak spatial correlations. In contrast, the developmentally essential high-amplitude events (high mode) showed strong spatial correlations. There were no clear spatial patterns in the early preterm, but clear frontal and parieto-occipital modules at term age. A significant fronto-occipital gradient was also seen in the development of the graph measure clustering coefficient. Strikingly, no bimodality was found in the fMRI recordings of the fullterm babies, suggesting that early EEG activity and fMRI signal reflect different mechanisms of spatial coordination. The results are compatible with the idea that early developing human brain exhibits intermittent long-range spatial connections that likely provide the endogenous guidance for early activity-dependent development of brain networks.
The patterns and mechanisms of spontaneous brain activity, often referred to as resting-state or intrinsic brain activity, have recently gained a lot of interest. Much of the work in this field has been based on studying the slow, spatially correlated fluctuations in the blood oxygen level-dependent (BOLD) signals in the functional magnetic resonance imaging (fMRI) signals (Biswal et al. 1995), while subjects are apparently resting (Fox and Raichle 2007; Zhang and Raichle 2010). It is now known that the brain is hierarchically organized into large-scale resting-state networks (RSNs) across functionally integrated, anatomically remote brain regions (Greicius et al. 2003; Fransson 2005). Prior studies have also shown that the fMRI signals recorded from sleeping human newborns, including the early premature infants, exhibit interhemispherically spanning RSNs that are in many ways comparable with those seen in the older subjects (Fransson et al. 2007; Doria et al. 2010; Smyser et al. 2010). According to the assumption of a causal link from neuronal activity to fMRI signal (Logothetis 2008), the networks observed in resting-state fMRI (rs-fMRI) recordings are commonly considered to reflect spatially defined networks of neuronal activity. In support of this idea, several recent studies in adults have shown how the power of electroencephalography (EEG) or magnetoencephalography (MEG) signals at specific frequency bands exhibits robust spatial correlations with resemblance to the RSNs described from rs-fMRI signals (Mantini et al. 2007; Pasquale et al. 2010; Brookes et al. 2011; Hipp et al. 2012).
In this context, it is intriguing that recent work in developmental neurophysiology has provided convincing data on how early neuronal activity is fundamentally different from the one seen in the more mature brain (Blankenship and Feller 2010; Hanganu-Opatz 2010; Colonnese and Khazipov 2012). The electrical activity of early brain networks is characterized by 2 alternating modes of activity: The self-organizing (Kilb et al. 2011), locally generated spontaneous activity transients (also known as [a.k.a.] SAT or bursts in human infants; Vanhatalo, Palva, et al. 2005; Jennekens et al. 2012) and the periods of relative silence between them (Blankenship and Feller 2010; Kilb et al. 2011). This type of intermittent network activity is believed to provide the crucial endogenous driver for brain wiring at the developmental period (mainly last trimester in pregnancy in humans) before the cortical networks become modulated by the exogenous stimuli, the sensory input (Hanganu-Opatz 2010; Colonnese and Khazipov 2012). The immaturity of anatomical brain pathways (Kostović and Judas 2010) and the characteristic cellular level network mechanisms underlying early brain activity are such that they have been considered to be unable to sustain long-range functional connections seen in RSNs using fMRI signals (Colonnese and Khazipov 2012).
Our present study was set out to examine the striking discrepancy in the current literature where the cellular level mechanisms of early electrical brain activity support the notion of an absence of RSNs in the newborn brain on one hand, while recent fMRI studies report a presence of large-scale RSNs in the hemodynamic activity in human newborns on the other hand. In the following, we will use the term network as a functional concept referring to a set of brain regions (fMRI) or EEG signals that display some degree of functional connectivity as measured with the given mathematical paradigms. We will use the more specific term electric resting-state network (eRSN) referring to the measures of functional connectivity, or spatial amplitude correlations, between the EEG-derived signals presented in this work.
Materials and Methods
For additional details of and schematic illustration of signal analysis procedures, see Supplementary Material.
EEG Data, Subjects and Data Acquisition and Preparation
Eleven healthy human preterm babies (conceptional age 28.3–33.8 weeks) and 10 healthy human fullterm babies (conceptional age 37.6–43.5 weeks) were recruited. The preterm babies were recruited from a larger, ongoing research project that develops EEG methodology and examines the central nervous system development in this age group. These babies were selected on the basis that they had normal serial brain ultrasound examinations, they did not have any diagnosed abnormalities or chronic illnesses, no acute inflammations, as well as no ongoing drug treatments, that would be known to potentially confound the EEG activity. The EEG was recorded either in the neonatal intensive care unit (preterm babies) or in the department of children's clinical neurophysiology (fullterm babies), Children's Hospital, Helsinki University Central Hospital, Finland. The EEG signals were collected during sleep at 256 Hz using a NicOne (Cardinal Healthcare/Natus, USA) or Cognitrace (ANT-Neuro, Germany) EEG systems and EEG caps (sintered Ag/AgCl electrodes; Waveguard, ANT-Neuro), which included 28–64 channels with positioning according to the international 10-20 standard. For further details of the recording multichannel newborn EEG, see Vanhatalo et al. (2008) and Stjerna et al. (2012) and http://www.nemo-europe.com/en/educational-tools.php/. A 5-min artifact free segment during quiet sleep (i.e., background activity trace alternant or trace discontinue; see Fransson et al. 2012 for more details) was selected from each subject by an EEG expert (S.V.). Each segment underwent a preprocessing stage consisting of bandpass filtering at 0.5–20 Hz, offline remontaging into a reference-free Laplacian current source density (CSD) montage with 25 channels (minimal smoothing, γ = 10−8), and finally exporting into European data format for further processing with custom-scripted software in the MATLAB environment. The study was approved by the Ethics Committee of the Hospital for Children and Adolescents, Helsinki University Central Hospital.
Preprocessing of the EEG Signals
Since reliable source reconstruction models do not yet exist for neonates, we decided to work at the signal space after transforming the EEG data into a Laplacian CSD montage, which was done in BESA software based on the spherical spline interpolation of scalp potentials with minimal smoothing (a.k.a. reference-free montage). Minimal smoothing was supported by our recent demonstration of high spatial specificity of neonatal scalp EEG signals (Odabaee et al. 2013). To further alleviate spatial smearing of the scalp-recorded EEG signals by volume conduction, we decided to implement the recently introduced method of orthogonalization (Hipp et al. 2012). This technique takes advantage of the orthogonal projection of a signal on another to alleviate the common instantaneous amplitude components (see Supplementary Material for further details).
Since higher frequencies in the neonatal EEG signal are mainly nested in the multifrequency SAT (or bursts) events (Vanhatalo, Palva, et al. 2005), the fluctuations at higher frequency band amplitudes do indirectly depict the occurrence of SAT events (Tokariev et al. 2012). We extracted these event-level activities using Hilbert transform (see Supplementary Material for further details). The infraslow components are an inherent property of SATs and can be recorded using Full-band EEG methodology (FbEEG; Vanhatalo et al. 2005a); however, they are shown to link via the tight nesting mechanism (a.k.a. phase-amplitude coupling) to the higher frequency oscillations (Vanhatalo, Palva, et al. 2005). This coupling makes it possible to use higher frequency amplitude envelope as the temporally more accurate proxy marker of SAT events, which is useful in spatial analyses of the present kind that would be easily confounded by any unavoidable artifacts in the infraslow signal components (for further details, see Vanhatalo et al. 2005b).
fMRI Data, Subjects and Data Acquisition and Preparation
MRI data were acquired at 1.5 T (Philips Intera, 6-channel receive only coil) in 18 newborns that were scanned during natural sleep. Further, all subjects were born via a planned cesarean section at full term with a normal Apgar score (Fransson et al. 2009). No visible sign of brain abnormality could be observed from the anatomical MR scans, and detailed information regarding age, age at scanning, and head circumference is given in Fransson et al. (2009). Written parental consent was obtained prior to MR examination, and the study was approved by the local ethics committee in Stockholm. BOLD-sensitive echo-planar images of the newborn brain (time repetition/time echo/flip = 2000 ms/50 ms/80°, matrix size = 64 × 64, field-of-view = 180 × 180 mm, 20 axial slices, slice thickness = 4.5 mm, spatial resolution 2.8 × 2.8 × 4.5 mm3) were acquired during 10 min (300 image volumes).
Preprocessing of fMRI Data
All fMRI data preprocessing was carried out using the SPM5 (Wellcome Department of Imaging Neuroscience, London, UK) software package running in MATLAB and in-house developed routines. rs-fMRI volumes were realigned and normalized to an MR neonatal image template (Kazemi et al. 2007). We a priori set the maximum allowable amount of estimated movement during fMRI scanning not to exceed 1.0 mm (translational) and 1.5° (rotational movement), which resulted in that only 200 contiguous echo planar imaging (EPI) volumes for each newborn were used for further analysis. The normalized EPI volumes were resampled to 3 × 3 × 3 mm3 and spatially smoothed using a Gaussian filter with the full-width at half-maximum of 6 mm. Subsequently, rs-fMRI networks in the newborn brain were sampled by defining 39 regions of interest (ROIs) in brain regions that span auditory, visual, attention, default, saliency, and subcortical networks. The fMRI signal intensity time course was extracted from each ROI (for a complete list of networks and coordinates for all ROIs see Supplementary Table 1 in Fransson et al. 2012) defined as a sphere with a radius of 6 mm. Lastly, spurious nonneuronal signal sources were removed from the fMRI signal intensity time courses by linear regression. Covariates of no interest included the global signal averaged over the whole brain, signals averaged from the cerebrovasulcar space (ROIs located in the lateral and fourth ventricles), as well as signals sampled from deep white brain matter. Additionally, 6 motion parameters were included in the regression model as well as the time-shifted version of all covariates to account for temporal derivatives of spurious signal variance.
Spatial Correlation Analysis
Spatial correlations between brain areas have been conventionally studied by computing linear correlation (Pearson correlation coefficient) between 2 signals. While this procedure requires normal distribution of time series to work properly, the band amplitude fluctuation (BAF) time series are by nature highly skewed (Fig. 1). To examine the bivariate amplitude dependencies, we decided to rescale one of the compared amplitudes by quantizing them into 20 quantiles (i.e., 5% of data points in each quantile). Re-scaled scatter plots revealed a robust amplitude relationship that was different in lower and higher amplitude ranges. These amplitude relationships were quantified by computing slopes of linear regression between the 2 signals (one signal quantized and the other signal original amplitudes) separately in the low-amplitude range (low mode: 5–50%) and the high-amplitude range (high mode: 70–95%). The lowest and highest quantiles (0–5% and 95–100%) were removed to exclude possible very high-amplitude artifacts or near zero values. As shown in Figure 2, a linear regression slope can be reasonably computed for these ranges. To test whether the slopes between real signals are significantly higher than the chance level (as could be observed between any 2 signals), we employed a multistep, hypothesis-driven statistical procedure.
We employed surrogate data to test significance of the rescaled scatter plots (Fig. 1) and proceeded to estimating the group-level significance as follows:
Individual-level analysis: Within each mode (low and high modes), we computed the regression slopes between BAF signals of channel i(CHi (t, f), i = 1, … ,CH) and the orthogonalized BAFs of all other channels to channel , where CH is the number of EEG channels or fMRI ROIs (25 and 39, respectively). In addition, we created 49 time-shifted surrogates of each orthogonalized BAF , and used them to compute 49 surrogate slope values to form a null distribution. The regression slope in the original signals was taken as significant if it was higher than 95th percentile of the null distribution. The connectivity matrices of each individual (the matrices containing the original pair-wise slope values at each mode; see Fig. 1) were then thresholded with this 95% significance level.
Group-level analysis: The thresholded individual connectivity matrices were converted into binary matrices and pooled together by using binomial testing at the 5% probability of success. This resulted in a group-level matrix where each cell shows the binomial probability of that cell being significant. We then used false detection rate (FDR) correction at 1% level to control for the multiple comparison procedure [CH × (CHi− 1)/2 simultaneous binomial tests], which gave us group-level binary mask showing that signal pairs are significantly correlated at the group level (see Fig. 1C for the schematic block diagram of the procedure). This mask was used to threshold the grand average of individual connectivity matrices. The comparison of connectivity matrices between groups (modes or ages; see Figs 3 and 4) was performed by using pair-wise comparison (Wilcoxon), with FDR correction at 1% level leading to another binary matrix of significant differences.
Statistical Testing with Surrogates
To evaluate the significances through hypotheses testing, we generated surrogate data that share statistical characteristics with the original data. Several methods have been previously used to generate surrogate datasets, and the main challenge in our situation is to minimize the chance of false positives coming from using surrogates that are too different from the original signals. To maximally preserve the temporal structure in each surrogate realization, we generated 49 surrogates obtained by time shifting of one signal in the pair by increments of 6 s (i.e., 1:50 of the whole 300-s signal). The schematic procedure of generating surrogates is presented in Supplementary Figure 1. For each electrode pair, we then calculated the connectivity measures between the first signal and all time-shifted versions of the second signal, which yielded a null distribution for statistical testing. We used the Mantel test for assessing global similarity of connectivity matrices between the 2 age groups (preterm vs. fullterm) and functional modes (low vs. high). See Supplementary Material for further details.
Graph Theoretical Measures of eRSN Networks
During the course of our study, it became apparent that eRSNs in the neonatal brain do not form such spatially distinct and statistically independent constellations as has been reported in previous fMRI studies. While this is compatible with the idea of spatially overlapping, temporally mixing and variable network configurations now well established in the electric activity of adult brains (Stam et al. 2009; Stam and van Straaten 2012), it also calls for employing metrics suitable to capture such dynamics. To this end, we used graph theoretical metrics that were recently developed for quantifying features of complex networks, such as those measured by EEG (Bullmore and Sporns 2009; Stam et al. 2009; Stam and van Straaten 2012). In computing the weighted graph measures, we took each EEG time series as the nodes, and the values of linear regression slopes were taken as the link strength. This approach has been recently employed in a number of studies in older children and adults (Stam et al. 2009; Boersma et al. 2011; Tsiaras et al. 2011; Hipp et al. 2012; Smit et al. 2012; Boersma et al. 2013). While the use of a raw EEG signal as the node has been recently challenged (cf. Rubinov and Sporns 2010) due to spatial smearing of at least adult scalp EEG, it is notable that our recent work on spatial characteristics of neonatal EEG indicated a high spatial specificity of neonatal scalp EEG (Odabaee et al. 2013).
For each node, we calculated the weighted clustering coefficient (CC) that expresses the weighted likelihood that neighbors of a given node are connected. Averaging CC over all nodes (i.e., EEG signals) resulted in an average, whole brain, weighted CC that reflects the level of local organization of the network. We also examined average weighted CC separately from the pre- and postcentral brain regions to see potential network differences between these otherwise differentially developing brain regions (cf. Fransson et al. 2012; Kostović et al. 2012). Global network communication capacity was assessed by computing the commonly examined weighted characteristic path length (PLw), that is, a path that minimizes the sum of the inverse link weights. Finally, we assessed network modularity, a key feature in early network development, by computing Optimal Community Structure and Modularity (OCSM; Newman 2006; Rubinov and Sporns 2010), which quantifies the degree to which the network under investigation may be subdivided into clearly delineated subnetworks. Finally, we also assessed whether small world-like properties (S = CC/PL; see Boersma et al. 2011) may be seen to emerge at this early age. All graph analyses were computed using openly available Brain Connectivity Toolbox (Rubinov and Sporns 2010).
Before computing group-level statistics or comparisons between groups, all graph measures were normalized at an individual level by dividing the absolute graph value by the mean of 50 values obtained from random shuffling of the given network (i.e., 50 surrogates; see also Boersma et al. 2011). Notably, normalization of graph measures made it possible to compare networks between individuals, between modes and across age groups, the latter of which would have been compromised by the substantial developmental differences in the gross cranial anatomy of preterm and term babies.
BAF Correlations Depend on Instantaneous Amplitudes
Following the design of adult EEG studies (e.g., Mantini et al. 2007; Brookes et al. 2011), we first examined BAFs at distinct frequency bands as obtained by bandpass filtering followed by Hilbert transformation. This was done for BAFs within a broad frequency band (3–15 Hz), as well as by using narrower frequency bands (3–8 and 8–15 Hz) that were previously shown to represent relatively independent oscillations within SAT activity (cf. Tokariev et al. 2012), the main component of high amplitudes in our analysis. Frequency spectra of BAF signals (see Supplementary Methods and Tokariev et al. 2012) suggested that it can be down-sampled at 4 Hz while still capturing its full time-domain behavior.
The common method in analyzing functional connectivity within or between different RSNs is based on computing pair-wise linear (Pearson product–moment) correlation coefficients between 2 fMRI time series or EEG/MEG power envelopes (a.k.a. seed-ROI based correlation analysis), followed by group statistics using Fisher r-to-z-transformation (Dijk et al. 2010). This procedure builds on the assumption that the samples are comparable, independent and belong to a bivariate normal distribution. None of these is the case with neonatal BAF signals where distributions are highly asymmetric with long tails toward higher amplitudes (Fig. 1B), resembling the Rayleigh distributions typically seen in power time series. Formal testing with Kolmogorov–Smirnov also confirmed nonnormality of BAF distributions (K–S 0.51, P < 0.01), while the rs-fMRI time series exhibited reasonably normal distribution (K–S 0.13; P = 0.09). The heavy clustering of low-amplitude values seen in the EEG-derived BAF scatter plots (Fig. 1B) implies an unavoidable bias in linear correlation analyses. A straightforward technical rejection of high-amplitude values as “outliers” was precluded by the prior neurophysiological knowledge that the high values represent genuine and developmentally crucial brain activity, the SAT events (Vanhatalo, Palva, et al. 2005; Tokariev et al. 2012). We therefore tested the possibility that spatial relationships between BAF signals might be nonlinear across the amplitude scale. To test this hypothesis, we quantized the amplitude in each BAF signal pair into 20 quantiles (5% cumulative percentiles each). Inspection of such rescaled scatter plots (Fig. 1B) suggests that the linear relationship between BAF amplitudes is negligible at the lower amplitude range (hereafter called “low mode”: percentiles 5–50%), while there is a steep relationship at the upper amplitude range (hereafter called “high mode”; percentiles 70–95%). In full agreement with the prior neurophysiological literature (Vanhatalo and Kaila 2006; Blankenship and Feller 2010; Hanganu-Opatz 2010; Colonnese and Khazipov 2012), this finding implies that functional network connectivity in the newborn brain is expressed as 2 different modes of neuronal communication with significantly different degree of amplitude dependance. A comparison of all available BAF signal pairs (25 channels, i.e., 300 comparisons in each baby) demonstrated that bimodality in spatial correlations of this kind seemed to be a nearly global feature (Figs 3 and 4). Strikingly, such bimodality was not observed in the fMRI data, which showed a remarkably monotonic, linear correlation across all amplitudes (Fig. 1B) and between all rs-fMRI time series (Fig. 3).
To test whether the observed correlations are statistically significant, we employed a conservative testing approach using surrogate signals generated by repeated (n = 49) time shifting of the same data. This procedure allows testing the null hypothesis that any BAF or fMRI signal with the given temporal behavior would give rise to comparable relationships between amplitudes, and group-level testing using binomial statistics yielded connectivity matrices where each pair-wise relationship could be tested individually.
Comparison of Low and High Modes
Our results revealed a striking difference between low and high modes of EEG-based connectivity (Figs 2–4). A detailed, pair-wise comparison of connectivity matrices showed that the levels of amplitude relationship were significantly higher in almost every signal pair during high mode activity (Figs 3 and 4). Global assessment of similarity between connectivity matrix patterning indicated an overall change from the low to high modes in the preterm infants (Mantel statistics −0.32; P = 0.99), whereas the overall pattern was stable across modes in the fullterm babies (Mantel statistics 0.55; P < 0.01). These observations are fully compatible with the existing idea (Blankenship and Feller 2010; Hanganu-Opatz 2010; Kilb et al. 2011; Colonnese and Khazipov 2012) that network communication in the preterm brain predominately relies on the bursts (or SATs) that are represented as the high mode in our analysis.
Development of Bimodality
The early development of network activity is thought to begin with the developmentally unique, intermittent activity patterns (i.e., the high mode) that will be gradually replaced by the continuous EEG oscillations (Vanhatalo and Kaila 2006). Conceivably, this developmental trajectory should also be reflected in the spatial EEG correlations when the fetus develops from the preterm to fullterm age. Assessing next EEG recordings from 10 fullterm newborns gave us a possibility to study developmental window that is characterized by a massive growth of long-range brain connections (Kostović and Judas 2010), hence providing the physical framework to support large-scale neuronal interactions.
We found that the difference between the low and high modes is markedly reduced when the babies approach fullterm age (Fig. 3); however, one-third of the signal pairs (34%; at 3–15 Hz) still showed a significantly higher amplitude relationship during the high mode. Such developmental reduction of bimodality is fully compatible with the idea that the nearly silent periods (low mode) of early brain activity will become replaced by the more organized ongoing spontaneous activity near term age (Vanhatalo and Kaila 2006; André et al. 2010). A significant change in global patterning of connectivity was only seen in the low activity mode (Mantel test −0.14, P = 0.93; high mode 0.57, P < 0.01), which is compatible with the idea that the spontaneous brain activity during low mode (i.e., inter-SATs) undergoes a significant qualitative, spatial organization toward fullterm age (Tolonen et al. 2007; Jennekens et al. 2012).
Spatial Topology of Correlated EEG Activity
We next explored the spatial topology of the pair-wise signal correlations in more detail using a graph theoretical approach (Rubinov and Sporns 2010; Telesford et al. 2011). This aimed to characterize differences in spatial network organization between low and high activity modes, and how they develop during the last months of pregnancy. To this end, the nodes in the network graphs were represented by the electrodes, and the links were defined by the strength association between the slope of pair-wise linear relationships. The graph metrics were rendered comparable between age groups by normalizing them using surrogate graphs, and we studied first the network metrics within the dominant 3- to 8-Hz frequency band (André et al. 2010; Palmu et al. 2010; Jennekens et al. 2012; Tokariev et al. 2012).
As a first measure of local spatial eRSN organization, we compared the normalized weighted CCs between age groups and found CC during high mode to be significantly increased during development (preterm: 1.02 ± 0.01 (SD) vs. fullterm: 1.04 ± 0.02; P = 0.001). In addition, fullterm babies showed a significant increase in CC when they shifted from the low to high activity modes (1.02 ± 0.01 vs. 1.04 ± 0.02, P < 0.01). The normalized characteristic path length, however, was not significantly changed during development or between activity modes. Functional modularity was measured with the metric OCSM to assess the presence of clearly delineated subnetworks within the brain (Rubinov and Sporns 2010). Despite of the substantial difference in connectivity levels between activity modes (Fig. 3), we found modularity to be comparable between the low and high modes within both age groups (preterm: Low mode 0.41 ± 0.04 vs. high mode 0.58 ± 0.10, P = 0.19; fullterm low mode 0.70 ± 0.09 vs. high mode 0.84 ± 0.07, P = 0.29). However, there was a significant developmental increase in modularity in both activity modes (low mode: Preterm 0.41 ± 0.04 vs. fullterm 0.70 ± 0.09, P = 0.01; high mode: 0.58 ± 0.10 vs. 0.84 ± 0.07, P = 0.04, Wilcoxon). Finally, we searched for the evidence of small-worldness as a proxy to suggest a shift from random toward more ordered networks (Rubinov and Sporns 2010). While there was no developmental change, small-worldness showed a trend-level increase in the fullterm brain when shifting from the low to high activity modes (P = 0.07).
While the above analyses are all global network measures, there is ample anatomical and functional evidence indicating a developmental fronto-occipital gradient (Brockmann et al. 2011; Fransson et al. 2012; Kostović et al. 2012). Hence, we wanted to examine potential topological differences in network clustering (CC) between pre- and post-central regions, and we found that precentral CC increases while postcentral CC decreases during late fetal development in both activity modes (Fig. 5). A precentrally emphasized CC was also seen in the pair-wise comparisons (high mode: Both preterm P < 0.04 and fullterm P < 0.01; low mode: Fullterm P < 0.01; Wilcoxon test). Finally, the precentral CC was significantly increased in both age groups when they shifted from the low to high activity modes (P < 0.01, Wilcoxon test). These findings together are strikingly compatible with the notion of distinct developmental trajectories, and in particular, a protracted developmental window of frontal brain areas in humans (Kostović and Judas 2010).
eRSN at Different Oscillatory Frequencies
Prior studies in adults have shown that long-range correlations in the electric brain activity may take place at a wide range of oscillatory frequencies, but they exhibit frequency-specific spatial patterning (Mantini et al. 2007; Palva et al. 2010; Pasquale et al. 2010; Brookes et al. 2011; Hipp et al. 2012). Prior work on neonatal EEG has reported spatial differences and developmental changes (Tolonen et al. 2007; André et al. 2010; Jennekens et al. 2012) in the frequency composition of early bursts/SATs. Therefore, we examined whether eRSN could be frequency specific by repeating all-above analyses after first filtering the EEG signal into 2 frequency bands (3–8 and 8–15 Hz) that cover functionally independent oscillatory entities within the SATs (André et al. 2010; Tokariev et al. 2012). We found that a functional bimodality can be observed at all frequencies in a qualitatively comparable manner (Figs 3 and 4).
A detailed network analysis with graph metrics, however, revealed frequency-related differences. In contrast to the 3- to 8-Hz frequency band, the higher 8- to 15-Hz frequency band showed no differences in the mean clustering (CC) between activity modes or between brain areas. Compared with the 3- to 8-Hz frequency band, however, modularity at 8–15 Hz was more dynamic, and it was found to increase from the low to high modes in both age groups (fullterm 0.67 ± 0.08 vs. 1.41 ± 0.13, P < 0.01; preterm 0.34 ± 0.05 vs. 0.85 ± 0.15, P < 0.03). There was also a significant developmental increase in modularity at both low and high activity modes (P = 0.01 and 0.02, respectively). These shifts from relatively random to modular structures consisting of anterior and posterior network entities can be clearly observed in the weighted graph images in Figure 5.
Absence of Bimodality in the Newborn rs-fMRI data
Earlier studies have reported consistent and robust RSN in newborns using rs-fMRI (Doria et al. 2010; Smyser et al. 2010; Fransson et al. 2012), which as a method implicitly assumes a linear correlation between 2 signals. In our present analysis procedure, such linear correlation over the whole amplitude range means no difference between higher and lower amplitudes, that is, “unimodality” in the amplitude relationship (see also the rescaled scatter plots in Fig. 1B). Since we are not aware of prior studies exploring how spatial fMRI correlations depend on instantaneous amplitudes, it is possible that bimodality of this kind has only been ignored due to the widely used analysis procedures. In light of our EEG observations above, we hence wanted to see if the network correlations in the newborn fMRI time series are genuinely linear, or could they potentially exhibit such bimodality that was seen in the electric brain activity. As expected from the inspection of scatter plots of raw fMRI time series, the rescaled fMRI time series showed a clear linear relationship and no evidence of inflection of the regression line at the upper amplitude range akin to that seen with BAF signals (Fig. 1B). Moreover, a statistical comparison between lower and higher amplitude quantiles (5–50% and 70–95%) confirmed the visual observation: There were very few (0.9%; 7 of 741 comparisons) pair-wise differences in regression slopes between lower and higher amplitude ranges in each signal pair, constituting a strong case against bimodality in the amplitude relationships between newborn rs-fMRI time series.
The current findings demonstrate the presence and development of spatially large-scale eRSN in the early electric activity of the newborn brain during a time window of the active growth of long-range anatomical connections (Dudink et al. 2008; Huang et al. 2009; Kostović and Judas 2010). Most importantly, these networks exhibit a robust temporal bimodality with alternating episodes of low and high degrees of correlations. The dichotomic behavior in temporal domain declines toward term age when the brain activity is known to develop into more stable, continuous oscillations (Vanhatalo and Kaila 2006). Our findings suggest 2 main implications: First, the observed bimodality in spatial connectivity supports the notion that the early developmental network events guide the development of long-range connectivity in an activity-dependent manner. Secondly, the absence of comparable bimodality in the rs-fMRI time series suggests at least partly distinct underlying mechanisms for the early networks observed using EEG and fMRI modalities.
Early Brain Network Activity
All structures of the immature nervous system are known to produce endogenous, self-organizing network activity, which is characterized by alternation between periods of high activity interrupted by relative inactivity (Vanhatalo and Kaila 2006; Blankenship and Feller 2010; Hanganu-Opatz 2010; Colonnese and Khazipov 2012). This network behavior is driven by an interplay of early depolarizing (or even excitatory; Ben-Ari et al. 2012) γ-aminobutyric acid (GABA)ergic interneuronal transmission, and synaptic as well as extrasynaptic network mechanisms located within the cortex and the developmentally transient subplate zone (Blankenship and Feller 2010; Kilb et al. 2011). Early network events generated by neuronal networks themselves or triggered by activity in the sensory systems are considered to be of crucial importance for the activity-dependent survival of young neurons (Hanganu-Opatz 2010; Kilb et al. 2011), and they are believed to provide the essential, activity-dependent guidance for early brain wiring that in human primarily takes place in utero (Hanganu-Opatz 2010; Kilb et al. 2011).
Our present work is fully compatible with the idea that brain establishes long-range connectivity during brief events of endogenous network activity (high modes). The apparent random connectivity in the early preterm baby (Fig. 5) is in full agreement with the absence of direct anatomical connectivity between remote cortical locations at this age, whereas the observed increase in modularity, or spatial organization toward fullterm age, coincides nicely with the histological growth of long-range connections (Kostović and Judas 2010). It would be attractive to assign any of the findings to a special pathway structure, such as callosal, thalamo-cortical, or subplate connections. However, the prolonged, spatially varying and overlapping developmental time windows of all these pathways during the age span from our preterm to fullterm babies (Kostović and Judaš 2006, 2010; Petanjek et al. 2011; Kostović et al. 2012) make it markedly ambiguous to speculate on relative contributions of each individual network structure in the spontaneous EEG activity. Obtaining such spatially and structurally more specific information may be possible, however, by exploring task-dependent activities (e.g., Vanhatalo and Lauronen 2006; Milh et al. 2007; Vanhatalo et al. 2009; Stjerna et al. 2012; Mahmoudzadeh et al. 2013). Intriguingly in this context, recent modeling studies have shown how modular network structure emerges optimally when anatomical growth is guided by functional network communication (Stam et al. 2010), which was observed in our study during high mode activity in the human infants.
Due to challenges in recording dense array EEG in early preterms, there are no comparable prior studies on early electric networking in humans. However, several recent works using rodent models are strikingly compatible with our present findings. Imaging propagation patterns of neuronal activity in whole-brain slices has shown that the earliest network events spread widely and with low spatial specificity across the cortical anlage (Lischalk et al. 2009; Conhaim et al. 2011), being strikingly compatible with our observation of apparently random global networking during high activity mode in the early preterm infants (Fig. 3). These animal experiments are also indirectly challenging the idea of neuronally driven, temporally stable, and spatially distinct networks that are reported from the early human rs-fMRI (Doria et al. 2010; Smyser et al. 2010). Moreover, early network events in young rat pups develop differentially in posterior and frontal areas, and frontal areas were recently shown to exhibit a robust networking at developmental stages that compare with human fullterm babies (Brockmann et al. 2011). These findings provide an intriguing translational parallel and a tentative cellular level mechanistic explanation to our observations of the fronto-occipital gradient in the developing network clustering, as well as the developmental emergence of frontally prominent networks in the fullterm babies.
Difference Between Electric Activity and Blood Oxygenation
Clearly, cortical activity sampled in EEG and fMRI recordings differs not only at the spatial and temporal resolution achieved, but also in the neurophysiological origin of the signal measured. Most previous comparisons between EEG and fMRI signals have investigated relationships of fMRI signals to different spectral components of the EEG (e.g., Tagliazucchi et al. 2012) or the local field potential activity (e.g., Magri et al. 2012). There is, however, very little information available on the relationships of underlying neuronal mechanisms associated with these modalities at earliest ages (cf. Colonnese et al. 2008). The robust bimodality in the electric activity and the virtual absence of such behavior in the blood oxygenation signals seems apparently conflicting in the context of the widely assumed causal mechanism of electrical brain activity, being the driving mechanism for fMRI signals in the brain (Logothetis 2008). One possible technical caveat in the direct comparison of EEG and fMRI signals is the markedly lower sampling rate of fMRI, which in theory could result in ignorance of shorter brain events. To the best of our understanding, however, the dynamics in the rs-fMRI (Niazy et al. 2011) and in the preterm-evoked fMRI responses (cf. Arichi et al. 2012) are slow enough to be captured with our rs-fMRI sampling frequency. The duration of high activity epochs (Fig. 2) is also long enough to trigger fMRI responses (Hu et al. 1997; Janz et al. 1997; Fransson et al. 1999) that should such take place during early endogenous network events. Since fMRI responses are known to rapidly deteriorate with decreasing interstimulus intervals (Boynton et al. 1996; Fransson et al. 1998), the most likely explanation for the seeming discrepancy between our eRSN and rs-fMRI findings may be that individual responses in the oxygenation signal are flattened by the absence of sufficient recovery time between activity bouts. This may be especially true in the newborn brain, which is shown to have a remarkably slow and variable hemodynamic response function after neuronal activation (Colonnese et al. 2008; Kusaka et al. 2011; Arichi et al. 2012). The differential effect of inhibitory versus excitatory neurotransmission on fMRI signals (Buzsáki et al. 2007) adds further complexity to reconciling fMRI and EEG activities in neonates, because (1) the early network activity is based on excitatory GABAergic transmission (Ben-Ari et al. 2012), and (2) the basal metabolic rate in the associated developing neuronal structures may be high enough to mask any metabolic demands set by intermittent synchronized network events. Due to the fact that rs-fMRI data from preterms were not included in this study, we cannot exclude the possibility that the hemodynamic signal in preterms would show a bimodal behavior. However, we think that this scenario is unlikely given that previous research in rodents (Colonnese et al. 2008) has shown how BOLD response exhibits a systematic decline in latency and an increase in amplitude during the early development. Such developmental transition from slower and smaller to faster and larger BOLD responses speaks against the idea that bimodality could be present in the preterm rs-fMRI data without being seen in the fullterm data. Taken these together, we find it plausible to suggest that essential dynamics of early network activity in the human brain may not be directly reflected in the rs-fMRI signals. Indeed, this conclusion and our findings in general are fully compatible with the recently proposed framework, where resting-state spatial correlations in the fMRI activity are viewed in the context of serving developmental circuit formation rather than reflecting the presence of higher cognitive operations (Colonnese and Khazipov 2012).
Transient States in the Brain Activity
Several recent studies on adults have reported distinct, temporally delineated brain states in the context of RSN. The fMRI amplitude fluctuations were reported to follow bimodal networking (Petridou et al. 2013) bearing resemblance with our EEG observations; however, those temporal dynamics were an order of magnitude slower than what we found in the newborn EEG. Recent studies on adult brain activity have demonstrated multistability in distinct network oscillations (Heitmann et al. 2012), and spatial networks driven by rapid “microstates” (Britz et al. 2010; Musso et al. 2010; Yuan et al. 2012), at time scales that are considerably shorter than the high mode epochs in our newborns. While the idea of temporally limited functional modes/states is thus gaining wider attention in neuroscience, the unique cellular networks and transmitter mechanisms of early developing brain call for a need to interpret newborn findings in their own specific context. Speculatively, this implies that the mechanisms driving transience in developing brain activity are searched from cellular structures and molecular mechanisms that are unique to the early developmental period.
The present observations open a window to extending the assessment of early brain function from the conventional phenomenological EEG reading (André et al. 2010) to a quantifiable assessment of global network integrity. A large number of clinical conditions, ranging from early prematurity to intraventricular bleeding, stroke, and asphyxia, are known to associate with later neurocognitive disabilities (Mwaniki et al. 2012), and the ontogenetic mechanisms of those disabilities entail early disordered network function that readily escapes capture by neuroimaging. Further longitudinal studies may establish links from early brain insult to disorganized eRSN mechanisms as well as further to compromised neurocognitive development. Harnessing such network measures for clinical work (Fox and Greicius 2010) may open novel avenues for constructing early functional biomarkers (cf. Ment et al. 2009) with promise for both guiding early individualized therapies and for greatly expediting studies on early clinical interventions.
A.O., P.F., and S.V. designed and performed research, analysed data and wrote the paper. S.V., P.F., and M.M. collected the data.
This work was supported by Helsinki University Hospital and Juselius Foundation. S.V. received support from the European Community's Seventh Framework Programme (FP7-PEOPLE-2009-IOF, grant agreement no. 254235), P.F. was supported by the Swedish Research Agency, and A.O. received financial support from the Perinatal Research Centre of Prof. Paul Colditz, UQCCR.
We thank Mr Anton Tokariev for providing us with the 3D scalp image for eRSN visualization. In addition, we would like to thank Drs Ulrika Åden, Mats Blennow, and Hugo Lagercrantz for their continuous interest and support. Conflict of Interest: None declared.