Breakdown of Whole-brain Dynamics in Preterm-born Children.

The brain operates at a critical point that is balanced between order and disorder. Even during rest, unstable periods of random behavior are interspersed with stable periods of balanced activity patterns that support optimal information processing. Being born preterm may cause deviations from this normal pattern of development. We compared 33 extremely preterm (EPT) children born at < 27 weeks of gestation and 28 full-term controls. Two approaches were adopted in both groups, when they were 10 years of age, using structural and functional brain magnetic resonance imaging data. The first was using a novel intrinsic ignition analysis to study the ability of the areas of the brain to propagate neural activity. The second was a whole-brain Hopf model, to define the level of stability, desynchronization, or criticality of the brain. EPT-born children exhibited fewer intrinsic ignition events than controls; nodes were related to less sophisticated aspects of cognitive control, and there was a different hierarchy pattern in the propagation of information and suboptimal synchronicity and criticality. The largest differences were found in brain nodes belonging to the rich-club architecture. These results provide important insights into the neural substrates underlying brain reorganization and neurodevelopmental impairments related to prematurity.

Introduction

A healthy brain operates close to a critical point, which has been identified as the boundary between order and disorder. Order has been described as a synchronized cluster between a group of brain areas and disorder as desynchronization between the remaining brain areas (Deco and Jirsa 2012; Deco et al. 2017a). This ensures that the information processing function operates at an optimal level (Deco et al. 2017a). Even when it is at rest, the brain is still balancing the competing activities of different network assemblies so that it can respond flexibly and rapidly to stimuli (Deco et al. 2009; Moutard et al. 2015). Any disturbance to these dynamics, due to genetic or environmental factors, can have a significant impact on healthy brain functioning. Preterm birth affects 11.1% of all live births worldwide (Blencowe et al. 2013) and is one of the most common clinical risk factors for neurodevelopmental disorders (Twilhaar et al. 2018), as it has a significant effect on structural (Ball et al. 2014; van den Heuvel et al. 2015) and functional brain connectomes (Fransson et al. 2011; Cao et al. 2017). Being born very preterm, at less than 32 weeks of gestation, and without major neurological impairment can potentially lead to altered spontaneous cortical activity. This may affect the optimal flow of information and neural synchronization (Doesburg et al. 2011b).

The late effects of being born preterm can be seen during childhood, as it affects the organization of the structural and functional networks of the child’s brain. Children born preterm have been reported to have more segregated and less integrated structural brains than their term-born counterparts, (Thompson et al. 2016). The hierarchical organization of their brains is disrupted, and this has important repercussions on the higher order areas (Lax et al. 2013). This leads to them having a less efficient and robust brain (Fischi-Gomez et al. 2016), which may contribute to the cognitive difficulties reported in preterm populations (Twilhaar et al. 2018). However, there is little information regarding the spatiotemporal organization of functional networks in preterm children. Studies using magnetoencephalography (MEG) have shown disrupted functional network organization (Ye et al. 2016) and alterations in oscillatory neural synchrony affecting short-memory performance in very preterm children (Moiseev et al. 2015). Despite that, no studies have been carried out that have used magnetic resonance imaging (MRI) data to delve into the effects that extremely preterm (EPT) birth, being born before 27 weeks of gestation, may have on the propagation of neural activity and criticality of the brain.

We applied two model-based approaches that allowed us to explore resting-brain functional organization on both a local and global basis. The first approach used a local-to-global framework that identified local ignitory events that propagated neural activity in the brain and led to global integration (Deco and Kringelbach 2017; Deco et al. 2017b). The second approach used a whole-brain Hopf model, which examined whether the brain operated synchronously, noisily, or under a critical regime (Deco et al. 2017a). To do this, we used structural, functional, and diffusion MRI (dMRI) data. Our hypothesis was that when EPT-born children were compared with healthy term-born controls they would demonstrate reduced ignitory events that affected the hierarchy of information processing, together with changes in neural synchronization and criticality during rest.

Materials and Methods

Population

The study population was part of the Extremely Preterm Infants in Sweden Study cohort that comprises perinatal data and long-term follow-up data (Serenius et al. 2016). The participants in this study were recruited from a previously described cohort of 111 preterm children born at less than 27 weeks of gestation in Stockholm County, Sweden, between 2004 and 2007. All the participants had successfully undergone an MRI scan at term-equivalent age (Padilla et al. 2015). The children were followed up and rescanned when they were a mean of 10 years ±2 months, and 66 EPT children and 46 term-born controls were included in this study (Supplementary Fig. 1). EPT infants with congenital infections and malformations were excluded. In addition, we did not include EPT infants with periventricular leukomalacia or intraventricular hemorrhage (grade III–IV), focal brain lesions (cysts and malformations), persistent ventricular dilatation, or moderate or severe white matter abnormalities qualitatively defined by MRI examinations at term-equivalent age (Skiold et al. 2010). The exclusion criteria for the term-born children were the same as the ones we used for the EPT sample, with the addition of any history of birth complications. The final sample included 33 EPT and 28 term-born children. The regional ethical review board in Stockholm granted ethical approval for the study, and written, informed consent was obtained from the parents of the children.

