Neonatal brain dynamic functional connectivity: impact of preterm birth and association with early childhood neurodevelopment

Brain functional dynamics have been linked to emotion and cognition in mature individuals, where alterations are associated with mental ill-health and neurodevelopmental conditions (such as autism spectrum disorder). Although reliable resting-state networks have been consistently identified in neonates, little is known about the early development of dynamic brain functional connectivity and whether it is linked to later neurodevelopmental outcomes in childhood. In this study we characterised dynamic functional connectivity with functional magnetic resonance imaging (fMRI) in the first few weeks of postnatal life in term-born ( n = 324) and preterm-born ( n = 66) individuals and evaluated whether early dynamic functional connectivity: i) changes with age in the neonatal period ii) is altered by preterm birth and iii) is associated with neurodevelopmental and behavioural outcomes at 18 months. Global brain dynamics in preterm-born infants were atypical when compared with term-born controls, and associated with atypical social, sensory, and repetitive behaviours measured by the Quantitative Checklist for Autism in Toddlers (Q-CHAT) scores at 18 months of age. On a modular scale, we identified six transient states of neonatal dynamic functional connectivity: three whole-brain synchronisation states and three regionally constrained states occupying occipital, sensorimotor, and frontal regions. Modular characteristics of these brain states were correlated with postmenstrual age and postnatal days at scan. Preterm-born infants had increased occurrence of frontal and occipital states. Higher neonatal sensorimotor synchronisation was associated with lower motor and language outcome scores at 18 months. Lower frequency of occurrence of whole-brain synchronisation states and higher frequency of occurrence of the sensorimotor state were associated with higher Q-CHAT scores at 18 months. Overall, we show that a dynamic landscape of brain connectivity is already established by the time of birth in the human brain. This landscape is altered by preterm birth and its profile is linked to neurodevelopmental outcomes in toddlerhood.


