Functional dissection of prenatal drug effects on baby brain and behavioral development.

Prenatal drug exposure (PDE) is known to affect fetal brain development with documented long-term consequences. Most studies of PDE effects on the brain are based on animal models. In this study, based on a large sample of 133 human neonates and leveraging a novel linear mixed-effect model designed for intersubject variability analyses, we studied the effects of six prenatally exposed drugs (i.e., nicotine, alcohol, selective serotonin reuptake inhibitor, marijuana, cocaine, and opioids) on neonatal whole-brain functional organization and compared them with five other critical nondrug variables (i.e., gestational age at birth/scan, sex, birth weight, and maternal depression). The behavioral implications were also examined. Magnitude-wise, through summing across individual drug effects, our results highlighted ~5% of whole-brain functional connections (FCs) affected by PDE, which was highly comparable with the combined effects of the five nond rug variables. Spatially, the detected PDE effects featured drug-specific patterns with a common bias in higher-order brain regions/networks. Regarding brain-behavioral relationships, the detected connections showing significant drug effects also demonstrated significant correlations with 3-month behavioral outcomes. Further mediation analyses supported a mediation role of the detected brain FCs between PDE status and cognitive/language outcomes. Our findings of widespread, and spatially biased PDE effect patterns coupled with significant behavioral implications may hopefully stimulate more human-based studies into effects of PDE on long-term developmental outcomes.

INTRODUCTION

Drug use among pregnant women is on the rise in the United States with an estimated 5–6% prevalence for illegal drugs and 16–17% for legal drugs (Substance Abuse and Mental Health Services Administration [SAaMHSA], 2013). Prenatal drug exposures (PDEs) affect the developing fetus through a variety of mechanisms including maternal physiology, placental function, and direct disruptions of endogenous neural signaling of the fetal brain (Lester & Padbury, 2009). Consistently, adverse effects of PDE on developmental outcomes in different domains, including learning, memory, attention, emotion regulation, and executive functions are frequently reported (Ross, Graham, Money, & Stanwood, 2015). However, human studies on the brain basis of PDE effects are scarce and most existing ones focus on later developmental periods thus are often complicated by other adverse postnatal risk factors accompanying PDE (Crandall et al., 2007; Ren, Malanga, Tabit, & Kosofsky, 2004; Thadani, 2002; Thompson, Levitt, & Stanwood, 2009; Tronick & Beeghly, 1999). Therefore, human studies of PDE effects during the neonatal period, a time more proximal to the drug exposures and less affected by postnatal factors, represents an urgent research priority to better understand PDE‐related neurobehavioral sequelae.

Although scarce, existing human studies of PDE effects point toward significant disruptions in distributed brain regions and circuits (Derauf, Kekatpure, Neyzi, Lester, & Kosofsky, 2009; Morie, Crowley, Mayes, & Potenza, 2019; Ross et al., 2015). For example, our previous studies on prenatal cocaine exposure reported reduced prefrontal volume (Grewen et al., 2014), aberrant functional connectivity (FC) between the medial prefrontal cortex and the amygdala (Salzwedel et al., 2015), as well as disrupted thalamus connectivity to different cortical targets (Salzwedel, Grewen, Goldman, & Gao, 2016). In a smaller subset of neonates with prenatal marijuana exposure we also reported disrupted amygdala, insula, and striatal FC (Grewen, Salzwedel, & Gao, 2015). While informative, each of the existing studies examined a single drug‐type based on a case–control design and only focused on prehypothesized regions of interest. These studies are unlikely to fully unveil the complex and intertwined effects of different drug exposures on the whole‐brain functional system. Moreover, the applied case–control design is unable to compare the potential PDE effects with other critical nondrug participant characteristics (PCs, e.g., age, birth weight, sex, maternal traits), which are known to be related to brain and behavioral development (Gilmore, Knickmeyer, & Gao, 2018). Therefore, a whole‐brain system‐wide study that simultaneously models the potential effects of multiple drugs as well as other nondrug PCs is urgently needed to improve our overall understanding of the relationships between PDE, the developing brain, and behavior.

In this study, leveraging a cohort of 133 newborns including 75 neonates with PDEs to 1 or more of 6 drugs (i.e., nicotine, alcohol, selective serotonin reuptake inhibitors [SSRIs], opioids, cocaine, and marijuana) and 58 drug‐free control neonates with 2‐week magnetic resonance imaging (MRI) and 3‐month behavioral assessment, we sought to systematically characterize different drugs' effects on the developing brain and behavior. FC measures based on resting‐state functional MRI (rsfMRI) (Biswal, Yetkin, Haughton, & Hyde, 1995) scans were used as the primary brain measures. In addition to the six drug variables, nondrug PCs including gestational age at birth, gestational age at scan, birth weight, sex, and maternal depression scales were also included and similarly assessed/compared. A novel linear mixed‐effect model (LME) specifically designed to handle intersubject variability (ISV) measures was employed to detect the multivariate relationships among drug/nondrug PCs, brain FC, and behavioral outcomes (Chen, Taylor, Shin, Reynolds, & Cox, 2017). This model was designed to answer the question of how ISV in one domain relates to ISV in the other (Seghier & Price, 2018). Given the well‐documented effects of PDE on the developing fetal brain in both human (Behnke et al., 2013) and animal models (Castaldo et al., 2010; Thompson et al., 2009), we hypothesize wide‐spread PDE effects on neonatal FC. Based on the reported dominance of higher‐order brain function disruptions associated with PDE (Ross et al., 2015), we expected a bias of PDE‐affected connections within high‐order functional regions/networks. Finally, with the brain lying at the interface between PDE and behavioral outcomes, we further hypothesized that FCs affected by PDE would mediate the relationship between PDE status and 3‐month behavioral outcomes. Consistent with these expectations, three main findings were reported in this study. First, magnitude‐wise, the estimated overall effects of PDE (through summarizing across individual drug effects) on newborn FC were wide spread (account for ~5% of whole‐brain connections) and highly comparable with the combined effects of five nondrug variables including sex, gestational age at birth, birth weight, gestational age at scan, and maternal depression. Second, spatially, although drug‐specific heat maps of effects were observed, one common pattern of biased distribution among higher‐order functional networks was observed when comparing PDE effects with nondrug PC effects. Finally, significant behavioral correlations were observed between the detected drug‐related connections and 3‐month cognitive scores. These findings greatly improved our understanding of different drugs' effects on newborn brain FC and may prompt more longitudinal studies of PDE effects on the developing brain and behavior.

MATERIALS AND METHODS

Participants