MRI Data Acquisition

All participants were scanned at 10 ± 2 years and were asked to keep their eyes closed during the scans. MRI data were acquired using a Sigma HDx 3 T MR scanner (GE Healthcare). The MRI protocol included a sagittal 3D-T1 weighted with a BRAVO SPGR sequence: time to inversion = 400 ms, field of vision = 240 × 240 mm2; flip angle = 12°; voxel size 1 × 0.938 × 0.938 mm3; slice thickness = 1.0 mm. Resting-state functional data were acquired, namely: GE-EPI, time repetition(TR)/time echo (TE) = 2000/30 ms; flip angle = 70°; voxel size 3.0 × 3.0 × 3.5 mm3 with full-brain coverage; 300 volumes. dMRI (twice-refocused SE-EPI, TR/TE = 6300/72 ms; voxel size 2.3 × 2.3 × 2.3 mm3; 60 slices for full-brain coverage; ASSET factor 2) was acquired with a multi-shell acquisition with 4 b = 0, 30 b = 1000 s/mm2, 60 b = 2500 s/mm2, including an additional 2 b = 0 volumes with reversed phase-encoding polarity.

Functional MRI Preprocessing and Motion Censoring

Image preprocessing was carried out using FSL software, version 5.0.5 (FMRIB). All the raw resting-state functional MRI (fMRI) data and outputs from each preprocessing step were examined visually. Extreme motion data were eliminated from the beginning. The preprocessing included three elements. The first was volume realignment with the MCFLIRT intra-modal motion correction tool (FMRIB) (Jenkinson et al. 2002). A distortion correction with field map, namely reversed phase-encoding from the dMRI, was applied. The second was co-registration of functional data to the high-resolution structural image using boundary-based registration (Greve and Fischl 2009). The third was registration from high-resolution structural data to the study-specific template using the FSL FNIRT nonlinear image registration tool (FMRIB). We calculated outlier volumes and root mean square frame displacement for each participant using the FSL motion outliers tool (FMRIB) (Dipasquale et al. 2017). We also used that tool to calculate DVARS, where D was the temporal derivative of time courses and VARS referred to the root mean square variance over voxels. The number of outliers (>75th percentile + 1.5 interquartile range) was also recorded. We used a frame displacement threshold of 0.2 mm to identify potentially motion-contaminated scans (Power et al. 2015). Runs that had more than 15% of outlier volumes in the frame displacement or DVARS, with more than 0.2 mm of motion displacement, were flagged up as being potentially contaminated by motion and were excluded from the study. We did not use scrubbing in our data, since synthetic data and temporal discontinuities could be introduced by data interpolation and censoring, respectively, disrupting the temporal correlation of the signal (Caballero-Gaudes and Reynolds 2017). Subsequently, we used aggressive ICA-AROMA, which provides ICA-based automatic removal of motion artifacts, to identify and remove residual motion-related artifacts (Pruim et al. 2015). Components identified as head motion were removed from the signal by using linear regression, namely non-aggressive denoising. Mean white matter and cerebrospinal fluid signals were computed by averaging the values within the masks that were derived from the segmentation of each subject’s structural image. These masks were eroded to prevent the inclusion of gray matter signals via partial-volume effects (Ciric et al. 2017). Lastly, a band-pass temporal filter with a cut-off frequency of 0.01–0.07 Hz was used.

dMRI Preprocessing

Prior to any processing, visual inspection of all the raw diffusion data was carried out to detect excessive motion and/or signal dropouts. We estimated and corrected susceptibility and eddy current–induced image distortions and performed motion correction using the eddy tool in FSL (Andersson and Sotiropoulos 2016) and the reversed phase-encoding data set. The intra-volume movements were successfully corrected using the slice-to-volume correction routine in the eddy tool (Andersson et al. 2017). FSL’s higher-order diffusion model bedpostx (FMRIB) (Jbabdi et al. 2012) was used to estimate crossing fibers within each voxel of the brain, and probabilistic whole-brain fiber tractography was performed using FSL’s probtrackx (FMRIB) (Behrens et al. 2007) (see details in the Supplementary Material).

Brain Network Construction

First, both diffusion and functional data were registered to the T1-weighted structural image, then to the population template. The resulting transformations were concatenated, inversed, and used to warp the automatic anatomical labeling atlas with 90 anatomical regions (Tzourio-Mazoyer et al. 2002) from the population template space to the individual native space using nearest-neighbor interpolation.

