Functional Connectivity-Derived Optimal Gestational-Age Cut Points for Fetal Brain Network Maturity

The architecture of the human connectome changes with brain maturation. Pivotal to understanding these changes is defining developmental periods when transitions in network topology occur. Here, using 110 resting-state functional connectivity data sets from healthy fetuses between 19 and 40 gestational weeks, we estimated optimal gestational-age (GA) cut points for dichotomizing fetuses into ‘young’ and ‘old’ groups based on global network features. We computed the small-world index, normalized clustering and path length, global and local efficiency, and modularity from connectivity matrices comprised 200 regions and their corresponding pairwise connectivity. We modeled the effect of GA at scan on each metric using separate repeated-measures generalized estimating equations. Our modeling strategy involved stratifying fetuses into ‘young’ and ‘old’ based on the scan occurring before or after a selected GA (i.e., 28 to 33). We then used the quasi-likelihood independence criterion statistic to compare model fit between ‘old’ and ‘young’ cohorts and determine optimal cut points for each graph metric. Trends for all metrics, except for global efficiency, decreased with increasing gestational age. Optimal cut points fell within 30–31 weeks for all metrics coinciding with developmental events that include a shift from endogenous neuronal activity to sensory-driven cortical patterns.


Introduction
Recent advances in brain imaging have enabled researchers to investigate in vivo human brain development. Resting-state functional connectivity MRI (RS-fMRI), specifically, has allowed us to noninvasively probe developing functional connections across multiple brain networks simultaneously in fetuses [1,2]. Understanding neural circuitry formation in healthy fetuses is critical in recognizing early signs of functional alterations caused by genetic or environmental insults. The timing of developmental events can help differentiate aberrant from neurotypical processes. Thus, perinatal RS-fMRI aims to detect when functional networks emerge (e.g., appearance of visual and sensorimotor RSNs) and to identify when networks mature (e.g., become more similar to newborn architecture) [3][4][5][6].
While fetal imaging is a nascent field, findings over the past decade have begun to shed light on patterns of neural circuitry formation: connectivity follows a mediolateral, posteroanterior timing [4,7,8], occipital and sensorimotor circuits develop earlier than parietal connections [4], lateralization appears in superior temporal cortices in utero [3], and intraand interhemispheric connectivity strength increase with maturity [9]. Age-related changes are not restricted to localized or regional connections. Systems-wide changes in network integration and segregation associated with brain maturation have been described using graph-theoretic techniques, a quantitative framework for describing complex networks such as the brain [10,11]. Studies have shown the emergence of small-world network topology as early as the second trimester and reduced segregation with advancing GA based on decreasing modularity, normalized clustering, and local efficiency [5,12,13]. While trends have been described, it is unclear when the transition to mature network patterns occurs.
When studying brain development, individuals are commonly grouped into 'young' and 'old' cohorts [13][14][15]. In some cases, this boundary is clear, such as when comparing the brains of adults versus children. Functional connectivity-based cut points are less defined during gestation, for example, when comparing brain maturity at mid-and latefetal stages. The goal of our study was to estimate optimal gestational-age (GA) cut points, or transition periods, where we can subdivide fetal groups based on age-associated changes in network topology and perform meaningful between-group comparisons. In this study, cut points were defined for global functional network metrics that we have recently described; these metrics are normalized clustering and path length, global and local efficiency, modularity, and the small-world index, metrics we have previously described in this fetal cohort [5]. We selected GAs between 28 and 33 weeks as potential cut points based on previous studies reporting transitional electrophysiologic events [16,17] and restingstate studies that demonstrated significant connectivity changes around this period [4,9]. We hypothesized that optimal cut points would likely fall between 30 and 32 weeks coinciding with the dissolution of the transient fetal compartments and consolidation of thalamocortical connections.

Participants
A total of 110 data sets from fetuses of 95 pregnant women with healthy pregnancies were included in the study. RS-fMRI data were collected as part of a prospective, longitudinal study at Children's National in Washington DC investigating pre-and postnatal brain development in the setting of complex congenital heart disease. We previously described global network topology in this cohort in a recently published study [5]. Conventional T2 for all fetuses revealed structurally normal brains. Pregnant women with known psychiatric/metabolic/genetic disorders, complicated pregnancies (i.e., preeclampsia and gestational diabetes), multiple pregnancies, alcohol and tobacco use, maternal medications, and contraindications to MRI were excluded from the study. For fetuses, those with dysmorphic features on antenatal ultrasound, chromosomal abnormalities by amniocentesis, presentation after 28 weeks' gestational age, and evidence of congenital infections were excluded.

