Altered resting-state whole-brain functional networks of neonates with intrauterine growth restriction

The feasibility to use functional MRI (fMRI) during natural sleep to assess low-frequency basal brain activity fluctuations in human neonates has been demonstrated, although its potential to characterise pathologies of prenatal origin has not yet been exploited. In the present study, we used intrauterine growth restriction (IUGR) as a model of altered neurodevelopment due to prenatal condition to show the suitability of brain networks to characterise functional brain organisation at neonatal age. Particularly, we analysed resting-state fMRI signal of 20 neonates with IUGR and 13 controls, obtaining whole-brain functional networks based on correlations of BOLD signal in 90 grey matter regions of an anatomical atlas (AAL). Characterisation of the networks obtained with graph theoretical features showed increased network infrastructure and raw efficiencies but reduced efficiency after normalisation, demonstrating hyper-connected but sub-optimally organised IUGR functional brain networks. Significant association of network features with neurobehavioral scores was also found. Further assessment of spatiotemporal dynamics displayed alterations into features associated to frontal, cingulate and lingual cortices. These findings show the capacity of functional brain networks to characterise brain reorganisation from an early age, and their potential to develop biomarkers of altered neurodevelopment. in individual functional brain IUGR theory We further characterised functional spatiotemporal dynamics and assessed network nodes with altered temporal IUGR but also on occipital, parietal and sub-cortical regions. Nodal strength of RSN long-term effects of this The specific connectivity fingerprint could underlie the increased prevalence of IUGR subjects developing disorders of neural development some the links of IUGR with neurodevelopmental individualised biomarkers of the long-term prognosis of infants with IUGR in order to be able to clinically intervene. Although multiple comparisons problems could preclude interpretations of direct partial correlations of network features with NBAS, the ordinal regression of functional network features and NBAS severity score is especially relevant as it is robust to this problem. These results are in line with previous findings reporting an association of network features of structural brain networks with neurodevelopment in IUGR (Batalle Batalle In addition, MRI in fetal functional brain early biomarker We combination of multi-modal features of brain networks improve the assessment of the risk of neurodevelopmental problems of prenatal origin since an early age, being its translation to the clinical practice in the mid-term horizon.


M A N U S C R I P T
A C C E P T E D ACCEPTED MANUSCRIPT Batalle et al. 2 Title: Altered resting-state whole-brain functional networks of neonates with intrauterine growth restriction

INTRODUCTION
Intrauterine growth restriction (IUGR) affects 5-10% of all pregnancies in developed countries and it is a major public health issue, being associated with short-and long-term neurodevelopmental and cognitive dysfunctions (Arcangeli, Thilaganathan, Hooper, Khan, & Bhide, 2012;Baschat, 2013;Løhaugen et al., 2013). The characterisation of underlying brain alterations supporting these dysfunctions and the prediction of the subset of the population with a higher risk of altered neurodevelopmental outcomes are among the challenges of modern fetal medicine and paediatrics. Magnetic resonance imaging (MRI) has been used to characterise structural brain alterations underlying neurodevelopmental dysfunctions of subjects with IUGR at different stages of development, starting in-utero (Egaña-Ugrinovic, Sanz-Cortes, Sanz-Cortes et al., 2013), persisting at neonatal and early infancy (De Bie et al., 2011;Dubois et al., 2008;Esteban et al., 2010;Lodygensky et al., 2008;Padilla et al., 2011;Tolsa et al., 2004) and at adolescence (Martinussen et al., 2009;Skranes et al., 2005). In the recent years, the knowledge of structural brain organisation has significantly advanced with the assessment of the macroscopic circuitry of connections of the brain with structural brain networks obtained from MRI (Hagmann, 2005;Sporns, Tononi, & Kotter, 2005). Importantly, graph theoretical features have been used to characterise brain networks (Bullmore & Sporns, 2009), allowing to comprehensibly describe with a few network features the underlying brain connectivity organisation. This approach has been demonstrated to be useful to characterise a wide-range of pathologies and conditions that affect brain connectivity (Bassett & Bullmore, 2009). Based on anatomical and diffusion MRI, this technique has been promising in the study of IUGR, allowing to demonstrate alterations in the structural brain network organisation and its association with altered neurodevelopment in one- year-old infants (Batalle et al., 2012;Batalle et al., 2013), school-age infants (Fischi-Gomez et al., 2014), and in an animal model of long-term IUGR (Batalle et al., 2014). However, it remains unknown if there is brain reorganisation at a functional level in this population, and if it can be detected at neonatal age.

M A N U S C R I P T
A C C E P T E D ACCEPTED MANUSCRIPT Batalle et al. 5 Since the seminal study of Biswal et al. (Biswal, Yetkin, Haughton, & Hyde, 1995), the potential of low-frequency components of resting-state functional MRI (rs-fMRI) to obtain whole-brain functional brain networks based on partial correlations of blood oxygen level-dependent (BOLD) signal (Salvador et al., 2005) has been demonstrated. Several studies have demonstrated the feasibility to use rs-fMRI to characterise the functional organisation of the healthy neonatal brain, opening the opportunity to characterise also the alterations in brain organisation due to prenatal conditions, such as IUGR. Using independent component analysis (ICA), the emergence of synchronised spontaneous low-frequency rs-fMRI BOLD signals exhibiting resting state networks (RSN) has been demonstrated in term and preterm infants during light sedation and natural sleep (Fransson et al., 2009;Fransson et al., 2007). Both ICA and seed-based correlation approaches have also been used in longitudinal studies showing the emergence of connections partially or completely matching several RSN during neonatal development (Doria et al., 2010;Gao et al., 2009;Lin et al., 2008;Smyser et al., 2010). However, studies considering whole-brain functional brain networks of the neonatal brain are scarce in the literature. Neonatal networks composed of selected regions of interest (ROIs) were studied by Gao et al. (2009), while voxel-wise networks obtained in a normalised space were obtained by Fransson et al. (Fransson, Aden, Blennow, & Lagercrantz, 2011), showing the presence of cortical hubs and sub-networks associated with these hubs. Finally, Gao et al. (2011) studied the normal evolution of ROI-based functional brain networks from neonatal age to two years of age and its resilience to random attacks, and recently van den Heuvel et al. (2015) studied the evolution of both structural and functional connectivity during preterm brain development.
In the present study we used partial correlations of rs-fMRI BOLD signals averaged into 90 regions of an anatomical brain atlas (Tzourio-Mazoyer et al., 2002) in 13 controls and 20 subjects with IUGR scanned around 44 weeks equivalent post menstrual age (PMA). Using the whole-brain functional networks obtained we characterised alterations in the individual functional brain connectivity of neonates with IUGR using graph theory features. We further characterised functional spatiotemporal dynamics and assessed network nodes with altered temporal characteristics.

M A N U S C R I P T
A C C E P T E D ACCEPTED MANUSCRIPT Batalle et al. 6 Finally, the association of individual network features with neonatal neurobehavioral outcomes was also assessed.

M A N U S C R I P T
A C C E P T E D

Participants, neurobehavioral assessment and MRI acquisition
The infants of the study were part of a larger prospective research program in IUGR involving fetal assessment and short-and long-term postnatal follow-up at Hospital Clínic (Barcelona, Spain).
The local Ethics Committee approved the study protocol, and written informed consent was obtained from the parents or legal guardians of all the participants (CEIC: 2012/7715). The original sample of the study included a sample of 45 pregnancies with 30 late-onset IUGR and 15 control fetuses. Late-onset IUGR was defined as those fetuses with estimated fetal weight below the 10 th centile according to local reference standards (Figueras et al., 2008) confirmed at birth and delivered after 34 weeks of pregnancy. Control subjects were sampled from general pregnant population and defined as fetuses with fetal estimated weight between 10 th and 90 th centile confirmed at birth. Infants with chromosomal, genetic or structural defects and signs of intrauterine infection or neonatal onset sepsis were excluded from the study. Neonatal data was prospectively recorded including gestational age (GA), birth weight, gender, Apgar at 5 min, umbilical artery pH, neonatal complications and maternal smoking status during pregnancy.
MRI was performed around one month corrected age during natural sleep after feeding the baby. An expert neuroradiologist analysed anatomical images and those infants without overt brain lesions were further analysed. In addition, all acquired images were visually inspected for apparent artefacts and subjects excluded accordingly. Thus, ten of 30 cases and two of the 15 controls were excluded from the study due to white matter (WM) lesions, awakening or excessive movement during acquisition, obtaining a final sample of 13 controls (5 males) and 20 subjects with IUGR (13 males). Characteristics of the population are described in Table 1. Particular care was taken in order to ensure neonatal welfare during the MR acquisition. A pulseoximetry probe was placed around the baby's wrist to monitor oxygen saturation levels throughout the scan, and acoustic noise was minimised with the use of neonatal ear muffs (MiniMuffs ® Neonatal Noise Attenuators, Natus Medical Incorporated, USA). The infant was swaddled with one or two infant sheets before being placed within a vacuum immobiliser, air was removed from the bag and the infant was contained M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT Batalle et al. 8 within a rigid cradle that is shaped to its body, effectively swaddling the infant. The lighting in the scanner room was reduced to aid the infant's sleep but was kept at a level that allows safely monitoring.
Neonatal neurobehavioral performance was assessed at neonatal age with the Neonatal Behavioral Assessment Scale (NBAS) (Nugent & Brazelton, 2000), which evaluates cortical and subcortical functions in 35 items grouped into 6 clusters: habituation, social-interactive, organisation of state, regulation of state, autonomous nervous system and attention (Sagiv et al., 2008). Cluster scores were transformed to z-scores according to a standard population (Costas Moragas, Fornieles Deu, Botet Mussons, Boatella Costa, & de Caceres Zurita, 2007;Sagiv et al., 2008) and defined as abnormal if they have a z-score below minus one. NBAS severity score was defined as the number of abnormal NBAS clusters for each subject. Due to the low successful rate of a valid habituation cluster (12 out of 33), it was excluded from the calculation of NBAS severity score.
Therefore, severity score was only calculated on 29 out of 33 subjects that had valid estimation of social-interactive, organisation of the state, regulation of the state, autonomous nervous system and attention clusters, in a range from 0 to 5.

Pre-processing and network extraction
All anatomical T2-weighted volumes were first skull-striped using BET (Smith, 2002). Brain tissue was segmented into WM, grey matter (GM) and cerebrospinal fluid (CSF) with unified segmentation model (Ashburner & Friston, 2005) available with SPM8. The default adult templates were replaced with specific neonatal tissue probability maps (Shi et al., 2011). AAL atlas (Tzourio-Mazoyer et al., 2002) recently adapted to neonatal population in a T2-weighted template (Shi et al., 2011) was used to parcellate each subject into 90 cortical and sub-cortical regions (Supplementary  (Warfield et al., 2002) was used to obtain an elastic transformation matching the template with each subject's T2 volume. AAL labels were propagated to each subject using this elastic transformation with discrete labelling preserved by nearest neighbour interpolation.
Image pre-processing of BOLD images was mainly performed with SPM8 package. First, intravolume time differences between slices were corrected. Inter-volume geometric displacements were corrected using a six-parameter rigid transformation for each acquired volume and spatially smoothed using a Gaussian kernel with 2 mm full width at half maximum. Correction of head motion effects in the signal was performed by regressing out the 6-parameter head motion profiles previously estimated in each voxel across all the acquisition time. Root mean square frame displacement (FD-RMS) between each volume was estimated using average of rotation and translation parameters differences (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012) using matrix RMS formulation , as available in FSL5.0 Motion Outliers tool. Number of outliers in each subject (>75 th percentile + 1.5 × InterQuartile Range) was also recorded. GM/WM/CSF segmentation and AAL ROI parcellation obtained in T2 anatomical volume were registered to an averaged BOLD volume using an affine transformation. The representative averaged time series corresponding to each ROI were estimated for those voxels belonging to the GM mask and band pass filtered (0.01 -0.15 Hz). This "broad band" has been suggested to be the more reliable for graph theory analysis (Braun et al., 2011). Network edges were calculated as the partial correlation Batalle et al. 10 coefficients obtained between each pair of ROI averaged signal excluding the effects of the signal of the other 88 ROIs, obtaining a 90 x 90 partial correlation matrix for each subject. Pearson coefficients were transformed according to Fisher z-transformation and negative correlation coefficients were excluded, i.e., thresholding the networks with ρ>0. The correlation between each pair of ROIs excluding the effect of the signal common to the rest of ROIs was therefore inferred to be proportional to the connectivity of each given pair of ROIs, obtaining weighted matrices that represent the raw connectivity of RSN of each subject, as shown in Figure 1A. Note that only positive correlations were included

Network normalisation
In order to disentangle network infrastructure from network organisation, i.e., evaluate the organisation of networks independently of their average strength and density (cost), three different approaches were followed. First, different average strength among subjects was neutralised by means of normalisation of each subject's brain network by its total energy (Batalle et al., 2014), that is, given a network C with weights ‫ݓ‬ , for each pair of nodes i, j, we defined a normalised version of the network C with normalised weights defined as ‫ݓ‬ , = ‫ݓ‬ , ∑ ‫ݓ‬ , ∀, ⁄ . Secondly, differences of network density (cost) were neutralised by means of a cost-corrected analysis of network features (Achard & Bullmore, 2007). Following this approach, a binary network was created at different network costs, selecting for each cost-value x the connections with a strongest weight that yield to a cost x. The cost range was limited by the minimum value of network density for any subject included into the analysis, which is the maximum network cost where is possible to fairly compare among all the subjects. Finally, the effect of differences in strength and density were neutralised at the same time by means of a combination of both methods: instead of obtaining a binary network at each network cost, as it was performed in the second approach, the weighted network obtained at each network cost was considered in this third approach. Given that weighted features are much related to the strength of a network (weighted cost), to obtain pure organisational descriptors, the weighted networks obtained at each network cost were normalised as described in the strength normalisation approach. Particularly, for a given network ‫ܥ‬ሺ݅, ݆ሻ = ‫ݓ‬ , where ‫ݓ‬ , Batalle et al. 11 represents the weight of the connection between nodes i and j, we defined a cost-corrected network ‫ܥ‬ ௪ ሺ݅, ݆, ‫ݔ‬ሻ for each value of cost ‫ݔ‬ . Note that ‫ܥ‬ ௪ ሺ݅, ݆, ‫ݔ‬ሻ = ‫ݓ‬ , ሺ‫ݔ‬ሻ where ‫ݓ‬ , ሺ‫ݔ‬ሻ = ‫ݓ‬ , if the link between node i and j belongs to the subset of strongest connections that ensure a network density of ‫ݔ‬ and ‫ݓ‬ , ሺ‫ݔ‬ሻ = 0 otherwise. In this context, normalised cost-corrected network at cost ‫ݔ‬ was defined as ‫ܥ‬ ௪ ሺ݅, ݆, ‫ݔ‬ሻ = ‫ݓ‬ , ሺ‫ݔ‬ሻ ∑ ‫ݓ‬ , ሺ‫ݔ‬ሻ ∀, ⁄ .

Network analysis
Graph theory features allow summarising infrastructure and organisation of a brain network represented as an adjacency matrix (binary or weighted). Global functioning of each network was assessed by its infrastructure (average strength), integration (global efficiency) and segregation (local efficiency). Regional characteristics were evaluated by means of nodal strength, assessing the total connectivity of a node in a given network, and nodal efficiency, measuring the efficiency of the sub-network associated to a given node. Nodes with a high nodal efficiency indicate a high tolerance of the network to the elimination of the given node, which is associated to a high clustering of the neighbourhood of this node (Achard & Bullmore, 2007). Formulation and calculation of the graph theory features used to assess each network was based on the definitions and code compiled by Rubinov and Sporns (Rubinov & Sporns, 2009).

Dynamic functional connectivity
Dynamic functional connectivity (DFC) has been defined as the functional connectivity over a sliding time window (Sakoglu et al., 2010). Recently, phase synchronisation has been used to solve the resolution/reliability trade-off of windowing the signal in functional MRI time-series by means of conversion of the real signal into a complex analytic version, with techniques such as Hilbert transform. In the present study, we used similar techniques to study the temporal dynamics of the Note that as ∆߮ ሺ‫ݐ‬ሻ is constrained between 0 and ߨ, ‫ݓ‬ ு ሺ‫ݐ‬ሻ will be constrained between 0 and 1, being ‫ݓ‬ ு ሺ‫ݐ‬ሻ = 1 when the phase of the two regions is identical.
Using this approach, a weighted connectivity matrix based on the phase similarity was constructed at each instant of time, allowing the extraction of temporal-dependent global and regional network features that assess dynamic functional connectivity organisation.
However, given that the subjects are not receiving a specific stimulus that would allow direct correspondence with temporal points and specific brain state, features obtained at each time point are not easily comparable among subjects. Hence, we used Kuramoto order parameter to reorder temporal axis according to the relative level of global synchronisation of each subject, from its lowest to its highest level. Therefore, although after reording we still cannot give interpretation to each time point, overall they can now be interpreted as a percentage of time points showing the M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT Batalle et al. 13 same percentage of synchronisation relative to each subject, this way it is possible to reorder the temporal features obtained for each subject and can be compared as they are obtained based on the same relative level of synchronisation.

Statistical analysis
Comparisons among groups were performed by general linear models (GLM) with gender and smoking status of the mother as a co-factor and GA and PMA at MRI as co-variables. Significance was declared at p<0.05 (uncorrected). Regional alterations were shown using BrainNet viewer (Xia, Wang, & He, 2013

Resting-state networks in IUGR neonates
Infrastructure of the raw weighted partial correlation functional brain networks obtained ( Figure   1B-D) was assessed by comparing weighted graph theoretical features among groups by means of GLM. Importantly, this analysis showed significantly increased values of IUGR average strength (p=0.013), suggesting an increased pattern of weighted connectivity in IUGR networks. As expected of more strongly connected networks, weighted measures of global (p=0.015) and local efficiency (p=0.028) were also increased in IUGR. Network density, however, did not significantly differ between cases and controls.
Analysis of the pure organisational components of brain networks is strongly influenced by differences in the network infrastructure among subjects (Ginestet, Nichols, Bullmore, & Simmons, 2011). In order to assess the effectiveness of the organisation between groups independently of average strength, the effect of different values of average strength among subjects was neutralised M A N U S C R I P T

A C C E P T E D ACCEPTED MANUSCRIPT
Batalle et al. 14 by means of normalisation of each subject's brain network, as previously described in section 2.4.
As a result, significantly reduced normalised local efficiency was observed in IUGR group (p=0.003, Figure 1E). Although non significantly different between IUGR and control group, differences among subjects in network density (cost) could also be influencing the analysis of network organisation when using global network features. In this case, this effect was neutralised by means of a costcorrected analysis (Achard & Bullmore, 2007) also significantly reduced for global and local efficiencies (p=0.013 and p=0.016 respectively). All three different approaches of normalisation yielded the same conclusion: although IUGR had an increased RSN infrastructure characterised by increased average strength yielding to increased raw efficiencies, its organisation was sub-optimal when compared with controls.
In addition to the analysis of global network features, in order to find the regional components that have a higher influence in the network reorganisation observed in IUGR group, differences among groups were assessed using a GLM for each nodal feature using a similar approach as previously used for global network features. Particularly, nodal strength and nodal efficiency were further assessed in its raw and normalised versions (Figure 2). We observed that nodal features of raw networks showed a pattern of alterations in IUGR for both nodal efficiency and strength in frontal areas, but also on occipital, parietal and sub-cortical regions. Nodal strength and efficiency M A N U S C R I P T

Batalle et al. 15
of normalised networks showed a reduced number of alterations in IUGR, although there were still present some significant differences in subnetworks belonging to frontal, temporal and occipital cortices and sub-cortical regions such as amygdala and hippocampus. Note, however, that the nodal alterations reported must be considered exploratory, given that the differences observed did not withstand a false discovery rate (FDR) correction (Benjamini, Krieger, & Yekutieli, 2006).

Dynamic functional connectivity
In order to assess if there were any notable differences in the DFC of IUGR, we first used Kuramoto order parameter R to analyse the global level of synchronisation of each subject's brain at each instant of time ( Figure 3A), being zero for complete asynchrony and 1 for full synchronisation. As the subject' acquisition was not under specific stimuli, rather than analysing the results in a temporal basis, the average value over time was compared. A tendency towards significance of having an increased Kuramoto order parameter was observed in IUGR group (p=0.055), showing that, overall, subjects with IUGR tend to present a more synchronised brain than their control counterparts. The Kuramoto order parameter was further used to sort the assessed, finding a set of nodes with statistically significant differences in IUGR across time ( Figure   3E). Observing the percentage of time that each nodal efficiency is significantly different in IUGR allowed highlighting four regional features altered in IUGR more than two standard deviations

DISCUSSION
Characterisation of the brain changes underlying neurodevelopmental problems in IUGR is a current challenge in modern fetal and paediatric medicine (Ment, Hirtz, & Huppi, 2009). A better understanding of the pathophysiology of this condition is essential to start developing early biomarkers to detect the infants at high risk of having altered neurodevelopmental problems.
Importantly, it has been shown that early individualised interventions significantly improves IUGR neurobehavioral performance at short-and mid-term (Als et al., 2012;McAnulty et al., 2013).
However, given the high prevalence of IUGR and the economic cost of individualised care units, selecting those IUGR infants with a higher risk is essential to appropriately advise parents and optimally use clinical resources. With this long-term goal in mind, in the present study we investigated functional brain networks in IUGR.
The results obtained showed a very specific pattern of alterations in IUGR whole-brain RSN, characterised by a hyper-connectivity of their raw networks. However, when assessing the pure organisational features by means of three different normalisation procedures, we observe a suboptimal organisation in IUGR characterised by decreased global and local efficiency. Further analysis of nodal features showed a spread pattern of regional alterations in IUGR, however, not strong enough to withstand a FDR correction. Spatio-temporal analysis of the signal supported previous results, revealing a tendency of increased overall synchronisation in IUGR, and a set of nodal features that might have an important role in the reorganisation of functional brain networks obtained. These results, together with previously reported alterations in IUGR structural brain networks (Batalle et al., 2012) support the hypothesis that IUGR is a condition strongly associated with brain organisation, as has been suggested to happen with mood and psychiatric disorders.
Albeit it is important to note that establishing a causality with the evidence available is still premature, association of network features with neurobehavioral performance is reported in this and other studies (Batalle et al., 2012;Batalle et al., 2013), postulating IUGR as a candidate to be a brain-network disorder (Rubinov & Bullmore, 2013).

M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT
Interestingly, in a previous study of structural brain networks in one-year-old population, reduced fractional anisotropy (FA) weighted global and local efficiency was found in IUGR (Batalle et al., 2012). Intuitively one would expect to obtain reduced efficiencies in functional brain networks due to a weaker structural connectivity, but in contrary, significantly increased raw efficiencies were found in IUGR neonates. A possible explanation for this discord could be the difference of age in the population under study, as it is a critical period in terms of brain connectivity changes . Although the constraining of functional connectivity to the structural substrate has been demonstrated (Deco et al., 2013;Honey et al., 2009;van den Heuvel et al., 2015), functional and structural brain networks might also be capturing different aspects of brain organisation, and they do not necessarily need to behave in the same manner (Cabral, Hugues, Kringelbach, & Deco, 2012;. Note also that in a previous study on healthy infants, Gao et al. (2011) showed that, from neonatal age to one year of age, more connections presented a reduction of connectivity than those presenting an increase. Thus we cannot discard a possible retard in maturation as an explanation of the general increase in IUGR average strength.
IUGR has been suggested to be a risk factor of developing disorders such as attention deficit hyperactive disorder (ADHD) (Heinonen et al., 2010;Linnet et al., 2006), autism spectrum disorders (ASD) (Gardener, Spiegelman, & Buka, 2011;Moore, Kneitel, Walker, Gilbert, & Xing, 2012) and schizophrenia (Nielsen et al., 2013). Although there are some contradictory reports, generally ASD has been characterised by having functional hyper-connectivity of salience and default mode network (Menon, 2013). However, whole-brain RSN have been reported to show reduced raw local efficiency but increased global efficiency after cost-correction (Rudie et al., 2012). In adult schizophrenia patients, reduced global and local efficiency after cost-correction of RSN has been reported (Liu et al., 2008), but also decreased average functional connectivity (Lynall et al., 2010).
RSN of children with ADHD have been reported to have increased local efficiency for several costs, although raw weighted measures have not been described (Wang et al., 2009). Overall, after comparison with the body of literature we note a unique pattern of alterations in neonatal IUGR brain networks. This pattern of alterations must be confirmed at a later age, being the follow-up of IUGR population a crucial aspect for the characterisation of the evolution of the alterations and M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT Batalle et al. 19 long-term effects of this condition. The specific connectivity fingerprint trajectories that could underlie the increased prevalence of IUGR subjects developing disorders of neural development are not clear, and will be the likely focus of future prospective studies. We hypothesise that some of the links of IUGR with neurodevelopmental disorders could be partially associated with dysfunctions in some specific regional sub-networks produced by brain reorganisation. It is worthwhile to note that although partial correlation networks and temporal dynamics yield to same global functioning conclusions (IUGR has a hyper-connected, hyper-synchronised brain but with a sub-optimal organisation), the regional results obtained by each of these techniques showed different altered regions, suggesting that they assess complementary features of brain organisation, although both highlighted the role of frontal areas. Particularly, analysis of temporal dynamics allowed to obtain a reduced set of regions altered, comprised by frontal, cingulate and lingual cortices. In agreement with the results obtained, several frontal regions with altered structural brain network features have been previously reported in IUGR (Batalle et al., 2012), while evidence of alterations in frontalposterior networks in ASD has been reported in several studies (Maximo, Cadena, & Kana, 2014).
Interestingly, altered nodal efficiency of left lingual gyrus has been specifically reported in structural brain networks of IUGR (Batalle et al., 2012;Batalle et al., 2013) and has also been shown to be decreased in RSN of ADHD patients (Wang et al., 2009). Concerning cingulate cortical areas, FA measures of this area have been strongly associated with altered neurobehavioral performance in a long-term rabbit model of IUGR (Illa et al., 2013). In addition, functional alterations of cingulate areas have been reported in different studies of ASD (Maximo et al., 2014), schizophrenia (Lynall et al., 2010) and ADHD (De La Fuente, Xia, Branch, & Li, 2013).
Finding alterations in the brain network features associated with pathology is by itself relevant for the characterisation of its pathophysiology. In the case of IUGR, the results suggest that brain reorganisation previously demonstrated in the structural substrate is also present at a functional level since neonatal age. This significantly improves the knowledge of this condition, serving as a potential physiological basis for the neurobehavioral alterations reported in these infants (Figueras et al., 2009). Note that given the heterogeneity of the aetiology and progression of this condition, big samples are needed to find differences in neonatal neurobehavior, being even more crucial to M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT Batalle et al. 20 find individualised biomarkers of the long-term prognosis of infants with IUGR in order to be able to clinically intervene. Although multiple comparisons problems could preclude interpretations of direct partial correlations of network features with NBAS, the ordinal regression of functional network features and NBAS severity score is especially relevant as it is robust to this problem. These results are in line with previous findings reporting an association of network features of structural brain networks with neurodevelopment in IUGR (Batalle et al., 2012;Batalle et al., 2013). In addition, recent reports showing the feasibility to assess functional MRI in fetal period reinforce the potential to use functional brain networks as an early biomarker (Thomason et al., 2013). We are confident that combination of multi-modal features of brain networks will significantly improve the assessment of the risk of neurodevelopmental problems of prenatal origin since an early age, being its translation to the clinical practice in the mid-term horizon.
Finally, there are several issues of the study that must be noted. First, we would like to note that the MRI acquisition was performed during natural sleep. Previous reports have associated deepness of sleep with functional connectivity (Horovitz et al., 2009;Larson-Prior et al., 2009); however, although the effects of sedation and sleep in functional connectivity are not fully understood, significant differences have not been found between sedated and non-sedated infants using ICA and seed-based correlation approaches (Doria et al., 2010;Fransson et al., 2011;Fransson et al., 2009). Importantly, neonatal sleep has been suggested to be mainly in active sleep (Biagioni et al., 2005), minimising the possibility that different sleep deepness could be partially explaining some of the results obtained. Regarding the comparability of the results obtained, two previous studies use ROI-parcellation to study whole brain functional brain networks in neonates (Gao et al., 2011;van den Heuvel et al., 2015). Franson et al. (2011) also assessed whole-brain networks in a neonatal population, however, it was performed voxel-wise in a normalised space, obtaining very large networks (4966-by-4966 elements). Albeit the use of voxel-wise networks has some advantages, constraining the networks obtained to an anatomical brain atlas allows a more comprehensible and manageable comparison among studies, especially given the broad use of AAL atlas in the literature. However, regional parcellation of the neonatal brain is also a critical issue and could be a source of bias. This was alleviated by the use of T2-weighted anatomical M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT Batalle et al. 21 volumes which improve WM contrast in neonates (Williams et al., 2005), and by the use of a specific neonatal atlas (Shi et al., 2011). Note that although we correct for gender in the statistical model, in our population the percentage of males is much higher in IUGR group. Previous studies have described gender-related differences in brain connectivity (Ingalhalikar et al., 2014;Tian, Wang, Yan, & He, 2011), although the specific effects of gender differences on functional brain connectivity at such an early age are still largely unknown. We acknowledge the relatively small sample size available for this study. In order to assess the reliability of the results obtained, we applied a resampling approach allowing us to reassess the differences between cases and controls for global network features. After applying bootstrap with 10000 samples we observed that most of the differences reported are still statistically significant after resampling (Supplementary Table 2).
Although the effects associated with IUGR appeared to be consistent when using different approaches, statistical power is limited by this reduced sample size. Further effort to study this very challenging population is needed to confirm the results here presented. With respect of motion during scan, the estimation of motion for both groups was comparable with previously reported for infants and adolescents (Pruim, Mennes, Buitelaar, & Beckmann, 2015), and importantly, we didn't find any significant differences in motion between cases and controls (Table 1), minimising the risk that any of the results reported could be related to different levels of motion in the two groups under study. We acknowledge that the robustness to motion of the approach used to analyse the data, i.e., constructing networks by means of partial correlations within the averaged signal among ROIs, needs to be further assessed (Power, Schlaggar, & Petersen, 2015). However, it should inherently remove signals present throughout the brain, therefore be robust to jerk motion typical in neonatal acquisitions. As the signal of the 88 remaining ROIs is used as a confounder in each of the partial correlations, any signal global to all the brain will be neglected as a possible source of pair-wise correlation. From a specific graph theoretical point of view, partial correlations have been shown to be an effective way to remove head motion effects (Yan, Craddock, He, & Milham, 2013). In addition, the effect of low motion, will also be further compensated by voxel-wise regressing out of the signal the 6 rigid-body features previously used to correct inter-volume motion. We further assessed the effect of motion by adding average FD-RMS and number of frames with high motion M A N U S C R I P T A C C E P T E D ACCEPTED MANUSCRIPT Batalle et al. 22 (outliers) as a covariate in the GLM assessing group effects in main network features, without finding any relevant difference with reported results (Supplementary Table 3). Furthermore, we also tested censoring of frames with high motion and regressing them out from correlations when computing connectivity matrices (Supplementary Material), also confirming the main findings reported, suggesting that the motion during acquisition did not have a relevant effect in the graph theory characteristics here reported. Finally, since haemodynamic adaptation in IUGR occurs with blood flow redistribution preferentially to the brain, i.e. the brain sparing effect, we could not discard that part of the changes observed could be related with changes in brain vasculature. Since most of IUGR babies included in our study did not have brain vasodilation, as assessed by Middle Cerebral Artery pulsatility index, we could expect that differences would not be related with changes in neurovasculature. However, early changes in brain blood redistribution could also affect pulsatility index (Garcia-Canadilla et al., 2014), not allowing us to discard this effect in our population.

CONCLUSIONS
In conclusion, the results presented show the feasibility of using functional brain networks at neonatal age to characterise alterations of prenatal origin. Using IUGR as a model of prenatal condition allowed finding a unique pattern of alterations in the functional brain network organisation, associated with neurobehavioral scores. Overall, the observed functional reorganisation in IUGR neonates could be a potential substrate of altered neurodevelopment in infants with IUGR, and together with previous findings, postulate IUGR condition as a possible brain-network disorder.
Association of network features with neurobehavior since neonatal age opens the opportunity to develop early image biomarkers of altered neurodevelopment, a clinical chance to improve the management of a condition that affects up to 10% of the population.   Kuramoto order parameter reordered according to relative value for each subject, from its lowest to its highest value (B). Global (C) and local efficiency (D) at each time point reordered according

M A N U S C R I P T
A C C E P T E D ACCEPTED MANUSCRIPT