A structural connectivity matrix for each participant was created, by calculating the connections (edges) between each of the 90 areas (nodes) using probabilistic tractography as described previously. The resulting structural connectivity matrices were then averaged across the groups. The list of the regions and abbreviations are provided in Supplementary Table S1. For functional data we used the extracted and averaged time courses (time-series length 10 min) of 90 cortical and subcortical brain areas, defined by the automatic anatomical labeling atlas, to compute subject-specific functional connectivity matrices. The functional connectivity matrices were computed with MATLAB and using Pearson’s linear correlation coefficient.

Intrinsic Ignition

We used the intrinsic ignition method that has previously been described in adults (Deco and Kringelbach 2017; Deco et al. 2017b; Deco et al. 2019). Intrinsic ignition corresponds to naturally occurring intrinsic perturbations in the brain. This method quantifies the capability of a given local node, namely a neuron or brain area, to propagate neural activity to other regions. This then leads to different levels of integration and determines the capacity of the whole network to become interconnected and exchange information. Each node drives or “ignites” global integration each time its local activity passes an amplitude threshold of one standard deviation (Tagliazucchi et al. 2012). The global integration of the functional connectivity matrix is then measured at that specific time point. More specifically, for a given brain region, we average across the events a measure of the integration elicited at time t relative to events. We defined the ignition event–driven integration of a given brain area as the averaged elicited integration during a time window of four TRs, which was determined by the time that takes the integration to return to basal values. (Deco and Kringelbach 2017) (Fig. 1). Integration is computed by measuring the area under the curve of the largest connected component of the functional connectivity matrix. This is computed at a range of thresholds, from 0 to 1 (Saenger et al. 2017). The larger this area is, the higher the integration of the network. It should be noted that integration is only measured when a node passes the activity threshold, and this means that the resulting ignition is paired with driving events. The integration values are normalized between 0 and 1, where 1 represents a fully connected network and 0 denotes a fully disconnected network. The mean of all these paired integration events is called ignition, while the standard deviation refers to the variations across the ignition events. This means that the standard deviation indicates the variability of the local events that drive the different levels of integration (Deco et al. 2017b). Defining the variability of the ignition events across time makes it possible to characterize the degree of hierarchy of information processing (Deco and Kringelbach 2017). The regions with high standard deviations, and therefore high variability, correspond to the regions that are more computationally relevant, because they play a central role in broadcasting information that is located at the top of the hierarchy. In contrast, regions with a low standard deviation, and therefore low variability, are more likely to be related to local processing, which is located at the bottom at the hierarchy. The entire process was repeated for each subject in both groups, which generated subject-specific ignition and hierarchy vectors.

Figure 1

Overview of the Methodology. (a) Intrinsic ignition. The blood oxygen level–dependent (BOLD) neuroimaging signal of a previously defined brain region (A, green signal) can be treated as spiking data by applying a threshold method (see Materials and Methods) to define an ignition event. For each driving event the activity of the rest of the network is measured (B, read area) in the time window (gray rectangle). For each ignition event this gives rise to a connectivity matrix where integration can be measured. (b) Whole-brain computational model. SC, structural connectivity; FC, functional connectivity. SC data were computed from dMRI-based tractography. The anatomical labeling atlas was used to parcellate the dMRI data in 90 anatomical regions. Then, a structural connectivity matrix (90 × 90) for each participant was created. The functional data from fMRI-BOLD activity was also used to create an FC matrix (empirical) following similar aforementioned steps (see Materials and Methods). Based on the SC matrix, a computational model is built to construct a simulated FC matrix (FC model) that is fitted to the empirical FC matrix, using different model parameter combinations. The optimal fit corresponds to a minimal Euclidian distance. (c) Hopf bifurcation model. The local dynamics of each node is modeled by using a normal form of a Hopf bifurcation. The nodes can behave synchronously (oscillations), asynchronously (noise), or critically (noise + oscillations), and these characteristics are represented by a bifurcation alpha value (bifurcation parameter, BP) with positive, negative, and values around 0, respectively. The simulated signal looks like the empirical data when the nodes are at the border between noisy and oscillatory behavior (bifurcation point).

Whole-brain Model

The whole-brain model was based on the average structural connectivity matrix (empirical) of each group and was parcellated in 90 regions, namely nodes. The local nodes were coupled through the underlying structural connectivity matrix, which contains the fiber density between all pairs of brain areas. Based on this matrix, a Hopf whole-brain model was built in order to simulate the resting activity of each of the 90 brain nodes. The simulated functional connectivity, known as the functional connectivity model, was then fitted to the empirical functional connectivity matrix (Fig. 1).

Hopf Bifurcation Model