Preprocessing of Resting-State Data
Fetal RS-fMRI data were preprocessed to attenuate the effects of noise on the measured blood oxygen level dependent (BOLD) signals. These steps were previously described in [5,18]. EPI images were slice time corrected, followed by removal of the first four RS volumes to allow for magnetic gradients to stabilize. The data were then oriented to radiologic orientation and despiked. We also performed bias-field correction to address image intensity nonuniformities due to variations in the magnetic field [19]. Next, we corrected for head motion using an algorithm validated in fetuses and newborns [20,21].
Then, RS images underwent intensity normalization to a global mode of 1000 [22]. This was followed by smoothing using an isotropic 4.5 mm full-width half-maximum Gaussian kernel. Band-pass filtering, retaining signals in the range 0.01 Hz-0.1 Hz, and nuisance regression were then simultaneously performed. Along with regression, censoring of high motion frames and volumes with a high number of voxel intensity outliers was performed [23,24]. High motion was defined as frame-by-frame translational and rotational motion >1 mm and >1.5 • [7,25], respectively. Frames where more than 10% of voxels had intensities deviating from the voxel time series' median absolute deviation were also censored from the time series. Regressors included in the general linear model were white matter/CSF signals [26][27][28], linearly detrended rigid motion parameters, and their temporal derivatives [29]. Residual BOLD signals from this regression were analyzed.

Graph Construction
The fetal functional connectome was formed from 200 regions of interest (ROIs, or nodes) (see Supplementary Figure in [5]) defined using a widely used functional clustering technique [30]. In this approach, a normalized-cut spectral algorithm was used to group voxels into nonoverlapping, functionally homogenous ROIs. Clustering was performed using the temporal correlation between voxel BOLD signals, followed by a 2-level clustering approach to create group-level parcellations. As previously described [5], the BOLD signals for each ROI were measured by averaging signals from high-quality voxels that make up each region [31]. Then, all possible pairwise Pearson correlations for all ROIs were computed yielding a 200 × 200 matrix (or 19,900 correlations). The correlation, r, between an ROI pair is its functional connectivity. Only significant positive connections (p FDR < 0.05) were included [32]. In addition, for reliable graph estimates, thresholds were applied such that nodes of individual graphs were at least 95% connected and had average degree k > 2 * ln (number of nodes) [33,34]. Resulting graphs had an edge density range between 0.10 and 0.51. Global network metrics were computed from each fetus' 200 × 200 binary, undirected graphs at a 0.01 interval. Metrics were averaged across the full density range [35].

