Atypical neuromagnetic resting activity associated with thalamic volume and cognitive outcome in very preterm children.

Children born very preterm, even in the absence of overt brain injury or major impairment, are at increased risk of cognitive difficulties. This risk is associated with developmental disruptions of the thalamocortical system during critical periods while in the neonatal intensive care unit. The thalamus is an important structure that not only relays sensory information but acts as a hub for integration of cortical activity which regulates cortical power across a range of frequencies. In this study, we investigate the association between atypical power at rest in children born very preterm at school age using magnetoencephalography (MEG), neurocognitive function and structural alterations related to the thalamus using MRI. Our results indicate that children born extremely preterm have higher power at slow frequencies (delta and theta) and lower power at faster frequencies (alpha and beta), compared to controls born full-term. A similar pattern of spectral power was found to be associated with poorer neurocognitive outcomes, as well as with normalized T1 intensity and the volume of the thalamus. Overall, this study provides evidence regarding relations between structural alterations related to very preterm birth, atypical oscillatory power at rest and neurocognitive difficulties at school-age children born very preterm.

Introduction

Worldwide, the incidence of preterm birth with a very low gestational age (VLGA, born ≤ 32 weeks of GA) is 1 in 100, and extremely low gestational age (ELGA, born ≤ 28 weeks of GA) is 1 in 200, with a mortality rate of 5–10% and > 10%, respectively (Torchin et al., 2015). Even in the absence of overt brain lesions, very preterm infants are at higher risk of cognitive, behavioral and motor problems (Aarnoudse-Moens et al., 2009, Anderson, 2014, Moore et al., 2012) compared to healthy full-term children. Risk factors in very preterm infants (ELGA & VLGA) include post-natal infection (Ranger et al., 2015, Zwicker et al., 2016), as well as exposure to numerous daily invasive procedures such as blood draws and line insertions, during weeks to months in the neonatal intensive care unit (NICU) (Roofthooft et al., 2014, Vinall and Grunau, 2014). The disruption of normal intrauterine maturation by preterm birth, neonatal clinical factors, and pain-related stress alters ongoing brain development during critical periods involving prolific axonal growth, dendritic sprouting, and synapse formation (Mrzljak et al., 1992), that leads to abnormal brain morphology and activity (Kapellou et al., 2006, Schneider et al., 2018, Smith et al., 2011). Thus, to ensure proper extrauterine environment and adequate preterm interventions, it is of vital importance to understand how structural and functional preterm alterations linked to environmental stressors are associated with cognitive and behavioral outcomes in children born very preterm.

During the neonatal period in very preterm infants born 24-32 weeks of gestation, thalamic afferents are synapsing in the subplate, a transient structure that generates endogenous activity and is critical to the formation of long term thalamocortical circuitry (Kostovic and Judas, 2002), gradually synapse onto cortical neurons and become sensory driven. Disruption in this developmental process might cause neural apoptosis (Anand et al., 2007, Dührsen et al., 2013) and altered development of thalamocortical axons (Dean et al., 2013, Molnar and Rutherford, 2013). Moreover, very preterm infants present abnormal myelination and fractional anisotropy in gray and white matter areas (Dubois et al., 2008, Eaton-Rosen et al., 2015), decreased cortical and subcortical gray matter volume (Boardman et al., 2006, Chau et al., 2019), thalamocortical alterations (Ball et al., 2012, Cai et al., 2017) and pain-induced volume reduction of white matter and subcortical gray matter (Brummelte et al., 2012, Duerden et al., 2018). These structural alterations persist into school age (Hohmeister et al., 2010, Lax et al., 2013) and adulthood (Menegaux et al., 2017, Nosarti et al., 2014) and are associated with lower IQ (Breeman et al., 2017, Nosarti et al., 2014).

The thalamus acts as a sensory relay by directing peripherical information to the cortex and as an integrative hub for cortical representations by mediating cortico-cortical communication (Hwang et al., 2017, Jones, 1998, Malekmohammadi et al., 2015, Sherman and Guillery, 2013). It has been associated with numerous cognitive functions (Jones, 2012, Sherman, 2016, Theyel et al., 2010) and is affected in many disorders (Giraldo-Chica et al., 2018, Llinás et al., 1999, Llinás et al., 2005, Ribary et al., 2019, Tona et al., 2014). Functionally, the thalamus regulates cortical power at different frequencies and their inter-relationships during cognition and perception (Ribary et al., 2017). Neurological conditions are often associated with the thalamus slowing alpha oscillatory activity to theta (Li et al., 2015, Ribary et al., 2019, Tarapore et al., 2013), and in very preterm children at school age decreased alpha power and increased gamma have been reported (Cepeda et al., 2007, Doesburg et al., 2010, Doesburg et al., 2011), as well as atypical neurophysiological connectivity (Doesburg et al., 2013, Kozhemiako et al., 2019b, Moiseev et al., 2015b).

