Identifying Emergent Mesoscopic–Macroscopic Functional Brain Network Dynamics in Infants at Term-Equivalent Age with Electric Source Neuroimaging

Abstract

Aim:

To identify and characterize the functional brain networks at the time when the brain is yet to develop higher order functions in term-born and preterm infants at term-equivalent age.

Introduction:

Although functional magnetic resonance imaging (fMRI) data have revealed the existence of spatially structured resting-state brain activity in infants, the temporal information of fMRI data limits the characterization of fast timescale brain oscillations. In this study, we use infants' high-density electroencephalography (EEG) to characterize spatiotemporal and spectral functional organizations of brain network dynamics.

Methods:

We used source-reconstructed EEG and graph theoretical analyses in 100 infants (84 preterm, 16 term born) to identify the rich-club topological organization, temporal dynamics, and spectral fingerprints of dynamic functional brain networks.

Results:

Five dynamic functional brain networks are identified, which have rich-club topological organizations, distinctive spectral fingerprints (in the delta and low-alpha frequency), and scale-invariant temporal dynamics (<0.1 Hz): The default mode, primary sensory-limbic system, thalamo-frontal, thalamo-sensorimotor, and visual-limbic system. The temporal dynamics of these networks are correlated in a hierarchically leading–following organization, showing that infant brain networks arise from long-range synchronization of band-limited cortical oscillation based on interacting fast- and slow-coherent cortical oscillations.

Conclusion:

Dynamic functional brain networks do not solely depend on the maturation of cognitive networks; instead, the brain network dynamics exist in infants at term age well before the childhood and adulthood, and hence, it offers a quantitative measurement of neurotypical development in infants. Clinical Trial Registration Number: ACTRN12615000591550.

Impact statement

Our work offers novel functional insights into the brain network characterization in infants, providing a new functional basis for future deployable prognostication approaches.

Introduction

In humans, the development of neuronal systems in late gestation determines the organization of specialized neural functions that underpin the large-scale functional organization of the brain (Grayson and Fair, 2017; Tau and Peterson, 2009). Infants born very preterm and preterm are at high risk of neurodevelopmental disabilities in the brain higher order functions, such as learning, cognition, and motor behaviors in childhood (Marlow et al., 2005; Moore et al., 2012). Structural magnetic resonance imaging (MRI) and clinical assessments have not provided a basis for sufficiently accurate early diagnosis (Ball et al., 2015). As a measure of large-scale functional organization of resting-state (RS) brain activity, functional connectivity MRI (fcMRI) and electroencephalography (EEG) have revealed richly structured RS networks (RSNs) (Biswal et al., 1995; Mehrkanoon et al., 2014a,c), which fluctuate over time and form dynamic functional connectome (dFC) in adults (Allen et al., 2014; Boersma et al., 2011; Mehrkanoon et al., 2014a, 2016).

The EEG in infants born between ∼30 and 42 weeks of gestation has shown the emergence of low-frequency cortical oscillations (Arichi et al., 2017; Myers et al., 2012), which may exhibit an oscillatory activity of large-scale cortical correlation structure measured by scalp-EEG FC (Meijer et al., 2014; Tokariev et al., 2016). The neuronal oscillations, with their specific frequencies, open the gates of communications between cortical regions to establish patterns of local and large-scale functional brain network interactions (Hipp et al., 2012), which underpin frequency-specific dFC (Hipp et al., 2012; Mehrkanoon et al., 2014a). Thus, FC measurements based on the time and frequency components of interacting cortical oscillations may provide new insights into the evolutionary macroscopic neuronal circuits and their interaction patterns in the infant brain.

Despite numerous RS-fcMRI and EEG studies on the infant brain, to the best of our knowledge, there is no systematic description on the spatial organization, spectral fingerprints, and temporal dynamics of local and large-scale cortical activity in the infant brain. The infant scalp-EEG-based FC does not show details of the gray matter (GM)-based cortical networks (Kaminska et al., 2017; Meijer et al., 2014; Schumacher et al., 2015), due to an incomplete anatomical localization of RSNs and the absence of spatial topologies, spectral components, and temporal dynamics of RSNs at term age (Fransson et al., 2013; Mohammadi-Nejad et al., 2018; Roche-Labarbe et al., 2008; Routier et al., 2017; Tokariev et al., 2016). Because the brain and neuronal circuits are highly plastic in infancy, it is thought that the implementation of neurorehabilitation treatments is likely to minimize the emergence of neurodevelopmental disabilities (Sutton, 2013; Lee, 2018).

The EEG-based dFC may be useful both as a screening tool of developmental aberration as well to assess the effectiveness of neurorehabilitation treatments in real time. Such process has the potential to allow optimized and individualized treatments to be developed. To characterize the functional organization of macroscopic cortical activity in infants at term age, we infer source-space spatial topologies, spectral fingerprints, and temporal dynamics of RS-FCs that overcome current limitations in fcMRI and scalp-EEG topography.

Materials and Methods

Participants

The Human Research Ethics Committee of the University of Queensland and Royal Brisbane Women's Hospital approved the protocol (HREC/15/QRCH/5). Parents of participants provided informed consent. One hundred healthy infants were included in this study: 84 were preterm-born infants (born at 30–32 weeks of gestation), and studied at term-equivalent age (40–42 weeks of gestational age), and 16 were term-born infants.

High-density EEG data acquisition

We recorded high-density EEG from 100 infants in eyes-closed active-sleep condition. The active-sleep condition was defined on the basis of visual inspections of behavior, including rapid eye movements, frowns, smiles and sucking, body and limb movements, and irregular respiration coincident with low-voltage EEG. A single 40-min session was selected to acquire EEG by a 61-channel appropriately sized neonatal EEG cap arranged according to the international 10–20 system (Fig. 1A). The experimental setup included EGI HydroCel GSN 130 Geodesic Sensor Nets and an EGI Net-Amps-300 amplifier (EGI; Philips) at sampling rate of $f_{s} =$ 250 Hz. The first 5 min of artifact-free EEG data were discarded to avoid any transient brain activity (i.e., non-RS or awake condition).

FIG. 1.

Construction of large-scale brain network dynamics in infants at term age. (A) EEG acquisition, artifact correction, and beamforming-based source estimation of cortical oscillations using (B) construction of MNI space volume conduction model of the head including 90 ROIs in the gray matter volume of the infant brain. (C) Network derivation using spatial ICA of subject-level PCA-driven time–frequency functional networks, yielding (D) spatial architectures, spectral fingerprints, and temporal dynamics of the networks, and (E) z-score statistics of the group-level power spectra of cortical oscillations in the 90 ROIs. EEG, electroencephalography; ICA, independent component analysis; MNI, Montreal Neurological Institute; PCA, principal component analysis; ROI, region of interest. Color images are available online.

For each infant, we selected a 4-min epoch from the annotated active-sleep period. We also ensured the 4-min epoch did not contain high-voltage chin movement artifact. All EEG data were referenced to Cz and band-pass filtered (0.5–30 Hz). The InfoMax independent component analysis (ICA) algorithm was then used to remove cardiac, ocular, and muscular artifacts (Mehrkanoon et al., 2016, 2014a,b,c) (Fig. 1A, middle panel). The artifact-free EEG was then used for the EEG source reconstruction (Fig. 1A, right panel, see below). To ensure that developmental stage was similar in the infants born preterm and studied at term equivalent and the infants born at term, we statistically compared the group-level power spectral density (PSD) of the 4-min EEG epoch across the two groups. For each infant, we estimated the Morlet wavelet PSD. The group-level PSD was obtained by averaging the PSDs across infants in each group (see the Results section).

Construction of the volume conduction and leadfield model of the infant head

The volume conduction model of the head represents the EEG forward problem that estimates the electric potentials on the scalp for a known configuration of electric dipole sources, provided that the physical properties of the head tissues such as conductivities are known (Mehrkanoon et al., 2014c). A publicly available neonatal normalized T2-weighted MR image model that consists of the tissue priors of scalp, skull, cerebrospinal fluid (CSF), GM, and white matter of the infant brain anatomical template was used in this work (Shi et al., 2011). Scalp EEG electrodes were first co-registered with this publicly available neonatal MRI model using a nine-parameter transformation of the EGI HydroCel GSN 130 Geodesic Sensor Nets such that the electrodes were placed according to the 10–20 system (Fig. 1B) (Mehrkanoon et al., 2014c).

