The analysis of resting-state functional connectivity (RS-FC) constitutes a promising tool as it provides complementary information to structural imaging related to brain physiology. Indeed, since its discovery in 1995 (Biswal et al., 1995) resting state functional MRI (rsfMRI) studies have provided new insights in the understanding of brain architecture and cerebral development (Power et al., 2010; Smyser and Neil, 2015; Gao et al., 2017; Grayson and Fair, 2017; Keunen et al., 2017; Zhang et al., 2019). Smyser et al. (2013) demonstrated the feasibility of using rsfMRI to explore the alterations in resting state networks (RSN) associated with preterm birth and white matter injury. Alterations of the default mode and ventral attention networks at birth, are associated with behavioral inhibition at age of two years, (Sylvester et al., 2018) which suggests early alterations of the RSN present a correlation with clinical manifestations, and opens the opportunity of early diagnostics and treatment. Additionally, neonatal RSN are consistently identifiable and present with high similarities to older populations (Fransson et al., 2007, 2009; Gao et al., 2009). RS-FC is based on low frequency regional fluctuations (<0.1 Hz) in the Blood-Oxygen-Level-Dependent (BOLD) (Ogawa et al., 1990, 1993) signal while the participant is not performing any task, a useful feature when evaluating neonates (Smyser and Neil, 2015). RSN signal is very stable across subjects (Lee et al., 2013), but vulnerable to several artifacts such as head-motion (Maknojia et al., 2019), susceptibility distortions and or white matter (WM) and cerebrospinal fluid (CSF) signals (Jo et al., 2013; Power et al., 2014). Robust rsfMRI data processing is key to reduce the nuisance effects of the non-neural signals in the data to identify reliable resting state activity (Lund et al., 2006; Giove et al., 2009). Its clinical potential and implementation present several methodological challenges that need to be addressed before considering its use to develop a new generation of biomarkers. For this reason, straightforward to use and open-source tools for the neonatal rsfMRI data processing need to be readily available. Tools for mature brains already exist to process rsfMRI data, but analyzing the neonatal brain presents challenges that need to be addressed with new approaches (Smyser and Neil, 2015). There are several straightforward rsfMRI data processing pipelines developed for adults such as Conn toolbox (Whitfield-Gabrieli and Nieto-Castanon, 2012),fmriprep (Esteban et al., 2019), the Human Connectome Pipeline (HCP) (Glasser et al., 2013), the Resting-State Analysis Toolkit (REST) (Rubinov and Sporns, 2010), or the Connectome Computation System (CCS) (Xu et al., 2015), however, those are not adapted to the newborn brain which presents additional challenges, such as different contrast due to myelination (Enguix et al., 2018). T2-weighted images are usually needed for tissue segmentation in place of T1-weighted images. Further, varying brain sizes between subjects (Smyser and Neil, 2015) makes adult skull stripping less robust on the neonatal brain. Additionally, different age specific atlases and tissue probability maps are required for accurate segmentations, common space normalization and seed-based analysis.

To the best of our knowledge the only existing open-access pipeline to process neonatal rsfMRI data is the one developed by the developing Human Connectome Project (dHCP) (Fitzgibbon et al., 2020). While this pipeline has proven to provide excellent results with the dHCP data, its implementation on smaller or clinical datasets remains challenging, as it requires large datasets for independent component (IC) denoising. Furthermore, the dHCP pipeline can be difficult to set up for cohorts acquired at other centers, because the pipeline was developed/optimized from the dHCP database specifically. For example, the dHCP denoising step is based on spatial independent component analysis (sICA), which separates independent correlating signals that can be classified as neural or non-neural signal. This denoising technique has been shown to provide superior results in adults and infants when the dimensionality is accurately set (Griffanti et al., 2017; Alfaro-Almagro et al., 2018). However, the identified signals need to be classified as neural signal or structured noise, which in most cases is performed manually and a difficult process to automate. To overcome this limitation, the dHCP pipeline uses a machine learning approach (ICA-based Xnoiseifier) (Salimi-Khorshidi et al., 2014) to classify the independent components as neural signals or noise. The machine learning algorithm requires a minimum of 35 manually labeled subjects to be trained, which is not always possible in smaller cohorts and requires specialists to manually classify the independent components (Fitzgibbon et al., 2020).