We explored the difference in the underlying global and local dynamics further by using a whole-brain Hopf bifurcation model with a whole-brain fiber tractography-based structural backbone network (Deco et al. 2017a). The activity of the brain networks is characterized by a mix between order and disorder. Order is expressed as cluster synchronization between a group of brain areas, and disorder is expressed as desynchronization between the remaining brain areas. Furthermore, those states of partial synchronization are metastable. This means that the brain is constantly switching between different states of partial synchronization, for example the group of areas that show synchronized changes across time. This state, which lies on the boundary between order and disorder, is called the critical state, and a system that is in this state is referred to as being “at criticality”. When the system is at the point of transition between the two different states, this is known as the “bifurcation point”.

In this context, the Hopf bifurcation model can be used to simulate the resting activity of each of the 90 nodes derived from the anatomical labeling atlas parcellation, as previously described in this paper. In short, this model can be used to identify how nodes behave locally in terms of bifurcation. The nodes can behave synchronously, asynchronously, or critically, and these characteristics are represented by a bifurcation “alpha” value with positive, negative, and values around 0, respectively (Fig. 1). The resting brain dynamics (functional connectivity model) that emerge from the coupled structural nodes (empirical, structural connectivity) simulated by the Hopf model use two parameters to scale and approximate simulated empirical activity. The first is a global coupling strength scaling factor (G) that determines the impact of the structural connectivity in brain dynamics; a higher optimal G means that modeled functional connectivity is more constrained by the structural connectivity and a lower optimal G means that structural information becomes less relevant. The second is the local (nodal) bifurcation parameter alpha. This local critical behavior is approximated by using a gradient-descent algorithm that reduces the error between the modeled and empirical power spectrum in each node. Once this error is lower than 10%, then the final bifurcation value for each node is fixed. Those values are then considered to be optimized, because they represent the optimal power spectrum in the empirical data (Donnelly-Kehoe et al. 2019). This means that they can be used to compare the dynamic local (Saenger et al. 2017) and global (Jobst et al. 2017) differences in brain activity between the groups. More information about this procedure can be found in the Supplementary Material and in previously published studies (Deco et al. 2017a; Donnelly-Kehoe et al. 2019; Saenger et al. 2017).

We used two metrics to measure the level of agreement between the modeled and empirical data. The Kolmogorov–Smirnov distance was used to reflect the similarity of intrinsic fluctuations (dynamic functional fitting), while the structural fitting measured the Pearson correlation between the modeled and empirical functional connectivity matrices. Given that the relevance of dynamic rather than static information was prioritized in the current study, the most important score is the Kolmogorov–Smirnov distance, which is better as it gets closer to 0, as it portrays dynamic information. (Donnelly-Kehoe et al. 2019). Further details of this method have previously been described (Deco et al. 2017a; Jobst et al. 2017; Saenger et al. 2017) and are also provided in the whole-brain modeling section of the Supplementary Material.

Statistical Analysis

Data were analyzed using SPSS version 20 (IBM Corp.). The independent samples t-test and Wilcoxon rank-sum test (equivalent to a Mann–Whitney U test) were applied, where appropriate, to assess group differences, and two-sided P < 0.05 was considered statistically significant. When it came to the bifurcation analyses, we applied a permutation test to investigate statistical differences between the distributions of bifurcation values between both groups. This was because the distribution of bifurcation values is highly non-normal and can present binomial characteristics (Saenger et al. 2017). The first step was to compute three distribution statistics for the bifurcation values for each of the two groups. “Kurtosis” captured the shape of the distribution, “second-order raw moment” captured the data dispersion from zero, and finally the “median” separated the lower values from the higher half of the distribution. We computed the absolute difference (delta) for each of these three distribution statistics, so that we could quantitatively represent the differences in the distribution characteristics between the term and preterm groups. Once the observed absolute differences were computed for the three statistics, the bifurcation values of both of the distributions were pooled together. This process was repeated 10 000× to finally compute the significance, which was when the proportion of delta random was larger than the observed delta for each statistic (Saenger et al. 2017). Finally, we tested the relationship between the ignition and bifurcation distributions by using Spearman’s correlation test.

Results

Characteristics of the Population

This study comprised 33 children born EPT and 28 children born full term with a mean age of 10.14 ± 0.81 versus 9.98 ± 0.90 years, respectively. As expected, there were significant differences between the two groups in terms of gestational age and birth weight. The two groups did not differ in terms of motion parameters, namely translation and rotation (Table 1). The characteristics of the two cohorts are summarized in Table 1, and the neonatal morbidity of the preterm group is shown in the Supplementary Table S2.

Table 1

Characteristics of the children born EPT and full term

Characteristic	EPT	Term	Statistic (P)	EPT	Term	Statistic (P)
	Ignition analysis Value ± SD			Hopf model analysis Value ± SD
	N = 33	N = 26		N = 28	N = 21