Introduction
A fundamental feature of the brain's activity is the dynamic properties of its functional connectivity, with continuous shifting between different connectivity profiles or 'states'. These dynamic properties are linked to processes such as language 1,2 , cognition 3-6 , and motor function 7,8 . Importantly, altered brain dynamics have also been linked to the clinical features and/or cognitive dysfunction in neurodevelopmental conditions such as schizophrenia 9 , attention deficit hyperactivity disorder (ADHD) 10 and autism spectrum disorder (ASD) 11 . Individuals with ASD, for example, have been reported to switch between different connectivity profiles more directly, whereas typically developing individuals switch between those same brain states via an intermediate connectivity profile 11 . However, although it is increasingly appreciated that neurodevelopmental conditions likely have their origins in the perinatal period, little is known about the brain's dynamic properties at this critical juncture. Moreover, it is also not clear whether its dynamic characteristics are associated with later childhood neurodevelopment, and in particular, which of these characteristics signal a higher likelihood of later neurodevelopmental difficulties.
Conventional studies of "static" functional brain connectivity in early life have shown that the spatial representation of resting state networks (RSNs) appears to be relatively mature and adult-like even soon after birth [12][13][14] . It has also been confirmed that these RSNs may be disrupted by perinatal exposures, such as preterm birth 12,15-17 , which increase the likelihood of neurodevelopmental impairments 18 . However, despite the evidence for the importance of dynamics in older groups, the dynamics of these networks in early life remain to be described.
In this study we applied state-of-the-art techniques to evaluate functional Magnetic Resonance Imaging (fMRI) in a cohort of term-born (n=324) and preterm-born babies (n=66) scanned at term equivalent age from of the developing Human Connectome Project (dHCP), the largest publicly available population-based dataset of the healthy new-born brain 19 . We used two methods that tap into dynamic brain function. First, we characterised global dynamics using Kuramoto Order Parameter (KOP) 20 based measures, namely mean synchronisation and metastability 21,22 as global measures of brain synchronicity and flexibility 21,23 . Second, we characterised modular dynamics. That is, we identified subnetworks involved in temporal 'states', i.e., paroxysmal modes representing synchronisation of the brain, using Leading Eigenvector Analysis (LEiDA) which is a time-resolved metric 24,25 . We assessed whether neonatal brain state features (fractional occupancy, dwell times, state mean synchronisation, and state metastability) and state transition probabilities were associated with postmenstrual age (PMA) at scan, postnatal days (PND) at scan and preterm-birth; and correlate with neurodevelopmental outcomes at 18 months measured using the Bayley Scales of Infant and Toddler Development, 3 rd Edition (Bayley-III) 26 , and atypical social, sensory and repetitive behaviours measured by the Quantitative Checklist for Autism in Toddlers (Q-CHAT) 27 . This allowed us to test three hypotheses: 1) that neonatal brain dynamics quickly develop with age at scan; 2) that preterm birth alters the typical pattern of functional brain dynamics; and 3) that neonatal brain dynamics are linked to neurodevelopmental and behavioural outcome measures at age 18 months.  transitions from intermediate whole-brain synchronisation state C to the occipital state decreased with PMA at scan, while transitions in the opposite direction (from occipital state to synchronisation state C) were increased with PND at scan.  Atypical modular brain state features in preterm born neonates. Compared to term-born participants, preterm-born babies had shorter dwell times for the global state A (t = -4.6; p < 0.001, Figure 4A); decreased fractional occupancy for the global state A (t = -4.1; p < 0.001); and increased fractional occupancy for global state B (t = 3.4; p = 0.001), occipital state (t = 2.2; p = 0.030), and frontoparietal state (t = 2.7; p = 0.009) ( Figure 4B). Preterm birth was also associated with lower mean synchronisation of global state A (t = -5.1; p < 0.001), global state B (t = -3.4; p = 0.001), global state C (t = -2.8; p = 0.005), occipital state (t = -2.2; p = 0.031), and frontoparietal state (t = -2.6; p = 0.010)) ( Figure 4C); and reduced metastability for the global state A (t = -2.5; p = 0.014) and frontoparietal state (t = -4.3; p < 0.001) ( Figure 4D).
Effect of preterm birth on brain dynamics and brain state transition probability. Preterm birth was associated with an increased transition towards an occipital connectivity profile, i.e., an increased transition probability from global state A to C (t = 4.7; p = < 0.001) and from global state C to occipital (t = 3.1; p = 0.002); as well as a reduction in the probability staying (dwelling) in global state A (t = -4.4; p < 0.001); see Figure 4F. A similar profile of results was obtained when assessing how brain dynamic features changed with gestational age at birth which essentially captures preterm and term birth in a linear way, see Supplementary Figure S2. Significant changes in transition probabilities associated with preterm-birth and their relation to associations with PMA and PND at scan are summarised in Figure 4H.  Neonatal brain states and early childhood neurodevelopment. Lastly, we investigated associations between brain states dynamics after birth and neurodevelopmental outcomes at 18 months of age measured with Bayley-III and Q-CHAT scores. Higher mean synchronisation of the sensorimotor state in neonates was associated with poorer performance on both Bayley's cognitive and motor scores assessed at 18 months of age (t = -3.2, p = 0.002; and t = -3.4, p = 0.002 respectively; Figure 5A). Higher Q-CHAT scores at 18 months of age were associated with reduced mean synchronisation (t = -2.7, p = 0.008) and higher fractional occupancy of sensorimotor state and reduced fractional occupancy of global state A (t = 2.5, p = 0.014; and t = -2.6, p = 0.010 respectively; Figure 5B).  A similar pattern of regional brain dynamics was obtained with an alternative parcellation scheme. The six brain states obtained with the alternative M-CRIB neonatal atlas were compatible with the ones obtained for the AAL atlas (Supplementary Figure S3 with three global synchronisation states and three more regionally constrained. States 1, 2 and 3 obtained with the M-CRIB atlas were consistent with global synchronisation states obtained with the AAL atlas. The M-CRIB's State 4 featured some of the structures present in occipital state, State 5 featured similar structures to the sensorimotor state, and State 6 was concordant with FP state obtained from AAL. Moreover, for the analyses using the M-CRIB atlas in term-born babies (Supplementary Figure S4), PMA was also associated with increased fractional occupancy (t = 2.9; p = 0.002) and mean synchronisation (t = 3.0; p = 0.003) for State 5 (~SM); and PND was associated with reduced dwell times (t = -2.5; p = 0.012) for State 5 (~SM). Transitions from States 6 (~FP) to State 2 (~Glb.) also increased with PND (t = 3.3; p = 0.002). Pre-term birth was also associated with similar changes in state metrics for M-CRIB atlas (Supplementary Figure S5) with increased fractional occupancy for State 4 (~Occ.) (t = 2.3; p = 0.020) and State 6 (~FP) (t = 3.0; p = 0.002); and reduced mean synchronisation for State 6 (~FP) (t = -2.2; p = 0.027). See Supplementary Figure S6 for a similar analysis with GA.
Association with developmental outcomes in cognitive (t = -4.0; p < 0.001) and motor (t = -3.4; p < 0.001) components of Bayley-III showed negative associations with mean synchronisation for State 5 (~SM) and higher Q-CHAT scores were associated with reduced fractional occupancy for State 1 (t = -2.6; p = 0.008), see Supplementary Figure S7 for a summary. A detailed description of results using M-CRIB atlas for parcellation is available in Supplementary Materials.