Given the head tissue priors, the Matlab toolbox FieldTrip was used to construct the realistic volume conduction model of the infant head by numerically solving the forward model (Dang and Ng, 2011). The resulting leadfield matrix $H$ is a mathematical operator that linearly transforms electric fields generated at N dipole sources in the brain into scalp potentials observed in EEG signals: $Y (t) = \sum_{i = 1}^{3} H_{i} S_{i} (t) + e (t), i = x, y, z$ , where $Y (t) \in ℛ^{n_{e} \times T}$ denotes the EEG signals (i.e., n_e sensors/electrodes and $T$ data points). Note that the artifact-free EEG was used in this work (see above). The $H_{i = x} \in ℛ^{n_{e} \times N}$ denotes the leadfield matrix in the x orientation, and hence, $H_{y}$ and $H_{z}$ are the leadfield matrices in the y and z orientations, respectively, $S_{i} (t) \in ℛ^{N \times T}$ is the dipole source, $e (t) \in ℛ^{n_{e} \times T}$ is the measurement noise, which is assumed to have the Gaussian distribution with zero-mean and diagonal covariance matrix $\sum_{e e}$ , that is, $e \sim N (0, \sum_{e e})$ . The conductivity values of the infant's head-component tissues were set to 0.43 $S ∕ m (s c a l p)$ , 0.2 $S ∕ m (s k u l l)$ , 1.79 $S ∕ m (C S F)$ , and 0.3 $S ∕ m (b r a i n)$ (Roche-Labarbe et al., 2008).

The electric dipole locations were constrained by the infant GM volume with respect to the predefined Automated Anatomical Labeling neonatal brain atlas that comprises 90 regions of interest (ROIs; Fig. 1B, middle panel) (Shi et al., 2011). This process resulted in $N =$ 4882 dipole locations in the 90 ROIs, each with 3 orientations in the x, y, and z axis—each dipole source represents 5 mm $^{3}$ in volume. The leadfield array $H \in ℛ^{(62 \times 4, 882) \times 3}$ for electric source imaging was used in each of the 100 infants in the subsequent analyses.

Electric source imaging of neonatal EEG

We used the linearly constrained minimum variance (LCMV) beamforming algorithm to estimate the magnitude of the nth dipole source at each orientation (Veen et al., 1997): ${\hat{S}}_{n} (t) = {[\frac{C^{- 1} H_{n}}{H_{n}^{T} C^{- 1} H_{n}}]}^{T} Y (t)$ , where $C = ℰ [Y Y^{T}]$ denotes the covariance matrix of sensor-space averaged-reference EEG signals and T is the matrix transposition. Note that the LCMV beamformer is robust to partial correlation between sources when the interdependencies of periodic cortical activity are preserved (Hipp et al., 2012). From 4882 seed voxels, we selected the seed voxels that are in the Euclidian center of each ROI in the left and right hemispheres. This process resulted in 90 seed-voxel time series that represent the large-scale cortical activity in 90 ROIs (Hipp et al., 2012; Mehrkanoon et al., 2014c).

Derivation of functional brain networks

To quantify the coherent spontaneous activity of spatially distributed cortical activity in each neonatal brain, the functional connectivity was measured by the imaginary part of the complex-valued time–frequency coherence (i.e., non-zero-lag synchronization) between all pairs of 90 ROI signals [Eq. (2)]. Taking the imaginary part of the coherence minimizes the occurrence of spurious correlation among the nearby seed voxels due to the volume conduction property (Mehrkanoon et al., 2013, 2014a; Nolte et al., 2004). For N = 90 ROIs, the maximum number of network edges is given by $M = \frac{N \times (N - 1)}{2} = 4005$ at each time point and frequency for the kth infant, yielding the functional connectivity array $γ_{k} (t, f) \in ℛ^{T \times F \times M}$ : $\begin{matrix} γ (t, f) = [\begin{matrix} Γ_{(1, 2)} (t, f) & Γ_{(1, 3)} (t, f) & Γ_{(1, 4)} (t, f) & \dots & Γ_{(1, 90)} (t, f) \\ † & Γ_{(2, 3)} (t, f) & Γ_{(2, 4)} (t, f) & \dots & Γ_{(2, 90)} (t, f) \\ † & † & Γ_{(3, 4)} (t, f) & \dots & Γ_{(3, 90)} (t, f) \\ ⋮ & ⋮ & ⋱ & ⋮ & ⋮ \\ † & † & \dots & † & Γ_{(89, 90)} (t, f) \end{matrix}], \\ t & = 1, \dots, T, \\ f & = 1, \dots, F, \end{matrix}$ (1)

where $Γ_{i, j} (t, f)$ is the imaginary part of the complex-valued time–frequency coherence between the ROI signals $x_{i} (t)$ and $x_{j} (t)$ computed by $\begin{matrix} Γ_{i j} (t, f) & = ℑ (\frac{S \{P_{i j} (t, f)\}}{\sqrt{S \{P_{i i} (t, f)\} S \{P_{j j} (t, f)\}}}) i, j = 1, \dots, 90, i \neq j, \\ P_{i j} (t, f) & = X_{i} (t, f) {\tilde{X}}_{j} (t, f), P_{i i} (t, f) = {|X_{i} (t, f)|}^{2} \\ X_{i} (t, f) & = x_{i} (t) * ψ (t, f), \\ = x_{i} (t) * {(a^{2} π)}^{- 1 ∕ 4} e x p (- \frac{t^{2}}{2 a^{2}} - i 2 π f t), \end{matrix}$ (2)

where $ℑ (.)$ denotes the imaginary operator, $P_{i j} (t, f)$ and $P_{i i} (t, f)$ [or $P_{j j} (t, f)$ ] are, respectively, the power and cross-spectral densities, $\tilde{X} (t, f)$ is the complex conjugate of $X (t, f),$ the asterisk $*$ is the convolution operator, $ψ (t, f)$ is the complex-valued Morlet wavelet, a is the wavelet scale, $i = \sqrt{- 1}$ , and f is the wavelet center frequency (Mehrkanoon et al., 2013). The scale and frequency are proportionally related by $a = \frac{6}{2 π f}$ (Mehrkanoon et al., 2013), and we chose a broad frequency spectrum ranging from 1 to 30 Hz by increment of 0.5. The smoothing operator $S {.}$ used in this study is defined by the Gaussian kernel $K (t, f) = e x p (- (\frac{t^{2}}{2 σ_{t}^{2}} + \frac{f^{2}}{2 σ_{f}^{2}}))$ , where $σ_{t} = 0.66 s e c$ and $σ_{f} = 1.32$ Hz denote the time and frequency spreads of the kernel.

The smoothing process was implemented by convolving the kernel $K (t, f)$ with the cross- and power-spectral densities to improve the reliability of the coherency estimate (Mehrkanoon et al., 2013). The symbol $†$ denotes array entries that, due to the undirected nature of Eq. (1), are symmetric across the main diagonal and these are not considered further. Since Eq. (1) is a lossless decomposition of the ROI signals, the elements in the array $γ_{k} (t, f)$ will be mutually correlated and contain considerable redundancies. That is, the spatiotemporal patterns and spectral contents of the brain FC structure are interrelated, generating folded spatiotemporal and spectral dynamics (Mehrkanoon et al., 2014a). To unfold the complex dynamics of the brain FC structure in each neonatal brain, multilinear principal component analysis (mPCA) was used throughout the eigen-decomposition of the covariance matrix of the time–frequency array $γ_{k} (t, f)$ .

Since $γ_{k} (t, f)$ has three factors of time (t), frequency (f), and network edges (M = 4005), we transformed $γ_{k} (t, f)$ into the two-dimensional array $Θ_{k} (t f, l) \in ℛ^{[T \times F] \times M}$ as (time–frequency)-by-edges. The subject-level RS networks were then obtained from the eigen-decomposition of the covariance matrix of $Θ_{k} (t f, l)$ as follows (Mehrkanoon et al., 2014a): $\begin{matrix} C_{Θ Θ_{k}} & = ℰ [Θ_{k}^{T} Θ_{k}] = V_{k} D_{k} V_{k}^{T}, \\ Y_{k} (m, t f) & = V_{{(1 : m, :)}_{k}}^{T} Θ_{k}^{T} \in ℛ^{m \times t f}, l = 1, 2, \dots, M, \end{matrix}$ (3)

where T denotes the matrix transposition, k is the subject index, $V_{k} \in ℛ^{M \times M}$ is the orthogonal matrix whose eigenvector columns (or modes) correspond to the edges of the RSNs for the kth neonatal brain, D_k is the diagonal matrix representing the eigenvalues of $C_{Θ Θ_{k}}$ such that $d i a g (D_{k}) = {d_{1_{k}} > d_{2_{k}} > \dots > d_{M_{k}}}$ shows the variance of each evident RS mode. That is, each $M \times 1$ eigenvector contains the contributions of every ROI pair combinations (or the edges of the network) to that specific eigenmode and hence reflects the spatial topology of each network. Notably, $M = 4005$ , and hence each $M \times 1$ eigenvector can be easily reshaped into a 90 $\times$ 90 FC matrix. Note that the sum of the variance explained by the first five (m = 5 of 4005) eigenvalues was 96%. Therefore, we selected the first five eigenvectors (or eigenmodes) as being sufficient to span the phase-space of RS dynamics for the subject-level network derivation [Eq. (3)].