The study cohort (N = 133) consisted of healthy drug‐free control infants (CTR; n = 58) and infants with PDE (n = 75), including; cocaine, marijuana, alcohol, nicotine, SSRIs (setraline, cipramil, or fluoxetine), and opioids (heroin, oxycodone, oxycontin, methadone and the mixed agonist/antagonist, buprenorphine). Pregnant women were recruited in the third trimester of pregnancy. Primary recruitment sites for the drug group were local residential and outpatient treatment programs for women with perinatal substance abuse and their children. In addition, we recruited CTR and PDE mothers from Chatham, Orange, Durham, Alamance, and Wake County Health Department obstetric clinics, the University of North Carolina hospital low‐income obstetrics clinic, and flyers, local advertisements, and Craigslist. At enrollment, mothers were required to be between 18 and 44 years of age and free from (a) chronic medical or psychiatric disease, (b) untreated current clinical depression or anxiety disorder, and (c) language barrier that might prevent informed consent. Several subjects with opiate abuse were treated with methadone or buprenorphine maintenance treatment for part of their pregnancies. The infants took part in the imaging experiment during the neonatal period (2–6 weeks of age) and participated in behavioral assessments at approximately 3 months of age. All infants were required to be living with biological mothers at time of testing. Infants were excluded for multiple reasons, including gestational birthweight <4.75 lbs, delivered at <32 weeks or >42 weeks gestation, history of mechanical ventilation or surgery of any kind, >24 hours in neonatal intensive care unit, or chronic illness of any kind.

Mother–infant dyads were characterized and compared on the following criteria: sex, gestational age at birth, gestational age at scan, birth weight, PDE, infant cognitive, language and motor functions, and maternal depression at the time of scan. In addition, maternal body mass index (MBMI) and education (MEDU) were available in different subsamples (n = 123 for MBMI and n = 110 for MEDU). MBMI did not differ between PDE and control groups while MEDU was higher in the control group (Table S1). Therefore, post hoc analyses were conducted for MEDU in the subsample to test the robustness of our primary results against this variable. All participants were tested for prenatal drug use using interviews, medical record review, and postnatal urine toxicology at study visits. PDE status was based on three criteria: (a) maternal self‐report with Time Line Follow Back interview (Robinson, Sobell, Sobell, & Leo, 2014) conducted in third trimester and again at neonatal MRI visit; (b) response to a questionnaire about maternal substance use done at 3 months; and (c) medical record queries of prenatal urine toxicology. Maternal self‐report or positive urine toxicology for drug use qualified the mother–infant dyad for PDE status. Maternal depression was indexed by score on the Edinburgh Postnatal Depression Scale (EPDS) at time of scan (Murray & Carothers, 1990). Infant behavior was assessed using the Bayley III Scales of Infant and Toddler Development (Bayley, 2006) for cognitive, language, and motor development. Research staff that conducted the Bayley assessments was blinded to drug‐exposure status. Group means (PDE vs. CTR) for each PC were compared using students t tests. See Appendix S1 for a detailed description of the Bayley III. This study was approved by the University of North Carolina at Chapel Hill and Cedars‐Sinai Biomedical Institutional Review Boards.

Imaging

Each participant was fed, swaddled, and fitted with ear protection prior to imaging. Subjects were sleeping during image acquisition. Head position was secured in the scanner using a vacuum‐fixation device. Vital signs (heart rate, SaO2) were monitored continuously by a nurse throughout the examination. Data were collected using two scanners: (a) 3‐T head‐only Siemens Allegra with circular polarization head coil and (b) 3‐T Siemens Tim Trio with 32‐channel head coil. T1‐weighted structural images were collected using a 3D magnetization‐prepared rapid acquisition gradient echo pulse sequences: repetition time (TR), 1,820 ms; echo time (TE), 3.75 ms; inversion time, 1,100 ms; flip angle, 7°; 144 slices; voxel size, 1 mm3. Resting‐state fMRI data were acquired using a T2*‐weighted echo planar imaging pulse sequence: TR, 2 s; TE, 32 ms; 33 slices; voxel size, 4 mm3; number of volumes, 150 (5 min).

Image preprocessing

Functional data were preprocessed using the FMRIB Software Library (FSL v5.0.8) (Jenkinson, Beckmann, Behrens, Woolrich, & Smith, 2012), Analysis of Functional NeuroImages suite (AFNI v16.0.10 February 25, 2016) (Cox, 1996), and MATLAB (R2018a). Steps included discarding initial volumes (n = 3), motion correction, motion censoring (Power, Barnes, Snyder, Schlaggar, & Petersen, 2012; i.e., frame‐wise displacement [FD] >0.3 and <3 continuous volumes], interpolation of censored time points using an autoregressive model (MATLAB gapfill), and band‐pass filtering (0.01–0.08 Hz). Subjects with less than 90 volumes postscrubbing were excluded (n = 13). Confound regression was used to reduce distance dependence (Ciric et al., 2017). Specifically, the confound regression strategy consisted of motion censoring plus a 32‐parameter nuisance signal model; 8 regressors (i.e., eroded white matter [WM], eroded cerebral spinal fluid [CSF]), and six parameters corresponding to rigid‐body motion correction), their derivatives, quadratic terms, and squares of derivatives. Nuisance signals were also band‐pass filtered to prevent frequency‐dependent mismatch (Hallquist, Hwang, & Luna, 2013). Nuisance regression was performed via linear regression (AFNI 3dTproject). Censored time points were ignored so as not to influence the fit and then excised from the data. Finally, the global signal was extracted using a whole‐brain mask, excluding the eroded WM and CSF regions, and then regressed from the data.

Spatial normalization was achieved using a combination of within‐ and between‐subject transformations. The University of North Carolina (UNC) neonate template was used for co‐registration (Shi et al., 2011). Specifically, within‐subject functional‐to‐anatomical alignment was achieved using rigid‐body registration (FSL flirt). The inverse of this transformation was used to align the tissue‐specific (i.e., WM, CSF, gray matter [GM]) masks for confound signal extraction. Subjects with very poor structural images (i.e., excessive “banding”) due to in‐scanner motion were excluded (N = 5). Between subject anatomical‐to‐standard alignment was achieved using nonlinear warping (FSL fnirt). Functional‐to‐anatomical and anatomical‐to‐standard transformations were then combined and applied to the preprocessed functional data. Finally, the data were spatially smoothed using a Gaussian filter (FWHM = 6 mm).

FC measures