Discussion
We used global and modular fMRI signal analysis tools to investigate the characteristics of dynamic functional connectivity in a large sample of term and preterm born neonates. This is the first study to characterise the fundamental characteristics of the repertoire of brain states and their dynamics this early in human development. We found that preterm birth disrupted brain dynamics and that the profile of brain dynamics in early postnatal life is associated with a range of early childhood neurodevelopmental and behavioural outcomes at 18 months of age.
Global dynamics. We found that global dynamic features remained relatively stable in early postnatal development (37-44 weeks PMA at scan) in a term-born population, although lower metastability was observed with more postnatal days at scan, suggesting that ex-utero life experience reduces brain dynamic flexibility after birth, and promotes more stable connectivity patterns. This is consistent with our observation that preterm born infants scanned at term equivalent age, and thus with greater exposure to postnatal life experience at the time of scan, had lower metastability than term born children. However, preterm babies also had lower mean synchronisation at term-equivalent age suggesting a unique pattern of alteration in brain dynamics associated with preterm birth which is independent of the extent of ex-utero life exposure (PND at scan). Lower metastability has been previously associated with impairments in cognitive flexibility of mature individuals after traumatic brain injury 6 . However, the negative association between duration of postnatal life and metastability observed here is unlikely to relate to cognitive flexibility, but rather reflect a refinement of network dynamics driven by primary sensory stimulation associated with ex-utero life experience. Metastability seems to be additionally reduced by term equivalent age in preterm born babies, which may be consistent with theories suggesting perinatal stress atypically accelerates brain maturation, with potentially negative long-term impact 30 .
Our study extends prior work which has observed that, later in childhood, very preterm born infants have suboptimal neural synchrony and altered global dynamic connectivity patterns when scanned later in childhood 31 . Our metrics indicate that the impact of preterm birth on brain dynamics begins much earlier in life. Such differences in preterm-born babies may have their roots in alterations in the framework of structural and functional networks reported to accompany preterm birth. For example, previous studies have shown that preterm birth is associated with altered brain structure 32,33 , global functional architecture 15,34 , and structural network changes in the neonatal period 35 which continue to be present at school-age 36 and into young adulthood 37,38 .
Global dynamic features were linked to later behaviour: lower mean synchronisation in the neonatal brain was associated with higher Q-CHAT scores at 18 months. Although high scores on the Q-CHAT indicate more autistic traits, in this study, the Q-CHAT captured a continuum of social, sensory, and repetitive behaviours across the normal distribution 27,39 . Thus, while global metrics might usefully signpost the trajectory of foundational dynamic steps of human brain development, their association with Q-CHAT should not be overinterpreted as our study was not a study of ASD and we did not examine children at an age where neurodevelopmental diagnoses begin to be formalised. Modular brain state profile in term-born neonates. Summary metrics of global dynamics are likely themselves to be underpinned by much more complex activity. Therefore, we explored the emergence and behaviour of modular brain states in the early postnatal period. We characterised six transient states in the newborn brain at term equivalent age. Amongst the three states that showed widespread concordance (global states A-C), the first encompassed nearly all of the analysed cortical regions, the second showed a higher contribution of all sensory regions (auditory, sensorimotor, and visual) and the third state showed a higher contribution of visual and frontoparietal cortices. The other three states had a more restricted/regional span with distinct contributions from sensorimotor, visual and frontoparietal regions. The majority of these states thus encompassed primary sensory networks which are already known to mature earlier than higher order networks 12,15,40 . This adds confidence to the results reported in studies of dynamic FC in neonates around birth. For example, Ma et al. (2020), described dynamic functional connectivity with four brain states that encompassed default-mode, dorsal attention, auditory, sensorimotor, and visual networks in 37 term neonates 41 . Here we establish a series of six brain states and describe associations with age at scan and postnatal days in a larger cohort with 324 term neonates.
We observed that in newborn infants, whole-brain synchrony state A had the highest mean synchronisation, as well as largest fractional occupancy and dwell times; thus suggesting that newborn infants spend a large amount of their time in a state of global phase synchrony which is a similar to the previously described dominant pattern of whole-brain synchronisation seen with both fMRI 41,42 and EEG 43 . Together with prior studies, our work supports the idea that large scale activity plays a crucial role in early brain development. Our results also align with the concept that this activity could support large scale cortical network formation and may foster the associated long-range connections, which are known to subsequently mature during the first postnatal year 44,45 .
Age-related changes in modular brain states. We observed a positive correlation of occupancy and mean synchronisation within the sensorimotor cortices with increasing PMA at scan. This supports existing evidence that this system is relatively mature, both in function and structure in comparison to other systems at birth. This may be the product of significant short-range functional reorganisation during the last foetal trimester 46,47 to support the functional specificity needed to respond to sensory information coming from feet, hands and mouth 48,49 .
We also found that increased PMA at scan, was associated with an increased probability of transitioning into frontoparietal synchronicity states comprising the anterior part of the Default Mode Network (DMN) 50 . This is concordant with other evidence that, although this system is relatively immature at birth, it undergoes significant changes postnatally with increasing recruitment of frontoparietal areas into the network 12,15 . Our work suggests that this system is recruited more and more with age in the postnatal period.
There was also a significant effect of PND at scan on an increasing probability of transition from occipital to whole-brain synchronisation states. This finding is in line with key stages of neurodevelopment. Specifically, changes in transitions from an occipital cortex profile may help the maturation of visuomotor abilities and sensory integration in early infancy in the environment outside the womb 51 . Finally, we observed increased transitions from sensorimotor to frontoparietal structures with increasing duration of postnatal life, perhaps reflecting the emerging of functional maturation of frontoparietal systems which coincides with high interneuron migrations to theses regions 52 , relative to the already mature sensorimotor systems 47 .
Effect of preterm birth on modular brain states. Preterm birth has been associated with a higher likelihood of atypical neurodevelopment 18 including a greater rate of autism diagnosis 53-58 . Previous studies from our group have shown preterm born infants have alterations in their functional architecture 12,15,34 . Here, we extend this work to report that preterm birth also has an impact on dynamic functional connectivity, increasing fractional occupancy of occipital and frontal states and increased transitions from global to occipital state; and decreased dwelling for whole-brain high synchronisation state. Only one other study has investigated the effect of prematurity on dynamic functional connectivity 41 . They reported significantly shorter mean dwell times in a state with stronger connectivity between sensorimotor and auditory cortices and significantly higher mean dwell times for a global state 41 . Direct comparison with our findings is challenging however, given the higher resolution of our study in terms of a larger number of ROIs included in our study but also the power available from our larger sample. We observed multiple global states in neonates with diverse contributions from occipital and frontoparietal regions and a negative effect of prematurity on dwell times in state A. In summary, preterm-birth shifts dynamic functional connectivity towards occipital and frontoparietal synchrony profiles and suppresses wholebrain synchronisation modes. Preterm-birth is known to impact on cognition throughout the lifespan 59,60 . Our results raise the possibility, that alterations in brain dynamic functional connectivity that are present soon after birth have functional consequences.
Association with early childhood neurodevelopment. Our work and others consistently recognise the neonatal period as a key time for sensorimotor cortical development 46,48,49 and subsequent transition to higher order network connectivity 12,15 . We extended this observation here, to show that altered brain dynamics contributes to both general developmental outcomes and more specific social, sensory, and repetitive behaviours at age 18 months. Lower levels of synchrony in the sensorimotor state around term were positively correlated with better cognitive and motor outcomes (Bayley-III) at 18 months. This association is evident in the analyses with both the AAL and M-CRIB atlases. However, when there was lower fractional occupancy of the high whole-brain synchronisation state A and increased fractional occupancy of the sensorimotor state around birth, there were more atypical social, sensory, and repetitive traits present at 18 months, as captured by the Q-CHAT in AAL atlas.
Thus, our work indicates that the link between brain dynamics and autistic traits is not only limited to state transitions in adulthood 11,61 but may comprise alterations in state occupancy and overall synchronicity as well as transitions established during early development. One possible explanation for the association between the sensorimotor cortex dynamics and later social, sensory and repetitive traits, is that altered brain dynamics in the neonatal period may predispose individuals to exhibit unusual responses to sensory stimuli 62 . The positive correlation of higher fractional occupancy of the sensorimotor state and a higher Q-CHAT score could also represent an overreliance on that particular network during the neonatal period -which may impact upon the development of higher order networks.
We emphasise that we did not follow-up the children with higher Q-CHAT scores beyond 18 months, thus our work did not evaluate predictors of a confirmed diagnosis of ASD, nor is our study about diagnosed autism. However, children who go on to receive a diagnosis of autism and those with a broader phenotype may already show emerging traits of the condition. This was the basis for development of the Q-CHAT 27,39 . There is also accumulating evidence across different modalities that infants who go on to receive a diagnosis of autism have differences in their neurobiology and physiology from as young as 6-9 months 63-68 Thus our work links emerging social cognitive profile relevant to ASD, also to dynamic functional connectivity at birth, especially within sensory networks. This correlation between higher scores on an instrument, which captures early features relevant to ASD (though not necessary diagnostic) to the dynamics of sensory systems, is in agreement with the importance of sensory processes throughout the lifespan in individuals who have an ASD diagnosis. Longitudinal studies are clearly needed, but our work adds to the evidence for a neurodiverse trajectory of sensory systems to autistic features across the lifespan. Sensory differences are among the first features to signal a diagnosis of ASD 66,69 and persist, remaining core to the diagnosis 70 . For example, we have previously reported that neonates with an increased familial likelihood of ASD had higher regional homogeneity in the sensorimotor cortex 71 , which fits with the higher fractional occupancy we record in this region. Studies of diagnosed individuals have reported atypical activation of the motor cortex 72 , which likely affects the translation of visual input into motor understanding, with a potential impact on social interaction. These functional brain differences in sensory systems are thought to arise from alterations in excitation-inhibition pathways, especially GABA neurotransmission. Evidence includes a tight relationship between sensory processing and differences in sensory cortex GABA levels 73 and we have reported direct experimental proof in adults, that visual processing differences in ASD are GABA-dependent 74 .

