Functional and diffusion MRI reveal the neurophysiological basis of neonates’ noxious-stimulus evoked brain activity

Understanding the neurophysiology underlying neonatal responses to noxious stimulation is central to improving early life pain management. In this neonatal multimodal MRI study, we use resting-state and diffusion MRI to investigate inter-individual variability in noxious-stimulus evoked brain activity. We observe that cerebral haemodynamic responses to experimental noxious stimulation can be predicted from separately acquired resting-state brain activity (n = 18). Applying this prediction model to independent Developing Human Connectome Project data (n = 215), we identify negative associations between predicted noxious-stimulus evoked responses and white matter mean diffusivity. These associations are subsequently confirmed in the original noxious stimulation paradigm dataset, validating the prediction model. Here, we observe that noxious-stimulus evoked brain activity in healthy neonates is coupled to resting-state activity and white matter microstructure, that neural features can be used to predict responses to noxious stimulation, and that the dHCP dataset could be utilised for future exploratory research of early life pain system neurophysiology.

evoked responses (HRF goodness-of-fit z-scores, Supplementary Figure 1) and resting-state network timeseries (total number of outlier timepoints and percentage of unique outlier timepoints, Supplementary Figure 2). The associations are displayed in Supplementary Figure 3.
In all cases, we performed Pearson correlation tests with statistical significance assessed in FSL's PALM 2 using a two-tailed permutation test with 10,000 permutations. We did not correct for multiple testing, and no tests displayed significance or trends towards significance in the uncorrected p-values. The Pearson correlation coefficients and associated p-values are displayed in Supplementary Figure 3. These results suggest that the observed variability in noxious-evoked response amplitudes (i) do not reflect noxious-evoked response data quality related to the degree of (mis)alignment when fitting the HRF, or (ii) do not correspond to resting-state network timeseries data quality, and by extension, the quality of the estimated resting-state network timeseries amplitudes.

Resting-state coupling: global effects, non-uniformity, and nociception-relevance
In our core analyses, we tested three sets of predictors for their capacity to predict neonates' noxious-evoked response amplitudes: nine resting-state networks, three resting-state imaging confounds, and six clinical variables. The multivariate resting-state network amplitudes could predict neonates' noxious-evoked response amplitudes with statistically significant accuracy (main text Table 2). The resting-state imaging confound multivariate model, which included head motion as well as CSF and white matter signal amplitudes, was not predictive, demonstrating that the resting-state network predictions could not be attributed to head motion or artefactual global signals. Finally, the multivariate clinical variables model (which included postmenstrual, gestational, and postnatal age, as well as birth weight, sex, and total brain volume) was not predictive, suggesting these biologically interesting variables similarly could not account for the resting-state network predictions. Here, we examine each of these 18 predictors individually to gain further insight into our prediction results.

Univariate correlation analysis
Performing univariate correlations between each individual resting-state network's amplitudes and the noxious-evoked response amplitudes (adjusted for resting-state and noxious-evoked response imaging confounds, see Methods), eight of nine resting-state network amplitudes exhibited positive correlations with noxious-evoked response amplitudes (Supplementary Figure 4a) -the executive control network had a negligible negative correlation coefficient (r=-0.03). This large proportion of positive associations suggests a near global effect that may be attributable to brain maturity, such that larger evoked response amplitudes and larger resting-state network amplitudes are associated with a more functionally mature brain (see Discussion in main text). Statistical significance for each correlation was assessed using permutation testing in FSL's PALM using 10,000 permutations and FWER-correction for multiple testing across the 18 univariate tests. Of the 18 variables, only two resting-state network correlations were statistically significant: somatomotor network (SMN; r=0.73, p= 0.011) and the occipital pole visual network (VNop; r= 0.77, p=0.003). Interestingly the role of visual brain regions in pain and nociception is also reported in the adult literature -see the recent adult pain prediction paper by Spisak and colleagues 3 and the discussion and further reference therein. No individual resting-state imaging confound or clinical variable exhibited trends towards significance, even without corrections for multiple testing.