FC measures were derived using an functional parcellation‐based atlas (UNC‐CEDARS INFANT) (Shi, Salzwedel, Lin, Gilmore, & Gao, 2018). Specifically, average time series from all seed region (n = 222) were used to construct a Fisher's Z‐transformed correlation matrix (Biswal et al., 1995; Taylor & Saad, 2013) for each subject. Note, the original atlas contains 223 regions; however, one region was lost when the atlas was converted to the resolution of the functional data (4 mm3). For each FC measure (i.e., cell in the FC matrix), the effects of scanner and subject motion (residual FD [rFD] and number of volumes scrubbed [NVOLS]) were further minimized using post hoc regression. The resulting residual FC measures were then used in the ISV analysis.

ISV analysis

ISV measures

The goal of this study was to characterize the relationships between ISV of different PC variables and the ISV of brain FC measures. Specifically, for each measurement we computed pair‐wise ISV using the Euclidean distance metric where higher distances between subjects equate to higher ISV. For FC measures (i.e., cells in the correlation matrix) and continuous PCs (e.g., birth weight, gestational age at birth and scan, maternal depression, and motion [rFD + NVOLS) this resulted in continuous and positive difference measures for each subject pair. Categorical measures (i.e., sex, drug exposure, and scanner) were similar; however, the resulting difference measures were binary; 0 for subject pairs that share the same label for the variable of interest (e.g., male–male or female–female, with cocaine exposure‐ with cocaine exposure) and 1 for subject pairs that differ on the status of the variable of interest (e.g., male–female, with cocaine exposure‐without cocaine exposure) for each subject pair.

The rational to use the ISV model rather than conventional models based on raw individual FC measures is twofold. First, there is an increasing interest in neuroscience research to directly study ISV (Gao et al., 2014; Mueller et al., 2013). Through directly modeling ISV as the dependent variable, our results may more specifically reflect how individual differences in certain drug or nondrug variables influence individual differences in brain FC. Second and more importantly, the selected ISV model may better utilize the latent information embedded in our data set to allow more sensitive detection of drug and/or nondrug effects. Specifically, given N subjects in question, for each subject, by calculating his/her pair‐wise FC difference with every other subjects in the cohort (i.e., ISV), the N − 1 difference vector would explicitly reflect the gradient brain differences between the subject in question and the rest in the cohort, likely related to different degrees of differences in certain PCs (i.e., drug and nondrug variables in this study). Therefore, by directly modeling ISV as the dependent variable and pair‐wise differences in drug and nondrug PCs as independent variables, we have a set of much larger (i.e., total number of observations increase from N to [N * N − 1]/2) and richer information reflecting the fine gradient differences between every pair of subjects. However, given the highly related nature between the ISV model and the raw FC‐based models, we also expect correlated estimates of effects between the two models. To test this hypothesis, parallel general linear models (GLMs) using either ISV or raw FC values as the dependent variable, were carried out. Overall, with the proposed ISV model, we aim to more explicitly and better answer (i.e., with richer information thus potentially higher sensitivity) the question of how individual differences in drug and nondrug variables affect individual differences in brain FC.

Modeling and related statistics

Pair‐wise intersubject differences in FC and different PC variables were modeled using a novel LME model with crossed‐random effects (CRE) (Chen et al., 2016, 2017). Briefly, this technique leverages random‐effect terms to account for the dependence in the ISV data structure (i.e., one of the subjects is shared across several subject pairs that include this subject). Note, as explicitly pointed out in our original paper presenting this model (Chen et al., 2016, 2017), the (N * N − 1)/2 pair‐wise ISV values are not totally independent since every subject is repeated in N − 1 observations. Therefore, special attention has to be given to properly model the degree of freedom to control false positives. Specifically, our LEM model properly controls for false positives by adjusting the canonical LME statistics using the true sample size (i.e., number of unique subjects) rather than the number subject pairs to calculate the degree of freedom. Like GLM, LME‐CRE produces effect estimates, SEs, t statistics, and p values in the multivariate sense and the specific model we used was:where ISV (FC) is the FC ISV measure for a given connection across subject pair X and Y (i.e., the dependent variable), ISV (PC) represents pari‐wise differences in PC variable i between Subjects X and Y (i.e., independent variables), and (1 | Subject) represents the random effects added to all subject pairs that share Subject or Subject. Significance was assessed using the false discovery rate (FDR) correction across connections and PC terms (FDR; q = 0.05). To compare our ISV models with the conventional model using raw FC values, results based on a similar general linear regression model (GLM) were also obtained: FC ‐ PC1 + PC2 + ⋯ + PC. Where FC represents the raw FC value for individual subjects and PC represents the raw PC variable for individual subjects.

Spatial characterizations and specificity analyses

After the detection of FCs showing significant PC effects, we computed the number of significant effects (i.e., degree) within each region for each PC variable (e.g., cocaine status) separately. These values were then projected onto adult space to create PC‐specific “heat maps” showing the whole‐brain distribution of effects (BrainNet) (Xia, Wang, & He, 2013). At network level, the Cole‐Anticevic Brain‐wide Network Partition (CAB‐NP) (Spronk et al., 2018) was used to generate network‐level summary measures. Specifically, regions within the neonatal UNC‐INFANT functional atlas where assigned to networks using a winner‐takes‐all‐approach based on spatial overlap with the CAB‐NP (Figure S1 and Table S2). For every region we were thus able to generate a network assignment. The CAB‐NP is linked to the Human Connectome Project—Connectome Workbench software and data formats, therefore we used a series of commands to convert the CAB‐NP into a standard volumetric format for the spatial computations across atlases (see Appendix S1 for more details).

We used Fisher's exact tests to determine the specificity of effects at the network‐level. Fisher's exact test determines if there is a nonrandom association between two categorical variables, summarized as an odds ratio (OR), 95% confidence interval (CI), and corresponding p value. Specifically, we created contingency tables detailing the number of significant effects for each PC or grouping variable (i.e., drug or nondrug) and each network or network grouping variable (i.e., primary or high order). For individual PCs and network assessments there were many separate tests therefore significance was assessed using FDR correction (q = 0.05) across all tests.

Behavioral analyses

For brain–behavior analyses (i.e., FC–BEHAVIOR), we repeated the ISV model within connections showing a significant PC–FC relationship. The behavioral distance measure (Euclidean distance in cognitive, language, or motor composite scores) was the only explanatory variable of interest. Separate models were implemented for each behavioral measure and then the results were pooled within the previously identified PC–FC connections to establish an overall significance threshold (q = 0.05, FDR corrected). We also evaluated the multivariate performance of the identified brain–behavior connections in explaining the variance of each 3‐month behavioral outcome measure at the level of the residual FC measures using GLM. The model performance (adjust R 2) was compared across different groupings of identified connections (i.e., drug, nondrug, and drug + nondrug effects).

The ISV model was again used to assess the direct relationships between drug and nondrug PCs with behavioral measures. The following variables were included in the PC–BEHAVIOR model: sex, gestational age at birth, birth weight, maternal depression (EPDS), and the drug status, to detect significant drug or nondrug PC–behavior relationships. To test for mediation effects, we generated canonical variables representing the combined effects of either drug or nondrug‐related connections detected previously. Briefly, a canonical model (MATLAB canoncorr; FC1 + FC2 + ⋯ + FC–BEHAVIOR) was generated within each behavioral domain and for each set of identified brain–behavior connections (drug or nondrug). The sample canonical coefficients corresponding to the FCs were then used to compute ISV across subjects, resulting in a single variable (i.e., Canonical FCISV) for drug and nondrug‐related connections, respectively. These canonical FC variables were then incorporated in the corresponding original PC‐behavioral model showing significance (i.e., BEHAVIOR–PC + Canonical FCISV), separately, to test for mediation effects. Significant mediation was detected if the original significant PC–behavior relationship became nonsignificant after inclusion of the corresponding canonical variable.

RESULTS

Study cohort; PCs, and behavioral outcome measures

PCs and behavioral outcome measures for the study cohort (n = 133) are summarized in Table . The distribution of drug exposures revealed a predominance of polydrug effects; 62 out of 75 (i.e., 82.7%) neonates with PDE were exposed to more than 1 drug type and 30 unique polydrug profiles were identified (Table S3). For the six drugs characterized, the number of exposures varied (in descending order); nicotine > cocaine > marijuana > alcohol > opioids > SSRIs. In addition to the six drug PCs, five nondrug PCs (sex, gestational age at birth, birth weight, gestational age at scan, and maternal depression scales [EPDS]) were also available in the full sample and examined similarly. Full domain (cognitive, language, and motor) behavioral outcome measures at 3 months of age, indexed by performance on the Bayley Scales of Infant and Toddler Development III (Bayley, 2006), were also available in a subset of subjects (n = 80). Consistent with previous studies(Ross et al., 2015), significant group (CTR vs. PDE, p < .05) differences were observed for multiple PCs and two of the three behavioral measures. Specifically, the PDE group was on average born marginally earlier (gestational age at birth, p = .089), lighter (birth weight, p = .006) and the PDE mothers were on average more depressed (EPDS, p = .004). The groups did not significantly differ in sex distribution, race, and age at scan. For 3‐month behavioral outcomes, neonates with PDE scored significantly lower in the language (p = .018) and motor domains (p = .012) and marginally lower on cognitive measures (p = .069).

TABLE 1

Participant characteristics and behavioral outcome measures

	PDE (n = 75)				CTR (n = 58)				CTR versus PDE		ALL (n = 133, 67 F)
	Mean	SD	MIN	MAX	Mean	SD	MIN	MAX	t (131)	p	Mean	SD	MIN	MAX
Gestational age at scan (days)	305.44	13.16	275	345	306.43	11.40	285	330	0.46	.650	305.87	12.38	275	345
Gestational age at birth (days)	275.56	10.28	234	290	278.38	8.16	251	290	1.71	.089	276.79	9.49	234	290
Birth weight (pounds)	7.07	1.11	4.75	10.63	7.59	0.99	5.375	10.5	2.79	.006	7.29	1.09	4.75	10.63
Maternal depression (EPDS)	5.74	5.65	0	25	3.35	3.04	0	13	‐2.91	.004	4.70	4.82	0	25

Note: Values in italic indicates p < 0.1 and values in bold indicates p < 0.05.

Abbreviations: CTR, drug‐free control infants; EPDS, Edinburgh Postnatal Depression Scale; PDE, prenatal drug exposure; SSRI, selective serotonin reuptake inhibitor.

The magnitude of PDE effects on whole‐brain FC in comparison with those associated with nondrug PCs

The relationships between drug/nondrug PCs and neonatal brain FC were characterized using a LME with ISV of FC as the primary dependent variable and the ISV of 11 PCs (i.e., the 6 drugs plus 5 nondrug PCs) as primary independent variables. Scanner and motion parameters (i.e., rFD and the number of volumes scrubbed) were also included as control variables. Sex and drug PCs were modeled as categorical variables while the other four nondrug PCs as continuous variables. By modeling all 11 PCs simultaneously, we aimed to detect drug‐specific effects on brain FC while controlling for effects from other drug and nondrug PCs. The neonatal brain was divided into 222 functionally homogenous regions based on our previously generated neonate‐specific functional atlas (Shi, Salzwedel, Lin, Gilmore, & Gao, 2017), yielding a 222*222 FC matrix for each subject (Figure 1a, “upper triangle” showed the group mean). Among the 24,531 unique connections assessed covering the whole brain, 2,553 or 10.4% showed one or more significant effects (p < .05 after FDR correction) related to at least 1 PC examined (Figure 1a, “lower triangle”), resulting in 2,724 total significant effects. Minimal effects from scanner and motion parameters were observed (total number of significant effects for all three control variables was 133, details in Figure S2 and Table S4). Comparing with the conventional GLM model using raw FC measures as the dependent variable rather than ISV, highly correlated estimates were observed for all 11 PC variables of interest (Figure S3). However, only 27.3% (i.e., 745 out of 2,724) of the observed significant effects in the ISV model also show uncorrected p < .05 significance in the conventional GLE model, confirming higher sensitivity from the proposed ISV model.

FIGURE 1

Degree of PDE effects on whole‐brain functional organization. (a) FC matrix depicting group mean correlation in the upper triangle and detected significant ISV effects (FDR corrected; q = .05) based on multivariate LME‐CRE model in lower triangle. Connections sorted by region‐to‐network‐level affiliation. (b) Ranked proportion of effects in different categories: single drug > single nondrug > combined > multiple nondrugs > multiple drugs. (c) Ranked proportion of effects across individual drug and nondrug PC variables. (d) Distribution of hyper/hypo‐ or positive/negative associations based on individual FC measures rather than ISV for each set of PC‐related effects. CRE, crossed‐random effects; FC, Functional connectivity; FDR, false discovery rate; ISV, intersubject variability; LME, linear mixed‐effect model; PC, participant characteristics; PDE, prenatal drug exposure

Several features regarding the specificity, proportion, and distribution of these effects are noteworthy (Figure 1b,c). First, the effects were highly specific as evidenced by the dominance of singular PC associations (n = 2,390, 93.6%), likely partly resulting from the simultaneous modeling of all 11 PCs in the same model. Second, when summarized across individual PC effects to obtain an overall estimate of drug‐related effects (n = 1,206, 47.2%) and nondrug‐related (n = 1,184, 46.4%) effects, they were highly comparable indicating a similar overall degree of impact. The remaining connections (n = 163, 6.4%) demonstrated multi‐PC associations; multiple drugs (n = 35, 1.4%), multiple nondrugs (n = 41, 1.6%), and combined effects (n = 87, 3.4%). Third, the distribution of effects across PCs varied (in descending order); sex (n = 375, 13.8%) > gestational age at scan (n = 352, 12.9%) > SSRIs (n = 279, 10.2%) > nicotine (n = 271, 10%) > alcohol (n = 250, 9.2%) > maternal depression (n = 241, 8.9%) > gestational age at birth (n = 220, 8.4%) > opiates (n = 222, 8.2%) > marijuana (n = 199, 7.3%) > birth weight (n = 159, 5.8%) > cocaine (n = 147, 5.4%). Lastly, when assessed using individual FC measures rather than ISV, there were relatively balanced hyperassociation/hypoassociation and positive/negative correlations for categorical and continuous PCs, respectively (Figure 1d), indicating connection specific patterns with no dominant hyperconnectivity/hypoconnectivity or positive/negative correlations across all detected effects.

Spatial distribution of PDE and nondrug PC effects on brain FC

After quantifying the overall degree of impacts for drug and nondrug PCs, we next aimed to elucidate the spatial distribution of the observed effects at the regional (Shi et al., 2017) and network (Spronk et al., 2018) levels. First, the distribution of all effects (n = 2,724) encompassing the 2,553 identified connections were examined at the regional level by computing the number of significant connections (i.e., degree) associated with each region for each drug or nondrug PC, separately. These values were then visualized on whole brain surface, creating PC‐specific “heat maps” of the spatial distribution of effects. The heat maps for all six drug effects along with glass–brain representations of the individual connections are shown in Figure 2a (nondrug PCs are presented in Figure S4). Notably, the heat maps were highly distinct for each drug type, featuring high levels of concentration for the following drug—region pairs: SSRI—visual, orbital/medial/lateral frontal; nicotine—medial/lateral frontal; alcohol—sensory and lateral frontal; opioids—middle frontal and angular gyrus; marijuana—medial/lateral parietal, sensorimotor, and orbital/lateral frontal; cocaine—supplementary motor area and medial/lateral frontal.

FIGURE 2

Spatial distribution of PDE effects and their relation to different functional networks. (a) Drug‐specific heat maps depicting regional degree and corresponding glass–brain connection models. Warmer colors and larger nodes represent regions with higher degree or number of significant effects, in the surface and glass brain plots, respectively. (b) Drug versus nondrug composite heat maps. Warmer colors depict regions with more drug‐related effects while cooler colors depict regions with more nondrug‐related effects. (c) Summary of significant network‐level effects grouped by drug and nondrug PCs. A, alcohol; AB, gestational age at birth; AS, gestational age at scan; AUD, auditory; BW, birthweight; C, cocaine; CO, cingulo‐opercula; DA, dorsal attention; DM, default mode; E, maternal depression indexed by Edinburgh Postnatal Depression Scale value (EPDS); FP, fronto‐parietal; M, marijuana; N, nicotine; LA, language; O, opioids; OA, orbital affective; PC, participant characteristics; PDE, prenatal drug exposure; PMM, posterior multimodal; S, selective serotonin reuptake inhibitor (SSRI); Sex, infant sex; SM, sensorimotor; VIS1, primary visual, VIS2, secondary visual; VMM, ventral multimodal

Further comparison of overall drug versus nondrug effects (i.e., degree difference, Figure 2b), revealed a striking pattern that featured low levels of drug effects in certain primary sensory regions (e.g., right sensorimotor cortex and bilateral visual areas) but high levels of expression in association brain areas (particularly bilateral medial/lateral frontal regions and left orbital frontal area). When examined against 12 functional networks (Ji et al., 2019), broadly categorized as either primary (i.e., primary visual, secondary visual, sensorimotor, and auditory) or high‐order (i.e., cingulo‐opercular, dorsal attention, language, frontoparietal, default‐mode, posterior multimodal, ventral multimodal, and orbital affective) Fisher's exact tests showed an overall interaction effect of PC group (drug and nondrug) and network (primary and high order); the odds of drug‐related effects overlapping with high‐order networks was significantly greater (OR: 1.27, 95% CI: 1.12–1.44, p = 1.68E−04), confirming there was preferential localization of drug effects to higher‐order networks. When further decomposed into their network affiliations, 20 significant nonrandom associations were detected (Figure 2c; also see Figure S5 and Table S5). Of these, 13 were drug‐related and only 2 out of them involved primary networks (i.e., SSRI—primary visual and marijuana—sensorimotor). However, for nondrug PCs, the detected associations were more balanced (n = 7, 4 involving primary networks).

Behavioral implications of PDE and nondrug PC‐related FCs

Having demonstrated the relationships between drug/nondrug PCs and brain FC, we sought to further characterize the relationships between the detected connections and 3‐month behavioral outcomes (cognitive, language, and motor scores) in order to evaluate their behavioral significance. To this aim, individual FC‐behavior models were constructed for each behavioral measure and each connection showing significant drug and/or nondrug FC‐PC effects. Overall, 85 significant (q = 0.05, FDR‐corrected) relationships (cognitive: n = 21, language: n = 31, and motor: n = 33) were detected encompassing 73 individual FCs (Figure 3a). Of these, 42 were drug related (single drug: n = 40, multiple drug: n = 2) and 26 were nondrug related (single nondrug: n = 25, multiple nondrug: n = 1). Five connections (4 motor and 1 language) showed combined effects. The degree of overlap between FC–PC and FC–behavior effects varied with SSRIs showing the most pronounced correspondence (Figure S6).

FIGURE 3

Behavioral correlations of the detected functional connections showing either PDE or nondrug PC effects. (a) Glass‐brain models depicting significant (FDR corrected; q = 0.05) FC–BEHAVIOR (cognitive, language, and motor) relationships within identified functional connections showing either significant PDE or nondrug PC effects. (b) Multivariate prediction of behavioral outcome measures using different grouping schema; drug‐related connections (DRUG), nondrug‐related connections (NONDRUG), and all connections (DRUG + NONDRUG). Solid line represents least‐squares fit. Dashed‐lines correspond to 95% confidence bounds. Dots denote individual subjects (n = 80). X‐axis represents actual cognitive, language, and motor scores while Y‐axis represents fitted values based on multivariate regression using corresponding groups of connections (i.e., drug related, nondrug related, and combined). FDR, false discovery rate; PC, participant characteristics; PDE, prenatal drug exposure

Multiple regression analyses were conducted to quantify the explained variance of behavioral outcomes associated with drug‐related connections, nondrug‐related connections, and all significant connections (i.e., drug + nondrug), respectively. FC was significantly (p < .001) predictive of all three behavior scores (Figure 3b ) but the explained variance (i.e., adjusted R 2) differed. Specifically, the drug‐related connections accounted for 42.0% (cognitive), 50% (language), and 31% (motor) compared to 28% (cognitive), 29.0% (language), and 34% (motor) for nondrug connections, respectively. Thus, drug‐related FCs accounted for 1.5–1.7 more adjusted variance for cognitive and language outcomes compared to nondrug‐related connections, while the two were effectively equivalent for motor score correlations.

Finally, having shown significant relationships between PCs–brain and brain–behavior, the relationships between PCs and behavior were further characterized to determine whether drug and nondrug PCs was directly related to 3‐month behavioral outcomes and whether such direct relationships were mediated by brain FC. First, eight significant (q = 0.05, FDR‐corrected across PCs and domains) positive associations were detected for the three behavioral measures, six of which were related to drug PCs and two to nondrug PCs (Table 2). Specifically, maternal depression (p = .001) and SSRIs (p = .001) were significantly associated with cognitive scores (note, cocaine was marginally significant at p = .042). Nicotine (p < .001) and SSRIs (p = .004) were significantly associated with language scores, while cocaine (p < .001), nicotine (p < .001), SSRIs (p < .001), and gestational age at birth (p < .001) were significantly associated with motor scores. For mediation analysis, we included a canonically derived drug‐ or nondrug‐related composite variable (calculated based on the detected connections in the corresponding FC–behavior models) in the PC–behavioral models, and then reassessed their relationships. Our results showed that all four significant drug PC‐cognitive composite/language score relationships became insignificant (p > .05) after including the drug‐related variable in the model (the marginally significant cocaine–cognitive composite score relationship also became insignificant, p = .806). However, the two drug PC motor‐related relationships remained significant. For nondrug PC–behavior relationships, the birth age–motor relationship became insignificant while the maternal depression–cognitive relationship became less significant (p value changed from .001 to .019).

TABLE 2

Direct PCs–behaviora relationships and mediation analyses of the detected brain functional connections

	Cognitive composite				Mediation		Language composite				Mediation		Motor composite				Mediation
	e	SE	t	p	DRUG	NONDRUG	e	SE	t	p	DRUG	NONDRUG	e	SE	t	p	DRUG	NONDRUG
Birth age	0.00	0.02	0.08	.936			0.01	0.02	0.57	.568			0.13	0.03	4.53	.000	0.023	0.575
Sex	0.03	0.28	0.11	.915			−0.17	0.20	−0.84	.401			−0.24	0.31	−0.78	.436
Birth weight	−0.27	0.23	−1.18	.241			−0.28	0.17	−1.66	.100			0.25	0.26	0.97	.336
EPDS	0.19	0.06	3.33	.001	0.008	0.019	0.02	0.04	0.42	.679			0.06	0.07	0.86	.392
Cocaine	0.67	0.32	2.06	.042	0.806	0.028	0.02	0.23	0.07	.944			2.05	0.36	5.61	.000	0.000	0.000
Nicotine	0.37	0.29	1.30	.198			0.76	0.21	3.69	.000	0.127	0.001	1.20	0.32	3.71	.000	0.000	0.000
Alcohol	−0.33	0.42	−0.80	.425			−0.04	0.30	−0.12	.906			−0.45	0.47	−0.96	.338
Marijuana	0.10	0.38	0.27	.786			−0.40	0.27	−1.48	0.142			−0.34	0.43	−0.79	0.434
SSRIs	2.35	0.69	3.42	.001	0.270	0.005	1.47	0.50	2.94	.004	0.106	0.006	3.66	0.79	4.65	.000	0.409	0.000
Opiates	0.38	0.58	0.66	.509			–0.23	0.42	−0.54	.588			−0.15	0.66	−0.22	.823

Note: Bold p values indicate significance (FDR corrected, q = 0.05).

Abbreviations: CTR, drug‐free control infants; E, estimate; EPDS, Edinburgh Postnatal Depression Scale; FDR, false discovery rate; PDE, prenatal drug exposure; SSRI, selective serotonin reuptake inhibitor.

Bayley Scales of Infant and Toddler Development III.

DISCUSSION

PDE has long been a serious public health concern which will only grow amid the recent surge in opioids use and the increasing availability accompanying legalization of marijuana in the United States and other countries. In this study, we used an ISV model to determine the potential impacts of different drugs on neonatal brain FC and compared them with those of nondrug PCs. The behavioral implications and the potential mediating role of the detected brain FCs on drug–behavior relationships were also explored. We showed that 10.4% of the whole‐brain FCs were significantly associated with at least one PC, with the six drug PCs and the five nondrug PCs each accounting for roughly half of these connections (i.e., ~5% for each). Spatially, different drugs disproportionally affected the high‐order functional networks compared to nondrug PCs with each drug demonstrating unique whole‐brain heat maps. Lastly, significant relationships between the detected FCs and 3‐month behavioral outcomes were detected which, in part, mediated the direct relationships between PDE status and behavioral outcomes, particularly in cognitive and language domains. Compared with conventional GLM models based on raw FC values, the proposed ISV model provided highly correlated estimates but with much improved sensitivity (Figure S3), underscoring the value of directly modeling the enriched intersubject differences to dissect the individual drug and nondrug PC effects. Since results from the ISV model were used as our primary findings, the following discussions on related PC effects and associations should bear this model in mind.

To the best of our knowledge, there are currently no existing reports on the overall impact of PDE on the developing brain's functional system. In this study, we simultaneously modeled all six drugs in the same LME model and estimated their individual effects on ISV of brain FC. Through post hoc summary of individual drug effects, our findings pointed to ~5% of FCs significantly associated with one or more drug type. Importantly, by including both drug and nondrug PCs in the same model, we were also able to compare these drug‐related effects with those of nondrug PCs. To our surprise, the six drug PCs were associated with slightly more connections than the five nondrug PCs combined. Although we expected widespread PDE effects, this finding is still striking and puts the overall effect of PDE on par with the collective impacts of the five critical nondrug PCs, which are known to be associated with early brain development. Specifically, sex, birth weight, and gestational age at birth and scan were documented to be among the most significant predictors of structural brain development in a large normative sample of neonates (Knickmeyer et al., 2017). Particularly, sex has been shown to exert a comparable degree of effects (i.e., 0.51% of connections compared to 1.53% in this study) in children and youth based on similar resting‐state fMRI measures (Satterthwaite et al., 2015). Gestational age at birth was linked to FC of a left‐hemisphere prelanguage region in 30‐week‐old fetuses (Thomason et al., 2013), as well as thalamocortical connections in preterm babies scanned at term age (Rogers et al., 2017). Maternal depression and stress/anxiety have been reported to alter amygdala FC in newborns (Qiu et al., 2015) and structural brain growth in school‐aged children (Lebel et al., 2016). Overall, the detection of significant nondrug PC effects on neonatal brain FC is consistent with previous reports, but perhaps the more important and intriguing finding was the documentation of significant PDE effects that were of similar scale to those of the five nondrug PCs combined. This finding strengthens the ever‐growing concerns of PDE and underscores the potential need to better understand its route of actions with the hope to one day develop effective early intervention strategies.

In descending order, the three legal drugs (SSRIs, nicotine, and alcohol) were most impactful, followed by opioids, marijuana, and cocaine. This finding is in line with the reports of adverse effects and birth defects related to nicotine (e.g., still birth) and alcohol exposures (e.g., fetal alcohol syndrome), which are often reportedly more severe than those from exposures to other drugs (Ross et al., 2015). However, the wide‐spread SSRI effects were less expected. The use of antidepressants medication during pregnancy is thought to be relatively safe (Sujan, Öberg, Quinn, & D'Onofrio, 2019), but the topic is heavily debated (Talati & Weissman, 2019), with some reports showing increased risks for speech and language disorders in childhood (Brown et al., 2016), mood disorders in early adolescence (Brown et al., 2016), and other nonpsychiatric conditions (Gingrich et al., 2017). The topic is further complicated by the fact that maternal depression, the treatment target of SSRI, itself is linked to adverse neurodevelopmental outcomes (Lebel et al., 2016; Posner et al., 2016; Qiu et al., 2015; Rifkin‐Graboi et al., 2013; Wonch et al., 2016). Given this consideration, we included maternal depression scales in the same model and showed that prenatal SSRI exposure was related to neonatal FC after covarying for maternal depression, supporting independent effects of SSRI beyond the related maternal trait. More importantly, our behavioral analysis further showed that SSRI was the only drug that showed significant associations between drug status and all three 3‐month behavioral outcomes (Table 2), paralleled by the highest percentage of overlap between PC–FC and FC–behavioral relationships (Figure S6). If validated and further linked to adverse long‐term outcomes, these findings could add another consideration to the optimal ways of treatment of depression in pregnant women.

Among the three illicit drugs, opioids had the highest proportion of effects (0.91% of connections) followed by marijuana (0.81%) and cocaine (0.60%). Consistent with this finding, illicit opiate use during pregnancy is linked to multiple problems including premature labor, preeclampsia, intrauterine growth retardation and/or death (Bashore, Ketchum, Staisch, Barrett, & Zimmermann, 1981; Hulse, Milne, English, & Holman, 1998; Kaltenbach, Berghella, Finnegan, & Finnegan, 1998). Furthermore, opiate exposed neonates are typically lower in birthweight, have smaller head circumference, and more likely experience neonatal abstinence syndrome compared to other illicit drugs (Binder & Vavrinkova, 2008). Together, these findings represent an urgent call for more studies into the underlying mechanisms and potential remedies for prenatal opioid exposure. In line with this recognition, there has been a recent request for application from National Institutes of Health (NIH) calling for a nation‐wide multicenter study of early brain development with an emphasis on the delineation of effects related to prenatal opioid exposure (the Healthy Brain and Child Development study [HEALthy BCD]).

Besides degree, the spatial distribution of individual drug effects could shed additional light on their potential action pathways to behavioral outcomes. To this aim, we derived drug‐specific heat maps quantifying the number of affected connections associated with each region across the whole brain for each drug type. As hypothesized, drug‐related heat maps exhibited a common pattern consisting of distinct foci in high‐order brain areas which is highly in line with the general theme of high‐order functional deficits associated with PDE (Ross et al., 2015). However, the individual drug‐specific heat maps may reveal more intriguing brain bases of PDE effects on behavior. For example, the nicotine heat map featured an exaggerated concentration in bilateral medial and right lateral prefrontal regions, which is consistent with the reported links between prenatal nicotine exposure and higher likelihood for attention deficit hyperactivity disorder (Dong et al., 2018; Schmitz et al., 2006), externalizing behaviors (Tiesler & Heinrich, 2014), and executive function (Piper & Corbett, 2012) deficits. Similarly, a hallmark behavioral deficit associated with prenatal cocaine exposure relates to arousal and emotional dysregulation (Dennis, Bendersky, Ramsay, & Lewis, 2006; Li et al., 2009), which was again consistent with the cocaine heat map featuring focal enhancement of effects within bilateral anterior cingulate and left middle frontal areas regions known to be involved in such regulations (Etkin, Egner, & Kalisch, 2011; Wager, Davidson, Hughes, Lindquist, & O'chsner, 2009). For opioids, the right angular gyrus and left middle frontal gyrus stood out as two brain areas with the highest concentration of effects. As both areas are multimodal association areas involved in a range of functions including language, number processing, spatial attention, theory of mind for the angular gyrus (Seghier, 2013) and language, executive control for the left middle frontal gyrus (Sierpowska et al., 2018), this finding agrees with previous reports of nonspecific disruptions of attention, cognitive, and regulatory functions associated with prenatal opioids exposure (Nygaard, Slinning, Moe, & Walhovd, 2016). Exposures to alcohol and marijuana prenatally have been both associated with a range of attention, memory, and executive control problems (Coles, 2011; Fried & Smith, 2001; Goldschmidt, Day, & Richardson, 2000), which is in part consistent with our observed more scattered distribution of effects in their corresponding heat maps. Overall, the observed drug‐specific heat maps agree well with the most commonly reported behavioral deficits associated with the corresponding drugs. Although detailed analysis of each drug‐specific heat map with respect to the corresponding drug‐related cellular/molecular action pathways is beyond the scope of this study, the current set of “signature” heat maps for different drugs could serve as a starting point to dig into the relationships between the “hot” drug‐related brain regions and corresponding long‐term behavioral outcomes, with the hope to one day derive imaging‐based biomarkers for the prediction of behavioral outcomes in this at‐risk population.

Indeed, the detected FCs associated with PDE show significant correlations with 3‐month behavioral outcomes. Although the percentage of connections showing such behavioral correlations were relatively small (73 of 2,553 or 2.86%; 42 drug, 26 nondrug, and 5 combined), they combined to explain a substantial amount of variance for the 3‐month behavioral outcomes: 42/28% (cognitive), 50/29% (language), 31/34% (motor) for drug/nondrug PC‐related FCs, respectively. These numbers reflect a bias of effects from drug‐related FCs on cognitive and language performances compared to motor scores, which is also 1.5 ~ 1.7 times higher than those from nondrug‐related connections. This is again consistent with the preferential association of drug effects with high‐order functional networks. When directly assessing the potential mediating role of these connections on drug–behavior relationships, this bias was again evident; the five drug–cognitive/language relationships all became insignificant while the two motor‐related relationships remained significant suggesting limited mediation (Table 2). Taken together, through a three‐step mediation analysis, this study provides strong evidence supporting an action pathway of PDE on later developmental outcomes through the developing brain's functional connectivity.

There are several limitations in this study. First, drug exposure in this sample is complex and the reported individual drug effects are based on a statistical model which includes all drug and nondrug variables simultaneously to control for their effects. Future studies with a large‐enough sample that encompasses enough single‐drug users (e.g., the NIH HEALthy BCD study) may be needed to validate the current drug‐specific findings. However, since polydrug use is the norm, rather than the exception, inclusion of infants with prenatal exposure to more than one drug is a more naturally valid representation of the exposed population. Moreover, the drug information was largely qualitative (i.e., 1 or 0) limiting our capability for detecting dose‐dependent effects. The lack of longitudinal follow‐up beyond 3 months of age represents another limitation of the current study. Given the very dynamic nature of postnatal brain and behavioral development (Gao, Lin, Grewen, & Gilmore, 2016; Gilmore et al., 2018), there are likely age‐dependent changes in the reported degree and distribution of PDE effects throughout infancy and beyond. A third limitation relates to the potential interaction effects between drugs and nondrug PCs as well as between different drugs. Given the multitude of PCs in this study, it was not feasible to assess all pair‐wise interactions in our model. However, we did perform post hoc interaction analyses between sex and PDE status (i.e., 0 for no PDE and 1 for any PDE) and our results (Figure S7) showed dominant negative interactions indicating differential PDE effects in male and female infants, which is consistent with previous reports of sex‐dependent drug effects (Dow‐Edwards et al., 2014; Terasaki, Gomez, & Schwarz, 2016). Regarding polydrug exposure, we further coded drug exposure differences in a semi‐continuous manner (i.e., from 0 to 5 corresponding to the number of drug exposure differences) and repeated our ISV‐based model with FC as the dependent variable and semi‐continuous drug exposure differences as the independent variable. Our results (Figure S8) demonstrated dominant positive correlations suggesting larger differences in the number of drug exposures were associated with greater brain FC differences. Although rudimentary, this trend supports the idea that more severe effects can be associated with polydrug compared to single drug exposures. Limited sample measures, particularly in the area of socioeconomic status (SES), which are known to affect brain/behavioral development (Gao et al., 2014) and may interact with PDE during this process, represents a fourth limitation. With this in mind, we repeated the main analysis in the subset of neonates (n = 110) whose MEDU information was known. Here, two models with and without MEDU were examined, thus allowing us to evaluate the robustness of our primary results against the potential effects of MEDU. Our results showed a high degree of consistency in the beta estimates of the 2,724 reported significant effects (nearly perfect correlation with r ~ .99, Figure S9a) between the two models, indicating there were likely minimal effects of MEDU in our primary findings. However, there were indeed a small number of connections showing significant MEDU effects within the reported effects (Figure S9b), suggesting that SES effects are in fact detectable, thus future studies with larger sample sizes involving this and other related measure (e.g., income, family environment, etc.) are needed to further elucidate the potential effects of SES factors. Finally, two scanners were used for image acquisition and we used post hoc regression and covariate modeling to control for scanner‐related effects. Future studies based on single scanner acquisitions could provide independent validation of the current findings.

In conclusion, findings from this study revealed three major features of PDE effects on overall infant brain and behavioral development. First, PDE affected ~5% of whole‐brain FCs, which was comparable to the combined effects of five critical nondrug PC variables (sex, gestational age at birth/scan, birth weight, and maternal depression scales). Interestingly, the number of connections affected by licit drugs (i.e., SSRI, nicotine, and alcohol) were greater than illicit drugs (i.e., opioids, marijuana, and cocaine). Second, the distribution of PDE effects concentrated on high‐order functional networks rather than primary ones, which is consistent with the numerous reports of high‐order functional disruptions associated with PDE (Ross et al., 2015). Beyond this shared pattern, our heat map analysis revealed drug‐specific heat maps, which could potentially link drug‐specific action pathways to corresponding developmental outcomes. Indeed, the detected FCs associated with PDE showed significant correlations with 3‐month behavioral outcomes and collectively accounted for more variance of cognitive and language outcomes compared to motor outcomes. Consistently, our results support a significant mediation role of these brain FCs on the relationship between PDE and 3‐month cognitive/language outcomes. Overall, these findings provide the first overall assessment of the impacts of PDE on neonatal functional brain and behavioral development. If validated, these findings could provide critical guidance on future studies of the brain basis of PDE effects and potentially inform early intervention designs.

APPENDIX S1. Supporting information

TABLE S1. Maternal body‐mass‐index (MBMI) and education (MEDU) summary measures

TABLE S2. functional parcellation, AAL, and network‐level affiliations.

TABLE S3. Prenatal drug exposure (PDE) profiles

TABLE S4. Scanner and motion related summary measures

TABLE S5. Summary of network‐level effects

FIGURE S1. Functional parcellation network‐level correspondence

FIGURE S2. Scanner and motion effects

FIGURE S3. Comparison of the estimated effects of drug and non‐drug PCs on newborn brain functional connectivity based on the proposed ISV model and conventional GLM model using raw functional connectivity values

FIGURE S4. Spatial distribution of non‐drug PC effects on brain functional connections

FIGURE S6. Correspondence between functional connectivity (FC) ~ BEHAVIOR and FC ~ participant characteristics (PC) effects

FIGURE S7. Prenatal drug exposure (PDE)‐by‐sex interactions

FIGURE S8. Poly‐drug effects

FIGURE S9. Maternal education (MEDU) impact

Click here for additional data file.