Strengths and limitations.
We studied state-of-the-art infant fMRI acquired with a dedicated neonatal multiband EPI pipeline which features high temporal resolution (TR = 392 ms) but acknowledge this has relatively low signal-to-noise ratio in the deep grey matter structures, including thalamus, basal ganglia, and brainstem 75 . We chose a LEiDA processing pipeline, as it provides time-resolved metrics of brain states and has been successfully applied to study brain dynamics in adults 24,25 . A majority of prior studies of dynamic brain functional connectivity have adopted sliding window approaches, which necessitates the choice of arbitrary parameters in the analysis, such as window and step sizes. There is little consensus on these parameters as the use of small windows can magnify spurious variations and large windows can soften sharper changes in brain dynamics [76][77][78] . Thus, the choice of a timeresolved approach like LEiDA is a strength in our methodology, though other techniques such as Hidden Markov Models (HMMs) are also available 79 . Our analysis treated brain state occurrence in an independent fashion, i.e., without memory. Future studies could benefit from developing a metric that considers how brain states are impacted by previously occurring states. Our results are also potentially limited by the use of the AAL atlas to define our regions of interest. While this atlas has been adapted for neonatal use 28 , structurally defined atlases could potentially poorly correlate with functional boundaries in the brain, particularly during the neonatal period 80,81 . A possible solution for future work could be the adaptation and use of multi-modal generated surface atlases for neonatal fMRI studies 82 . However, we reproduced our pipeline using an alternative parcellation scheme (the M-CRIB atlas), obtaining similar results to those featured by the neonatal AAL atlas.
The dHCP cohort focused on characterising typical development, hence this study features an unbalanced number of term vs preterm infants. Preterm-born babies included were predominantly moderate or late preterm, and mostly "healthy", with no incidental findings of clinical significance. While this may be better representative of general preterm population (80% of preterm infants are born moderate or later preterm 83 ) we cannot extrapolate our results to very or extremely preterm born babies, or to those with significant white matter damage. In our analysis, we also didn't consider the effect of socio-demographic factors and social deprivation in early neurodevelopmental outcomes beyond using the Index of Multiple Deprivation scores as a covariate. However, multiples studies have shown that family psychosocial and socio-demographic factors have a significant impact on brain development in childhood; with factors such as maternal stress, depression, low education, maternal immigration status, maternal age greater than 35 years, paternal age over 38 years and low household income all being linked to poorer developmental outcomes [84][85][86][87][88] . Thus, further studies could benefit from evaluating links between these wider socio-demographic markers and brain dynamics in early childhood 89 .