Graph Analysis
We used the Brain Connectivity Toolbox (https://sites.google.com/site/bctnet/, accessed on 2 June 2021) to compute normalized (i.e., relative to 100 reference random graphs) values of clustering coefficient (γ) and path length (λ), global efficiency (GE), local efficiency (LE), modularity (Q), and the small-world index (σ) [36,37]. Global metrics related to network segregation and integration were computed. Clustering, local efficiency, and modularity were used to describe network segregation. Normalized path length and global efficiency were used to describe integration. Briefly, the clustering coefficient describes the tendency of neighbors of a node to also be linked to each other. Local efficiency describes the ease of communication among neighbors of one node when that node is removed [38]. Modularity captures the tendency of regions to form densely connected subnetworks while being sparsely linked to other clusters [39]. Characteristic path length is the average shortest distance between two nodes in a network. Global efficiency, the inverse of path length, determines the ability to transmit information across the network. The small-world index describes the balance between segregation and integration commonly found in complex networks such as the brain. A network is considered small-world when σ > 1. For a detailed, mathematical description of these metrics, see [36].

Statistical Analysis
We modeled the effect of GA at scan on connectivity metrics using separate repeatedmeasures generalized estimating equations (GEE). A repeated-measures approach was utilized for GEE models, given that 15 patients had two scans. In fetuses with two scans, the mean interval between scans was 7.38 ± 2.81 (mean ± SD) weeks. The interval range was 3-12.86 weeks. Our modeling strategy took three approaches. First, separate overall GEE models per metric were computed with estimated GA slopes at select GA cut points (i.e., 28-33 weeks) calculated. Second, the first approach was repeated, with the participants stratified by gender. Third, the first approach was again repeated with the fetuses grouped into 'young' or 'old' cohorts based on each GA cut point (e.g., 28 weeks vs. 29+ weeks, ≤ 29 weeks vs. 30+ weeks, etc.), with models stratified by cohort. Lastly, we used the quasi-likelihood independence criterion (QICu) statistic to compare model fit between 'old' and 'young' cohorts. The QICu is analogous to the Akaike Information Criterion (AIC), and models with a smaller statistic are preferred. Had our data been binary, traditional cutoff analyses would be based on a binary outcome with subsequent use of the Youden's Index or some other combined measure of sensitivity or specificity. Given that our connectivity outcomes are continuous, we defaulted to the QICu for determining cutoff. The QICu is generally used to distinguish between regression models with a working correlation structure. For consistency, we used an identical correlation structure (exchangeable) across models. Model fit indices such as the -2 loglikelihood (-2LL) utilized in logistic regression are traditionally relied upon to build out or contract models by reducing or increasing parameters. Unlike the -2LL function for binary outcomes data (e.g., logistic regression), the QICu does not utilize a chi-square function that easily translates into statistically significant testing of the addition or subtraction of parameters. There is no currently available way to statistically test the QICu. For this exercise, we are not model building in a traditional way, i.e., adding or subtracting covariates. Each model across each individual GA cut point was identical to avoid issues of parsimony and fitting-we did not want any covariate effects to essentially affect the GA-connectivity relationship. Our model fitting is essentially static. GA cut points were then assessed following graphical depiction (see Figure 1). For each approach, to test the robustness of our GEE connectivity models, GA estimates were evaluated using a bootstrap approach with 1000 iterations (sampling rate = 75%) from which 95% CIs were derived. Confidence intervals were derived using the 2.5 th and 97.5 th percentiles generated from the bootstrapped estimates distribution. For each connectivity metric, we then utilized a graphical approach to compare QICu fit indices for 'young' and 'old' cohorts at each GA cut point by plotting both groups against each other. All analyses were performed using SAS (ver. 9.4, Cary, NC, USA).

Results
A total of 110 fetal resting-state scans from 95 healthy fetuses (49 females, 46 males) between 19.14 and 39.71 gestational weeks (median: 34.93; 25, 75 IQR: 31.29, 36.57) were analyzed. All fetuses in the study were eventually born full term, with a median GA of 39.57 (25,75 IQR: 39.71,40.29). After rigorous preprocessing, an average of 108 ± 17 volumes were available for analysis. Of these, we only included 80 volumes (i.e., 4 min scan duration) for each participant to reduce bias introduced by variability in time series length. Nevertheless, for this cohort, global metrics for partial time series data closely estimated metrics computed from all available time points (Supplementary Figure S1). Average and maximum frame-by-frame head displacement [24] did not correlate with GA at scan or any of the global network metrics. Demographic and motion findings were previously reported in [5].

Optimal GA Cut Point
For each GA cut point, we then stratified the cohort into 'young' or 'old' groups based on values below or above the GA threshold of interest (Table 2) and subsequently separately modeled the GA on metric association. The proportion of fetuses allocated to the 'younger' cohort based on GA ranged from 11.8% (GA = 28) to 30.9% (GA = 33). Across all metrics, among fetuses classified to the 'younger' cohort, the percent change in estimated GA slope ranged from a low of -9.84% (Q) to a high of -39.21% (σ). Lastly, we present model performance QICu metrics (Table 3, Figure 1) for each metric stratified by younger and older cohorts for each GA cut point. For each metric, at each GA cut point, we compared younger and older QICu scores using a modified Youden index to determine the optimal cut point. Optimal cut points fell within the range of 30-31 weeks for all metrics. Supplementary Figure S2 shows global metrics relative to this cut point.

Discussion
We statistically defined optimal gestational-age cut points that can be used to dichotomize human fetuses into 'young' and 'old' subgroups based on resting-state functional connectivity. Most network measures' trends, except for global efficiency, decreased with advancing gestational age. Statistically, optimal GA cut points for small-world index, normalized clustering and path length, global and local efficiency, and modularity were at 30-31 weeks, coinciding with transitional physiologic events that have been shown to be critical to the development of fetal neural circuitry.
That network organization evolves during gestation is not surprising, given how rapidly the brain changes in the third trimester [40]. For example, reduced network segregation with advancing fetal gestational age, demonstrated by decreasing modularity, normalized clustering, and local efficiency, has previously been reported [5,9]. In our previous work, we focused on identifying trends in global network metrics across gestation; here, we focused on identifying transition points where we could meaningfully split fetal cohorts when examining age-related changes in the fetal connectome. The late second to early third trimester of gestation marks a transition in the neural circuitry of the fetal brain, from immature, autonomously generated neuronal bursts to a more sensory-driven pattern [16,41]. The neurobiological underpinnings of the identified inflection points in the network metrics are undetermined, but we speculate that the reported changes in the fetal connectome may be related to maturing electrical patterns. Before 30-32 weeks, electrical activity in the brain is centered around the subplate [42][43][44]. During the midfetal period, or around 15 weeks GA, transient fetal compartments essential to the formation of the fetal neural circuitry have already been formed. The subplate, specifically, serves as a waiting area where postmigratory neurons temporarily connect with subcortical afferents [45]. The establishment of these neural circuits has been linked to the emergence of spontaneous activity transients (SATs), oscillatory neuronal activity that predates evoked brain responses, on electroencephalography recordings in extremely premature infants [46,47]. At around 28-30 weeks, with subcortical afferents synapsing with the cortical plate and refinement of these developing circuits, there is eventual dissolution of the subplate. Electrical activity gradually shifts to the more permanent thalamocortical (subcortical-cortical plate) circuitry, and SATs are replaced by complex cortical EEG patterns [46,48]. The histologic and electrical events described coincide with the identified cut points providing a potential substrate for the transitions seen in the fetal connectome at around 30-31 weeks.
The link between physiologic events and age-related changes in the connectome is currently speculative, but a previous study by Jakab and colleagues [4], where they examined resting-state connectivity in 32 healthy fetuses, also emphasized the relevance of this period. They showed that functional connectivity rapidly increased between 24 and 31 weeks' gestation and gradually stabilized thereafter, likely due to the third-trimester consolidation of thalamocortical connections.
It is important to point out the limitations of our study. First, we had few late second-to early third-trimester fetuses in our group. A uniformly distributed sample across gestational weeks would enable more robust statistical testing of optimal cut points. Nevertheless, we based our transition points on more than 100 resting-state data sets, a substantially large sample that we think provide reasonable baseline cut points for future investigations. Second, the coupling between the BOLD signal and neural activity, especially in utero, is poorly understood [49]. More studies are needed to elucidate the neural mechanisms underlying the relationship between gestational age, BOLD signal changes, and global network topology. A methodological consideration that bears mention is node selection. Here, we opted to use a functional clustering algorithm to define brain regions. Studies have suggested that global and intermediate metrics are stable across different parcellation strategies, but future studies that systematically test the influence of atlas choice on cut points would help refine our current estimates. In-scanner fetal head motion remains a challenging issue in MRI, and excessive head motion led to the exclusion of a substantial number of otherwise healthy fetuses from our study cohort. To minimize the effect of motion on functional connectivity, we performed a rigorous visual quality assessment, used motion correction algorithms optimized for fetuses, and censored high motion volumes from the hemodynamic time series. Of note, and as previously reported, we did not observe any association between the metrics reported in this study and head motion and gestational age at scan [5]. Lastly, here, we only focused on identifying transition points in global network features which may be related to brain maturation. We did not address whether global network metrics changed in a linear fashion throughout gestation [5,9], an important avenue to explore in future studies.

Conclusions
We identified GA cut points based on global network topology that could potentially be used for future analyses comparing 'young' and 'old' fetuses. The period between 30 and 31 weeks coincides with a period in neural circuitry formation where spontaneous activity of early circuit transitions to sensory-driven activity. Additional studies are needed to understand the mechanisms underlying the association between gestational age, neural circuitry changes, and global network topology.