Part 1: Relating noxious-evoked response amplitude to resting-state activity
Subject information
We recruited healthy neonates from the John Radcliffe Hospital postnatal ward (Oxford University Hospital NHS Trust). Neonates were considered healthy if they were inpatients that never required admission to the neonatal unit, had no history of congenital conditions or neurological problems, and were clinically stable at time of study. Written informed consent was obtained from parents prior to the study. Ethical approval was obtained from an NHS Research Ethics Committee (National Research Ethics Service, REC reference: 12/SC/0447), and research was conducted in accordance with Good Clinical Practice guidelines and the Declaration of Helsinki.
Twenty one neonates had noxious-evoked and resting-state data collected. Subjects were excluded from analysis if (i) either scan session was not fully completed, in order to remove inter-subject variability in data quality related to scan length, and (ii) the vertex of the cerebral cortex left the scan field of view for more than 5% of the scan session, in order to ensure reliable data in this functionally relevant brain region. Using these criteria, three subjects were excluded resulting in a final sample of n=18 neonates. Demographic details for this sample are displayed in Table 3.
Experimental design
Scans occurred in the Wellcome Centre for Integrative Neuroimaging (Oxford, UK). Neonates were fed and swaddled, fitted with ear protection, then placed on a vacuum-positioning mattress with padding around the head. Heart rate and blood oxygen saturation were monitored throughout scanning, but were not of sufficient quality for analysis. The noxious stimulation paradigm was an event-related design in which a mild non-skin-breaking noxious stimulus, the 128 mN sharp-touch pinprick (PinPrick Stimulator, MRC Systems), was applied to the dorsum of the left foot: ten trials, 1 s per trial, 25 s minimum inter-stimulus interval10. Stimuli were applied only when the neonates were naturally still to minimise the influence of motion artefacts. For all other scan types, neonates lay passively in the scanner. No sedatives were used at any stage.
MRI data acquisition
All data were collected on a 3T Siemens Prisma with an adult 32-channel receive coil using the following scan parameters. Structural: T2-weighted, TSE (factor 11), 150° flip angle, TE=89 ms, TR=14,740 ms, GRAPPA 3, 192 × 192 in-plane matrix size, 126 slices, 1 mm isotropic voxels, acquisition time (TA)=2 mins 13 s. Fieldmap: gradient echo, 2DFT readout, dual echo TE1/TE2=4.92/7.38 ms, TR=550 ms, 46° flip angle, 90 × 90 in-plane matrix size, 56 slices, 2 mm isotropic voxels, TA=1 min 40 s. Resting-state and noxious stimulation fMRI: T2* BOLD-weighted, gradient echo, EPI readout, 70° flip angle, TE=50 ms64, TR=1,300 ms, multiband 465,66, 90 × 90 in-plane matrix size, 56 slices, 2 mm isotropic voxels, AP phase encode direction. Resting-state TA=10 mins 50 s (500 volumes), mean noxious stimulation TA=6 min approx. (277 volumes approx.). Diffusion MRI: T2 diffusion-weighted, spin echo, EPI readout, 90° flip angle, TE=73 ms, TR=2,900 ms, multiband 3, 102 × 102 in-plane matrix size, 60 slices, 1.75 mm isotropic voxels, AP phase encode direction, multishell (b=500, 1000, 2000 s/mm2), 143 directions uniformly distributed over the whole sphere, TA=8 mins. Phase-reversed b0 images were collected to derived a spin-echo fieldmap for distortion correction of the diffusion data. The entire MRI acquisition protocol had a nominal duration of 40 mins.
MRI data preprocessing
All MRI data were preprocessed using Developing Human Connectome Project (dHCP) pipelines. The T2 structural data were processed using the MIRTK Draw-EM neonatal pipeline67,68. Noxious stimulation and resting-state fMRI data were preprocessed using the dHCP fMRI pipeline47,69. Data were motion and distortion corrected using FSL’s EDDY70,71, which included slice-to-volume motion correction72 and susceptibility-by-movement distortion correction73. Noxious stimulation data were high-pass temporally filtered at 0.01 Hz, and resting-state data at 0.005 Hz. Data were denoised using FSL’s FIX74,75, low-pass spatially filtered with a 3 mm FWHM filter using FSL’s SUSAN76, and scaled to a common global spatiotemporal median. For spatial normalisation, data were registered from functional to structural space using BBR77 with FSL’s FLIRT78,79, then from structural space to the 40-week standard template80 using ANTs’s SyN81.
Diffusion data were analysed using the dHCP dMRI pipeline35,82. Phase-reversed fieldmaps were processed using FSL’s TOPUP83,84. Data were corrected for motion, distortion, and eddy currents using FSL’s EDDY, which included outlier detection and replacement85, slice-to-volume motion correction, and susceptibility-by-movement distortion correction. Spatial normalisation followed the same sequence of registrations as the functional data.
Noxious-evoked response amplitudes
Noxious-evoked response maps were generated using subject-level voxelwise GLM analysis in FSL’s FEAT86, fitting the term-neonate double-gamma HRF44,69. A group average t-statistic map was generated using the 18 subjects’ regression parameter maps. Regression parameter maps were used as subject-level response maps; the group-average t-statistic map was used as the group-level response map. The group-level response map was regressed onto each subject-level response map producing a spatial regression coefficient, constituting each subject’s overall noxious-evoked response amplitude (Supplementary Figure 10 Step 1). The influence of HRF goodness-of-fit on noxious-evoked response amplitudes was assessed in Supplementary Information (Supplementary Figures 1 and 3).
To localise the noxious-evoked activity, we generated a thresholded group average activity map using group-level voxelwise GLM analysis in FSL’s Randomise87. Statistical significance was assessed using permutation testing with 10,000 permutations, variance smoothing (6 mm FWHM kernel) due to limited degrees of freedom88, cluster-based thresholding with z=3.1 (p=0.001) cluster-defining threshold89, and a FWER-corrected cluster p-value of p=0.05. Regions of activity were identified using the probabilistic Harvard-Oxford Cortical Structural Atlas90.
To gain insight into the processes underlying the noxious-evoked activity, we assessed expression of several adult pain signature templates and Neurosynth meta-analysis association test map templates by performing whole-brain correlations between adult template maps and subject-level response maps (Supplementary Figure 10 Step 2). For each template, n=18 correlation coefficients were generated, and group average template expression (average correlation) was assessed using two-tailed t-tests with statistical significance assessed non-parametrically in FSL’s PALM87 using 10,000 permutations.
The adult pain signatures tested were the Neurologic Pain Signature (NPS) and Social Rejection Pain signature, where the social rejection template was used as a negative control, as per the original NPS adult study37. Association test maps were used for all Neurosynth terms, as these maps display brain regions that are preferentially related to the term-of-interest. The primary Neurosynth term-of-interest was “pain”; “visual” was used as a Neurosynth negative control; a subset of pertinent pain dimension was assessed, including the sensory-discriminative term “nociceptive”, the cognitive terms “arousal”, “salience”, and “attention”, and the emotional term “unpleasant”. Due to the whole-brain nature of the expression correlations and the default thresholded nature of all adult template maps included, we used the thresholded group average activity map generated from the neonates’ noxious-evoked responses (thresholding performed in Randomise as described above) as the positive control. While the correlation between this thresholded noxious-evoked response map and the neonates’ response maps is likely inflated due to the circular analysis, the correlation strength of this positive control sets a useful upper limit reference.
Finally, inter-subject variability in noxious-evoked response amplitudes (regression parameters) and adult template correspondences (correlation coefficients) were assessed for all functional templates using two-tailed Pearson correlation tests with statistical significance assessed non-parametrically in PALM using 10,000 permutations.
Noxious stimulation imaging confounds
In all analyses using noxious-evoked response, amplitudes were adjusted for (i) mean head motion: mean framewise displacement, (ii) stimulus-correlated head motion: multiple correlation coefficient between the predicted BOLD response (stimulus timeseries convolved with the HRF) and the 24 head motion timeseries (estimated during preprocessing), and (iii) CSF signal amplitude: the mean regression coefficient within the CSF ROI, intended to capture residual cardiac pulsatility. Details on ROI construction are provided in Supplementary Information (Supplementary Figure 8).
Resting-state network amplitudes
To define a robust set of resting-state networks (RSNs), RSNs identified in the noxious stimulation paradigm dataset (n=18) were compared to those identified in a dHCP dataset previously produced as part of the dHCP38 (Supplementary Figure 10 Step 3). Robust networks were defined as those replicated across datasets to ensure networks were not dataset-specific. This dHCP dataset included 242 healthy term-aged neonates: mean GA at birth = 38.6 weeks; mean PMA at scan = 40.4 weeks; 112 females38. The RSN analysis performed on our noxious stimulation paradigm dataset was matched to that of the dHCP dataset38,47. In brief, probabilistic functional mode (PFM) analysis using FSL’s PROFUMO91,92 was run with a pre-specified dimensionality of 25, using the term-neonate double-gamma HRF44,69 as the temporal prior. PROFUMO’s Bayesian model complexity penalties eliminated modes unsupported by the data, thus returning a number of group-level modes less than the pre-specified dimensionality. The data-determined dimensionality for the noxious stimulation paradigm dataset was 11, nine of which corresponded to the dHCP dataset resting-state networks assessed using spatial correlation followed by visual confirmation. Due to the larger sample size, the dHCP RSN maps had greater SNR and were thus used throughout our analyses as the nine RSN template maps.
These RSN template maps were spatially regressed onto each neonate’s resting-state data using multiple regression, resulting in network timeseries. Timeseries amplitudes were quantified as the median absolute deviation (MAD), due to MAD’s increased robustness to outliers compared to the more commonly used standard deviation (Supplementary Figure 10 Step 4). The association between network timeseries outliers and noxious-evoked response amplitudes was assessed in Supplementary Information (Supplementary Figures 2 and 3).
Resting-state imaging confounds
In all analyses using RSN amplitudes, amplitudes were adjusted for (i) mean head motion: mean framewise displacement, (ii) CSF amplitude: timeseries MAD extracted using the CSF ROI, intended to capture residual cardiac pulsatility, and (iii) white matter amplitude: timeseries MAD extracted using the white matter ROI, intended to capture global signal amplitude. Details on ROI construction are provided in Supplementary Information (Supplementary Figure 8). These confounds were also directly tested for association with noxious-evoked response amplitudes.
Clinical variables
Six clinical variables were tested for association with noxious-evoked response amplitudes: postmenstrual age (PMA), gestational age (GA), postnatal age (PNA), birth weight (BW), total brain volume (TBV), and sex (Table 3). Age variables are defined according to the American Academy of Paediatrics93, and TBV was calculated from neonates’ structural MRI tissue segmentation outputs. Testing the clinical variables assessed whether associations between RSN amplitudes and noxious-evoked response amplitudes could be explained by these biologically interesting variables.
Predicting noxious-evoked response amplitudes
For prediction analyses (Supplementary Figure 10 Step 5), the primary responses to be predicted were the neonates’ overall noxious-evoked response amplitudes (Figure 1 scalar values). Three sets of predictors were tested for predictive capacity: nine RSN amplitudes, six clinical variables, and three resting-state imaging confounds. To further assess whether pain components of the overall noxious-evoked response amplitudes can themselves be predicted from resting-state data, we assessed RSN amplitude-based predictions for both the NPS and Neurosynth pain map response magnitudes. Adult pain signature response magnitudes were quantified using cosine similarity, which is equivalent to the Pearson correlation coefficient without mean-centring, thus retaining magnitude information.
For each set of predictors, a multivariate linear support vector regression (SVR) model was used to generate predictions using leave-one-out cross-validation (LOO-CV). The linear SVR model was fit using scikit-learn packages94. Noxious-evoked, NPS, and Neurosynth pain responses were all adjusted for the noxious-evoked response imaging confounds, and RSN amplitudes were adjusted for the resting-state imaging confounds using cross-validated confound regression (https://github.com/lukassnoek/MVCA)95. SVR parameters were: kernel=linear, loss=epsilon insensitive, epsilon=0.1, regularization=ridge, regularization strength={0.001, 0.01, 0.1, 1}. Regularization strength optimisation was performed using an initial LOO-CV grid search.
Prediction performance was assessed using root mean squared error (RMSE), sums-of-squares coefficient of determination (R2), and Spearman’s rank correlation coefficient (RSp). We performed one-tailed significance tests for these measures using permutation tests, running 1,000 permutations through the full prediction pipeline. The relationship between individual predictors and noxious-evoked response amplitudes was assessed in Supplementary Information (Supplementary Figure 4).
Part 2: Relating noxious-evoked response amplitude to white matter microstructure
Structure-function analysis using an exploratory-confirmatory approach
The neonates’ noxious-evoked response amplitudes were assessed for structure-function associations by analysing white matter microstructure. Due to lack of research into the microstructural basis of noxious-evoked responses in healthy neonates, an exploratory analysis was required. However, due to the small sample size of the noxious stimulation paradigm dataset (n=17, Table 3), appropriate corrections for multiple testing would prohibit identification of true positives.
We thus adopted a two-armed two-dataset exploratory-confirmatory analysis approach. In the initial exploratory arm, we used a large age-matched dHCP sample to test a range of microstructural features to identify candidate white matter tracts and diffusion model parameters (Supplementary Figure 10 Step 6). Structure-function relationships identified here were used to formulate data-driven hypotheses. In the confirmatory arm, these hypotheses were subsequently tested for validation in the noxious stimulation paradigm dataset (Supplementary Figure 10 Step 7). Cross-dataset consistencies in structure-function associations constitute initial indirect external validation for the resting-state prediction model, as these consistencies rely on the predicted response amplitudes in the dHCP dataset being similar in nature to the true response amplitudes in the noxious stimulation paradigm dataset.
dHCP dataset sample selection
Neonates were included in our dHCP sample if data were of reasonable quality and neonates were age-matched to our noxious stimulation paradigm dataset. The three data quality criteria were: (i) both fMRI and dMRI data passed dHCP QC pipelines47,82, (ii) both scan sessions completed fully, and (iii) the vertex of the cerebral cortex remained within the scan field of view for at least 95% of the scan session. The two age criteria were: (i) neonates were between 36–42 weeks for both GA and PMA, and (ii) neonates were scanned within the first 10 postnatal days. These criteria resulted in a sample of n=215 neonates. Key functional and diffusion acquisition parameters for this dHCP dataset are displayed in Table 4 to facilitate comparison with our noxious stimulation paradigm dataset.
dHCP dataset noxious-evoked response amplitudes
The dHCP dataset does not included noxious stimulation paradigm data. To analyse structure-function relationships relevant to noxious-evoked responses, the dHCP resting-state data were mapped to noxious-evoked response amplitudes using the RSN amplitude-based SVR prediction model described above. This model was trained on the noxious stimulation paradigm dataset (n=18) using the nine confound-adjusted RSN amplitudes as predictors and the confound-adjusted noxious-evoked response amplitudes as responses. In our dHCP dataset (n=215), predictors and confounds were extracted in an identical manner to the noxious stimulation paradigm dataset analysis, and the RSN amplitudes were used to generate predicted noxious-evoked response amplitudes.
White matter microstructural features
Our dHCP sample (n=215) was used to generate 16 bilateral white matter tracts using the “baby autoPtx” approach established as part of the dHCP dMRI preprocessing pipeline development35. In brief, FSL’s probabilistic multi-shell ball and zeppelins model96 is fit, then probabilistic tractography using FSL’s PROBTRACKX97,98 is run using pre-defined seed, target, and exclusion masks. At the time of analysis, 29 tracts were available, 13 unilateral and three bilateral. Unilateral tracts were fused to create bilateral tracts analogous to our bilateral resting-state networks, resulting in a total of 16 bilateral tracts: acoustic radiation, anterior thalamic radiation, cingulate gyrus part of the cingulum, parahippocampal part of the cingulum, corticospinal tract, forceps major, forceps minor, fornix, inferior fronto-occipital fasciculus, inferior longitudinal fasciculus, middle cerebellar peduncle, medial lemniscus, posterior thalamic radiation, superior longitudinal fasciculus, superior thalamic radiation, uncinate fasciculus. Normalised probability value results of each tract were group-averaged in standard space and thresholded at a probability of 0.01.
We used FSL’s DTIFIT to generate mean diffusivity (MD), fractional anisotropy (FA), and mean kurtosis (MK) parameter maps. We used the MD and FA maps from the b1000 shell, and MK maps from all three shells in both datasets. We thresholded each parameter map to remove noisy voxels with values outside the expected theoretical range, likely due to poor SNR or head motion: negative values for MD; values outside [0,1] for FA; values outside [0,3] for MK. Mean parameter values for each tract constituted the white matter microstructural features for our structure-function analyses. Comparisons of these microstructural features between the noxious stimulation paradigm and dHCP datasets is assessed in Supplementary Information (Supplementary Figure 6).
Identifying nociception-related structure-function associations
Exploratory arm analyses:
Using the dHCP sample (n=215), univariate correlations between predicted noxious-evoked response amplitudes and each microstructural feature was assessed using two-tailed Pearson correlations adjusted for three dMRI imaging confounds: mean head motion, number of noisy voxels outside expected theoretical ranges, and TBV. Adjustment for TBV mitigates macrostructural tissue density and partial volume effect confounds, which can bias dMRI microstructural parameters. Statistical significance was assessed using 10,000 permutations and FWER-corrected for multiple testing across all 48 tests (16 tracts x 3 parameters)99. While significant correlations are statistically valid, they are tentative due to the use of predicted noxious-evoked response amplitudes without a known ground truth.
Due to the presence of consistent structure-function correlation polarities across tracts identified in the dHCP dataset, variance across multiple tracts for a single dMRI parameter was pooled using principal component analysis (PCA). The first principal component across these tracts (e.g. MD PC1) was assessed for associations with noxious-evoked response amplitudes using two-tailed Pearson correlations, significance assessed using 10,000 permutations. In the dHCP dataset, these effect size and statistical significance measures are biased due to circularity in explanatory variable selection39, but the bias is restricted to the exploratory arm analyses.
Confirmatory arm analyses:
The two hypotheses formulated in the dHCP dataset (see Results) were tested in the noxious stimulation paradigm dataset. First, the polarity of the Pearson correlation coefficient between noxious-evoked response amplitudes and individual tracts were qualitatively assessed. While similar correlation polarities alone are not strong evidence, when considered in combination with other elements of evidence, this correlation polarity information provides useful complementary insight. Second, the correlation between noxious-evoked response amplitudes and microstructural features of specific subsets of tracts was assessed using one-tailed Pearson correlation tests, the directionality of the tailed test hypothesised by the nature of the exploratory arm results. Statistical significance was assessed in PALM using 10,000 permutations. The effects of both global and local structure-function association effects observed in both the dHCP and noxious stimulation paradigm datasets were further assessed in Supplementary Information (Supplementary Figures 7-8).
Data Availability
Source data for figures 1-6 are provided with the paper. The data that support the findings of this study are available from the corresponding author upon reasonable request.