Conclusions
In this study we evaluated global functional brain dynamics and transient brain states in the newborn brain. Our approach allowed us to define a set of six fundamental transient brain states, which are comprised by structures previously shown to be established in earlier phases of brain development. We have highlighted the impact brain maturation has on brain dynamics, as well as atypical patterns associated with preterm birth. Brain state dynamics at birth appear to be functionally relevant as they are correlated with a range of neurodevelopmental outcomes in early childhood. This encourages further work to understand their prognostic value and regulation to guide support and intervention where appropriate.

Methods
Participants. We analysed a total of 390 (out of 809) datasets from the dHCP (release 3) 19 acquired from both term (n = 324) and preterm born (gestational age (GA) at birth < 37 weeks; n = 66) babies. Full inclusion and exclusion criteria and sample excluded at each step are detailed in Supplementary Figure S8. Term-born babies were born at median GA at birth of 40.14 weeks (IQR = 1.86 weeks) and scanned soon after birth (median postmenstrual age (PMA) at scan = 41.57, IQR = 2.43 weeks). Preterm-born babies were born at a median GA at birth of 33.36 weeks (IQR = 5.86 weeks) and scanned at term-equivalent age (median PMA at scan = 40.5 (IQR = 2.71)). Table 2 shows demographic data of the sample. The distribution of PMA at scan and GA at birth for the individuals included in this study is shown in Supplementary Figure S9. All children were invited to the Centre for the Developing Brain, St Thomas' Hospital, London, for neurodevelopmental evaluation by experienced paediatricians or psychologists at 18 months after expected delivery date. The Bayley Scales of Infant and Toddler Development, Third Edition (Bayley-III) 26 were used to assess general developmental outcomes across motor, language and cognitive domains in 305 individuals from the total population, comprising 257 infants born at term and 48 born preterm (higher scores indicate greater skills). The Quantitative Checklist for Autism in Toddlers (Q-CHAT) 27 at 18 months corrected age was available in 300 individuals (254 born at term and 46 born preterm), as a measure of atypical social, sensory and repetitive behaviours which occur as a continuum in the population 27 . Although higher Q-CHAT scores may indicate more threshold or subthreshold autistic traits, we emphasize our use of this instrument was to capture behaviours not tapped by the BSID-III, rather than to screen for ASD. The index of multiple deprivation (IMD) rank -a composite measure of geographical deprivation estimated from the address of the mother at the time of birth 90 -was obtained for every subject and included as a covariate for models aimed at assessing the relationship between neonatal brain dynamics and subsequent neurodevelopmental and behavioural outcomes.