Based on previous literature reporting associations between the thalamus and cortical power (Edgar et al., 2015, Lindgren et al., 1999, Ribary et al., 2017), in the present study, we hypothesized that reductions in thalamic volume present in children born very preterm will be associated with atypical cortical power measured using MEG. Given that power captured by MEG is mostly generated by pyramidal neurons in the cortex (Murakami and Okada, 2006), when investigating the association with thalamic volume and power, the cortical gray matter volume was included as a variable of interest. Previous studies indicated delayed neural pruning in children born preterm (Ment et al., 2009, Mürner-Lavanchy et al., 2014), and this could be explained by diminished thalamic regulatory activity. To this end, gray matter volume was included when investigating the relation with power and thalamic volume. Higher cortical gray matter and lower thalamic volume associated with atypical frequency-specific power could suggest that a functional reduction of the thalamus in controlling cortical oscillations due to reduced number of neurons or number of neural connections.

To further understand the effects of thalamic alterations, specifically, the possible reduced amount of neural connections in the thalamus, we explored the relation between thalamic myelination level and frequency-specific neurophysiological power at rest. Previous studies have suggested that T1 intensities reflect myelination (Koenig, 1991, Stüber et al., 2014). Gray matter intensity has been used as a biomarker for neurological diseases and aging (Kong et al., 2015, Norbom et al., 2018, Salat et al., 2009), and myelin alterations in the thalamus and cortical gray matter in very preterm infants has been reported (Counsell et al., 2002, Counsell et al., 2014, Eaton-Rosen et al., 2015). In the present study, mean thalamic normalized intensity was used to assess associations between thalamic myelination and neuromagnetic power while accounting for cortical myelinization reflected by mean cortical gray and white matter intensities. In addition, given the risk of neurocognitive deficits in the very preterm, we tested the relationship between cortical oscillatory power and cognitive, perceptual and behavioral outcomes, and with neonatal procedures such as gestational age (GA), illness severity, infections, morphine doses and skin-breaking procedures.

Methods

Participants

A total of 108 children participated in a resting state MEG recording. Data from 9 participants were discarded due to major brain injury (periventricular leukomalacia [PVL] or grade III-IV intraventricular hemorrhage [IVH] on neonatal ultrasound), autism spectrum disorder [ASD], and/or or excessive motion. The final sample size was 99 children with a mean age of 7.8 years: 23 were born ELGA (24 to 28 wks, age 7.7 ± 0.39, 10 girls), 36 VLGA (29–32 wks GA, age 7.7 ± 0.39, 24 girls), and 39 healthy full-term (40 wks GA, age 7.9 ± 1.02, 24 girls). Participants were born at British Columbia’s Women’s Hospital between 2000 and 2004.

Data were collected during two separate sessions, one for a MEG recording and neurocognitive assessments, and one for a MRI recording. Of the 99 children that participated in the first session, only 62 underwent MRI. After discarding scans with poor quality due motion induced blurriness and ringing, 51 participants remained, consisting of 13 ELGA, 24 VLGA and 13 full-term children. All children were recruited as part of a prospective longitudinal study on the effects of neonatal pain-related stress on neurodevelopment of children born very preterm (Grunau et al., 2007, Grunau et al., 2009). The current cohort is the largest cohort of very preterm children followed from birth with both MEG and MRI at school age, making the results valuable and unique despite the limited number of participants with MRI. Participants’ sociodemographic, neurocognitive, neonatal and structural characteristics are presented in Table 1.

Table 1

Study cohort characteristics. Significant differences between ELGA or VLGA and full-term *p ≤ 0.05 ,**p ≤ 0.001. Significant differences between ELGA and VLGA †p ≤ 0.05, †† p ≤ 0.001. Mother and father ethnicity categories: 1 = African/Afro-American, 2 = Hispanic, 3 = Asian, 4 = Caucasian, 7 = Native American. A complete study on the psychocognitive and social effects on this preterm cohort can be found in (Grunau and Whitfield, 2004).

	ELGA	VLGA	Term
Demographic measures (MEG) N = 99	mean (std)	mean (std)	mean (std)
Age	7.7 (0.38)	7.7 (0.39)	8 (1.02)
Males/Females	13/10	12/24	15/24
Mother age	40 (4.6)	41 (5.1)	42.3 (5.1)
Father age	43.3 (5.9)	41.8 (5) *	44.4 (5.3)
Mother years education	15.5 (2.5) *	15.7 (2.7) *	18.6 (4.6)
Father year of education	15.3 (3.8) *	15 (3) *	17.7 (4)
Number of children	2.3 (0.75)	2.4 (0.91)	2.3 (0.67)
Birth order	1.66 (0.86)	1.8(0.82)	1.8 (0.67)
	mode (%)	mode (%)	mode (%)