The time–frequency principal components (PCs) associated with these five eigenmodes were then obtained from the projection of the kth original functional connectivity array $Θ_{k} (t f, l)$ through five eigenmodes [Eq. (3)], resulting in $Y_{k} (m, t f)$ . The time–frequency projection array $Y_{k} (m, t f)$ contains the dynamics of RSNs in the kth infant brain. We focus on identifying between-infants (i.e., group-level) RSNs reliability and consistency. Thus, the five eigenmodes obtained from each infant were spatially concatenated in the matrix . A second PCA was computed to extract $p = 5$ group-level PCA modes from the matrix $V$ . The group-level (or spatial) ICA of these five group-level PCA modes was computed by the second-order blind identification method to derive the independent RSNs that are commonly present in 100 infants (Belouchrani et al., 1997; Huster et al., 2015; Mehrkanoon et al., 2014a) (Fig. 1C): $\begin{matrix} C_{V V} & = ℰ [V V^{T}] = U D U^{T}, \\ Q & = U_{(1 : p, :)}^{T} V \in ℛ^{p \times M}, \\ S_{I C A} & = W Q \in ℛ^{p \times M}, \end{matrix}$ (4)

where truncated $U \in ℛ^{100 m \times p}$ denotes the weighting matrix whose first five eigenvector columns contain the weights of the contribution of each infant's brain networks (or eigenmodes) to a common group-level RSN space. The truncated $d i a g (D) = {d_{1} > \dots > d_{p}}$ is the diagonal matrix representing the explained variance (or scaling factors) of each of the five group-level PCA modes across 100 infants, $Q \in ℛ^{p \times M}$ is a group-level PCA that maps the eigenmodes of individual infant's cortical activity onto a common space, and $S_{I C A} \in ℛ^{p \times M}$ denotes the $p = 5$ independent RSNs. The unmixing matrix $W \in ℛ^{p \times p}$ linearly decomposes the orthogonal group-level PCA modes $Q$ into five independent RSNs denoted by $S_{I C A}$ . Note that each row of the matrix $S_{I C A}$ contains all possible network edges ( $M = 4005$ ), and hence, $S_{I C A}$ was reorganized as five functional connectome matrices each with the size of $90 \times 90$ ; $S_{I C A} = [s_{1}, \dots, s_{5}]$ , where $s \in ℛ^{90 \times 90} .$ The group-level weighting matrix $U$ can be shown as vertically augmented individual's weighting matrices $[u_{1}, \dots, u_{100}] \in ℛ^{100 \times (m \times p)}$ . Therefore, the time–frequency characteristics of each of the five independent RSNs were obtained from mapping the time–frequency projection (or PCs) of the same kth infant, that is, $Y_{k} (m, t f)$ , onto spatial ICA space as follows: $Y_{k}^{I C A} (p, t f) = {W u}_{k}^{T} Y_{k} (m, t f) \in ℛ^{p \times t f},$ (5)

where $Y_{k}^{I C A} (p, t f)$ denotes the time–frequency characteristics of the five independent RSNs in the kth infant in the spatial ICA space. For each infant, these five time–frequency spectra were reorganized as $y_{i, k} (t, f), i = 1, \dots, 5$ , such that $y_{i, k} (t, f) \in ℛ^{T \times F}$ . We then used the time–frequency matrix $y_{i, k} (t, f)$ for subsequent group-level time and frequency analyses of network dynamics across 100 infants.

Graph theoretical analysis

To quantitatively evaluate the topological architecture of each of the five RSNs obtained from the spatial ICA, mathematical graph theory was used. First, each $90 \times 90$ RSN was thresholded with a connection density of 17% as is typically employed in brain network analyses (van den Heuvel and Sporns, 2011). Note that the amplitude and phase of the cortical oscillation signals that were obtained from the LCMV are the optimal estimation of source signals in the weighted least-squares beamforming sense. Therefore, we statistically analyzed the functional structure of these optimally estimated source signals using surrogate network statistics to infer the RSNs and the rich-club organization. Two sets of graph connectivity measurements, primary and complex connectivity metrics, were then used to quantify the functional properties of these thresholded five RSNs (see below).

Primary connectivity metrics

The set of graph metrics computed for each of the five thresholded RSNs comprises node degree (K, the number of connections linked to each node); clustering coefficient (C, a measure of the brain's ability in functional segregation or information processing) (Watts and Strogatz, 1998); betweenness centrality (a measure of a node hubness); global efficiency ( $E_{g l o b}$ , a measure of the brain's ability in functional integration or binding the information processed by the brain regions) (Watts and Strogatz, 1998); path length (L, which is related to a measure of the brain functional integration by averaging the shortest path lengths between all pairs of nodes in the network) (Watts and Strogatz, 1998); and assortativity coefficient (a relative measure of the resilience (or elasticity) of the network against sudden attacks to the network components, such as hub nodes and edges in the network) (Achard et al., 2006).

Complex connectivity analysis

Small-world topology

We identified the small-worldness of each of the five RSNs by computing the ratio of the normalized clustering coefficient $γ = C ∕ C_{r a n d o m}$ to the normalized path length $λ = L ∕ L_{r a n d o m}$ . The $C_{r a n d o m}$ and $L_{r a n d o m}$ , respectively, denote the average clustering coefficient and average path length of the surrogate networks constructed by randomizing all connections of the original RSN under analysis (10,000 realizations) while preserving the dimension and degree distribution of the original network. In a small-world network, we expect that the ratio $λ = L ∕ L_{r a n d o m}$ tends to 1, and the ratio $γ = C ∕ C_{r a n d o m}$ to be >1, and hence, the small-world index $σ = γ ∕ λ$ should be >1 (Watts and Strogatz, 1998). This is indicative of a dense local clustering (or cliquishness) between neighboring nodes yet a short path length between any distant pair of nodes due to the existence of relatively few long-range connections.

Rich-club organization

The centerpiece of the graph theoretical analysis of the functional connectomes in this article is the exploration of the rich-clubs of functional hub nodes in the infants' RSNs by computing the rich-club coefficient $ϕ (K)$ across a range of node degrees K. Rich-club organization implies that the hub nodes of a network are more strongly connected with each other than expected by their high degree K alone (Colizza et al., 2006; Schroeter et al., 2015). For each neonatal RSNs, a weighted rich-club coefficient $ϕ^{w} (K)$ , at each degree level of K, was computed as the ratio of the sum of the edge weights that connect a group of nodes (the club) of degree $\geq K$ to the sum of the strongest edge weights that such nodes could share with other nodes in the whole network (van den Heuvel and Sporns, 2011): $ϕ^{W} (K) = \frac{W_{> K}}{\sum_{l = 1}^{E_{> K}} W_{l}^{r a n k e d}}, l = 1, 2, \dots, E,$ (6)

where $W_{> K}$ denotes the sum of the weights of the number of edges between the members of the club with a degree $\geq K$ given by the strongest $E_{> K}$ edges of the network, where E is the total number of edges, and $W_{l}^{r a n k e d} \geq W_{l + 1}^{r a n k e d} \geq \dots \geq W_{E}^{r a n k e d}$ is a series representing the ranked weights on all network connections. A network will exhibit a rich-club topological organization if, for a range of degree K, the normalized weighted rich-club coefficient, $ϕ_{n o r m}^{W} (K) = \frac{ϕ^{W} (K)}{Φ_{n u l l}^{W} (K)}$ , is $> 1$ (Schroeter et al., 2015; van den Heuvel and Sporns, 2011). The $Φ_{n u l l}^{W} (K)$ is the rich-club effect assessed on the appropriate surrogate (null) model constructed by randomizing the network edges (10,000 realizations, see above). For each RSN, 10,000 randomized networks were first constructed. The weighted rich-club coefficient $ϕ_{n u l l}^{W} (K)$ was then computed for each level of degree K, for each randomized network. The average rich-club coefficient over the 10,000 randomized networks was then used as $Φ_{n u l l}^{W} (K)$ to compute the normalized weighted rich-club coefficient $ϕ_{n o r m}^{W} (K)$ .

Statistical analysis of rich-club organization

Permutation testing was used to assess statistical significance of the rich-club organization of each neonatal RSN. The null distribution of rich-club coefficients $ϕ_{n u l l}^{W} (K)$ was first obtained from the 10,000 random networks (see above). For each degree K for which $ϕ_{n o r m}^{W} (K) > 1$ , we tested whether the rich-club coefficient $ϕ_{n o r m}^{W} (K)$ is significantly greater than the 95% of the null $ϕ_{n u l l}^{W} (K)$ (van den Heuvel and Sporns, 2011). The null distribution was rank ordered and the 95% confidence interval determined nonparametrically from this rank order; (p $<$ 0.05, one-tailed test). The rich-club coefficients $ϕ_{n o r m}^{W} (K)$ above this threshold were considered statistically robust rich-club regimen.

Dynamics of RSNs in infants

The exploration of the spectral fingerprints and temporal dynamics of neonatal RSNs was performed by analyzing the time–frequency spectra [i.e., $y_{i, k} (t, f), i = 1, \dots, 5$ ] of each of the five RSNs for the kth infant [Eq. (5)]. The spectral fingerprints of a network reflect the frequency range at which the network of interest oscillates (Mehrkanoon et al., 2014a). The temporal dynamics of each network quantify the variability of the network of interest over time (Mehrkanoon et al., 2014a; Vidaurre et al., 2017).