MRI acquisition.
We evaluated fMRI scans obtained as part of the dHCP at the Evelina Newborn Imaging Centre, Evelina London Children's Hospital, using a 3 Tesla Philips Achieva system (Philips Medical Systems). Ethical approval was given by the UK National Research Ethics Authority (14/LO/1169), and written consent was obtained from all participating families prior to data collection. Scans were performed without sedation in a dedicated neonatal set-up with optimised transport system, positioning devices, hearing protection, and custom-built 32-channel receive head coil and acoustic hood 91 . Scans were supervised by a neonatal nurse or paediatrician who monitored heart rate, oxygen saturation and temperature throughout the duration of the scan. Blood-oxygen-level-dependent (BOLD) fMRI was acquired using a multi-slice echo planar imaging sequence with multiband excitation (factor 9) (repetition time (TR) = 392 ms, echo time (TE) = 38 ms, voxel size = 2.15 x 2.15 x 2.15 mm, flip angle = 34°, 45 slices, total time = 15 m 3 s, number of volumes = 2300) 19 . Anatomical images were acquired for brain morphometry and clinical reporting 19 . T1-weighted images had a reconstructed spatial resolution = 0.8 x 0.8 x 0.8 mm, field of view = 145 x 122 x 100 mm, TR = 4795 ms. T2-weighted images had a reconstructed spatial resolution = 0.8 x 0.8 x 0.8 mm, field of view = 145 x 145 x 108 mm, TR = 12 s, TE = 156 ms. fMRI datasets with excessive motion (more than 10% of motion outliers 15 ), or incidental MRI findings of clinical significance (radiology scores 4 or 5 which indicate major lesions within white matter, cortex, basal ganglia or cerebellum -as described in dHCP database 19 ) were excluded. In cases of twin/triplet scans only one infant was included (the one with least motion outliers during acquisition). Data processing. Individual fMRI datasets were pre-processed according to the dHCP dedicated neonatal pipeline (Fitzgibbon et al., 2020). Briefly, local distortion due to field inhomogeneity was corrected using topup 92 ; intra-and inter-volume motion correction; and associated dynamic distortions correction using rigid-body realignment and slice-to-volume eddy 93 . Residual motion, multiband acquisition, and cardiorespiratory artefacts were regressed out using FSL FIX 94 .
Head motion has been shown to produce spurious and systematic correlations in resting state fMRI 95 . In addition to specific processing steps within the dHCP pipeline implemented to minimise motion artifact on the BOLD signal 75 , we also evaluated motion magnitude for each dataset using framewise displacement (FD) 95 . To minimise the effects of motion (and/or any likely associated differences) on the determination of brain states, only participants with less than 10% motion outliers (defined as FD > 75 th centile + 1.5*IQR) were selected for the final analysed subsample. The total number of motion outliers was additionally used as a covariate of control in the analysis models described in the statistics section.
We segmented the T2-weighted volumes into 9 tissue types including white matter, grey matter, and cerebrospinal fluid with a dedicated neonatal tissue segmentation pipeline 96 . We parcellated each subject's T2-weighted volume in 90 cortical and subcortical parcels using the Anatomical Automated Labels (AAL) atlas 97 , mapped to the neonatal brain 28 , adapted and manually corrected into the dHCP high-resolution template 98 . We transformed the AAL atlas from template space into each subject's native space with a non-linear registration based on a diffeomorphic symmetric image normalisation method (SyN) 99 using T2-weighted contrast and the segmentation obtained previously. Grey matter segmentation and parcels were propagated from T2-weighted native space into each subject's fMRI space with a boundarybased linear registration available as part of the functional dHCP processing pipeline 75 . Average BOLD timeseries were then calculated for each of the 90 AAL parcels in their intersection with grey matter, deep grey matter, or basal ganglia segmentation masks, as appropriate. An alternative parcellation scheme was also used, following the same procedure but using 80 cortical regions from the M-CRIB atlas 100 , a neonatal adaptation of the Desikan-Killiany atlas 101 .
Analysing BOLD timeseries. We filtered the BOLD timeseries with a bandpass Butterworth second order filter in the range of 0.02-0.10 Hz 25 and obtained the phases ߮ j ‫)ݐ(‬ for each parcel in time with the Hilbert transform. For a given real signal ‫ݏ‬ ሺ ‫ݐ‬ ሻ , we built a complex signal z(t) 29,102 given by: In which represented a Hilbert transform applied to the real signal s(t) and is defined below, with p.v. consisting of Cauchy principle value 102 : The phases ߮ j ‫)ݐ(‬ for each parcel can be calculated directly from z(t): ݆ at a given time. In our study, each brain parcel (AAL region) is treated as an independent oscillator.
Once KOP was obtained for every time point, we obtained the mean KOP (synchrony aggregate over time, referred as "mean synchronisation"), KOP standard deviation (referred as "metastability") 21,22 . Mean KOP provides a broad measure of whole brain synchronicity, whereas metastability provides a measure of how synchronisation between different oscillators fluctuates over time, i.e., brain flexibility 6,21,23 .
The Leading Eigenvector Analysis (LEiDA). KOP analyses can provide insight on global dynamic properties over all oscillators (brain parcels); but cannot inform which specific brain structures might be involved in those changes. To evaluate such modular (local) properties we applied the LEiDA approach -which allowed us to investigate phase coherence in different sets of parcels. To do so we first calculated the phase difference between a parcel ݅ and a parcel ݆ at every instant (TR), using the cosine distance: This results in a symmetric dynamic functional connectivity matrix for each fMRI volume.
We then obtain a lower-dimensional representation with the LEiDA approach 24,25 , whereby the LEiDA vector corresponds to the first eigenvector of the decomposition of the matrix Δ ߮ ݅ ݆ ‫.)ݐ(‬ This method has been previously shown to reveal information on the community structure of networks and graphs 103 . Once the LEiDA vectors were obtained, we clustered them using K-Means 24,25 with the optimal ‫ܭ‬ (six) determined heuristically with the Calinski-Harabasz and Davies-Bouldin methods (Supplementary Figure S10.
Each cluster represents a set of LEiDAs, and we refer to each of these as a brain state. The dynamics of such states can be studied with three main metrics: fractional occupancy -which refers to the total proportion of time spent in a given state or probability of that state; dwell time -which consists of the average continuous time spent on each state; and Markovian probabilities of transitions between each state 24,25 . In addition, we also calculated values for mean synchronisation and metastability for each state by averaging those for the volumes belonging to each cluster.
Statistics. Firstly, we restricted our sample to the term-born individuals only (n = 324) and evaluated the effect of brain maturation and ex-utero experience (with PMA and PND at scan, respectively) in global and modular brain dynamics. Secondly, we evaluated the effects of prematurity in brain dynamics by studying the entire sample of 390 individuals. Prematurity was coded as a binary variable with 1 for preterm-born individuals (GA at birth less than 37 weeks) and 0 for term-born participants (GA at birth of 37 weeks or more). Finally, we evaluated the association of global and modular dynamic features with later neurodevelopmental outcomes at 18 months corrected age (n = 305 -Bayley-III; n = 300 -Q-CHAT). Modular dynamics: brain states. Firstly, we tested differences between the brain states defined in this study in terms of their mean synchronisation, metastability, fractional occupancy, and dwell times per subject via a type III ANOVA with Satterthwaite's method 104

Statistical significance and repeated measures.
We evaluated the statistical significance of each variable of interest with permutation tests with 10,000 repetitions for all GLMs. Pvalues are reported uncorrected, highlighting those surviving multiple comparison correction across states using Benjamini-Hochberg False Discovery Rate (FDR) method with α error at 5% 105 .
Data and software availability. The fMRI datasets and clinical data analysed in this study are available as part of the dHCP's third data release, which can be obtained from https://data.developingconnectome.org. Pre-processed BOLD timeseries data used in this study are available in https://dx.doi.org/ 10.5281/zenodo.7053984. Dynamic properties of the BOLD signal fluctuations were assessed with dynFC: CoDe-Neuro's Dynamic Functional Connectivity Tools, a set of scripts written in Python v3.7 (https://codeneuro.github.io/dynfc/), and supporting libraries Numpy, SciPy, Scikit-learn, pickle, h5py, pandas, os, sys, and feather.
Brain volume images were produced with BrainNet Viewer and Tools for NIfTI and ANALYZE images. All scripts used in this article's statistics and figures; and relevant instructions on how to run them, are available in https://github.com/CoDe-Neuro/neonatal_dfc.  L  o  n  g  i  t  u  d  i  n  a  l  A  n  a  l  y  s  i  s  o  f  N  e  u  r  a  l  N  e  t  w  o  r  k  D  e  v  e  l  o  p  m  e  n  t  i  n  P  r  e  t  e  r  m  I  n  f  a  n  t  s  .   C  e  r  e  b  r  a  l  C  o  r  t  e  x  2  0  ,  2  8  5  2  -2  8  6  2  (  2  0  1  0 T  h  e  d  e  v  e  l  o  p  i  n  g  H  u  m  a  n  C  o  n  n  e  c  t  o  m  e  P  r  o  j  e  c  t  (  d  H  C  P  )  a  u  t  o  m  a  t  e  d  r  e  s  t  i  n  g  s  t  a  t  e  f  u  n  c  t  i  o  n  a  l  p  r  o  c  e  s  s  i  n  g  f  r  a  m  e  w  o  r  k  f  o  r  n  e  w  b  o  r  n  i  n  f  a  n  t  s  .  N  e  u  r  o  I  m  a  g  e  2  2  3  ,  1  1  7  3  0  3  (  2  0  2  0  )  .   7  6  .  L  e  o  n  a  r  d  i  ,  N  .  &  V  a  n  D  e  V  i  l  l  e  ,  D  .  O  n  s  p  u  r  i  o  u  s  a  n  d  r  e  a  l  f  l  u  c  t  u  a  t  i  o  n  s  o  f  d  y  n  a  m  i  c  f  u  n  c  t  i  o  n  a T  h  e  l  i  n  k  s  b  e  t  w  e  e  n  p  r  e  n  a  t  a  l  s  t  r  e  s  s  a  n  d  o  f  f  s  p  r  i  n  g  d  e  v  e  l  o  p  m  e  n  t  a  n  d   p  s  y  c  h  o  p  a  t  h  o  l  o  g  y  :  d  i  s  e  n  t  a  n  g  l  i  n  g  e  n  v  i  r  o  n  m  e  n  t  a  l  a  n  d  i  n  h  e  r  i  t  e  d  i  n  f  l  u  e  n  c  e  s  .  P  s  y  c  h  o  l  .  M  e  d  .  4  0  ,   3  3  5  -3  4  5  (  2  0  1 T  z  o  u  r  i  o  -M  a  z  o  y  e  r  ,  N  .  e  t  a  l  .  A  u  t  o  m  a  t  e  d  A  n  a  t  o  m  i  c  a  l  L  a  b  e  l  i  n  g  o  f  A  c  t  i  v  a  t  i  o  n  s  i  n  S  P  M  U  s  i  n  g  a   M  a  c  r  o  s  c  o  p  i  c  A  n  a  t  o  m  i  c  a  l  P  a  r  c  e  l  l  a  t  i  o  n  o  f  t  h  e  M  N  I  M  R  I  S  i  n  g  l  e  -S  u  b  j  e  c  t  B  r  a  i  n  .  N  e  u  r  o  I  m  a  g  e  1  5  ,   2  7  3  -2  8  9  (  2  0  0