To overcome the aforementioned challenges, we developed NeoRS, with the goal of creating a robust open-source pipeline containing the necessary tools to preprocess rsfMRI data. The main advantages of NeoRS are it has been developed specifically for neonates, simple to implement and flexible to process different datasets. Additionally, it can process single subject data, utilizes parallelizable environment and includes visual quality control checkpoints at each step.

The data processing steps include T2-weighted image alignment to a common space, slice timing correction and segmentation, and rsfMRI procedures are slice timing, distortion correction using reversed phase encoding polarity acquisitions, alignment in a common space, motion correction, removal of nuisance confounds and noise compromising functional data interpretation. Further, simple resting state functional connectivity cross-correlations based on seed-to-voxel and seed-to-seed approaches are incorporated.

NeoRS is an open-source image processing pipeline dedicated to neonatal rsfMRI. It includes seed-based and seed-to-seed 1st level analysis. NeoRS has neonatal brain templates for term in 1mm and 3mm Talairach space, as well as a set of 31 seed regions defining seven common resting state networks. NeoRS relies on the open-source neuroimaging pipelines: SPM, FSL and AFNI and encompasses robust methods to segment, register and denoise neonatal rsfMRI data. Each of the various processing steps were evaluated separately. Output was carefully inspected to ensure the best quality products, by optimizing skull stripping, image registration, head motion and denoising related results. Additionally, functional connectivity of the motor, visual, default mode, language, dorsal attention, ventral attention and fronto-parietal networks using seed-based functional connectivity analysis were demonstrated.

Image registration to the atlas was meticulously inspected for every subject and no significant misalignments were observed for T2-weighted or rsfMRI images for affine registration. We compared single step to two-step registration for accuracy and computational times. The single step registration was chosen, as this process required less computational time and demonstrated no substantial variation when compared to the two-step counterpart.

In contrast with some of the most common adult pipelines, image registration is implemented prior to skull stripping, as alignment quality results were identical. Further, performing image registration preceding skull stripping produced vastly robust skull stripping results across varying brain sizes. This was not the case employing traditional skull stripping before image registration. Skull stripping of the T2-weighted images was found to be a crucial step when working with neonates because poor stripping lead to misclassification of segmented data and ultimately unreliable representation of RSN because of misregistrations. Additionally, image registration prior to skull stripping facilitated brain extraction without any user intervention. Furthermore, if the brain was previously aligned to an atlas, bet2 was implemented on the subject specific atlas aligned data instead of applying an atlas brain mask to avoid subtle geometric inaccuracies. This procedure is critical and needs to be properly assessed. For this reason, an output for the skull stripping is provided which contains the image of the non-skull stripped mask with an overlay of the contour of the skull stripped brain.

After skull stripping, tissue segmentation is performed using Mantis, without the need of any further intervention. Mantis is integrated into the NeoRS pipeline as it provided robust results for disparate brain sizes using only T2-weighted images. Contrary to the adult brain that uses T1-weighted images for brain segmentation, it is fundamental in neonates to provide the pipeline with T2-weighted images as the water/cholesterol ratio is reversed with respect to adults due to lack of myelination. Neonatal T2-weighted images present a better contrast between brain structures (McArdle et al., 1987). While Mantis needs to be installed to use the NeoRS pipeline, no additional setup steps are required as it is fully assimilated in NeoRS. The results showed the binary masks created from Mantis tissue probability maps are perfectly aligned with their corresponding structures for various subjects without any manual intervention.