Spectral fingerprints and temporal dynamics of RSNs

The spectral property and the temporal fluctuations of the networks were independently investigated by analyzing the marginal densities of the group-level time–frequency spectra $y_{i, k} (t, f)$ [Eq. (5)] (Mehrkanoon et al., 2014a). The fast timescale fluctuations of the ith RSN in the kth infant are given by $y_{i, k} (f) = \frac{1}{T} \sum_{t = 0}^{T} y_{i, k} (t, f)$ . The spectral fingerprints of the five RSNs were then computed by the grand average and standard deviation of $y_{i, k} (f)$ across 100 infants. We inferred the temporal dynamics of the ith RSN in the kth infant from the instantaneous mean of $y_{i, k} (t, f)$ as $x_{i, k} (t) = \frac{1}{F} \sum_{f = 1}^{F} y_{i, k} (t, f)$ . The temporal dynamics of the five RSNs were then obtained from the grand average and standard deviation of the instantaneous mean of $x_{i, k} (t)$ across 100 infants. Furthermore, the PSD of the temporal dynamics [ $x_{i, k} (t)$ ] of each RSN was computed to infer slow timescales (or modulations) of fast oscillatory activity of the RSN, revealing the envelope of the fast oscillatory activity of RS cortical network (Mehrkanoon et al., 2014a).

Long-range dependencies of RSN dynamics

The behavior of a complex dynamical system over time is sensitive dependence on an initial condition of the system. That is, the past behavior of the system influences its future behavior, indicating a long-range dependence property of the process that governs the system dynamics, and is related to fractality (Kobayashi et al., 1999). The Hurst exponent is a measure of the extent of long-range dependence of such process that is reflected by time series obtained from the system of interest. Here, for each infant, the behaviors of RS cortical networks are reflected by their affiliated temporal dynamics [ $x_{i, k} (t)$ , see above]. Therefore, we further assessed the temporal dynamics of the networks by computing the Hurst exponent of their temporal fluctuations as $H_{i, k} = \frac{l o g (R_{i, k} ∕ S_{i, k})}{l o g (T)}$ , where $R_{i, k} ∕ S_{i, k}$ denotes the so-called range statistics (Kobayashi et al., 1999). The $R_{i, k}$ is the difference between the maximum and minimum deviation from the mean, and $S_{i, k}$ is the standard deviation over total number of T samples for the ith network in the kth infant. The 0 $<$ $H_{i, k}$ $<$ 0.5 indicates short-range dependence, $H_{i, k} = 0.5$ indicates a random walk process that is uncorrelated, and 0.5 $<$ $H_{i, k}$ $<$ 1 indicates the existence of long-range dependencies or persistence (Kobayashi et al., 1999). The Hurst exponent of the temporal dynamics [ $x_{i, k} (t)$ ] of the ith RSN in the kth infant is given by $\frac{R_{i, k}}{S_{i, k}} = \frac{{max}_{1 \leq J \leq T} V_{i, k} - {min}_{1 \leq J \leq T} V_{i, k}}{\sqrt{\frac{1}{T} \sum_{t = 1}^{T} {(x_{i, k} (t) - μ_{x_{i, k}})}^{2}}}, i = 1, \dots, 5; k = 1, \dots, 100,$ (7)

where $V_{i, k} = \sum_{t_{0} = 1}^{J} x_{i, k} (t_{0}) - μ_{x_{i, k}}$ , $μ_{x_{i, k}} = \frac{1}{T} \sum_{t = 1}^{T} x_{i, k} (t)$ , $V_{i, k}$ denotes the deviation of the ith RSN temporal dynamics from its mean value $μ_{x_{i, k}}$ , and $J = T ∕ 4$ is a nonoverlapping sample size.

Temporal organization of interacting RSNs