Perinatal Data
Gestational age (weeks)	25.7 ± 0.9	39.9 ± 1.1	t − 52.8 (<0.001)	25.7 ± 0.9	39.8 ± 1.2	t − 45.9 (<0.001)
Range	(23.5–26.6)	(37.3–41.5)		(23.5–26,6)	(37.3–41.5)
Birth weight (g)	856.2 ± 173.4	3663 ± 421.0	t − 31.8 (<0.001)	854.6 ± 178.7	3600 ± 387.0	t − 31.8 (<0.001)
Range	[550–1161]	[2.8754100]		[550–1161]	[2875–4100]
Gender (boy/girl)	13/20	14/12	Fisher’s test (0.17)	11/17	9/12	Fisher’s test (0.17)
Age at MRI (years)	10.0 ± 0.8	9.9 ± 0.9	t 0.65 (0.52)	10.14 (0.82)	9.90 (0.92)	t 0.94 (0.34)
Range	9.0–11.5	8.0–11.8		9.0–11.5	8.0–11.8
Motion estimation after AROMA
Average FD-RMS mean SD	0.10 ± 0.05	0.08 ± 0.04	t 1.44 (0.15)	0.10 ± 0.06	0.11 ± 0.12	t − 0.41 (0.68)
Average DVARS-RMS mean SD	5.97 ± 0.84	5.84 ± 0.91	t 0.23 (0.81)	5.96 ± 0.98	5.99 ± 1.05	t − 0.10 (0.91)
Number of frame outliers FD-RMS	19.9 ± 10.12	19.52 ± 9.5	t 0.14 (0.88)	17.57 ± 9.39	19.75 ± 10.27	t − 0.76 (0.45)
Number of frame outliers DVARS	14.4 ± 6.96	14.4 ± 9.14	t 0.02 (0.98)	14.61 ± 5.92	17.14 (1247	t 0.94 (0.35)
Rotation X	0.11 ± 1.78	0.45 ± 1.61	t − 0.73 (0.46)	0.46 ± 1.71	−0.42 ± 1.76	t 1.0 (0.31)
Rotation Y	−0.26 ± 0.60	0.13 ± 0.54	t − 1.03 (0.30)	0.04 ± 0.57	−0.13 ± 0.54	t 0.34 (0.73)
Rotation Z	0.25 ± 1.53	−0.57 ± 1.69	t 1.87 (0.06)	0.30 ± 1.60	−0.60 ± 1.67	t 1.95 (0.06)
Translation X	−0.03 ± 0.19	0.02 ± 0.12	t − 1.29 (0.20)	−0.05 ± 0.25	0.02 ± 0.11	t − 0.54 (0.58)
Translation Y	−0.10 ± 0.48	0.003 ± 0.13	t − 1.09 (0.27)	−0.14 ± 0.53	0.02 ± 0.11	t − 1.10 (0.27)
Translation Z	−0.007 ± 0.22	0.009 ± 0.13	t − 0.31 (0.75)	−0.02 ± 0.29	0.02 ± 0.12	t − 0.26 (0.98)

FD, frame displacement; RMS, root mean square; DVARS, D is temporal derivative of time courses and VARS is RMS variance over voxels.

The mean value of the ignition events across the brain areas was significantly different between the two groups (preterm 0.708, term 0.710, Wilcoxon rank test, P < 0.0001, z = 4.35) (Fig. 2). When the top 10 nodes with the highest values were ranked, the two groups shared some of the regions with the highest ignition values, namely the orbito-frontal superior, left gyrus rectus, right superior temporal-pole, and occipital cortices. The regions with the highest ignition values in the preterm group included the bilateral anterior cingulate gyrus (affective regulation) and the globus pallidus (motor and non-motor behaviour). Meanwhile, the term group recorded the highest values in the left temporal inferior (high level of visual processing and object recognition), the right orbito-frontal middle, and the bilateral inferior cortices (cognitive processing of decision making) (Fig. 2). The preterm group had a smaller variability of ignition events across the brain regions compared with the term group (preterm 0.032, term 0.035, Wilcoxon rank test P = 0.01, z = −2.53) (Fig. 2). In order to define how important the brain nodes were for broadcasting information across the brain, we calculated the variability of the ignition events across the nodes in the two groups. The brain nodes with the highest variability in the preterm group included the right superior parietal cortex, the right gyrus rectus, the left precentral areas, the left middle temporal pole, and the bilateral superior occipital cortices (Fig. 2). In the term groups the highest variability nodes included the right parahippocampal gyrus, the left middle superior frontal gyrus, the left fusiform gyrus, the left precuneus, and the left paracentral cortices.

Figure 2

Ignition. (a) Ignition values across the brain and variability across events for the preterm (orange) and the term (green) groups. Standard error is represented at a shaded area. (b) Ranked (descending) ignition and variability values where top 10 regions (highlighted in gray within the plot) for both groups are described at the top right insert. Rendered brains show the top 10 regions with the highest values in ignition and variability across events.

Whole-brain Hopf Bifurcation Model

We simulated the resting-state fMRI blood oxygen level–dependent activity in each of the 90 nodes from the parcellation, explained above, by using the whole-brain Hopf model with the brain nodes coupled through the empirical structural connectivity matrix. The preterm group had a larger optimal G parameter when we used both the structural (0.175) and functional (0.45) fitting, compared with the term group (0.10, structural and 0.15, functional), which suggests that modeled preterm dynamics are more constrained by structural connectivity. This pattern can be observed in Figure 3, which shows that the ongoing simulated dynamics in the preterm group reached maximal similarity with the empirical dynamics at a higher coupling than in the term group. Coupling can be interpreted as a connectivity scaler, where higher coupling values for maximum similarity mean that more structural information is required to reach the optimal fit. This suggests that functional interactions were more directly related to the underlying structural connectivity and that this led to a more stringent and less dynamic brain function in the preterm group.

Figure 3

Fitting between empirical and simulated brain activity. Static (left) stands for the correlation between empirical and simulated functional connectivity by means of a Pearson correlation coefficient. Dynamic (right) represents the Kolmogorov–Smirnov distance between the distributions of ongoing dynamics. The closer simulated and empirical dynamics are the lower the distance. The optimal coupling values are described in the top right corner of each plot. Preterm group (orange) and term group (green). Solid lines represent mean fitting values, while shadowed areas represent the standard deviation over simulations (as many as subjects per group).

The second parameter of the model, the bifurcation parameter distribution, was then determined between the two groups. The median of the bifurcation parameters at the working point was −0.025 in the preterm group and 0.015 in the term group (Fig. 4). We used a permutation test to explore if these qualitative differences were reflected by robust statistical differences, and this found significant differences in the medians between the distributions (observed = 0.0301, P = 0.0001) and moment (observed = 0.0027, P = 0.0005) (Supplementary Fig. 2). These results indicate that the distribution in the term group tended to be more positive and significantly closer to the bifurcation point than the preterm group (Supplementary Fig. 2). In contrast, the bifurcation parameter distribution in the preterm group was significantly less centered around zero, with a more negative distribution in alpha values. The highest absolute differences between the two groups were found in the bilateral middle temporal gyri, right paracentral gyrus, bilateral orbital middle frontal gyrus, left globus pallidus, right inferior occipital gyrus, right Heschl’s gyrus, left temporal pole, and right superior temporal gyrus (Fig. 4 and c).

Figure 4

Whole-brain Hopf bifurcation model. (a) Bar plot of alpha bifurcation values across the brain for preterm (orange) and term (green) groups. The third blue plot shows the absolute difference (term and preterm) between alpha bifurcation values on each region. Asterisks show the top 10 regions with highest difference. (b) Ranked absolute differences between the term and preterm alpha vectors. The top 10 regions presenting the highest difference are marked with red and described in the top right insert and (c) displayed in the rendered brains

Correlation between the Highest Differences in Alpha Bifurcation and Ignition Values

Most of the regions with the highest differences in ignition values between the two groups also had the highest differences in the Hopf bifurcation analysis. Despite that, some nodes appeared exclusively in the ignition data, and these were the left inferior temporal gyrus, right gyrus rectus, right calcarine, and right inferior orbito-frontal gyrus. Meanwhile, others appeared exclusively in the Hopf data, namely the right paracentral gyrus, left globus pallidus, right Heschl’s gyrus, and right superior temporal gyrus. We found a positive significant correlation (Spearman’s rho = 0.207; P = 0.049) between the ignition and bifurcation distributions (Fig. 5).

Figure 5

Differences in bifurcation (alpha) and ignition values. (a) and (c) represent rendered brains showing the top 10 regions with the highest differences (delta) in alpha bifurcation and ignition values. (b) Plot representing the ranked delta in alpha (red) and ignition values (green), with the top 10 regions representing the highest differences described in the insert. (d) Scatter plot between delta alpha and delta ignition with a fitted linear regression (straight gray line) and 95% confidence bounds (dotted lines). Spearman’s rank correlation coefficient between these two variables is indicated by r.

Discussion

This study used the novel concept of intrinsic ignition to characterize the dynamical complexity of the propagation of the neural information in the brain of EPT-born children. We used a whole-brain Hopf bifurcation model to explore the local and global impact that EPT birth had on the brain, with regard to creating the critical states that provide the necessary conditions for cognitive processes. There were three main findings when we compared our 10-year-old EPT cohort with healthy term-born controls of the same age. The first was that EPT children had reduced intrinsic ignition events. The second was that we found a different hierarchy pattern in the processing of information. Thirdly, the EPT children experienced a suboptimal neural synchrony and criticality, with the largest differences between these subjects and the controls found in the brain nodes that belong to the rich-club architecture. The data collected by this study provide us with a unified framework that helps us to understand the neural substrates that underlie the brain reorganization and neurodevelopmental impairments associated with EPT birth.

Intrinsic Ignition Analysis

The EPT children had significantly lower intrinsic ignition events than the term group. As previously described, during the resting-state mode, brain regions tend to self-ignite, and they have the ability to start the propagation of neuronal activity to other regions over time, which drives global integration (Deco and Kringelbach 2017; Deco et al. 2017b). Despite the lack of comparable studies, our results are consistent with previous investigations that have shown poor structural integration in preterm children (Fischi-Gomez et al. 2016; Thompson et al. 2016), which could be related to reduced ignition events. This was in line with observations from a previous study using MEG during the resting state, which showed reduced spontaneous neuromagnetic activity in preterm children (Doesburg et al. 2011b). This could have been related to a reduced amount of spontaneous ignition events. During normal development, a progressive shift toward high-frequency oscillation continues throughout childhood. This means that reduced spontaneous activity, and therefore reduced ignition events, could be interpreted as a delay in the development of normal brain maturation, due to abnormal physiological (Kwon et al. 2014), neurostructural (Karolis et al. 2016; Munoz-Moreno et al. 2016), and functional developmental trajectories (Fransson et al. 2007; Smyser et al. 2016; Padilla et al. 2017b). If a brain area does not ignite well, it will have a limited ability to propagate neural activity to other areas of the brain, and this will reduce the integration and affect the efficiency of global neural communication. This could be critical for the generation of cognitive processes in this population (Ward 2003). It should be noted that these changes should be placed within a developmental framework, following a previous study that reported that the brains of individuals born very preterm displayed accelerated maturation from adolescence onwards (Karolis et al. 2017). Future studies are clearly needed to explore this issue in greater depth.

The regions with the highest ignition varied between the two groups. In the preterm group, the highest ignition regions included those previously reported as being involved in processing affective regulation, motor (planning, control, and behavior), and some complex aspects of cognitive function (semantic language and visual integration). In the term group, we found a more consistent pattern of high ignition in areas involved in the processing of more complex aspects of cognitive function, namely multisensory integration and working memory.

In the children born EPT, the hierarchy of information processing, based on the variability of the intrinsic ignition events, was predominantly driven by visual and sensory regions. This was different from the term group, where the fronto-temporal (higher-order regions) and associative areas (angular gyrus) were involved. A healthy brain follows a primary-to-higher order maturational sequence (Chomiak and Hu 2017; Friedrichs-Maeder et al. 2017), consistent with their behavioral correlates. This means that during development, higher-order processing areas, such as the fronto-temporal and precuneus regions, display increased maturational priority over the sensory–motor regions. The higher-order processing areas then become pivotal nodes for the regulation of neural activity and neurocognitive development (Kolskar et al. 2018). Our study found that EPT-born children exhibited a different pattern of intrinsic functional organization, where the brain nodes with the highest intrinsic ignition values were predominantly located in sensory–motor, limbic and paralimbic regions, and somatosensory association cortices. These results were in line with previous studies that showed how prematurity disturbed the typical hierarchical pattern of structural and functional cortical maturation in children (Lax et al. 2013; Fischi-Gomez et al. 2015). The brains of children born EPT have malfunctioning or replaced core components that transfer information from one level to another (Karolis et al. 2016). They may also be affected by jeopardized hierarchical organization, leading to functional alterations (Crossley et al. 2014). Importantly, aberrant hierarchical organization has been described in autism (Parr et al. 2018), which is a highly prevalent disorder among children born EPT (Padilla et al. 2017a).

Computational models enabled us to study the link between anatomical structure and resting-state network dynamics. We found that the EPT group showed significantly lower synchrony than the term controls. Our results are in line with previous studies that used MEG to show reduced oscillatory synchrony when very preterm children performed a task (Moiseev et al. 2015) and during their resting state (Doesburg et al. 2011a; Ye et al. 2016). The mechanisms underlying the reduced synchrony that we observed could have been related to atypical development of white matter in the EPT group. Functional networks are extensively reconfigured during late childhood and adolescence (Supekar et al. 2009; Grayson and Fair 2017). These phenomena involve progressive and regressive neurobiological changes in myelination, axonal diameter, remodeling of synapses (Hagmann et al. 2010), and neurochemical changes (Crews et al. 2007), which are reflected in the maturation of the frequency and synchronization of oscillatory activity (Uhlhaas et al. 2010). Extreme preterm birth may disrupt these processes, leading to abnormal neural synchrony (Moiseev et al. 2015; Ye et al. 2016), which has important implications for the cognitive development of this population (Doesburg et al. 2010).

The bifurcation parameter distribution in the preterm group was significantly less likely to be centered around zero than the term controls. Distributions close to zero are usually referred to as healthy (van Hartevelt et al. 2014; Saenger et al. 2017). A brain with nodes that work at the edge of criticality could be considered to be a flexible system, as this makes it possible for the brain to adapt rapidly and efficiently to external and internal demands (Deco et al. 2017a). Extreme prematurity may compromise this adaptability by disturbing the ability of the cortical networks to reconfigure themselves across time, with further detrimental effects on supporting processing demands (Deco et al. 2011). This was reflected in a higher optimal G in the EPT group, suggesting that the functional connectivity in the model was more constrained by the structural connectivity, as previously suggested in adults (van Hartevelt et al. 2014; Saenger et al. 2017). Parameter G is a global scaling factor, which means that the global conductivity parameter scales equally for all synaptic connections. In other words, the G parameter scales the underling white matter to fit the functional data, and a larger G parameter suggests a possible weakening of the underlying white matter, which in turn could be interpreted as a relative paucity of white matter resources (Karolis et al. 2016), reflecting weaker connectivity in the EPT group as previously demonstrated (Fischi-Gomez et al. 2016; Thompson et al. 2016). In contrast, proper fitting values were slightly higher in the EPT group (both static and dynamic fit; Fig. 3), which might be due to the larger data set in this group (33 compared with 26 subjects), allowing the model to reach higher similarity with empirical data. Future studies are clearly needed to explore this issue.

The regions with the largest bifurcation parameter differences between the two groups overlapped with those hubs described in the rich-club architecture, which is a highly central connected backbone that is thought to enable efficient network communication in the brain (Ball et al. 2014; Karolis et al. 2016). Overall, these nodes were asynchronous in the EPT group, pushing the brain further away from a critical regime where it works properly (van Hartevelt et al. 2014; Saenger et al. 2017). One of the most asynchronous regions was the globus pallidum, which is of interest in premature infants, as their rich-club architecture is preserved at the expense of disruptions in cortico-subcortical connectivity (Ball et al. 2014; Karolis et al. 2016). We could speculate that the set of nodes that are part of the dynamic repertoire might be conditioned by the brain reorganization that takes place when an infant is born premature. In turn, this could be compensatory (Brittain et al. 2014) or could be central to mediating a wide spectrum (Uhlhaas and Singer 2012) of cognitive deficits (Twilhaar et al. 2018), motor deficits (Bolk et al. 2018), and neuropsychiatric disorders, such as autism (Padilla et al. 2017a). Studies with larger samples are needed to further characterize the nature and extent of these changes.

The correlation between the differences in bifurcation parameter and ignition values between the preterm and term-born groups revealed the degree to which these two complementary dynamical properties are interrelated in the complex dynamics of the brain (Deco et al. 2017b). These findings indicate that positive bifurcation dynamics, which reflect stable oscillations, indicate higher ignition values, while noisy dynamics reflect lower ignition. These results should provide the grounds for future studies that carefully explore these relationships.

Strengths and Limitations

The strengths of this study included the methodology, which combined the intrinsic ignition method with a whole-brain computational model not previously used in children born EPT. In addition, the cohort we studied was well defined, the EPT group did not display major brain alterations on the MRI scans, and their mean gestational age at birth was lower than in previously published studies. Although the effect of motion artifacts was an issue, we undertook rigorous processing to minimize any effects this may have had on the data and excluded participants if the motion artifacts were too severe to correct. Although the neurodevelopmental data for the children included here are not available yet, further studies including these data would help to determine the significance of the results found in the EPT group in our study. In addition, we are aware of the valid concern that we cannot make causal inferences about model parameters and neurobiological signals using computational approaches for the data obtained during MRI brain scans. However, computational formulations offer a detailed and quantitatively precise characterization of processes that are not accessible with other methods, suggesting more sophisticated explanatory frameworks.

Conclusion

This study compared the functional organization of the resting brain in children born EPT by comparing them with term-born controls at 10 years of age. We found that the EPT brains were characterized by reduced intrinsic ignition events, a different hierarchy pattern in the processing of information and suboptimal neural synchrony and criticality. The largest differences between the two groups were found in functional brain nodes overlapping the rich-club architecture. These results provide important insights into the neural substrates that underlie brain reorganization and neurodevelopmental impairments related to extreme prematurity. They also identify areas of the brain that could be strengthened by tailor-made interventions, so that clinicians can optimize the functions and quality of life of vulnerable children following EPT birth.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.