Univariate prediction analysis
Despite the trend towards a global positive association between all resting-state network amplitudes and noxious-evoked response amplitudes, the association strengths across networks were not uniform but varied largely from network to network (Supplementary Figure 4a). To identify individual networks most predictive of neonates' noxious-evoked response amplitudes and to further assess the network-to-network variability in univariate correlations strengths with noxiousevoked response amplitudes, we examined univariate predictions. For completeness, we performed this univariate prediction analysis for all 18 predictors. Each univariate prediction model was trained in an identical manner to our multivariate prediction models described in the main text, using a linear support vector regression model, leave-one-out cross validation, and grid search optimisation of regularisation strength. Each of the nine univariate resting-state network amplitude prediction models included cross-validated confound regression using the three resting-state imaging confounds to adjust predictors, and all models included cross-validated confound regression using three noxious-evoked response imaging confounds to adjust responses. Prediction performance was assessed using sums-of-squares formulation of the coefficient of determination (R 2 ), and each univariate model performance was qualitatively compared to each other and to the prediction performance of the three multivariate prediction models ( Table 2 in Figure 4b). Of the nine resting-state networks, only the SMN exhibit notable predictive value achieving a prediction performance (R 2 = 0.60) comparable to our multivariate resting-state prediction model performance (R 2 = 0.62) (Supplementary Figure 4b). Thus, while eight of nine resting-state networks were positively correlated with noxious-evoked response amplitudes, notable network-to-network variability in positive association strength existed, and the most robustly associated and predictive network (SMN) was a network with high functional relevance to our noxious stimulation paradigm.
These results suggest that global network correlation polarities may reflect global brain maturity effects 4,5 , while local network-to-network variability in correlation strengths may reflect local maturational asynchronies 6,7 . Thus in order to understand neonatal noxious-evoked response amplitudes, both global characteristics not specific to neural processing of noxious stimulation (e.g. maturational state) and local nociception-related characteristics (e.g. resting-state activity of functionally relevant networks) are needed. While we did not find maturity-related variables such as age and brain volume to have much predictive value in this study, it is likely these variables would increase their predictive value when considering a wider age range, so we do not conclude that age and brain volume are irrelevant to understanding neonates' noxious-evoked response amplitudes.

Cross-dataset similarities in functional and structural features
Our structure-function association analyses assessed the relationship between neonates' noxiousevoked response amplitudes and three white matter microstructural properties: mean diffusivity (MD), fractional anisotropy (FA), and mean kurtosis (MK). For these analyses, we adopted a twopart two-dataset approach involving exploration of a range of structure-function associations in the dHCP dataset (n=215) followed by confirmation of findings in our noxious stimulation paradigm dataset (n=17). In the noxious stimulation paradigm dataset, measured noxious-evoked response amplitudes are available. In the dHCP dataset, noxious-evoked response data are not available, so noxious-evoked response amplitudes are predicted per neonate from their resting-state data using our resting-state prediction model.
The two datasets were collected on different scanners using different diffusion and functional scan protocols -see Table 3 in main text for cross-dataset comparison of key acquisition parameters. These acquisition differences can introduce undesirable variability in the data that could cause issues in identification of replicable structure-function associations between datasets. Additionally, the reliance on predicted noxious-evoked response amplitudes in the dHCP dataset will contribute further variability due to unavoidable limitations in the accuracy of the resting-state prediction model. Here we qualitatively compared the distributions of the functional and structural features to identify any major inconsistencies between datasets. To calculate distributions for the three diffusion parameters, an average parameter value was calculated for each of the 16 white matter tracts per neonate. The results are displayed in Supplementary Figure 5.
Despite the differences in functional acquisition protocols and the reliance on predicted noxiousevoked response amplitudes in the dHCP dataset, the frequency distributions of noxious-evoked response amplitudes were very well aligned between the two datasets (Supplementary Figure 5, left). For each diffusion parameter, there was reasonably high consistency in frequency distributions between datasets, with greatest alignment visible for MD. Importantly, it was the structure-function associations that we attempted to replicate between datasets, so it is within-dataset relations between MRI modalities that is of primary interest. The minor cross-dataset variability in these functional and structural feature distributions (Supplementary Figure 5) is likely to have minimal influence on subject-to-subject variability within a dataset, and is thus unlikely to detrimentally impact comparisons of structure-function associations between datasets.
The close correspondence of both the functional and structural features (Supplementary Figure 5) and the structure-function associations between datasets ( Figure 6 and Supplementary Figures 6-7) suggests that undesirable between-dataset inconsistencies has, at most, minor effects, and our resting-state prediction model generalises well to the age-matched dHCP data.