Because these five RSNs capture the dynamics of the neonatal brain states in a low-dimensional space, we sought to explore between-networks temporal interactions. Such interactions or transitions between the brain states that are captured by a set of networks are referred to as meta-state dynamics (Vidaurre et al., 2017). To identify the pattern of interactions between the five RSNs, we assess the cross-correlation between the temporal dynamics of the networks ${\bar{x}}_{i, k} (t)$ . The cross-correlation between the temporal dynamics of two networks m and $m'$ in the kth infant is given as a function of a time lag $τ$ by $\begin{matrix} r_{{(m, m')}_{k}} (τ) & = ℰ {{\bar{x}}_{m} (t) {\bar{x}}_{m'} (t + τ)}, m, m' = 1, \dots, 5, m \neq m', \\ τ = - T, \dots, T . \end{matrix}$ (8)

To extract the component of the cross-correlation functions that is not invariant to time reversal, we determine the asymmetric components of $r_{{(m, m')}_{k}} (τ)$ as follows: ${\tilde{r}}_{{(m, m')}_{k}} (τ) = \frac{r_{{(m, m')}_{k}} (τ) - r_{{(m, m')}_{k}}^{T} (τ)}{2} .$ (9)

The slope signs of the asymmetric component near the zero time lag, ${\tilde{r}}_{{(m, m')}_{k}} (0)$ , were used to represent the relative leader (positive sign) and follower (negative sign) positions of the networks m and $m'$ . Examining all pairs of networks in each infant allows a temporal rank order (sequence) to be obtained. Confidence intervals for the asymmetric component of the cross-correlation were estimated using a nonparametric permutation approach. To this end, 1000 surrogate time series were constructed by random permutation of the network time series. The asymmetric component of the cross-correlation functions [Eq. (8)] derived from these surrogate time series represents trivial non-zero fluctuations of the magnitude expected from data of that string length and amplitude distribution but with no further structure of interest. Two-sided 95% confidence intervals were estimated by rank ordering this null distribution and choosing the 25th and 975th values. This process was repeated independently for each of the five networks of all infants, and the signs of the grand average asymmetric component were used to determine the temporal sequence (or cycle) between the five robust RSNs. The time lags that the networks need to take their leading–following positions during network interactions were determined by the first time point at which the network crosses $t = 0$ sec. These time lags were then used to distinguish between the lagged expression of the networks, yielding the patterns of cyclic temporal interactions of the networks.

Results

No significant difference between the PSD of EEG across the two groups

We found no significant difference between the PSD of the EEG of the two groups in the low-delta (1–2 Hz [0.085–0.11 mV $^{2}$ /Hz]) [ $t (7) = 1.12$ , p = 0.15, one-tailed; Fig. 2A], high-delta (2–4.5 Hz [0.042–0.048 mV $^{2}$ /Hz]) [ $t (7) = 1.23$ , p = 0.13, one-tailed; Fig. 2B], and the theta and low-alpha (5–10 Hz [0.022–0.027 mV $^{2}$ /Hz]) [ $t (7) = 1.35$ , p = 0.1, one-tailed; Fig. 2C]. The spatial pattern of the PSD for each group shows frequency-specific distributions such that the low-delta band (1–2 Hz) is widely present in the frontal and posterior areas of the scalp in the two groups, whereas high-delta (2–4.5 Hz; Fig. 2B) and the theta and low-alpha (5–10 Hz) are bilaterally distributed across the central, prefrontal, temporal, and parietal regions (Fig. 2B, C). The result shows that there is no significant difference between the two group-level PSDs in the delta, theta, and low-alpha bands, suggesting the association of the intrinsic cortical activity with band-limited neuronal oscillations in infancy, which is not equivalent to the broadband neuronal oscillations in the adults (Mehrkanoon et al., 2014a, 2016). Therefore, we pooled the EEG of the two groups to infer the RSNs that are consistently present across 100 infants.

FIG. 2.

The group-level spatial patterns of the PSD in the term and pre-term born groups. (A) The spatial pattern of the PSD for each group shows frequency-specific distributions such that the low-delta band (1–2 Hz) is widely present in the frontal and posterior areas of the scalp. (B) Small cliques bilaterally distributed PSD on the central and prefrontal regions. (C) Large cliques bilaterally distributed across the temporal, parietal, and central regions. PSD, power spectral density. Color images are available online.

Dynamic RS brain networks inferred from EEG source estimates in infants

The spatial ICA of the concatenated subjects' mPCA of non-zero time-lagged time–frequency coherence between the cortical activity signals in these 90 ROIs revealed spatially, spectrally, and temporally resolved distinctive modes of functional brain networks that are consistently expressed in 100 infants (Jackknifing test, one-sided p $<$ 0.01; Fig. 1B, right panel, C, D). The PSD of the cortical oscillations in these 90 anatomical ROIs shows statistically significant oscillations ranging from the delta to low-alpha frequency ranges (1–10 Hz; one-sided p $<$ 0.01, z-score = 2.5), dominated by low-frequency ranges (1.5–5 Hz) in both hemispheres (Fig. 1E). The results suggest that EEG source analysis has the capacity to reveal the spectral fingerprints of macroscopic and large-scale neuronal population oscillation in developing cortical GM regions in infants.

Rich-club organizations reflect cortical network states in infants

Our network derivation approach revealed five spatially independent patterns of RSNs, each with distinctive rich-club curves associated with specific rich-club topological organization (Fig. 3). The rich-club of interconnected cortical hubs in these five RSNs, respectively, show the existence of the default mode network (DMN), primary sensory-limbic system network, thalamo-frontal network, thalamo-sensorimotor, and visual-limbic system networks (Fig. 3). The normalized weighted rich-club curve (red), as the ratio of the weighted rich-club curve [ $ϕ^{W} (K)$ , blue] to the rich-club of random network [ $ϕ_{r a n d}^{W} (K)$ , gray], is significantly >1 (10,000 permutation tests, one-sided p $<$ 0.01) for the five RSNs in a set of node-degree ranges: $27 \leq K \leq 35 (D M N)$ , $22 \leq K \leq 31$ (primary sensory-limbic system network), $24 \leq K \leq 34$ (thalamo-frontal network), $23 \leq K \leq$ 33 (thalamo-sensorimotor network), and $22 \leq K \leq$ 31 (visual-limbic system network; Fig. 3, left column). The statistically independent RSNs with distinctive rich-club organizations that are commonly present in 100 infants suggest that functional structure of the brain's typical ability in information processing—the dynamic brain networks—is already present in infants at term equivalent age. The integration of EEG source-reconstruction and RSNs in infants can be a reliable quantitative tool for monitoring brain development in clinical settings. The distinctive topological patterns that capture dynamic brain network states are described blow.

FIG. 3.

Rich-club coefficients (left column) and organization of RSNs (glass-brain panels). (A–E, left column) The measures of weighted [ $Φ^{W} (K)$ , blue], null [ $Φ_{r a n d}^{W} (K)$ , gray], and normalized [ $Φ_{n o r m}^{W} (K)$ , red] rich-club curves for the DMN [that $Φ_{n o r m}^{W} (K)$ significantly increases in a range of degrees $27 \leq K \leq 35$ ], primary sensory-limbic system network ( $22 \leq K \leq 31$ ), thalamo-frontal ( $24 \leq K \leq 34$ ), thalamo-sensorimotor ( $23 \leq K \leq 33$ ), and visual-limbic system ( $22 \leq K \leq 31$ ). (A–E, glass-brain panels) Rich-club of cortical hubs (gold) that are connected by the rich-edges (gold). The feeders (cyan) and peripheral edges (purple), respectively, connect the non-hub regions to cortical hubs and other non-hubs. DMN, default mode network; RSN, resting-state network. Color images are available online.

The default mode network

The rich-club organization of the DMN at the node-degree K = 31 spans the left and right cortical hubs (gold) that are connected with each other by the rich-edges (gold): right superior frontal lobe, homologous posterior cingulate cortices, left cuneus, and right precuneus in the basic visual processing system (Fig. 3A, glass-brain panels). The functional connections among the basal ganglia (pallidum and caudate—the limbic system), right superior frontal lobe, and posterior cingulate cortex are richly structured in the DMN throughout the homologous thalamus regions. The non-hub cortical regions of homologous occipital lobes, frontal mid and superior gyri, left supplementary motor area, and somatosensory cortex, have a broad range of degrees $1 \leq K \leq 30 (p u r p l e)$ , which are densely connected to cortical hubs throughout the feeders (cyan; Fig. 3A, glass-brain panels).

Primary sensory-limbic system network

The rich-club architecture of the second network at the degree K = 27 reveals the primary functional hubs, mainly interhemispherically and bilaterally, including the primary auditory cortex (Heschl), primary somatosensory cortex (left postcentral), left insula (involved in the motor control and cognitive functioning), inferior parietal lobe, right olfactory, Rolandic operculum, right hippocampus, amygdala, and basal ganglia (putamen and caudate; Fig. 3B, glass-brain panels). The topological organization of the feeders densely links the non-hub cortical regions of superior, middle, and inferior temporal, frontal, and occipital lobes ( $1 \leq K \leq 25$ ) with the rich-club network within and between the hemispheres throughout the limbic system, somatosensory and the primary auditory cortex.

Thalamo-frontal network

The rich-club of hub cliques in the third network spans the limbic system and frontal gyri including the homologous thalamus regions, posterior and anterior cingulate cortices, right hippocampus, superior, medial, and middle frontal orbits, and rectus at the degree K = 28 (Fig. 3C, glass-brain panels). The feeders link non-hub cortical regions of the right supplementary motor area, left hippocampus, parahippocampal, amygdala, homologous olfactory systems, basal ganglia (putamen, pallidum, and caudate) and occipital lobe (in a range of degrees $1 \leq K \leq 26$ ), with the cortical hub regions of homologous anterior cingulate cortices and orbitofrontal gyri in the network.

Thalamo-sensorimotor network

The spatial pattern of the rich-club of hubs in the fourth RSN reveals the sensorimotor-thalamic network spanning the thalamus, basal ganglia (putamen and caudate), supplementary motor area, primary motor cortex, somatosensory, superior/inferior frontal gyri at degree K = 27 (Fig. 3D, glass-brain panels). The feeders densely link the noncortical hubs of occipital lobes, superior, middle, and inferior frontoparietal lobes, and cingulate cortices, in a range of degrees $1 \leq K \leq 25$ , with the rich-club of cortical hubs. The temporal lobe and primary auditory cortices have minimal functional connections to the rich-club of hubs compared with the other non-hub regions.

Visual-limbic system network

The spatial pattern of the rich-club of hub cliques in the fifth RSN shows the visual-limbic system network spanning the superior/middle occipital gyri, cuneus, putamen, posterior cingulate cortex, parahippocampal gyrus, hippocampus, amygdala, and olfactory at K = 27 (Fig. 3E, glass-brain panels). The feeders, throughout the rich-club of occipital lobe and the limbic system, link the homologous frontal lobes and anterior/middle cingulate gyri and the right occipital lobe with the rich-club brain hubs in a range of degrees $1 \leq K \leq 25$ . Crucially, the rich-club of hubs is more prominent in the left visual cortex, cuneus, and putamen than the homologous regions in the left hemisphere.

Richly structured RSNs in infants are functionally vulnerable

Graph theoretical analysis of the networks revealed the coexistence of functional integration and segregation capacity, small-world topology, and functional vulnerability in each of the five spatially independent RSNs (Table 1). Although the global node-degree/strength of the DMN (31/25) is greater than the others confined in the interval of [27/13, 29/16], the visual-limbic system network shows a greater mean betweenness-centrality (282.14) than the others nested in a range of [128, 278]. All networks show high mean clustering coefficient and small short path length in the intervals of [65, 78] and [0.6, 1.07], respectively. These five RSNs exhibit small-world topology with the values significantly >1 (1000 network randomization tests, one-sided p $<$ 0.01), which are bounded in the interval of [1.15, 1.35] (Table 1).

Table 1.

Graph Theoretical Metrics of Resting-State Networks in Infants at Term Age

Graph metric	Default mode	Sensory-limbic	Thalamo-frontal	Thalamo-sensorimotor	Visual-limbic
Degree/strength^a	31/25	27.3/15.25	28.7/13.8	27.8/15.9	27/15.92
Clustering coefficient	76.66	68	76.01	65.17	77.54
Betweenness	128.9	277.9	201.6	245.35	282.14
Global efficiency^b	1.011	1.021	1.011	1.021	1.01
Assortativity	−0.47	−0.52	−0.35	−0.44	−0.45
Average path length	1.07	0.68	0.72	0.86	0.83
Small-world index	1.25	1.28	1.34	1.15	1.26

Measures represent mean + 1 SD.

Normalized values.

Interestingly, although the normalized values of the global efficiency measures of the networks are >1 (1000 network randomization tests), which indicate the existence of functional integration capacity in the infant brain at term age, the negative values of the assortativity coefficients of the networks (i.e., [−0.52, −0.35]) reveal the vulnerability (or weak resilience) of high-degree hubs in these five RSNs (Maslov and Sneppen, 2002). The results mean that neonatal cortical activity in the macroscopic large-scale level is functionally immature, suggesting the presence of neural plasticity in the large-scale functional organization of neonatal brain. The result offers a capacity to monitor brain health with respect to the functional resilience of brain function development in preterm infants.

Spectral fingerprints of RSNs in infants reflect the already established cortical network-dependent carrier frequency at term age

The spectral fingerprints of the five RSNs reveal the carrier frequency (or fast timescale) at which the spatially distributed large-scale cortical oscillations in the brain exhibit synchronous regimen over time (Hipp et al., 2012; Mehrkanoon et al., 2014a). In the DMN, a wide range of significant functional connectivities (0.28–0.57, one-sided p $<$ 0.01, z-score = 2.45) occur in the delta and low-theta bands (1–5.5 Hz) with small standard deviation across infants (red, Fig. 4A, left panel), whereas the primary sensory-limbic system network exhibits significant functional connectivity (0.08–0.28, one-sided p $<$ 0.01, z-score = 2.52) at multiple distinctive frequencies with the peaks at 1.5, 3, 4.5, and 11 Hz, which have a larger standard deviation across infants (green, Fig. 4B, left panel). The fast timescale of the thalamo-frontal network demonstrates significantly strong functional connectivity (0.15–0.4, one-sided p $<$ 0.01, z-score = 2.5) with three distinctive peaks at 2 and 3.5 Hz (the delta band) and 6 Hz in the theta band (blue, Fig. 4C, left panel).

FIG. 4.

Spectral fingerprints, temporal dynamics, and slow timescale oscillations of the five robust RSNs in 100 infants. (A–E, left column) Cortical activity synchronization profile of RSNs captured by their frequency spectra over the time interval of [0, 240] sec. (A–E, middle column) Instantaneous fluctuation profile of RSNs inferred from their temporal dynamics over the time interval of [0, 240] sec. (A–E, right column) Nested-dynamics profile of RSNs obtained from their envelopes of the temporal dynamics. Color images are available online.

The spectral fingerprints of the sensorimotor-thalamic and visual-limbic system networks, respectively, show significantly strong functional connectivity values of 0.38 (cyan) and 0.32 (pink; one-sided p $<$ 0.01, z-score = 2.4) at 1.6 Hz (low-delta band), followed by multiple peaks in the frequency intervals of [2.5–6] Hz (Fig. 4D), [2–3] Hz and 10 Hz (Fig. 4E). Interestingly, the fast timescales of the networks of primary sensory-limbic system (Fig. 4B), sensorimotor-thalamic (Fig. 4D), and the visual-limbic system (Fig. 4E) show the existence of a significant 10 Hz component in their spectral fingerprints (one-sided p $<$ 0.01, z-score = 2.4). The results suggest that macroscopic functional circuits between the neocortex and subcortical regions are already present at term age, offering a potential baseline to help determine the lack of functional cortical circuit development in very-preterm and preterm infants.

Infant brain networks exhibit richly organized temporal dynamics

The instantaneous fluctuations of the fast timescale of the RSNs are reflected by the temporal dynamics at slow timescales, unfolding the brain network dynamics (Mehrkanoon et al., 2014a). The temporal dynamics of the five RSNs show slow timescale oscillatory patterns of the connectivity strengths confined in [−0.3, 0.3] interval with multiple peaks between $20 s e c \leq t \leq 220 s e c$ (Fig. 4, middle column). The DMN largely fluctuates in the time interval of [100, 220] sec (red, Fig. 4A, middle column), and the temporal dynamics of other four RSNs show multiple fluctuations across a broader time interval of [10, 240] sec (green, blue, cyan, and pink), respectively (Fig. 4B–E, middle column). Together, all RSNs temporal dynamics show small standard deviation across infants ( $<$ 0.1; gray curves, Fig. 4B–E, middle column). The temporal dynamics (middle column) of the RSNs have a specific spectral density characteristic, showing 1/f-like distribution that is quite stationary across the five networks (Fig. 4A–E, right column). Note that the slow timescale of the dFC does not represent the low frequency of EEG, rather it shows the slow modulation of fast timescale RSNs in infants. Remarkably, temporal dynamics of the RSNs exhibit highest power at frequencies well below 0.16 Hz, which is the characteristic frequency of RS-fcMRI (Fig. 4, right column).

All RSNs show peaks at $f = 0.022$ and $f = 0.12$ Hz, respectively. Note that the timescales of the RSNs in the time domain are not limited by the 0.5 Hz high-pass filter applied to the raw EEG signals because the instantaneous mean of the time–frequency spectra of Eq. (5) reflects the envelope of the fast oscillatory activity and hence unfolds at much slower timescales. Crucially, the results suggest that neonatal brain dynamics are not constrained by the age of the brain, revealing a signature of healthy immature brain that exhibits an adult-like cortical network dynamics (Mehrkanoon et al., 2014a). This approach offers a monitoring tool to give functional insights into the cortical network dynamics in very preterm infants, providing a platform for clinical assessments. The characteristics of the temporal dynamics of RSNs are described in detail below.

Long-range correlations of temporal dynamics of networks reflect scale invariance structure of cortical oscillations in infancy

Because the human brain can be thought of as a complex dynamical system (Kobayashi et al., 1999; Mehrkanoon et al., 2014a), temporal activity of such a system has a sensitive dependence on the initial condition of the system (Glasner and Weiss, 1993), indicating that the past temporal activity of the brain influences its future behavior (Kobayashi et al., 1999). To quantify the extent of a long-range dependence property of the infant brain network dynamics, the Hurst exponent analysis of the temporal dynamics of the five RSNs that exhibit apparent power-law (or 1/f-like) scaling of the low-frequency activity was performed (Fig. 5). The average Hurst exponent of the temporal dynamics of the RSNs lies in the range of 0.85–0.9 (Fig. 5A), indicating that the dFCs in infants at term age have persistent, long-range temporal dependencies, a property that is far from a random process and is consistent with slow power-law decay.

FIG. 5.

Temporal organization of brain network dynamics in infants. (A) The Hurst exponent of each of the five RSNs. (B) Hierarchical clustering of the correlation matrix of the networks' temporal dynamics. (C) Representative temporally fast and slow network expressions at t = 0 sec and $t \neq 0$ sec inferred from the cross-correlation between the temporal dynamics of the DMN and the others. (D–H) The leading and following properties of the five RSNs were inferred from the analysis of positive and negative slopes of the asymmetric component of the cross-correlations between the networks' temporal dynamics in the time interval of [−40, 40] sec—the 95% confidence interval (dashed line) for trivial non-zero values due to finite sample length was estimated using a nonparametric permutation test. (I–J) The temporal organizations of RSN interactions in the fast (leading networks: DMN and thalamo-frontal, and following networks: primary sensory-limbic system, thalamo-sensorimotor, and visual-limbic system) and slow sequence of cycles (leading networks: thalamo-sensorimotor and visual-limbic, and following networks: DMN, primary sensory-limbic system, and thalamo-frontal). (K) A switch (or a swap) in the dynamics of leading and following networks among the temporally fast and slow network expressions with primary sensory-limbic system network being as a pass-through state. Color images are available online.

Temporal dynamics of RSNs reflect fast–slow sequence of cortical network states in infants

The correlation matrix of the RSNs temporal dynamics show four levels of hierarchically interconnected clusters (1000 permutation tests, one-sided p $<$ 0.01, z-score = 2.4; Fig. 5B). The temporal dynamics of the DMN and sensorimotor-thalamic system network form a cluster in level 1, that is, C1 = (1,4), followed by the three hierarchical levels C2 = (thalamo-frontal network, C1), C3 = (visual-limbic system network, C2), and C4 = (primary sensory-limbic system network, C3; Fig. 5B). Networks are modulated together by similar patterns of slow timescale (Fig. 3, right column); however, the dynamics in some networks are expressed earlier than others during each timecourse of co-expression (or co-modulation). This temporal asymmetry in networks expression can be observed by the symmetric (Fig. 5C) and asymmetric components of the cross-correlation between the temporal dynamics of RSNs (panels D–H). For example, the temporal fluctuation of the thalamo-sensorimotor network is strongly correlated with the others at zero time-lag correlation, slowly decaying toward zero over time delay (Fig. 5C)—similar findings pertain to all other networks. Note that, the positive/negative slopes of the temporal asymmetric cross-correlation between the networks at around t = 0 sec determine the leading–following position of the network under analysis.

Crucially, the existence of temporal asymmetry in RSNs expression shows fast and slow sequences of directional interactions between RSNs in cortical activity. Each of the two temporally fast and slow sequences of leading–following network interactions show two cycles: a three-network cycle and a four-network cycle. The three-network cycle consists of {(DMN, thalamo-frontal, primary sensory-limbic system); (primary sensory-limbic system, thalamo-sensorimotor, visual limbic)} (panel I) and {(DMN, thalamo-sensorimotor, primary sensory-limbic system); (primary sensory-limbic system, visual limbic, thalamo-sensorimotor)} (panel J). The four-network cycle includes {(DMN, thalamo-frontal, primary sensory-limbic system, thalamo-sensorimotor); (primary sensory-limbic system, thalamo-sensorimotor, visual limbic, thalamo-frontal)} (panel I) and {(DMN, visual limbic, thalamo-sensorimotor, primary sensory-limbic system)} (panel J).

As an example, the thalamo-frontal network follows the networks of default mode, primary sensory-limbic system, thalamo-sensorimotor, and visual limbic at the time delays of 33.8, 23.9, 23.3, and 21.3 sec, respectively (panel J). The DMN and thalamo-frontal network, with in-degree = 3, are the leading networks in the fast sequence of network expression (panel I), tending to follow the sensorimotor–thalamic system and visual-limbic system networks in the slow sequence of network expression (panel J), revealing the existence of a swap in the leading and following networks among the temporally fast and slow sequences of network expressions. The primary sensory-limbic system network that has equivalent scores of in-degree = 2 and out-degree = 2 among the fast and slow sequences of network expressions plays a pass-through network in the dynamical interplay between the other four networks (panel K).

Discussion

The analysis of RS-EEG source reconstruction in 100 infants at term age revealed five robust RSNs that reveal significant new information; these identified RSNs capture, for the first time in human infants' functional brain network studies, the dominant patterns of anatomically constrained spatially, spectrally, and temporally resolved dynamic synchronization between distributed cortical activity. These RSNs are constituted by a narrowband carrier frequency (1–8 Hz) that is nested by much slower and persistent fluctuations at $<$ 0.1 Hz, showing that the envelope of the networks carrier frequencies are modulated across a broad spectrum of very slow timescales. Together, the slow modulation of fast timescales of the networks provide a snapshot of the nested dynamics of RS brain network activity that is not obvious in functional MRI (fMRI)-derived networks. The networks have independent topological organizations, distinctive spectral fingerprints, and hierarchically correlated temporal dynamics with fast and slow sequences of network states, effectively consistent with dynamic repertoire and macroscopic RS cortical activity in the human adult brain (Vidaurre et al., 2017). Our findings of the DMN, primary functional systems, and thalamocortical networks are consistent with the immature DMN, dorsal attention network configuration, and primary functional system networks obtained from term-age infants' RS-functional/diffusion MRI data (Ball et al., 2014; Batalle et al., 2017; Cao et al., 2017; Doria et al., 2010; Fransson et al., 2011; Van Den Heuvel et al., 2015).

Functional connectivity between the thalamus and other brain regions (e.g., basal ganglia, primary sensory, and sensorimotor systems) in these five RSNs suggests that long-range functional connections are already established at birth. These functional connections regulate the cortico-cortical functional networks and high-order cortices that support salience, executive, integrative, and cognitive functions in infancy and childhood (Cao et al., 2017; Gilmore et al., 2018; Tau and Peterson, 2009). Coexistence of the small-world topology and vulnerability of these five RSNs suggests that the primary dynamical complexity of the brain functions has the necessary, but insufficient, conditions for information processing at term age; this is consistent with the topological principles of the RSNs identified by diffusion MRI and fMRI studies in infants (Ball et al., 2014; Cao et al., 2017; Doria et al., 2010; Fransson et al., 2011; Van Den Heuvel et al., 2015); however, such poor network vulnerability gradually decreases during the first 2 years of life (Cao et al., 2017; Gao et al., 2011).

These networks have distinctive spectral fingerprints with narrowband frequency bases: The DMN has a bell-shape peak within the delta-theta bands, and the primary sensory, thalamocortical, and limbic system networks have peaks in the narrowband delta-theta and low-alpha frequencies The delta-band spectral fingerprint of the thalamo-frontal and thalamo-sensorimotor networks shows that low-frequency macroscopic neuronal oscillations in the thalamus regulate sensorimotor and frontal cortex, consistent with previous reports (Timofeev and Chauvette, 2011). The fMRI findings have shown that the visual networks start (or function) at birth (Doria et al., 2010; Fransson et al., 2007). The coexistence of peaks in the delta and low-alpha bands (as opposed to only alpha or low-alpha band) in the visual-limbic system network suggests that the visual networks start immaturely at birth. Because alpha oscillations increase in power from 4 to 10 months (Cornelissen et al., 2015), the visual networks observed at birth cannot be seen as already established mature visual networks.

The neuronal mechanisms underpinning network dynamics (or brain state) are governed by a long-range correlation process. The existence of this nonrandom process in cortical activity at birth is proven by the Hurst exponents of the temporal dynamics of the RSNs. The values of the Hurst exponents are >0.5 (0.88–0.9), indicating long-range correlation in the underlying process of network dynamics, which is consistent with previous studies (Mehrkanoon et al., 2014a). The coexistence of long-range correlation in network dynamics, distinctive spectral and spatial patterns of RSNs in infants suggest that the brain's intrinsic functional organization underpins the establishment of dynamic functional networks in infancy. This is consistent with macroscopic cortical network dynamics that evolve on very slow timescale and have long-range persistence in the adult RSN dynamics (Linkenkaer-Hansen et al., 2001; Mehrkanoon et al., 2014a).

Networks with a wide range of peaks in the delta-theta bands (DMN and thalamo-frontal network) are typically expressed slightly earlier than others during each fast cycle of network expression (i.e., t = 0 sec). In contrast, networks with distant, but distinctive, peaks in the delta, theta, and low-alpha bands are expressed over slow cycle of network expression t $\approx$ 25 sec (e.g., thalamo-sensorimotor and visual-limbic networks).

This timing profile of network dynamics suggests an indirect brain state transition from spectrally and temporally ordered network state to those with distinctive peaks during each wave of network expression. This is consistent with the view of dynamic FC, which is the direct product of intrinsic brain electrical activity at distinctive frequencies in the adult brain (Tagliazucchi et al., 2012). The present findings suggest that infant brain network dynamics arise from long-range synchronization of band-limited cortical oscillations subjected to an interplay between temporally fast and slow sequences of coherent cortical oscillations. The timing profile of this interplay is not temporally symmetric among the networks, but rather the appearance of some networks such as the default mode and thalamo-frontal characteristically leads the primary functional limbic, thalamo-sensorimotor, and visual-limbic system networks (Fig. 5I).

The temporally fast and slow sequences of network interactions identify the leader and follower networks in the evolution of the dynamic repertoire of RS cortical activity; this provides key timing information by when the temporal formation and dissolution of the brain states are captured (Fig. 5B, I–K), reflecting immature meta-state dynamics and evolutionary traveling waves in the infant brain, consistent with the hierarchically organized cortical network dynamics that exhibit spontaneous travelling waves in the adult brain (Vidaurre et al., 2017). Furthermore, spectral fingerprints of the networks identify the neuronal oscillatory mechanisms underpinning the evolution of functional networks. Such functional network evolution as identified by thalamocortical FC (which is part of Fig. 3A, C, D) and cortico-limbic system network (Fig. 5B, E) reveals immature cognitive and memory network formation in infants, which has been shown to have the potential to predict neurocognitive outcomes in children born preterm (Ball et al., 2015). The fast timescale, temporal dynamics, and slow timescale (or modulation) of RSNs in infants offer a system neuroscience tool to help assess the development of functional brain networks. Particularly, our future work will investigate which of the spatial, spectral, and temporal information features of RSNs correlate with adverse neurodevelopmental outcomes.

In conclusion, this study of a large cohort of human infants has shown that integration of mPCA and spatial ICA algorithms in conjunction with time–frequency EEG source estimates provides a novel method for inferring patterns of robust dynamic functional brain connectivity. These richly organized dynamic brain networks captured the evolution of cortical oscillation synchronization at multiple frequencies. The characteristic temporal dynamics of the networks prove the existence of a scale-free property in the neonatal cortical activity in which the brain network states mutually correlate with each other to establish the functional organization of the brain in infants. The findings suggest that the functional organization of nested RS brain network dynamics does not solely depend on the maturation of cognitive networks; instead, the brain network dynamics exist in infants at term age well before the mature brain networks present in childhood and adulthood, and hence, the existence of both spontaneous cortical oscillations and network dynamics in infants at term age offer a quantitative measurement of neurotypical development in infants. Furthermore, this work offers a novel contribution to the assessment of brain network identification and assessments in infants by quantifying anatomically constrained cortical networks. This provides a new functional basis for future work to design an effective clinical diagnostic tool with the capacity for early accurate diagnosis of neurodevelopmental disability and prognostication of neurodevelopmental outcome in sick preterm and term babies.

Footnotes

Authors' Contributions

S.M. designed the research problem in collaboration with P.C. S.M. performed research, analyzed data, and interpreted the results. P.C. contributed to the clinical interpretation of the results. B.B. commented on the final version of the article. S.M., B.B., and P.C. wrote the article.

Acknowledgments

The preprint postings of this article can be found in the following servers: arXiv (2019/02/11, https://arxiv.org/abs/1902.03725) and bioRxiv (2019/02/11, https://www.biorxiv.org/content/10.1101/545616v2).

Author Disclosure Statement

No competing financial interests exist.

Funding Information

S.M. was supported by an Advance Queensland Fellowship Grant (AQRF06016-17RD2) funded by the Queensland State Government, Australia (2017–2019). Data collection was supported by the National Health and Medical Research Project grant APP1084032. B.B. was supported by Qatar Foundation grants NPRP6-885-2-364 and NPRP6-680-2-282.

References

Achard

, Salvador

, Whitcher

, et al. 2006. A resilient, low-frequency, small-world human brain functional network with highly connected association cortical hubs. J Neurosci, 26:63–72.

Allen

, Damaraju

, Plis

, et al. 2014. Tracking whole-brain connectivity dynamics in the resting state. Cereb Cortex, 24:663–676.

Arichi

, Whitehead

, Barone

, et al. 2017. Localization of spontaneous bursting neuronal activity in the preterm human brain with simultaneous EEG-fMRI. eLife, 6:e27814.

Ball

, Aljabar

, Zebari

, et al. 2014. Rich-club organization of the newborn human brain. PNAS USA, 111:7456–7461.

Ball

, Pazderova

, Chew

, et al. 2015. Thalamocortical connectivity predicts cognition in children born preterm. Cereb Cortex, 25:4310–4318.

Batalle

, Hughes

, Zhang

, et al. 2017. Early development of structural networks and the impact of prematurity on brain connectivity. NeuroImage, 149:379–392.

Belouchrani

, Abed-Meraim

, Cardoso

, et al. 1997. A blind source separation technique using second-order statistics. IEEE Trans Signal Process, 45:434–444.

Biswal

, Yetkin

, Haughton

, et al. 1995. Functional connectivity in the motor cortex of resting human brain using echo-planar MRI. Magn Reson Med, 34:537–541.

Boersma

, Smit

, de Bie

, et al. 2011. Network analysis of resting state EEG in the developing young brain: structure comes with maturation. Human Brain Mapp, 32:413–425.

10.

Cao

, He

, Dai

, et al. 2017. Early development of functional network segregation revealed by connectomic analysis of the preterm human brain. Cereb Cortex, 27:1949–1963.

11.

Colizza

, Flammini

, Serrano

, et al. 2006. Detecting rich-club ordering in complex networks. Nat Phys, 2:110–115.

12.

Cornelissen

, Kim

S-E

, Purdon

, et al. 2015. Age-dependent electroencephalogram (EEG) patterns during sevoflurane general anesthesia in infants. eLife, 4:e06513.

13.

Dang

, Ng

. 2011. Finite difference neuroelectric modeling software. J Neurosci Methods, 198:359–363.

14.

Doria

, Beckmann

, Arichi

, et al. 2010. Emergence of resting state networks in the preterm human brain. PNAS USA, 107:20015–20020.

15.

Fransson

, Åden

, Blennow

, et al. 2011. The functional architecture of the infant brain as revealed by resting-state fMRI. Cereb Cortex, 21:145–154.

16.

Fransson

, Metsäranta

, Blennow

, et al. 2013. Early development of spatial patterns of power-law frequency scaling in fMRI resting-state and EEG data in the newborn brain. Cereb Cortex, 23:638–646.

17.

Fransson

, Skiöld

, Horsch

, et al. 2007. Resting-state networks in the infant brain. PNAS USA, 104:15531–15536.

18.

Gao

, Gilmore

, Giovanello

, et al. 2011. Temporal and spatial evolution of brain network topology during the first two years of life. PLoS One, 6:e25278.

19.

Gilmore

, Knickmeyer

, and Gao

. 2018. Imaging structural and functional brain development in early childhood. Nat Rev Neurosci, 19:123–137.

20.

Glasner

, Weiss

. 1993. Sensitive dependence on initial conditions. Nonlinearity, 6:1067–1075.

21.

Grayson

, Fair

. 2017. Development of large-scale functional networks from birth to adulthood: a guide to the neuroimaging literature. NeuroImage, 160:15–31.

22.

Hipp

, Hawellek

, Corbetta

, et al. 2012. Large-scale cortical correlation structure of spontaneous oscillatory activity. Nat Neurosci, 15:884–890.

23.

Huster

, Plis

, Calhoun

. 2015. Group-level component analyses of EEG: validation and evaluation. Front Neurosci, 9:254–254.

24.

Kaminska

, Delattre

, Laschet

, et al. 2018. Cortical auditory-evoked responses in preterm neonates: revisited by spectral and temporal analyses. Cereb Cortex, 28:3429–3444.

25.

Kobayashi

, Misaki

, Nakagawa

, et al. 1999. Non-linear analysis of the sleep EEG. Psychiatry Clin Neurosci, 53:159–161.

26.

Lee

, Lee

. 2018. The effects of early-stage neurodevelopmental treatment on the growth of premature infants in neonatal intensive care unit. J Exerc Rehabil, 14:523–529.

27.

Linkenkaer-Hansen

, Nikouline

, Palva

, et al. 2001. Long-range temporal correlations and scaling behavior in human brain oscillations. J Neurosci, 21:1370–1377.

28.

Marlow

, Wolke

, Bracewell

, et al. 2005. Neurologic and developmental disability at six years of age after extremely preterm birth. N Engl J Med, 352:9–19.

29.

Maslov

, Sneppen

. 2002. Specificity and stability in topology of protein networks. Science, 296:910–913.

30.

Mehrkanoon

, Boonstra

, Breakspear

, et al. 2016. Upregulation of cortico-cerebellar functional connectivity after motor learning. NeuroImage, 128:252–263.

31.

Mehrkanoon

, Breakspear

, Boonstra

. 2014a. Low-dimensional dynamics of resting-state cortical activity. Brain Topogr, 27:338–352.

32.

Mehrkanoon

, Breakspear

, Boonstra

. 2014b. The reorganization of corticomuscular coherence during a transition between sensorimotor states. NeuroImage, 100:692–702.

33.

Mehrkanoon

, Breakspear

, Britz

, et al. 2014c. Intrinsic coupling modes in source-reconstructed intrinsic coupling modes in source-reconstructed electroencephalography. Brain Connect, 4:812–825.

34.

Mehrkanoon

, Breakspear

, Daffertshofer

, et al. 2013. Non-identical smoothing operators for estimating time-frequency interdependence in electrophysiological recordings. EURASIP J Adv Signal Process, 2013:1–16.

35.

Meijer

, Hermans

, Zwanenburg

, et al. 2014. Functional connectivity in preterm infants derived from EEG coherence analysis. Eur J Paediatr Neurol, 18:780–789.

36.

Mohammadi-Nejad

A-R

, Mahmoudzadeh

, Hassanpour

, et al. 2018. Neonatal brain resting-state functional connectivity imaging modalities. Photoacoustics, 10:1–19.

37.

Moore

, Hennessy

, Myles

, et al. 2012. Neurological and developmental outcome in extremely preterm children born in England in 1995 and 2006: the EPICure studies. BMJ, 345:e7961.

38.

Myers

, Grieve

, Izraelit

, et al. 2012. Developmental profiles of infant EEG: overlap with transient cortical circuits. Clin Neurophysiol, 123:1502–1511.

39.

Nolte

, Bai

, Wheaton

, et al. 2004. Identifying true brain interaction from EEG data using the imaginary part of coherency. Clin Neurophysiol, 115:2292–2307.

40.

Roche-Labarbe

, Aarabi

, Kongolo

, et al. 2008. High-resolution electroencephalography and source localization in neonates. Human Brain Mapp, 29:167–176.

41.

Routier

, Mahmoudzadeh

, Panzani

, et al. 2017. Plasticity of neonatal neuronal networks in very premature infants: source localization of temporal theta activity, the first endogenous neural biomarker, in temporoparietal areas. Human Brain Mapp, 38:2345–2358.

42.

Schroeter

, Charlesworth

, Kitzbichler

, et al. 2015. Emergence of rich-club topology and coordinated dynamics in development of hippocampal functional networks in vitro. J Neurosci, 35:5459–5470.

43.

Schumacher

, Stiris

, Larsson

. 2015. Effective connectivity in long-term EEG monitoring in preterm infants. Clin Neurophysiol, 126:2261–2268.

44.

Shi

, Yap

, Wu

, et al. 2011. Infant brain atlases from neonates to 1- and 2-year-olds. PLoS One, 6:e18746.

45.

Sutton

, Darmstadt

. 2013. Preterm birth and neurodevelopment: a review of outcomes and recommendations for early identification and cost-effective interventions. J Trop Pediatr, 59:258–265.

46.

Tagliazucchi

, von Wegner

, Morzelewski

, et al. 2012. Dynamic bold functional connectivity in humans and its electrophysiological correlates. Front Human Neurosci, 6:1–22.

47.

Tau

, Peterson

. 2009. Normal development of brain circuits. Neuropsychopharmacology, 35:147–168.

48.

Timofeev

, Chauvette

. 2011. Thalamocortical oscillations: local control of EEG slow waves. Curr Top Med Chem, 11:2457–2471.

49.

Tokariev

, Videman

, Palva

, et al. 2016. Functional brain connectivity develops rapidly around term age and changes between vigilance states in the human newborn. Cereb Cortex, 26:4540–4550.

50.

Van Den Heuvel

, Kersbergen

, De Reus

, et al. 2015. The neonatal connectome during preterm brain development. Cereb Cortex, 25:3000–3013.

51.

van den Heuvel

, Sporns

. 2011. Rich-club organization of the human connectome. J Neurosci, 31:15775–15786.

52.

Veen

BDV

, Drongelen

, Yuchtman

, et al. 1997. Localization of brain electrical activity via linearly constrained minimum variance spatial filtering. IEEE Trans Biomed Eng, 44:867–880.

53.

Vidaurre

, Smith

, Woolrich

. 2017. Brain network dynamics are hierarchically organized in time. PNAS USA, 114:12827–12832.

54.

Watts

, Strogatz

. 1998. Collective dynamics of ‘small-world’ networks. Nature, 393:440–442.