It is well known the GRE-EPI sequence for rsfMRI is prone to susceptibility-induced magnetic field inhomogeneity (Czervionke et al., 1988). The artifact primarily appears close to the extrema portions of the brain in the phase encoding direction (Andersson et al., 2018) and needs to be properly corrected. While several methods have been successfully used in these settings (Jezzard and Balaban, 1995; Cusack et al., 2003), NeoRS operates the standard topup/applytopup method. Two reversed phase encoding direction images are involved to correct for the deformations (Holland et al., 2010). The method is simple to implement in the acquisition protocol, provides high quality results and acquisition sequences require very short duration times.

Slice timing correction remains a controversial step when a very short repetition time (TR) is deployed (Parker and Razlighi, 2019). However, as it has been shown it can significantly improve z-scores (Parker and Razlighi, 2019) and as the aim was to make NeoRS work with the maximum number of datasets, it is included as an option. Slice timing correction can be deactivated for multi-band sequences with very low TR, as in this kind of low TR sequence, all slices in each volume are acquired very closer together (Glasser et al., 2013).

After data preprocessing, confounding signals and motion effects are removed. First, the framewise displacement threshold is defined as 0.25 mm, as performed by Smyser et al. on their neonatal study and shown to provide accurate results (Sylvester et al., 2018). Volumes with FD higher than 0.25 mm were removed from the time-series as described per Power et al. (2012). Motion censoring was performed prior to filtering to prevent spikes from passing through band-pass filtering, as this could introduce artifacts such as Gibb’s ringing and or skew correlation coefficients. Furthermore, extremely high motion acquisitions shouldn’t be taken into account as they can lead to inflated results. To do so, different metrics can be adopted, such as maximum framewise displacement, minimum number of low motion volumes or average framewise displacement. In NeoRS, acquisitions with average FD higher than 0.25 mm were introducing augmented correlations related to motion and improper denoising. Those acquisitions were completely removed. Therefore, the average framewise displacement was employed as the metric for exclusion.

To correct for nuisance variables NeoRS implements a traditional denoising strategy that performs global signal, white matter and cerebrospinal fluid signal regression, motion parameter regression and a band-pass filter. This simple approach provides robust results (Smyser et al., 2016), and doesn’t require manual intervention or large datasets for denoising purposes. In contrast to the aforementioned independent components denoising techniques, (Griffanti et al., 2017; Alfaro-Almagro et al., 2018) such as the one employed in the dHCP pipeline.

Finally, NeoRS incorporates seed-based functional connectivity analysis tools to assist the user in assessing initial results. Seed-based results across subjects showed patterns very similar to those observed in the literature for all the analyzed resting state networks (Fransson et al., 2009; Damaraju et al., 2010; Figure 11). Further data analysis can be carried out on the fully processed data, final_BOLD.nii, if desired.

NeoRS² is an open-source, straightforward to use rsfMRI data processing pipeline for the neonatal brain which relies on the open-source neuroimaging pipelines FSL, AFNI and SPM. NeoRS works with neuroimaging nifti format, BIDS folder structures and has been developed to work with different MRI vendors and diverse acquisition parameters with minimal user implication. After image processing with NeoRS, we observed resting state networks were in agreement with previously published studies at term age. Each processing step is easy to inspect to ensure consistent results through quality control checkpoint figures.

An open-source, rudimentary to use pipeline for neonatal resting state image processing will allow the community to process their data immediately after scanning sessions implementing a simple computational infrastructure. The democratization of rsfMRI processing will allow a higher number of centers to collaborate and process their datasets and optimistically bring the clinical biomarker application one step closer.

The original contributions presented in this study are included in the article/supplementary material, further inquiries can be directed to the corresponding authors.

VE developed NeoRS, processed and analyzed the data, and wrote the manuscript. JK provided guidance and helped comparing results with a non-open-source pipeline. DL contributed the writing and reviewing the first draft and provided guidance. JC-A provided guidance and reviewed the manuscript. GL supervised the project, contributed writing and reviewing all drafts, and provided guidance. All authors contributed to the manuscript revision, read, and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.