Structure-function associations: global effects, non-uniformity, and nociception-relevance
In our structure-function analysis, we defined 16 bilateral white matter tracts and assessed crosssubject variation in these tracts for three dMRI parameters: mean diffusivity (MD), fractional anisotropy (FA), and mean kurtosis (MK). We explored and identified nociception-related structurefunction associations in the dHCP dataset by correlating the neonates' predicted noxious-evoked response amplitudes with these 48 dMRI features (16 tracts x 3 parameters). The results of this exploratory analysis are presented in Figure 5 of the main text. We note two observations. First, there are global effects for both MD and FA: in all 16 tracts, there are negative associations between noxious-evoked response amplitudes and MD and positive associations between noxious-evoked response amplitudes and FA. Second, there is a local effect where statistically significant correlations were observed for the MD of five tracts, and these five tracts have high functional relevance to our noxious stimulation paradigm -see Discussion in main text. Furthermore, of these five tracts, the superior thalamic radiations and the corticospinal tracts form core structural connectivity for somatosensory and motor functions. This specific subset of functionally relevant tracts mirrors our identification of the resting-state SMN activity as most robustly associated with neonates' noxious-evoked response amplitudes (Supplementary Figure 4b). Thus, analogous to our resting-state functional coupling results above, our structure-function association results suggest both global and local effects are important features of the data. Here, we examine these effects further to gain additional insight into our structure-function association results.
To test for global structure-function effects in our noxious stimulation paradigm dataset, the Pearson correlation coefficients for each tract and diffusion parameter were generated using matched analysis to that in the dHCP dataset. The correlation results for both datasets are displayed together in Supplementary Figure 6a. As in the dHCP dataset, the noxious stimulation paradigm dataset exhibited global effects for both MD and FA: all 16 tracts exhibited negative correlations with MD, and 15 of 16 tracts exhibited positive associations with FA (the medial lemniscus exhibited a negative association). We attribute these global effects to brain maturity such that larger noxiousevoked response amplitudes and white matter tracts with smaller MD values and larger FA values are associated with a more functionally and structurally mature brain (see Discussion in main text).
Within both datasets, the structure-function association strengths are not uniform, but vary widely from tract to tract (Supplementary Figure 6a). In the dHCP dataset, the five statistically significant MD associations had the largest correlation values and were nociception-relevant. To test if this tract-to-tract variability in structure-function association strength was consistent between datasets, we correlated the association strength of individual tracts (each dataset's 16 Pearson r-values) between datasets for each diffusion parameter (Supplementary Figure 6b). The spatial autocorrelation between tracts means the observations within datasets are not independent, so we report effect sizes only and omit tests for statistical significance. Due to the presence of extreme values (e.g. medial lemniscus in FA), we exclude outliers from these assessments, where an outlier is defined as an observation with a Cook's Distance D>1 8 (Supplementary Figure 6b, yellow dots). Finally, we square the Pearson correlation to provide an estimate of shared variance (coefficient of variation). There was no meaningful similarity between datasets in tract-to-tract variability in association strength for FA (r 2 =0.0044) or MK (r 2 =0.019), suggesting the global positive associations for FA lack tract specificity. In contrast, there is a strong similarity between datasets for MD (r=0.62, r 2 =0.38), suggesting the global negative association for MD exhibits a degree of tract specificity i.e. the tracts with strongest associations between MD and noxious-evoked response amplitudes are relatively consistent between datasets.
Lastly, we examined the relative functional importance to noxious-evoked response amplitudes of the MD of the five white matter tracts identified in the dHCP dataset exploratory analysis. The five tracts were the superior and anterior thalamic radiations, corticospinal tract, uncinate fasciculus, and forceps minor (Supplementary Figure 6a, red box). As described in the main text, in our noxious stimulation paradigm dataset, the first principal component of MD across these five tracts (MD PC1) explained over 20% of the variance in noxious-evoked response amplitudes (Figure 6b; r 2 =0.21). We These results suggest that global correlation polarities of MD and FA may reflect global brain maturity effects 9 , while local tract-to-tract variability in MD correlation strengths may reflect local maturational asynchronies 7,10,11 . Thus, analogous to our resting-state functional coupling results above, in order to understand neonatal noxious-evoked response amplitudes, both global characteristics not specific to neural processing of noxious stimulation (e.g. maturational state) and local nociception-related characteristics (e.g. microstructure of functionally relevant white matter tracts) are needed. Future work investigating the microstructural basis of these global and local structure-function association effects would be a valuable route of enquiry, and would likely benefit from complementary insights from biophysical models such as NODDI 12 .

CSF and white matter regions-of-interest definition
When associating noxious-evoked response amplitudes with both resting-state fMRI features and white matter dMRI features, the functional data were adjusted for several imaging confounds. The noxious-evoked response amplitudes were adjusted for three noxious-evoked response imaging confounds (mean head motion, stimulus-correlated head motion, CSF amplitude), and the restingstate network amplitudes were adjusted for three resting-state imaging confounds (mean head motion, CSF amplitude, and white matter amplitude). The extraction of CSF and white matter amplitudes from the functional data required CSF and white matter ROI (region-of-interest) masks.
Due to the small size of the neonatal brain, partial volume contamination is problematic. The CSF and white matter masks used in our analyses were conservative to minimize grey matter contamination, and were defined using the dHCP neonatal template atlas 13 . The CSF ROI (Supplementary Figure 8 blue) is restricted predominantly to the fluid surrounding the brainstem and cerebellum and is a single region crossing the midline. Ventricular CSF was excluded from this ROI due to the small size of neonatal ventricles. The white matter ROI (Supplementary Figure 8 red) is restricted to the largest white matter regions in the cerebrum, consisting of a left and right half tightly localized to the centre of white matter regions. In the resting-state data, these ROIs were used to extract mean timeseries for both CSF and white matter. The amplitudes of these timeseries were quantified as the timeseries MAD values. In the noxious-evoked response data, the CSF ROI was used to extract the mean regression parameter from neonates' noxious-evoked response maps, which constituted the noxious-evoked response CSF amplitude. Figure 1

(main text). For all plots, the x-axis (time in seconds) runs from 0 to 25, where t = 0 s is the point of stimulus delivery, and the minimum inter-stimulus interval is 25 s. The y-axis is arbitrarily scaled to blue curve (HRF fit) peak equals one. Grey circles are individual timepoints from all 10 trials, the grey line is the trial average timeseries, and the blue line is the HRF fit to the data. Above each plot is the z-score quantifying the goodness-of-fit between the HRF (blue line) and the data (grey line). Source data are provided as a Source
Data file.    (1-7), each step summarised using five details (a-e), for the current study. Note in step 3, the dHCP dataset (n=242) analysis and results were generated prior to and independent of the current study in a previous dHCP research output 14 Figure 6 5. Response prediction a. Local b. rs-fMRI + task-fMRI c.