Mother ethnicity	4 (83%)	4 (85%)	4 (90%)
Father ethnicity	4 (75%)	4 (83%)	4 (90%)
Demographic measures (with MRI) N = 50	mean (std)	mean (std)	mean (std)
Age, years	7.7 (0.42)	7.7 (0.39)	7.7 (0.60)
Males/Females	8/5	5/19	5/8
Neurocognitive measures N = 98
Verbal Comprehension Index (WISC-IV)	93.8 (12.15) **†	103.1 (15.16)	108.7 (13.11)
Perceptual Reasoning Index (WISC-IV)	97.0 (14.47) **	103.4 (15.93) *	113.1 (13.31)
Working Memory Index (WISC-IV)	91.8 (10.79) **†	101.8 (11.51)	102.2 (11.26)
Processing Speed Index (WISC-IV)	89.6 (11.82) **†	98.5 (12.97)	104.8 (16.03)
Full-scale IQ (WISC-IV)	91.2 (11.15) **†	102.4 (14.31) *	109.9 (12.76)
Internalising Behavior (CBCL)	51.4 (10.94)	50.7 (10.70)	49.7 (10.88)
Externalizing Behavior (CBCL)	47.6 (10.23)	47.2 (10.58)	47.2 (10.55)
Behavioral Regulation Index (BRIEF)	51.3 (10.72)	51.8 (12.82)	50.2 (10.95)
Metacognition Index (BRIEF)	53.7 (12.98)	54.6 (14.00)	49.4 (10.48)
Visual Motor Integration (BEERY)	92.2 (9.88) *	95.4 (8.31) *	102.4 (12.44)
Visual Perception (BEERY)	97.7 (15.68) **	104.9 (14.53) *	113.8 (15.31)
Motor Coordination (BEERY)	89.5 (10.39) *	93.4 (9.35) *	97.7 (11.20)
Neonatal measures (with MRI) N = 36
Gestational Age, weeks	27.1 (1.36) **††	31(1.2) **	39.9(0.93)
Infection, yes/no	7/6	20/3	n/a
Number of Skin-Breaking Procedures	153.6 (68.68) ††	65.1 (50.87)	n/a
Morphine Dosage, mg/kg	2.7 (5.74) †	0.06 (0.18)	n/a
Score for Neonatal Acute Physiology	22.5 (9.73) ††	5.6 (7.5)	n/a
Structural measures (with MRI) N = 50
Thalamic mean volume, mm3	6,951 (6 9 7)**†	7,494 (7 4 5)	8,268 (6 0 3)
Cortical gray matter volume, mm3	536,606 (40,984)*	555,204 (57,458)	589,950 (6 0 3)
Intracranial volume, mm3	1,515,881 (122,911)	1,534,179 (142,595)	1,589,423 (58,071)
Thalamic mean intensity	88 (1.8)*	89.2 (1.6)	89.5 (1.5)
cortical gray matter mean intensity	88 (8.7)	91.54 (8.2)	88 (6.7)
cortical white matter mean intensity	108 (79)	114.73 (10.3)	110 (8.5)

MEG recording

Eyes open resting state MEG data were recorded for a total of two minutes using a CTF 151-Channel MEG system (CTF Systems; Coquitlam, Canada) in a magnetically shielded room with a sampling rate of 1200 Hz. Participants were in a supine position and instructed not to fall sleep and to gaze at a centrally presented fixation stimulus while minimizing eye movements and blinks. Continuous head localization was recorded by energizing three fiducial coils placed in the nasion and preauriculas. The head shape of the participants was recorded using a Polhemus Fastrak digitizer.

MRI recording

MRI scans were performed on a Siemens 1.5 Tesla Avanto (Berlin, Germany) using a 12-channel head coil. Each scan consisted of a 3D T1-weighted SPGR sequence 18msec/9.2msec/256/1 mm/0/256 × 256 (TR/TE/FOV/Thickness/Gap/Matrix).

Neonatal data

Neonatal data were collected from a daily chart review by an experienced research nurse during the neonatal period of the ELGA and VLGA participants as described previously (Grunau et al., 2009). In this study, we focused on gestational age (GA), if infection was positive, illness severity on day 1 (SNAP; (Richardson et al., 2001)), equivalent log-transformed cumulative morphine dose, and log-transformed number of skin-breaking procedures (pain), from birth to term-equivalent age.

Neurocognitive assessment

