Interaction of the salience network, ventral attention network, dorsal attention network and default mode network in neonates and early development of the bottom-up attention system

The salience network (SN), ventral attention network (VAN), dorsal attention network (DAN) and default mode network (DMN) have shown significant interactions and overlapping functions in bottom-up and top-down mechanisms of attention. In the present study, we tested if the SN, VAN, DAN and DMN connectivity can infer the gestational age (GA) at birth in a study group of 88 healthy neonates, scanned at 40 weeks of post-menstrual age, and with GA at birth ranging from 28 to 40 weeks. We also ascertained whether the connectivity within each of the SN, VAN, DAN and DMN was able to infer the average functional connectivity of the others. The ability to infer GA at birth or another network's connectivity was evaluated using a multivariate data-driven framework. The VAN, DAN and the DMN inferred the GA at birth (p < 0.05). The SN, DMN and VAN were able to infer the average connectivity of the other networks (p < 0.05). Mediation analysis between VAN’s and DAN’s inference on GA at birth found reciprocal transmittance of change with GA at birth of VAN’s and DAN’s connectivity (p < 0.05). Our findings suggest that the VAN has a prominent role in bottom-up salience detection in early infancy and that the role of the VAN and the SN may overlap in the bottom-up control of attention.


Introduction
To survive the changes and challenges of the external world, we need the ability to focus on the multiple sources of stimuli that constantly compete for our attention, the so-called "saliency detection". Saliency detection is a complex mechanism that requires a bottom-up mechanism for filtering stimuli standing out from a stream of sensory inputs and a higher order mechanism for the automatic attraction and consequent maintenance of attention on a specific task (Parr and Friston 2019). The salience network (SN) plays a dominant role in the detection of salient stimuli across multiple modalities (Crottaz-Herbette and Menon 2006). The role of reality filter is specifically attributed to the insula, which together with the anterior cingulate cortex, constitutes the main anatomical structures to which the SN is anchored (Menon and Uddin 2010). Saliency detection is also involved in the dynamic interaction of sensory and cognitive influences that control attention (Menon and Uddin 2010). Attention is controlled by two partially segregated networks: the ventral attention network (VAN) and dorsal attention network (DAN), forming a twofold attentional control system (Vossel et al.2014;Sridharan et al. 2007). The VAN includes the temporoparietal junction (TPJ) and ventral frontal cortex (VFC), is dominant in the right hemisphere and is generally activated when an unexpected event occurs and breaks one's attention from the current task (bottom-up processing). The DAN includes the intraparietal sulcus (IPS) and the frontal eye fields (FEF) of each hemisphere and shows sustained activation when focusing attention on an object and is thought to be responsible for goal-directed, top-down processing of attention (Corbetta and Shulman 2002). The complex interaction of bottom-up and top-down attention mechanism also requires the ability to disengage from the constant stream of our self-referential thoughts (the so-called "mind wondering"), a cognitive process attributed to the DMN (Raichle 2015;Buckner et al. 2008). This postulate is supported by the evidence that performance of a goal-directed, non-selfreferential task is accompanied by a decrease in activity in the DMN and a corresponding increase in activity in the DAN (Anderson et al. 2011;Esposito et al. 2018).
A number of fMRI studies have demonstrated that infants and even preterm newborns possess immature forms of many of the networks described in the adult (Doria et al. 2010;Smyser et al. 2010;Toulmin et al. 2015;Stoecklein et al. 2020). These studies show that higher order networks may be present, even if in a fragmented, immature form, even before term, as opposed to the previous concept by which the networks develop in parallel with the cognitive competences associated with stimulus-dependent thought.
Recent studies adopting machine learning methods, have shown that the functional connectivity data (expressed by different metrics) and volume data are able to predict prematurity both in cohorts of newborns (Smyser et al. 2016;Ball et al. 2016;Chiarelli et al. 2021) and adults (Shang et al. 2019) and to classify age in infants aged between 6 and 12 months (Pruett et al. 2015). These studies support the view that a functional architecture of the brain exists before birth and constantly evolves, especially during the first year of life. Changes in functional connectivity of prematurity also manifest in adults, which supports the hypothesis that neurocognitive disorders associated with preterm birth might represent a disease of brain connectivity (Lubsen et al. 2011).
Most studies based on the functional connectivity of newborns and infants have been based on a global analysis of brain networks. Few studies, however, have so far attempted to focus on the functional connectivity of highorder networks, whose structure and function is assumed to be immature in early life. The present study seeks to analyze the development of functional connectivity of the networks involved in the complex interaction of bottom-up and topdown attention systems and how they influence the functional connectivity of each other. We based our assumption on the evidence that the development of the insula begins early in the fetal period. This leads us to hypothesize that the salience system starts developing very early, probably interacting with the immature forms of the other networks that combine in the complex mechanism of attention.
First, we explored the association with GA at birth and the resting-state functional connectivity (rs-FC) extracted from the SN, the VAN, the DAN and the DMN.
Consequently, we examined the ability of each network to infer the mean FC of each other. Given the observation that the SN and the VAN have a different, but complementary activity of salience detection (Farrant and Uddin, 2015), we hypothesize that the SN and the VAN influence each other. Moreover, both the SN and the VAN may influence the DMN by decreasing its activity and functional connections.
Although the dorsal and ventral networks are specialized for distinct attentional subprocesses such us top-down controlled attentional selection and the bottom-up selection of relevant stimuli, it seems that the two networks work in concert to promote a specific attentional process and that top-down and bottom-up processing cannot uniquely be attributed to one system in isolation. The ventral and dorsal attention networks may, therefore, develop in a constant dynamic interaction, subsequently evolves into flexible attentional control (Shulman et al. 2003(Shulman et al. , 2007DiQuattro and Geng 2011;Vossel et al. 2014). This assumption led us to hypothesize that the ventral and dorsal attention networks affect the development of each other. To test this hypothesis, we employed a mediation analysis to assess whether the VAN acts as a mediator on the inference on GA at birth using DAN. A specular analysis was conducted for DAN acting as mediator on the inference on GA at birth using VAN. Our study was based on the same dataset of the recent study from Chiarelli et al. (2021) and adopted the same Machine Learning multivariate data-driven framework which allows to consider all the connections within the predicting network at once without any a-priori assumptions. However, we focused on pre-selected higher order networks rather than on exploring the effect on prematurity of different functional connectivity metrics extracted from a set of regions of interest (ROIs) covering the whole brain (Shi et al. 2011). We also further explored the ability of each network to infer the connectivity of each other.

Participants
The present study included a total of 88 healthy neonates with GA at birth from 28 to 40 weeks (mean = 33 weeks, SD = 3.7 weeks). 43/88 patients were female and 15/88 were born at term (> 37 weeks of GA at birth) (Table 1). Informed consent was obtained from the parents of all participants, and the experimental protocols were approved by the Institutional Review Committee. The neonates underwent standard clinical MRI examination within the 40th week of post-menstrual age (PMA, mean = 40.4; SD = 0.27).
Neonates, born before 37 weeks of GA, were selected based on the following exclusion criteria: The inclusion of premature neonates in the study aims at testing if the FC of the networks depends on the degree of maturation that the brain reaches before birth.
Neonates, born within or after 37 weeks of GA, were selected from a group of consecutive neonates without asphyxia. The neonates did not present signal abnormalities at standard MR sequences and had normal neurologic status at a 12-month clinical follow-up.

Data acquisition
MR imaging was performed with a 3 T whole-body system (Achieva 3.0 T X-Series) from Philips Healthcare (Best, Netherlands) using an eight-channel head-only receiver coil. In accordance with several previous fMRI studies in neonates (Ball et al. 2016;Stoecklein et al. 2020), participants were fed and then sedated with 0.05 mg oral Midazolam per kilogram of body weight immediately prior to clinical scans to minimize motion artefacts. Neonates were laid in the scanner in a supine position and swaddled in blankets. Molded foam was placed around the body to minimize head movement. Hearing protection was used through commercially available neonatal earmuffs (MiniMuffs; Natus Medical, San Carlos, California) and adapted ear-canal plugs. Heart rate and oxygen saturation were monitored during the MR imaging session by an intensive care neonatologist. Structural images used in this study were collected using the T1-weighted Turbo Field Echo (TFE) sagittal sequence (flip angle: 8°; TR: 9 ms; TE: 4.2 ms; voxel size: 1 × 1 × 1 mm 3 ; FOV: 200 × 200 × 150 mm 3 ) with a whole-body SAR below 0.2 W/kg. At the end of standard clinical MRI sequences, whole-brain functional images were additionally acquired using a T2*-weighted, echo-planar imaging (EPI), FFE axial sequence (flip angle: 90°; TR: 1555 ms; TE: 30 ms; voxel size: 2.5 × 2.5 × 3 mm 3 ; FOV: 180 × 180 × 75 mm 3 ; slice gap: 0 mm) with a whole-body SAR within 0.8 W/kg. The functional scan duration was 4 min and 15 s (255 s).

Preprocessing
The MR image pre-processing workflow is reported in Fig. 1.
BOLD connectivity was evaluated among 90 subcortical and cortical ROIs defined by the University of North Carolina (UNC) Infant Atlas (Shi et al. 2011). ROI names and associated labels are reported in the Supplementary Material.
To facilitate registration with the UNC atlas, an intermediate in-house brain template, built by averaging the infants T1-weighted anatomical images (Avants et al. 2010(Avants et al. , 2011 and by segmenting the brain of the average template by hand, was first registered to the atlas (Wang et al. 2005). The T1 in house template was constructed by using the Advanced Normalization Tools (ANTs, http:// stnava. github. io/ ANTs/) with default settings (Avants et al. 2010(Avants et al. , 2011. After warping each subject's anatomical image to the in-house template, inverse transformations into the in-house template and into each structural image were applied on the UNC Infant Atlas to successfully identify the ROIs in the original subject anatomical space. EPI T2*-weighted BOLD images, acquired at rest, were pre-processed according to a standard pipeline (Dolgin 2010), using AFNI (Cox 1996) and FSL (Jenkinson et al. 2012). The pipeline included: (i) slice time and motion correction using the 3dTshift and 3dvolreg functions; (ii) marking of motion outliers with the fsl_motion_outliers tool using DVARS metric and default setting for the definition of outliers (Jenkinson et al. 2012;Pruett et al. 2015); (iii) 4D image scaling using fslmaths; (iv) linear and quadratic temporal detrending using 3dDetrer (Churchill et al. 2012). The motion parameters (three translations, three rotations and motion outliers) were finally regressed out (without scrubbing) from the raw BOLD time series, employing the tool  (48) 15 (60) 0 (0) Birth weight at birth (g), mean (SD) 1821 (693) r = 0.88; p < 10 -3 1460 (328) 1938 (482) 3210 (423) APGAR score at birth, mean (SD) 6.4 (2.0) r = 0.14; p = n.s 6.3 (2.1) 7.4 (1.4) 5.3 (2.3) Number of fMRI volume deemed as outliers 7 (5) r = 0.04; p = n.s 7 (5) 6 (4) 7 (5) 3dREMLfit (Bright et al. 2017). Finally, the pre-processed BOLD images were co-registered into anatomical images and inverse transformations were applied to anatomical images to extract the average BOLD signals (expressed as relative signal change with time) from each ROI identified in the native BOLD spaces. Of note, all registrations were performed using the diffeomorphic registration method and the mutual information metric from ANTs (Avants et al. 2010(Avants et al. , 2011. Registered images were visually inspected by an expert neuroradiologist. Resting state functional connectivity (rsFC) matrices were built by evaluating pairwise associations of BOLD signals in the 90 ROIs considered also accounting for a global signal contribution. In particular, each normalized (z-scored) ROI's BOLD timecourse was regressed on each other normalized ROI's BOLD timecourse using the normalized average (among the 90 ROIs) BOLD signal as an additional independent variable within a general linear model (GLM) framework (Murphy and Fox 2017). To evaluate undirected connections for further analysis, the average between the lower and the upper diagonal portions of the connectivity matrices was computed.

ROIs selection within the networks of interest
We selected the ROIs for a given network based on previous research studies and SN-DMN-VAN-DAN were extracted from the compete rsFC matrix. For the SN, we selected the insula, anterior cingulate gyrus, amygdala and thalamus, both left and right for each region (Menon 2011;Menon and Uddin 2010). For the DMN, we selected the medial orbital cortex (as medial prefrontal cortex), the posterior cingulate gyrus, the pre-cuneus and the angular gyrus, both left and right for each region (Buckner and DiNicola 2019). For the VAN, we selected the medial frontal gyrus, inferior frontal gyrus, inferior parietal lobule, and superior temporal gyrus bilaterally. We selected these regions available on the neonatal brain atlas in use to include the ventral frontal cortex (VFC) and the temporo-parietal junction (TPJ), which are both recognized as nodes of the VAN (Corbetta and Shulman 2002;Vossel et al. 2014). For the DAN, we selected the supplementary motor area and the superior parietal gyrus bilaterally. We selected these regions available on the neonatal brain atlas in use to include the frontal eye filed (FEF) within the SMA and the intraparietal sulcus (IPS), which are both recognized as nodes of the DAN (Corbetta and Shulman 2002;Vossel et al. 2014).

Inference of gestational age (GA) at birth
A multivariate analysis framework was implemented to regress GA at birth on rsFC within different networks (Fig. 2). The number of independent connections were 28 for SN, DMN and VAN and 6 for DAN. To account for collinearity among ROI connections (Huopaniemi et al. 2009) a partial least square (PLS) regression was used (Wold et al. 2006), which reduces the predictors to a smaller set of uncorrelated components maximally associated with the dependent feature/s (Abdi et al. 2013). Of note, the learning process (fitting) of the PLS algorithm provides regression loadings that can be used to retrieve the weights (β-weights) of the original independent variables (Bishop 2006). In order to optimize the hyperparameter of the PLS (number of uncorrelated components) and to concurrently evaluate the out-oftraining-sample performance of the algorithm (generalization) a tenfold nested cross validation (nCV) was employed (Filzmoser et al. 2009). During the tenfold nCV, the number of components allowed during the hyperparameter optimization were constrained between a minimum of 1 and a maximum of 6. The expected β-weights of the PLS were finally computed by running a single analysis on the complete dataset using the rounded average number of components (i.e. the optimal number) delivered by the inner loops of the tenfold nCV analysis. β-Weights were transformed into z scores by dividing them by their expected standard deviation, which was computed relying on repeating the analysis on random shuffled outputs 10 6 times. As a control for possible effects of motion on the functional connectivity and activity results, the multivariate nCV regression was also performed on the variance of motion signals.

Inference of mean connectivity
The same multivariate analysis presented in the previous paragraph was implemented to test if each network could infer the average connectivity of the other three networks.

Mediation analysis
In the last set of statistical analysis, we conducted a mediation analysis (Mackinnon 2012) that was performed in order to determine whether the cross-validated inference on GA at birth using VAN acted as a mediator on the cross-validated inference on GA at birth using DAN. The rationale for the mediation analysis was to examine if the development of the VAN, assumed from the evidence of cross-validate inference on GA, can mediate the development of the DAN, assumed from the evidence of cross-validated inference on GA.
A second mediation analysis was performed in order to determine whether the cross-validated inference on GA at Fig. 2 Multivariate PLS analysis implemented to infer GA at birth or other networks average connectivity. The optimal number of PLS components, the multivariate β-weights and the inference perfor-mances were estimated through a tenfold nested cross-validation scheme. Connectivity matrices reported are from an exemplar subject birth using DAN acted as a mediator on the cross-validated inference on GA at Birth using VAN.
The mediation analysis was performed using the Sobel test (Preacher and Leonardelli 2010).

Data quality and motion artifacts
Thanks to sedation, motion artifacts were almost absent, therefore, only a small average number of 7 volumes per subject (SD of 5 volumes) was deemed as outliers by the algorithm (refer to Table 1) (Ciric et al. 2017). The number of motion outliers and the variance of the 6 motion DVARS signals showed no significant correlation with GA at birth (all r's < 0.1, all p's n.s.). Figure 3 reports the results of the multivariate framework in inferring GA at birth relying on connectivities within different networks.

Inference of gestational age at birth
The SN was not able to infer GA at birth using the multivariate framework implemented (p > 0.05).
The DMN was able to significantly infer GA at birth (r = 0.26; df = 86, p = 0.01, 2 PLS components). By looking at the statistical relevance (z-score) of the β-weights, the strongest positive effects on the inference of GA at birth (larger connectivities between DMN ROIs associated with older GA at Birth) was found for the connectivity of the right medial prefrontal cortex-right pre-cuneus (ORBmed-R-PCUN-R), left medial prefrontal cortex-right pre-cuneus (ORBmed-L-PCUN-L), right posterior cingulate-right angular gyrus (PCG-R-ANG-R), right posterior cingulateleft angular gyrus (PCG-R-ANG-L), left posterior cingulate-right angular gyrus (PCG-L-ANG-R).
The VAN could infer GA at birth (r = 0.29; df = 86, p = 7•10 -2 ) with an optimal number of PLS components of 5. The strongest effects on the GA predictability, were negative effects (lower connectivities between VAN ROIs associated with older GA at birth) found for the left medial frontal gyrus-inferior parietal lobule, right medial frontal gyrus-left inferior parietal lobule, right inferior frontal gyrus and left superior temporal gyrus.
The DAN could infer GA at birth (r = 0.25; df = 86, p = 0.017) with an optimal number of PLS components of 2. The strongest effects on the GA predictability, were found for the left and right superior parietal lobules.

Salience network
Results are displayed in Fig. 4a.
The SN was found to significantly infer the average connectivity of the VAN (r = 0.24; df = 86, p = 0.026, 2 PLS components). The strongest effect on the predictability of the VAN's mean connectivity, as shown on the beta-weight matrix, was found for the left insula-right insula connectivity. A negative effect was found for the right insula-right thalamus connectivity.
The SN also inferred the average connectivity of the DAN (r = 0.29; df = 86, p = 6•10 -3 , 6 PLS components). The strongest positive effect on the predictability of the DAN's mean connectivity, as shown on the beta-weight matrix, was found for the right anterior cingulum-right thalamus connectivity; the strongest negative effect was found for the left anterior cingulum-left thalamus connectivity).
All the regions of the SN did not predict the average connectivity of the DMN (p > 0.05).

Default mode network
Results are displayed in Fig. 4b.
The DMN was found to significantly infer the average connectivity of the SN (r = 0.29; df = 86, p = 7•10 -3 , 6 PLS components). The strongest effects on the predictability of the SN's mean connectivity, as shown on the beta-weight matrix, were negative effects found in the left angular gyrus-right pre-cuneus and right angular gyrus-left pre-cuneus.
DMN connections did not correctly infer the average connectivity of the VAN (p > 0.05).
The DMN was found to strongly infer the average connectivity of the DAN (r = 0.48; p < 10 -3 , 4 PLS components). The strongest positive effects on the predictability of the DAN's mean connectivity were found for the left and right angular gyri. The strongest negative effects were found for the left posterior cingulate gyrus-left prefrontal cortex and in the right posterior cingulate gyrus-right pre-cuneus.

Ventral attention network
Results are displayed in Fig. 5a.
The VAN was found to significantly infer the average connectivity of the SN (r = 0.34; df = 86, p = 10 -3 , 1 PLS component). The strongest positive effects on the predictability of the SN's mean connectivity were found for the connectivity of the medial frontal gyrus with the inferior parietal lobule of the contralateral side, bilaterally.
The VAN was found to significantly infer the average connectivity of the DMN (r = 0.22; df = 86, p = 0.044, 2 PLS components). The strongest positive effect on the predictability of the DMN's mean connectivity was found for the connectivity of the right inferior frontal gyrus with the right inferior parietal lobule.
The VAN strongly inferred the average connectivity of the DAN (r = 0.47; p = 4.8•10 -6 , 1 PLS component). The strongest positive effect on the predictability of the DAN's Fig. 4 Results of the multivariate framework to other network based on connectivities within different networks mean connectivity was found for the connectivity of the left and right inferior parietal lobule. The strongest negative effect was found for the connectivity of the left inferior frontal gyrus and left inferior parietal lobule.

Dorsal attention network
Results are displayed in Fig. 5b.
The DAN did not predict the average connectivity of the SN (p > 0.05) of the DMN (p > 0.05) and of the VAN (p > 0.05). Figure 6 reports the outcome of the mediation analyses.

Mediation analysis
The mediation analysis between VAN's and DAN's inference on GA at birth found that VAN's connectivity significantly modify the transmittance of change of the DAN's connectivity assumed form evidence of crossvalidated inference on GA at birth (VAN mediates DAN: Sobel Test t = 1.969, p = 0.048).
The mediation analysis between DAN's and VAN's inference on GA at birth found that DAN's connectivity significantly modify the transmittance of change of the VAN's connectivity assumed form evidence of crossvalidated inference on GA at birth (DAN mediates VAN: Sobel Test t = 1.910. p = 0.049).

Summary of findings
SN connectivity was found unable to infer the GA at birth. The VAN, DAN and the DMN inferred the GA at birth. The SN was able to infer the average connectivity of the VAN and DAN, but not that of the DMN. The DMN was able to infer the average connectivity of the SN and of the DAN, but not that of the VAN. The VAN was able to infer the average connectivity of all the other networks with the highest statistical significance in predicting the connectivity of the SN and of the DAN. The DAN was not able to predict the average connectivity of all the other networks. The mediation analysis showed that VAN's connectivity significantly modify the transmittance of change of the DAN's connectivity, as assumed from evidence of crossvalidated inference on GA. Mediation analysis also showed that DAN's connectivity significantly modify the transmittance of change of the VAN's connectivity, assumed from evidence of cross-validated inference on GA.

Inference of GA at birth
SN connectivity was found unable to infer the GA at birth. The SN is a large scale network of the heteromodal cortex involved in detecting and filtering salient stimuli, thus playing the role of "reality filter". A major function of the anterior insula (AI) node of the SN is the detection of behaviorally relevant stimuli (Crottaz-Herbette and Menon 2006;Eckert et al. 2009;Sridharan et al. 2008;Sterzer and Kleinschmidt 2010), while the anterior cingulate cortex (ACC) send strong motor output. The insula is the first cortical structure to differentiate beginning 6 weeks after conception and going under a complex process of maturation between the 13th and the 28th week of gestation, thus providing the structural basis for its hub role even before term. These findings may suggest that because the insula development starts early in fetal life, its connectivity, under the influence of early external stimulation, matures faster and earlier than other networks (Afif et al. 2007). The connectivity of the SN after birth might, therefore, be found at a stage that is too late to infer the GA at birth. In other words, the functional connectivity of the insula after birth might be too mature to variate accordingly with the gestational age.
The VAN and the DAN were found to significantly infer the GA at birth. The VAN and the DAN are known to play a complementary role in the complex mechanism of attention. Attention is a complex function that requires top-down sensitivity control and a bottom-up mechanism for filtering stimuli (Parr and Friston 2017). The bottom-up mechanism for filtering stimuli involved in the detection is attributed to the VAN (Corbetta and Shulman 2002;Vossel et al. 2014). The top-down mechanism (a process that regulates the relative signal strengths of the different information channels that compete for access to working memory) (Egeth and Yantis 1997) is attributed to the DAN. The VAN is composed of the temporoparietal junction (TPJ) and the ventral frontal cortex (VFC) and is thought to be lateralized to the right hemisphere of the brain (Corbetta and Shulman 2002;Corbetta et al. 2008). Our study found that the strongest effects of VAN predictability on the GA were negative effects (anti-correlation), which were found for the left medial frontal gyrus-inferior parietal lobule. A possible explanation might be that the VAN goes under early lateralization during infancy, which may account for the anti-correlation found in the left hemisphere. The process of lateralization might be still incomplete in early life (Vossel et al. 2014) and might, therefore, be able to reflect the degree of brain development. The DAN is a bilateral network comprising the intraparietal sulcus (IPS) and the frontal eye fields (FEF) of each hemisphere involved in top-down voluntary allocation of goal-driven attention. We found that the strongest effects of the DAN predictability on the GA were the left and right superior parietal lobules, this finding is coherent with the bilateralism of the network. The predictability of both the VAN and the DAN on the GA, supports the evidence by which the VAN operates with the dorsal attention network (DAN) forming a twofold largely interconnected attentional control system (Corbetta and Shulman 2002;Corbetta et al. 2008;Vossel et al. 2014).
The DMN was found to significantly infer the GA at birth. The DMN, a large-scale brain network primarily composed of the medial prefrontal cortex, posterior cingulate cortex/pre-cuneus and angular gyrus (Raichle 2015;Buckner et al. 2008), identified with the stream of self-referential thoughts, the so-called "resting state" (Andrews-Hanna et al. 2010;Raichle 2015). The mental state of stimulus-independent thoughts counteracts attention. This finding is supported by the evidence that the DMN is anticorrelated with the DAN (Fair et al. 2007). We found that the strongest effects of the DMN on the GA predictability was for the connectivity of the right medial prefrontal cortex-right pre-cuneus. These findings are in line with previous studies showing that the DMN is present, even if in a fragmented, immature form, even before term (Doria et al. 2010). More specifically, the cortical hubs in infants seem to be less present in the prefrontal cortex and in the pre-cuneus than in adults. The connectivity of these cortical hubs might therefore reflect the GA.

Machine learning and prediction of brain maturity
A number of studies have applied machine learning to predict brain maturity in pre-term infants and infants born at term. In a recent study by Chiarelli et al. (2021), the authors examined the effect of prematurity on measures of rs-FC, resting state functional connectivity nodal strength (rs-FCNS), fractional amplitude of low frequency fluctuations (fALFF) and regional volume in 90 ROIs, covering the whole brain, by performing region-based univariate analysis of each metric, to explore the association with GA at birth and the spatial consistency across metrics. The study confirms that prematurity is associated with a complex pattern of bidirectional alterations of regional functional connectivity brain volume and, and to a lesser extent, with modifications of fALFF. For rsFC, the results did not suggest strong focal effects and no ROI was significantly correlated with GA at birth after multiple comparison correction. We analyzed the same dataset of the recent study from Chiarelli et al. (2021), however, the goal of our study was different. Although we adopted the same Machine Learning multivariate data-driven framework using partial least square regression (PLS), which allows us to consider all the connections within the predicting network at once, as opposed to Chiarelli's study, we based our assumptions on an a-priori hypothesis based on the interaction of pre-selected networks. Moreover, Chiarelli's analyses were selectively conducted on positive rsFC correlations whereas in our study VAN's strongest prediction on GA were negative effects. Smyser et al. (2016) have applied a SVM-multivariate pattern analysis (MVPA) classification method and SVR (SVM regression) to infants' rs-fMRI data and developed a model to estimate an infant's GA at birth based upon rs-fMRI data collected at term equivalent PMA. Performing RSN-specific analyses revealed that higher order RSNs, such as DMN, cingulo-opercular, DAN and SN, contributed to successful categorization of pre-term infants and infants born at term. These findings are partially in line with our results since, in our study, the SN was not able to predict the GA. A possible explanation might be study subjects included very pre-term infants (23-28 weeks). We hypothesize that, since before the 28th week, the maturation of the insula is not completed, it might be capable of predicting brain maturity.
Ball et al. have combined high dimensional independent component analysis (ICA), representing the FC data into functional nodes and networks the with machine learning techniques to test the hypothesis that preterm birth results in specific alterations to FC. Although Ball did not divide the edges into networks, it is possible to compare Ball's edges to the statistical relevance of beta-weights on the inference of GA at birth. Interestingly, Ball's discriminatory edges comprehend the connections between the insula and anterior cingular region. However, Ball's edges do not contain the connection between the right prefrontal cortex and the pre-cuneus (left and right), the connection between the left medial prefrontal gyrus and left inferior parietal lobule and the connection between the left and right superior parietal lobules, representing the connections with the strongest effect on the inference on GA at birth, respectively, for the DMN, the VAN and DAN in our study. One possible explanation might be the inclusion of very pre-term infants (24-28 weeks) in Ball's cohort of patients. A second explanation might be that, using multiple iterations to select the discriminative edges in the network model, may induce the selection of isolated connections, which, within the architecture of a network, may instead show different effects on the inference of GA at birth. Representing brain connections through networks' architecture might better reflect the actual brain's connectivity structure (Rosenthal et al. 2018).
By adopting the SVM method, Pruett et al. (2015) were able to classify, above chance, 6-versus 12-month-old infants on FC data. The selection of functionally defined seeds enabled the identification of the cingulo-opercular, dorsal attention, fronto-parietal and salience network for 6-vs. 12-month old classification, thus proving that higher order networks are able to discriminate infants age during the first year of life. The SN might be able to discriminate between 6 and 12 month following the rapid development of the other structure other than the insula composing the SN (Gilmore et al. 2018).
Shang et al. adopted automated machine learning pipelines on gray matter volume and amplitude of low-frequency fluctuation (ALIFF) data to characterize multivariate brain patterns that separated full-term born from very preterm born adult subjects. Interestingly, the decrease of ALIFF in the right insula was able to discriminate between adults born at-term and adults born very pre-term, suggesting that, for very preterm infants, the insula might not have completed its process of maturation. The study highlights how the alteration of functional connectivity in preterm infants may persist through adulthood (Schafer et al. 2009).

Inference of average connectivity
The SN connections were found to significantly infer the average connectivity of the VAN (p = 0.026). Interestingly, also the VAN was found to strongly infer the mean functional connectivity of the SN (p = 10 -3 ). These findings seem to be in agreement with previous studies showing that the function of the VAN overlaps with that of the SN in salience detection (Dosenbach et al. 2008;Uddin 2015;Farrant and Uddin 2015). Uddin and Farrant proposed that the VAN and the SN have overlapping nodes in the region surrounding the VFC and anterior insula. They also speculated that in children, these two networks may be less segregated than in adults, and that bottom-up salience processes and attention to environmental stimuli may be over-represented in the child's brain (Farrant and Uddin 2015). The SN connections inferred the average connectivity of the DAN (p = 6•10 -3 ). A possible explanation to this finding is that the overlapping nodes between SN and VAN (Farrant and Uddin 2015) are part of the twofold, interconnected, attentional system composed by VAN and DAN, therefore the SN may infer DAN's connectivity.
The SN did not infer the average connectivity of the DMN (p > 0.05). The DMN, on the contrary, was found to significantly infer the average connectivity of the SN (p = 7•10 -3 ). This finding seems to oppose to the role of the SN in modulating the background activity of the DMN to elicit detection of salient stimuli and facilitate goal-directed behavior. However, since the function of the CEN is likely reduced in infancy (Gao et al. 2015), this finding might also indicate that the SN function in switching attention to goal-directed behavior matures with the development of the executive functions (Teffer and Semendeferi 2012). The VAN was found to significantly infer the average connectivity of the DMN (p = 0.044). This finding may suggest an interaction of the two networks during the development of cognitive and attention skills, unexplored by previous studies, with a functional interplay between the bottom-up attention system of the VAN and internally oriented processing attributed to the DMN (Kim 2010).
The VAN was able to strongly infer the average connectivity of the DAN (p = 4.8•10 -6 ). On the contrary, the DAN did not infer the average connectivity of any of the other networks, suggesting that it might be too early for the DAN to properly interact with other networks. A number of studies have shown that VAN and DAN are in competition during visual and verbal tasks (Anticevic et al. 2010;Todd et al. 2015;Matsuyoshi et al. 2010;Majerus et al. 2012) and that the higher the load of the task, the higher the activation of the DAN and the higher the deactivation of the VAN. Nevertheless, there is consistent evidence showing that VAN and DAN are specialized networks that work together, promoting a flexible attentional system in which top-down and bottom-up processing cannot be uniquely attributed to one single system in isolation (Shulman et al. 2003(Shulman et al. , 2007DiQuattro and Geng 2011;Vossel et al. 2014). Our findings lead us to hypothesize a conjoint development of the VAN and DAN with a prominent role of the VAN and, therefore, the bottom-up attention system during fetal life and early infancy. The two networks form an immature twofold attention system that matures with the growing ability of the topdown attention system to counteract the bottom-up system. Our hypothesis is further supported by a recent study by Suo et al. (2021), underpinning the stable and reliable anatomical connections of the DAN and the VAN via functional connectors, demonstrating that the functional interaction of the VAN and the DAN is supported by a solid anatomical structure.

Mediation analysis
The mediation analysis of the VAN's inference over the DAN, found that VAN's connectivity significantly modify the transmittance of change of the DAN's connectivity with GA at birth, assumed from evidence of cross-validated inference on GA at birth (VAN mediates DAN: Sobel Test t = 2.14, p = 0.032). A specular analysis of the DAN's inference over the VAN found that DAN's connectivity significantly modify the transmittance of change on the VAN's connectivity with GA at birth, assumed from evidence of cross-validated inference on GA at birth (DAN mediates VAN: Sobel Test t = 1.969, p = 0.048). This finding further supports the theory of a twofold, largely interconnected attentional control system (Corbetta and Shulman 2002;Corbetta et al. 2008;Vossel et al. 2014) and suggests that the bottom-up and top-down attention systems affect the development of each other. We hypothesize that maturation is crucial for the development of attention networks and prematurity may, therefore, substantially affect their architecture and function (Posner and Rothbart 2012;Nie et al. 2013;Wolf and Pfeiffer 2014;Wen et al. 2019).

Strengths and limitations
This is the first study performed on a cohort of neonates to show that the VAN may influence the connectivity of the SN, the DMN and the DAN. Interestingly, the SN influenced the VAN and the DAN, but not the DMN. We believe our findings suggest a prominent role of the VAN in the bottomup salience detection in early infancy and that the VAN and the SN may overlap in their roles of bottom-up control of attention. Although a small number of studies have adopted machine learning methods to assess if FC metrics were able to predict the GA at birth, this is the first study to focus on a small group of networks of higher order. The evidence that the SN is not able to predict GA at birth, as opposed to the other higher order networks we analyzed, is a novel finding that we thought to be the result of the early development of the insula, which goes under a complex process of maturation between the 13th and the 28th week of gestation (Afif et al. 2007). Such an early pattern of maturation would be less impaired by the effect of prematurity.
Our studies has some limitations. First, the atlas we used to select the seed regions for each network only identifies coarse regions. However, this limitation is provided by the absence of efficient brain anatomy segmentation tool in newborns. Second, the number of subjects was not very large considering that we studied four different networks using a multivariate analysis framework. A short BOLD acquisition time (4 min), driven by the limited time available for conducting a relatively non-clinical evaluation in a clinical environment. To reduce motion artifacts and maximize the efficacy of the standard clinical evaluation, newborns were also mildly sedated using Midazolam. Although Midazolam might have altered brain activity and hemodynamics, an effect of Midazolam should be observed in all subjects and should, therefore, not affect regression analyses that exploit subject-specific alterations. Another limitation of the present retrospective study is the minimally essential clinical information available. A limitation of the statistical analysis is that the results of the mediation analysis are uncorrected. Although infants with evident alterations at standard radiological assessment were excluded from the analysis and no relationship was found between the main available clinical variable, the APGAR score soon after birth, and the extent of prematurity, the presence of subtle clinical confounders cannot be definitely ruled out.
Future studies should replicate the present study with a larger cohort of patient with larger sets of clinical data available. Moreover, longer acquisition time should be taken into consideration together with a no sedation approach. Manual segmentations of ROI should be also performed for a more precise delineation of brain regions in neonates.

Conclusions
This is the first study performed with a cohort of neonates with a multivariate data-driven framework (i.e. machine learning framework) to suggest a prominent role of the VAN in the bottom-up salience detection in early infancy and that the VAN and the SN may overlap in their roles of bottom-up control of attention. We also found reciprocal influence of VAN and DAN on the development of each other network.
The SN was the only network in our analysis unable to infer the gestational age at birth.
This finding may indicate that the stage of maturation of the SN at birth might be too advanced to infer the GA.