On the day of the MEG scanning, IQ was assessed with the Wechsler Intelligence Scale for Children (WISC-IV; Wechsler, 2003), yielding standardized scores for verbal (VIQ), perceptual reasoning (IQ-PR), working (IQ-WM), processing speed (IQ-PSI) and full-scale (IQ-f). Visual-motor capabilities were assessed with Beery–Buktenica Developmental Test of Visual-Motor Integration, 5th Edn. (Beery et al., 2004), comprising the subscales visual-motor integration (BEERY-VMI), visual perception (BEERY-VP), and motor coordination (BEERY-MC). Behavior was assessed with the Child Behavior Checklist (CBCL; Achenbach and Rescorla 2001) questionnaire completed by a parent, measuring internalizing behavior (CBCL-INT and externalizing behavior (CBCL-EXT), and executive functions with the Behavior Rating Inventory of Executive Function (BRIEF; (Gioia et al., 2000) assessing behavioral regulation (BRIEF-BRI) and metacognition index (BRIEF-MCI).

MEG analysis

MEG sensor data were notch filtered at 60 Hz to remove line noise. Segments exceeding 5 mm displacement at any direction from the median head position were discarded. After visual inspection, segments with muscle artifacts were also discarded. Independent Component Analysis (ICA) was used to identify and reject components capturing eye and heart artifacts. The remaining artifact-free data was segmented into epochs of four seconds (trial number mean = 25, std = 4.5) and band-pass filtered at canonical frequency bands (delta:1–4 Hz, theta:4–8 Hz, alpha:8–12 Hz, beta:12–25 Hz, gamma:25–55 Hz). There were no group differences in motion or number of trials retained for analysis between groups. For source reconstruction, a forward model was computed using a single-shell head model (Nolte, 2003). For subjects without an MRI, the best match from a pool of child MRIs was used instead (Gohel et al., 2017). Source space activity was reconstructed on an 8 mm spaced grid using an LCMV beamformer with a regularized sensor covariance (Van Veen et al., 1997) which creates spatial filter weights that maximize activity from the target location while filtering out activity from elsewhere. Beamformer weights were estimated at each frequency band to optimize the spatial filtering (Moiseev et al., 2015a, Nunes et al., 2020). For each reconstructed source, relative power was obtained by dividing the absolute power within a frequency band by the sum of the absolute power in all the frequencies. Relative power is used as it is corrected for depth bias and it has been shown to be more stable across subjects (Suppiej et al., 2017). The relative power within a frequency represents its contribution to all the frequency spectrum of interest. Herein, power and relative power is used interchangeably. The analysis was performed using Fieldtrip toolbox (Oostenveld et al., 2011) and in house MATLAB scripts.

Morphometric measures

3D T1 MRIs were automatically processed using FreeSurfer (Fischl, 2012), in brief, steps included skull-stripping, Talairach transformation, intensity normalization, subcortical segmentation, tessellation of the gray matter/white matter boundary, topology correction and surface deformation to detect gray matter/white matter and gray matter/cerebrospinal fluid boundaries (Fischl et al., 1999). Some MRI scans were discarded due to bad segmentation caused by movement distortions, as explained in participants section.

For this study, the volumetric measures of interest were total cranial volume, cortical gray matter volume (CGV) and thalamic volume (TV). Total cranial volume was used for correcting subject differences in cranial size. CGV was used for accounting the cortical power produced by cortical neural mass and its relationship with the thalamic volume. TV was used as a measure of the neural mass that modulates cortical power with respect to the cortical neural mass.

The intensity measures of interest were the mean normalized intensity of the thalamus (TI), cortical gray matter (CGI) and cortical white matter (CWI). TI indicates the level of myelination in the thalamus. CGI indicates indicating myelin content and volume of dendrite/axons. CWI, sampled normal to the cortical layer at 2 mm under the cortical surface as per FreeSurfer, indicates the extent of myelination between the cortex and subcortex or the first layers cortical-cortical connections. Several studies report gray matter intensities contrasted with white matter intensities (Jefferson et al., 2015, Kong et al., 2015, Salat et al., 2009), while preserving the TI, to account for white matter intensity, it was included as a separate variable in the correlation analysis with power.

Statistical analyses

Partial Least Squares (PLS) and linear regression were used for statistical analysis. PLS is a multivariate technique used in neuroimaging (Krishnan et al., 2011, McIntosh and Lobaugh, 2004) based on singular value decomposition (available at: www.rotman-baycrest.on.ca/source/Pls.zip). It decomposes the data into latent variables (LV) composed of a left singular vector (left-SV), a singular value and a right-SV. Permutation is used to assess statistical significance by resampling without replacement the subjects’ group assignment and a p-value is obtained by counting the number of permutations where the permutation singular value exceeded the original singular value. PLS yields a single p-value for each LV, and accordingly does not require correction for multiple comparisons. Bootstrapping is used to measure the reliability of the features by resampling with replacement subjects within a group and taking the standard error (SE) of the bootstrapped left-SVs, then the original left-SV is divided by the bootstrap SE, and a bootstrap ratio is obtained. This ratio is presented as z-scores. The right-SV is termed the contrast, as it reflects the data-driven contribution of the groups or variables to the right-SV.

In this study, two types of PLS were used. First, a mean-centered PLS approach was used to test for differences between groups (ELGA, VLGA, and Term). Accordingly, the contrast reflects the contribution of a group, the z-scores reflect the reliability of a given feature (a frequency-specific relative power value in a source location) and the p-value indicates the confidence in rejecting the null hypothesis of no group differences represented by the contrast.

Second, a behavioral PLS approach was used assessing the correlation between brain features (frequency power at a location) and continuous variables (neurocognitive, neonatal, volumetric or intensity variables). In this case, the contrasts reflect the contribution of each continuous variable, the z-score the reliability of a feature correlating with the set of variables, and p-values indicate the confidence in rejecting the null where the set of variables are randomly correlated with the features. To increase statistical power, in all the behavioral PLS analyses all the subjects from all the groups with information on the variable of interest were included in the analysis.

To assess associations between neuroimaging measures and neurocognitive outcome at school age, 11 variables were included representing IQ, CBCL, BRIEF and BEERY subscales; in the intensity correlation, the included variables were thalamic intensity, cortical gray and cortical white matter mean intensities; for the volumetric PLS, the thalamic volume (TV) and cortical gray matter volume (CGV) ratio was obtained by dividing the former by the latter, to obtain the TV-CGV ratio. This was done to investigate associations between MEG power and TV not explained by CGV, as we hypothesized that decreased thalamic volume with increased cortical mass might better capture power changes across frequencies due to a reduction of thalamic cortical regulation relative to the amount of cortical neural mass. Cortical connections compensating deficits in thalamic connections could explain school-age very preterm delayed cortical thinning, a process where the brain optimizes communication by pruning unnecessary cortical synaptic connections (Mürner-Lavanchy et al., 2014).

Given that the TV-CGV ratio correlation with power does not tease out the sole contribution of TV or CGV, the TV and CGV were correlated with the so called ‘brain scores’ from the mean-centered PLS analysis. The brain scores are the dot product of the subjects by features matrix with the left-SV, and it expresses the covariance between the subjects’ features and SV. Essentially, the PLS ‘brain score’ for each subject indicates the degree to which that subject contributed to the pattern observed in the group contrast. A significant correlation between a volume variable and the brain scores would provide evidence of their contribution to the mean-centered PLS group differences.

Neonatal variables were not significantly correlated with power in the preterm group, and further analysis was conducted using linear regression to test if GA, sex, pain, morphine, number of days in ventilation, and illness (SNAP) were associated with TV-CGV ratio or with thalamic intensity, which in turn were associated with power and explained the association with negative neurocognitive outcomes and group differences, respectively.

Results

Atypical resting neurophysiological activity in ELGA children

A statistically significant group difference in resting spectral power was found with a p-value < 0.05 reflecting ELGA vs. VLGA and Term. As illustrated in Fig. 1, increased frontoparietal delta, frontal theta, and occipitotemporal gamma were observed in ELGA children relative to VLGA and Term (positive z-scores), whereas decreased in posterior areas in the alpha band and in frontal areas in the beta band were observed for the ELGA group (negative z-scores).

Fig. 1

Group differences in power across frequency bands. Plotted on the cortical surfaces are z-scores indicating the reliability of expressing the group contrast plotted on the top right. The contrast indicates ELGA > VLGA and term, thus red z-scores reflects higher power in the ELGA group. Statistical analysis was performed on an 8 mm spaced grid, for visualization purposes, the z-scores were interpolated onto cortical surfaces. N = 99. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Atypical resting activity associated with neurocognitive outcome

We tested for associations between power and a set of 11 neurocognitive measures and rendered a p-value < 0.05 and a data-driven contrast reflecting negative outcomes (CBCL, BRIEF) vs. positive outcomes (IQ, BERRY). Higher delta and occipitotemporal theta (positive z-scores), and lower alpha, beta and gamma (negative z-scores) were associated with higher negative outcomes and lower positive outcomes (Fig. 2). Overall, poorer outcome was associated with increased power at slower frequencies and decreased power at faster frequencies.

Fig. 2

Correlation between power and neurocognitive measures. Plotted on the cortical surfaces are z-scores indicating the reliability of expressing the correlation between the set of neurocognitive variables (explained in Methods section) and power. The bar graph contrast plotted on the top right indicates the correlation values for each of the variables and the whiskers indicate the 95% confidence interval. Confidence intervals not crossing the 0 line are considered significant. The contrast indicates negative > positive outcomes, thus higher power in red z-scores indicates higher negative outcomes and lower positive outcomes, i.e. lower or higher IQ, difficulties in executive functions and visual motor integration, conversely, lower powers represented in blue indicate lower negative and higher positive outcomes. Statistical analysis was performed on an 8 mm spaced grid, for visualization purposes, the z-scores were interpolated onto cortical surfaces. N = 98. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Atypical resting power associated with lower thalamic and cortical intensity

A significant association was found with a p-value < 0.05 between neurophysiological spectral power and mean intensity variables. Fig. 3. illustrate the results indicating that increased delta (frontoparietal), theta (frontal) and gamma (occipital) power, and wide-spread decrease in alpha and beta power was associated with lower thalamic intensity (TI) and increased cortical gray matter intensity (CGI) and cortical white matter intensity (CWI). Interestingly, the z-scores distribution resembled highly the PLS group power differences and we correlated both z-scores and obtained a correlation value of r = 0.7.

Fig. 3

Correlation between power and normalized mean intensity measures. Plotted on the cortical surfaces are z-scores indicating the reliability of expressing the correlation between the set of intensity variables (mean intensity in the thalamus (TI), cortical gray matter (CGI) and cortical white matter (CWI), explained in Methods section) and power. The bar chart contrast, plotted on the top right, indicates the correlation values for each of the variables, and the whiskers indicate the 95% confidence interval. The data-driven contrast indicates lower TI vs. higher CGI and CWI, thus higher power depicted by red z-scores imply smaller TI and higher CGI and CWI, and vice-versa with blue z-scores. Statistical analysis was performed on an 8 mm spaced grid, for visualization purposes, the z-scores were interpolated onto cortical surfaces. On the bottom right, the TI values are plotted with a box plot for each group, horizontal lines indicate significance between groups, * = p-val < 0.05. N = 50. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Associations between MEG power and volumetric measures

For PLS purposes, it would be equivalent to correlate power with TV and CGV as two separate variables. However, taking the TV-CGV ratio as a single variable allowed further investigation of the relations with cortical power by taking the highest PLS z-score percentile and visualizing their scatterplot and best fit and with neonatal variables by doing linear regression.

PLS identified a significant association between power and TV-CGV ratio with a p-value < 0.01. Decreased TV relative to the cortical volume of gray matter was associated with higher delta in general and theta in occipitotemporal areas, and lower alpha, beta and gamma (Fig. 4). The z-score distribution was very similar to the PLS association with neurocognitive measures, when correlating both z-scores we obtained a correlation value r = 0.9. A posteriori, to tease out if the ratio reflects the contribution of both variables or if it is driven by the TV while accounting for CGV, Pearson correlations with TV and CGV and the brain scores from the mean-centered PLS (group differences in power) were calculated. The correlation value for TV and brain scores was r = 0.4 and a p-val < 0.01, while CGV correlation was not significant, suggesting the TV as the main contributor.

Fig. 4

Correlation between power and the thalamic and cortical gray volumes ratio. Plotted on the cortical surfaces are z-scores indicating the reliability of expressing the correlation between the thalamic and cortical gray volumes (TV-CGV) ratio and power. The bar chart contrast, plotted on the top right, indicates the correlation value for the TV-CGV measure, and the whiskers indicate the 95% confidence interval. The contrast indicates > TV-CGV, thus higher TV-CGV is associated with higher power in red z-scores areas and blue z-scores with lower TV-CGV ratio . Statistical analysis was performed on an 8 mm spaced grid, for visualization purposes, the z-scores were interpolated onto cortical surfaces. On the bottom right, the TV-CGV values are plotted with a box plot for each group, horizontal lines indicate significance between groups, ** = p < 0.01. N = 50. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

We explored the association between the average power from the 2.5% locations with the strongest z-scores in each frequency band. For each frequency, a scatterplot and the line of best fit for each group was plotted and correlated the ratio with the average power for all the participants and within group (Fig. 5). The scatterplots and correlations indicate that the association of TV-CGV ratio and power is present across all frequencies and groups with the same effect direction, although, by combining all the groups higher statistical power is obtained.

Fig. 5

Correlation between average power from locations with the highest magnitude z-scores and the TV-CGV. Plotted on the cortical surfaces are the areas with the biggest absolute z-scores thresholded at 2.5%, with red to yellow indicating positive z-scores and blue to light blue indicating negative z-scores.. For delta and theta frequency bands, the z-scores were negatively skewed, whereas in the other bands were positively skewed with predominantly more negative z-scores. The 2.5% tail from the skewed side (plotted on the brain surfaces) at each frequency were selected and their power averaged. This tail represents the most stable features with a marked directionality (lower or higher between groups). The averaged power was correlated with the TV-CGV ratio, for each group and in total, and plotted in a scatterplot. The lines in the scatterplot represent the best linear fit, and in the box legend the correlation and significance are presented, with *= p < .05, **= p < .01, ***= p < .001. In all the groups and frequencies, the direction of the association (positive or negative) is common across groups, power in ELGA tends to correlate more, but overall, the association with all the participants is highest. N = 50. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Adverse neonatal experience associated with thalamic structure

Two linear regressions were calculated to predict either the TV-CGV ratio or the TI based on neonatal variables. The model for the TV-CGV was not significant, whereas the TI significance was p < 0.05, F(7,29) = 2.65, with R2 0.37. Beta and p-values for the predictor variables are reported in Table 2.

Table 2

Linear regression for thalamic mean normalized intensity. Predictor variables are gestational Age (GA), sex, GA*sex interaction, infection, illness severity (SNAP), cumulative morphine dose, number of skin breaking procedures (pain). P-values in red indicate statistical significance.

Model: TV-CGV
	Beta	Sig.
GA	−0.104	0.616
Sex	−14.14	0.030
Infection	−1.084	0.063
SNAP	−0.009	0.729
Morph	−2.948	0.008
Pain	−0.064	0.954
Sex*GA	0.487	0.027

Discussion

This study examined relationships between frequency-specific power, structural alterations associated with the thalamus, adverse neonatal variables, and neurocognitive difficulties at school age that characterize children born very preterm. We found that relative spectral power at rest in the extremely low gestational age (ELGA) group differed the most compared to the very low gestational age (VLGA) group and term control group, consistent with prior reports investigating absolute power (Kozhemiako et al., 2018), and likely related to prior observation of reduced alpha connectivity during task performance (Doesburg et al., 2010). Neurocognitive difficulties were associated with increased frontal and occipital delta, and occipital theta, as well as decreased power mostly involving frontal and occipital beta activity. Thalamic and cortical gray matter as well as cortical white matter intensities were associated with power and were very similar, in direction and spatially, to the group differences in power, differentiating ELGA from the VLGA and term groups, with a correlation between z-scores of r = 0.7. Similarly, the association between MEG power and the ratio thalamic/cortical volume revealed a spatial pattern very similar to that observed in the association between MEG power and neurocognitive difficulties. These similarities suggest that thalamic intensity, accounting for cortical gray and white matter intensity, is a prominent structural correlate of atypical power in children born at ELGA, and that thalamic volume, normalized by the total neural cortical mass, is a structural correlate of neurocognitive difficulties. Overall, this study provides novel evidence regarding the relationship between atypical spontaneous neurophysiological activity in preterms and structural thalamic alterations and negative outcomes at school age.

The central hypothesis of this study was based on reports indicating thalamic volume reduction and thalamocortical alterations in children born very preterm (McQuillen and Ferriero, 2006, Smyser et al., 2010), on studies indicating the role of the thalamus as a regulator of cortical power through low frequencies (FitzGerald et al., 2013, Malekmohammadi et al., 2015, Ribary, 2005), and on studies demonstrating the importance of the thalamus for higher-order cognitive processes (Mitchell, 2015, Saalmann, 2014, Sweeney-Reed et al., 2015). Based on this previous literature, we hypothesized that decreased thalamic volume might lead to frequency-specific alterations in spontaneous neurophysiological activity which would be associated with negative neurocognitive outcomes. Indeed, we found atypical power in the ELGA group marked by a relative increase at lower frequencies (delta and theta) together with a decrease at higher frequencies (alpha, beta, and partially gamma). This pattern of altered resting neurophysiological activity was also associated, in general, with more negative outcomes. These results are consistent with the phenomena of ‘alpha slowing’, which is characterized by decreased alpha power and increased theta and/or delta power (Garcés et al., 2013, Hughes and Crunelli, 2005, Llinás et al., 1999). This slowing has been associated with a decrease of inhibitory interneurons in the thalamus or a shift of inhibitory drive on the thalamus (Bhattacharya et al., 2011, Ribary et al., 2019). The results of the present study suggest that in ELGA children this may be structurally reflected with decreased gray matter intensity in the thalamus. This power imbalance seems to be already present in the neonatal period in preterm infants with a relatively higher contribution of low frequencies (Suppiej et al., 2017). Similarly, Doesburg et al. 2013, using a subset of the participants analyzed in the present study, found slowing of alpha peak, especially in the ELGA, in sensor-space during a visual-spatial memory task at age 7–8 years. In the present study, we show a source-space spatially-resolved characterization of atypical power, most distinctly evident in children born extremely preterm, and the association with structural thalamic deficits and neurocognitive difficulties.

In the analysis of associations between power and thalamic volume (TV) and intensity (TI), a priori, the cortical counterpart was included, i.e. cortical gray matter volume (CGV) or intensity (CGI), to account for power independent from thalamic alterations. Although decreased cortical thickness in preterms at school age has been reported (Lax et al., 2013, Ranger et al., 2013), very preterm birth is associated with neurodevelopmental problems reflected in slower neuronal optimization archived by pruning unnecessary synapses and strengthening necessary ones with myelination (Shaw et al., 2008). Longitudinal preterm studies have found that between 7 and 12 years of age there is decreased thinning of cortical gray matter volume indicating delayed cortical maturation (Ment et al., 2009, Mürner-Lavanchy et al., 2014). In our study, while CGV was significantly reduced in the ELGA, we found that the ratio TV-CGV was also significantly decreased, indicating a relatively larger CGV. In relation to the thalamus, we hypothesized that reduced thalamic regulation, either by a smaller number of thalamic neurons or pathways, might lead to cortical dysfacilitation were more cortical synapses are needed to compensate for the reduced role of the thalamus in integrating cortical processing. This hypothesis would also account for the increase of delta and theta seen in the ELGA, as more low-frequency integrative processing in the cortex would be needed, and, assuming cortical integration to be less effective than the thalamic, it would also explain the correlation between neurocognitive difficulties and higher low frequency power.

We expected to find a direct association between power and neonatal procedures, especially, pain (skin breaking procedures). However, at most, we found an indirect link with power by using neonatal variables as predictors of TV-CGV ratio or thalamic intensity (TI). While TV-CGV was not significantly associated with neonatal variables, TI reduction was predicted by morphine instead of pain, with a possible more pronounced effect on females. Gray matter intensity mostly reflects the amount of myelination or the amount of axonal mass in the cortical layer, while subcortical white matter intensity reflects cortical or subcortical axons. Thalamic gray matter intensity reflects the extent of myelinated axons inside the thalamus, and thus it is a direct structural measure of the thalamus. Alterations in myelination have been previously reported in preterm populations (Lubsen et al., 2011, Nossin-Manor et al., 2013), including in the thalamus (Eaton-Rosen et al., 2015), with early exposure to opioid affecting myelination (Sanchez et al., 2008). Nonetheless, the prediction of morphine by TI might also reflect its apoptotic effects in the thalamic neurons and microglia (Hu et al., 2002, McPherson and Grunau, 2014). In this study, we provide evidence that a higher cumulative dose of morphine is a predictor of TI which, in turn, is associated with atypical spontaneous neurophysiological activity in ELGA children. This is consistent with previous findings on the adverse effects of high dosing of morphine in relation to hippocampal (Duerden et al., 2016) and cerebellar (Zwicker et al., 2016) growth on MRI in the neonatal period in a different cohort, and cerebellar growth in the current cohort at school-age (Ranger et al., 2015).

One limitation of the present study is the limited sample size available for each group, imposed by the limited number of MRI scans available. It is necessary, however, to divide the very preterm children into ELGA and VLGA groups as prior research has demonstrated that alterations in spontaneous MEG activity are distinct for these populations (Kozhemiako et al., 2018). It was not possible to further split the groups by gender to study specific sex differences of preterm birth structural alterations and brain activity as we were under-powered for this analysis due to limited availability of MRI scans. We previously demonstrated cortical power and connectivity alterations associated with sex-specific preterm birth using an overlapping sample (Kozhemiako et al., 2019), and in the present study it is not possible to rule out the effects of sex differences between groups. Despite this limitation in sample size, the observations of the present study are novel and contribute to our understanding of structure-function relations as they pertain to school-age outcome in children born very preterm. Importantly, the data used in the present study is the largest dataset available, to our knowledge, that includes neonatal data together with school-age MEG and MRI, making it uniquely well positioned to address these questions.

In conclusion, our findings present evidence linking atypical spontaneous neurophysiological activity at rest in children born very preterm, especially ELGA, at school age with neurocognitive difficulties, structural deficits related with the thalamus, and morphine exposure. During the neonatal period following very preterm birth, the thalamocortical system is going through critical developmental windows lending heightened susceptibility to disruptions which may lead to alterations lasting until school age and possibly beyond.

Credit authorship contribution statement

Adonay S. Nunes: Conceptualization, Methodology, Formal analysis, Writing - original draft. Nataliia Kozhemiako: Methodology, Writing - review & editing. Evan Hutcheon: Data curation. Cecil Chau: Data curation. Urs Ribary: Funding acquisition, Resources, Writing - review & editing. Ruth E. Grunau: Supervision, Funding acquisition, Project administration, Writing - review & editing. Sam M. Doesburg: Supervision, Funding acquisition, Resources, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.