Network based statistics reveals trophic and neuroprotective effect of early high dose erythropoetin on brain connectivity in very preterm infants.

Periventricular white matter injury is common in very preterm infants and it is associated with long term neurodevelopmental impairments. While evidence supports the protective effects of erythropoetin (EPO) in preventing injury, we currently lack the complete understanding of how EPO affects the emergence and maturation of anatomical brain connectivity and function. In this case-control study, connectomic analysis based on diffusion MRI tractography was applied to evaluate the effect of early high-dose EPO in preterm infants. A whole brain, network-level analysis revealed a sub-network of anatomical brain connections in which connectivity strengths were significantly stronger in the EPO group. This distributed network comprised connections predominantly in the frontal and temporal lobe bilaterally, and the effect of EPO was focused on peripheral and feeder connections of the core structural connectivity network. EPO resulted in a globally increased clustering coefficient, higher global and average local efficiency, while higher strength and increased clustering was found for regions in the frontal lobe and cingulate gyrus. The connectivity network most affected by the EPO treatment showed a steeper increase graph theoretical measures with age compared to the placebo group. Our results demonstrate a weak but widespread effect of EPO on the structural connectivity network and a possible trophic effect of EPO reflected by increasing network segregation, predominantly in local connections.

RCT Entities:   Population Interventions Outcomes

Introduction

While the mortality rate of very preterm infants has decreased significantly within the last decade (Ruegger et al., 2014), long-term neurodevelopmental problems such as cognitive, motor and behavior impairments remain of concern (Latal, 2009; Saigal and Doyle, 2008). Periventricular white matter injury (PWMI) is the predominant form of perinatal brain injury in preterm infants PWMI is characterized by marked astrogliosis and microgliosis, and initially by a decrease in premyelinating oligodendrocytes (pre-OLs) (Haynes et al., 2003). PWMI is accompanied by neuronal/axonal deficits in the cerebral white matter, thalamus, basal ganglia, cerebral cortex, brainstem and cerebellum. This led to the term „encephalopathy of prematurity“, which describes a complex amalgam of primary destructive disease and secondary maturational and trophic disturbances (Volpe, 2009). Encephalopathy of prematurity has been linked to later neurodevelopmental impairments such as cerebral palsy, motor dysfunction, cognitive and behavioral problems as well as deficits in executive functions (Woodward et al., 2006; Woodward et al., 2011). The high incidence of neurodevelopmental impairments in preterm infants is reflected by the ongoing search for neuroprotective interventions that can prevent brain injury or enhance repair of the immature brain, with the ultimate goal of improving long-term motor and cognitive outcome (Gonzalez et al., 2009). Among several pharmacological candidates to prevent brain injury or improve its development in preterm infants, erythropoietin (EPO) has been shown to be among the most promising neuroprotective agents (Robertson et al., 2012). Protective effects of EPO important for reducing acute injury include decreased apoptosis, inflammation, excitotoxicity and glutamate toxicity (Juul, 2012; Shingo et al., 2001). Furthermore, EPO was shown to play a role in developmental mechanisms, as it stimulates neurogenesis (i.e. proliferation and differentiation of pre-oligodendrocytes (Sugawa et al., 2002)), angiogenesis and migration of regenerating neurons (Tsai et al., 2006). In 2012, recruitment of the first randomized, double-blind placebo-controlled, prospective multicentre trial (NCT00413946) in very preterm infants using high-dose EPO as a potential neuroprotective agent was completed (n = 450 preterm infants). In a subgroup of infants, MR imaging showed less PWMI on conventional MRI (Leuchter et al., 2014), and improved white matter microstructure assessed by tract based spatial statistics (O'Gorman et al., 2015).

Given the widespread brain changes observed in these prior studies and the clear relationship of global brain networks and brain function (van den Heuvel and Sporns, 2013), it is of interest to examine the global brain network connectivity in these infants. In recent years, due to improvements of MR imaging such as in diffusion MRI, it has become possible to map the global brain connectivity enabling us to study the global network organizational properties of the brain, the effect of neuroprotective interventions and their relation to brain development and pathologies (Hagmann et al., 2012). The complete structural or functional organization of the brain is represented by its connectome – the brain-network level architecture of all neuronal connections (Sporns, 2011). Numerous studies using dMRI proved that it is possible to map the developing connectome of infants with a granularity that is comparable to adult studies (Brown et al., 2014; Cao et al., 2017; Jakab et al., 2015; Tymofiyeva et al., 2012; van den Heuvel et al., 2015). We are only beginning to understand if neurodevelopment in congenital or acquired diseases can lead to a global impairment of the interactions between the neural elements, and whether this is reflected by impairments in the macro-scale connectome. A recent study using connectome analysis has demonstrated specific alterations in brain topology and structural organization in the brain network organization of high-risk preterm born children at school age and such altered regional brain connectivity was related to specific neurocognitive deficits in children born preterm (Fischi-Gomez et al., 2015; Fischi-Gomez et al., 2016). Hence, connectomic analysis might not only help to further understand the neural correlate of neurodevelopmental impairments in preterm infants but also to help understand the effects of EPO on global brain connectivity development. Our study is built on two hypotheses. First, we hypothesize that the neuroprotective effect manifests as changes in the brain's structural connections, which may only become apparent on the whole-brain network level and not locally. Connectomic analysis may therefore shed light on the effects of therapeutic agents not only on isolated brain regions, but also on how they affect the topography of brain connections (Bonilha et al., 2015; Crossley et al., 2017; Zeng et al., 2015). Second but most importantly, EPO treatment may steer the ongoing structural connectivity development into an alternative, more optimal trajectory, which would in turn be observed as global changes in the structural connectivity architecture compared to untreated subjects.

Materials and methods

Ethical approval was granted by the local ethical committee (KEK StV-36/04), and the study was approved by the Swiss drug surveillance unit (Swissmedic, 2005DR3179). The trial was registered at ClinicalTrials.gov (number NCT00413946).

Patient population

The inclusion criteria and selection process for the study have been described in detail previously (Leuchter et al., 2014; Natalucci et al., 2016; O'Gorman et al., 2015). In short, the preterm infants in this study represent a subgroup of infants enrolled in the randomized, double-blind placebo-controlled, prospective multicentre study “Does erythropoietin improve outcome in preterm infants” (NCT00413946). For this analysis, only infants with DTI performed in Zurich (n = 140) were included. Clinical co-morbidities such as sepsis (blood culture proven), NEC, PDA, and CLD (defined as oxygen requirement at corrected 36 weeks of gestation) were noted. Our database comprised the socio-economic status (SES) and corrected gestational age at the time of MRI (CGA).

Randomization, neuroprotective intervention and blinding

Two year neurodevelopmental outcomes and the treatment protocol have been has been previously published (Fauchere et al., 2008; Natalucci et al., 2016). Study medication was randomly assigned to each patient number in a 1:1 allocation using a computer-based random-number generator. Erythropoetin or an equivalent volume of normal saline (NaCl 0.9%) placebo was administered intravenously before 3 h of age after birth, at 12–18 and at 36–42 h after birth. A single dose consists of 25 μg (3000 IU) human erythropoietin per kilogram body weight dissolved in 1 ml sterile water. Hospital staff, parents, the outcome assessors and data analysts were kept blinded to the allocation.

Cerebral MRI protocol

Cerebral MRI was performed with a 3.0 T GE scanner (GE Medical Systems), using an 8-channel receive-only head coil. All infants were scanned under natural sleep using a vacuum mattress. Ear plugs (attenuation: 24 dB; Earsoft; Aearo) and Minimuffs (attenuation: 7 dB; Natus) were applied for noise protection. Oxygen saturation was monitored during scanning, and a neonatologist and a neonatal nurse were present during the MRI investigation.

We acquired T1-weighted MR images with a 3D fast spoiled gradient echo sequence (echo time = 2.6 ms, repetition time = 5.7 ms, inversion time = 750 ms, flip angle = 12°, voxel resolution 0.7 × 0.7 × 1.4 mm3), and T2-weighted images with a fast recovery fast spin echo sequence (echo time = 126 ms, repetition time = 6600 ms, resolution 0.7 × 0.7 × 1.5 mm3).

DTI was acquired using a pulsed gradient spin echo echo-planar imaging sequence with echo time = 77 ms, repetition time = 9 s, field of view = 18 cm, matrix = 128 × 128, slice thickness = 3 mm. The diffusion encoding scheme included 21 non-collinear gradient encoding directions with b = 1000 and four interleaved b = 0 images. Of the 140 infants scanned in Zurich, in 78 infants DTI data demonstrated motion artifacts (n = 78), 4 infants had with cystic lesions (n = 4) apparent on structural MRI and hence were excluded from further analysis, resulting in a final group size of 58 infants, of whom 24 were treated with EPO and 34 were treated with placebo.

Diffusion tensor image processing

DTI data were visually controlled for artifacts by observing each diffusion gradient volume and each image slice for signal dephasing or spin history artifacts caused by movements. If the newborn woke up or moved excessively during the DTI scan, data were discarded and the subject was excluded. Subjects were also excluded if >3 diffusion-weighting gradient volumes were corrupted by motion artifacts. Image frames and the corresponding entries in the b-matrix and b-value descriptor files were removed from further analysis if head movement caused extensive signal dropout throughout the brain in the given frame. A custom script (Digital Supplement) written in bash programming language for Linux was used to process the neonatal DTI images. Spurious image shifts originating from eddy currents were corrected with the eddy_correct command and the b-matrix was reoriented with the fdt_rotate_bvecs script in the FSL software. Next, image intensity shifts along the supero-inferior direction were corrected by normalizing the intensity of each axial slice to the corresponding slice in the motion-corrected, median diffusion-weighted image. Rician noise filtering of DTI was performed using a command line tool in Slicer 3D (Aja-Fernandez et al., 2009). FA and MD maps were calculated by running the dtifit command in FSL on the Rician noise filtered data. An overview of the image processing steps is provided in Fig. 1.

Fig. 1

Methodological overview. (a) image post processing and connectivity analysis steps, (b) thresholding procedure based on network density (cost).

Region of interest system

Anatomical parcellation of the brains was carried out by using the neonatal regions of interest (ROI) from UNC Infant atlas UNC (Shi et al., 2011). This comprised 90 cortical and subcortical areas that were propagated from the Automated Anatomical Labeling (AAL: (Tzourio-Mazoyer et al., 2002)) to the images of 33.4–42.1 gestational week neonates (for ROI nomenclature, see Table S1). We first created a custom fractional anisotropy template in the UNC space, which was used to transfer the UNC-AAL ROIs to the DTI space of each subject.

This FA template was constructed in three steps based on DTI data of a separate study population, consisting of 40 normally developing and term-born newborns imaged within 6 weeks after birth (CGA of the subjects: 42.5 ± 1.9 (39–48.7) weeks, 20 females and 20 males). First, we registered the B0 images of each neonate with the T2-weighted template of the UNC atlas using an affine registration followed by an elastic non-linear registration implemented in the NIFTIreg software (reg_f3d command, control grid size: 9 * 9 * 9 mm, weight of the bending energy penalty term: 0.05, gradient smoothing with a kernel of 4 mm). The transformations arising from this registration step were used to re-sample the Rician filtered FA images of each neonate to UNC template. An initial FA template was then created by averaging these FA images. Each subject's FA maps were co-registered to the initial FA template with NIFTIreg with an elastic deformation allowing for finer distortions (reg_f3d command, control grid size: 5 * 5 * 5 mm, weight of the bending energy penalty term: 0.005, gradient smoothing with a kernel of 4 mm). The FA template representative of the 40. week of development is available as digital supplement.

Structural connectivity network construction

Seed points for probabilistic diffusion tractography (Probabilistic index of connectivity method (Parker et al., 2003) were defined as the voxels with FA ≥ 0.1 within a whole-brain mask. Orientation density functions (ODFs) were estimated using fourth order Spherical Harmonics and a maximum of two local ODF maxima were set to be detected at each voxel and probability density function (PDF) profile was produced from the local ODF maxima. Fiber tracking was carried out on the voxel-wise PDF profile with the Euler interpolation method using 10 iterations per each seed point. Tracing stopped at any voxel whose FA was <0.2. The tractography script is available as digital supplement.

Two kinds of weighted, undirected structural connectivity networks were formed based on whole-brain probabilistic tractography: FAmean and SC. In both networks, nodes corresponded to the AAL ROIs in subject space. In the FAmean structural connectivity network, the strength of each edge was given by the mean fractional anisotropy (FA) value of all tractography streamlines connecting the regions at the two end-points.

Secondly, network edges were defined using the normalized structural connectivity strength, SC. Tractography streamline counts connecting any two nodes were normalized for the bias arising from the volumetric differences between the regions of interests corresponding to the nodes, and for the linear bias that arises from more distant brain regions showing more streamlines (Hagmann et al., 2008).where SC is the structural connectivity strength, D is the Euclidean distance between the node centre-points, S is the number of tractography streamlines connecting nodes i and j, V and V are the node volumes.

Graph theoretical analysis

Structural brain network analysis comprised connection (edge) selection with an iterative thresholding procedure. To avoid systematic differences in absolute number of edges, we used a thresholding procedure based on network costs, which ensures similar network topology across subjects by keeping the same number of strongest connections. First, we calculated the network cost for each subject without thresholding (actual network edges / all possible edges). We then selected the lowest possible network cost across the study cohort and used this as an upper value of the proportional threshold procedure (van den Heuvel et al., 2017), defined as Tmax. The structural connectivity (SC) networks were thresholded from a density of 1% to Tmax with the threshold_proportional command in the Brain Connectivity Toolbox (Rubinov and Sporns, 2010), and the same surviving edges were kept for both the FAmean and SC networks. For both networks at each cost-thresholded level, the following graph theoretical parameters were calculated using the Brain Connectivity Toolbox for Matlab.

Strength of region i:where wij is the connection weight between i and j, N is the entire network.

Betweenness centrality of region i (Freeman, 1979):where ρ is the number of shortest paths between h and j, and ρ(i) is the number of shortest paths between h and j that pass through i. N are all network nodes. h and j are the neighboring nodes of i.

Global efficiency of the weighted network N:where n is the total number of nodes in the network N, d is the shortest weighted path length between nodes i and j.

Weighted local efficiency is the efficiency calculated for a sub-network constructed from the 1st neighbors of the node i:where w, w are connection weights between nodes i,j and i,h, d is the shortest weighted path length between nodes i and h, k is the degree of node i.

Average local efficiency is the average of all local efficiencies across the network nodes.

Weighted clustering coefficient of the network (Onnela et al., 2005; Watts and Strogatz, 1998):where t is the weighted geometric mean of triangles around node i:

Core and periphery structure

The connectivity networks were subdivided into two non-overlapping groups of nodes, a core and a periphery based on the Kernighan-Lin algorithm for graph partitioning (Borgatti and Everett, 2000; Newman, 2006), implemented as the “core_periphery_dir” command in the Brain Connectivity Toolbox for Matlab. A common core for the study population was defined as the nodes that were in the core partition in >90% of the subjects. Core connections connected the core nodes, their first degree connections were defined as feeders, while the remaining non-zero connections were labeled as peripheral connections. For each subject, the coreness parameter was calculated that reflected the degree of separation between the core and peripheral nodes (Borgatti and Everett, 2000; Newman, 2006).

Statistical analysis

To test the effects of EPO treatment and neurodevelopment on the global and nodal graph theoretical parameters, a multivariate analysis of covariance in Matlab R2014 using the mancovan toolbox (Mathworks Inc., Mattick, USA) was performed. We included demeaned corrected gestational age as a confounder in the multivariate model. During the mass-multivariate testing of graph theoretical parameters across network nodes (brain regions), statistical significance was adjusted for multiple comparisons using the Benjamini-Hochberg procedure, implemented in the FDR toolbox in Matlab R2014. As networks were tested at various cost thresholds, only those results were accepted that survived FDR adjustment and consistently showed the same effect direction (positive or negative group difference) in 95% of all the analyzed network costs.

Hypothesis tests on structural connectivity networks were based on the network-based statistics (NBS) method described by Zalesky et al. (2010), which was implemented in the NBS Toolbox for Matlab R2014. NBS avoids the problem of multiple comparisons during mass univariate tests on connectivity networks by estimating statistical significance for subsets of mutually connected network nodes in topological rather than physical space (Reess et al., 2016). The NBS method comprises four steps. First, the t-statistic for each individual edge in the connectivity network is calculated. Second, a primary component-forming threshold (p < .05, uncorrected) is applied to identify edges displaying differences in connectivity strength. Third, subthreshold edges are assessed for mutual connections forming clusters in topological space that may point towards the existence of non-chance clusters. Fourth, a test with 5000 random permutations is applied to compute statistical significance for all previously identified network component. As NBS uses permutation test to build up the sample distribution, it can be applied to smaller study groups without assuming normality. The final hypothesis test is then carried out for the empirically determined components by comparing their extent with the proportion of permutations yielding a component with equal or greater size, correcting for the family-wise error rate at cluster level with p < .05. In our tests, a primary threshold t-score of 2 was chosen. Next, a post hoc analysis was performed using linear regression in SPSS V22 (IBM, Armonk, New York) to select any further demographic or clinical parameters that may have an effect on the structural connectivity values averaged over the networks from the NBS analysis.

We presented the results of the NBS as three-dimensional graph visualizations, which represented below p-threshold connection pairs surviving multiple comparison correction. The brain networks were visualized with the BrainNet Viewer for Matlab R2014 (Xia et al., 2013).

Results

Patients

58 preterm infants with mean (SD) corrected gestational age at birth 29.75 (1.44) weeks and at scanning of 41.1 (2.09) weeks) were included into the connectomic analysis. Table 1 shows the demographic and clinical details of the two groups. No significant differences were found in CGA at scanning or in any other neonatal or clinical details between the two treatment groups. Subjects with cystic lesions were not included in the current study. Based on the white matter injury scoring system of Woodward et al. (2006), 54 (93.2%) neonates displayed normal white matter and only 4 (6.8%) had abnormal white matter.

Table 1

Demographic and clinical parameters of the infants. Differences of clinical characteristics between groups were evaluated with the chi-square test or the Fisher exact-test as appropriate for the categorical variables, and with Student's t-test or Mann-Whitney test as appropriate for the continuous variables. SD – standard deviation, IQR – inter-quartile range (25th–75th percentiles).

Characteristic	Reported results	EPO group (N = 24)	Placebo group (N = 34)	P-value
Gestational age (weeks)	Mean ± SD	30.17 ± 1.44	29.5 ± 1.44	0.0952
Birth weight (g)	Mean ± SD	1337 ± 332	1192 ± 10	0.0982
Z-score of birth weight (SD)	Mean ± SD	−0.27 ± 0.71	−0.19 ± 0.95	0.722
Gestational age at MRI (weeks)	Mean ± SD	40.93 (2.09)	41.33 (2.12)	0.492
Gender (female)	N (%)	8 (33%)	14 (41%)	0.611
Retinopathy of prematurity	N (%)	3 (13%)	4 (12%)	0.771
Sepsis+	N (%)	3 (13%)	5 (15%)	0.881
NEC	N (%)	0 (0%)	1 (3%)	0.581
Mechanical ventilation	N (%)	8 (33%)	13 (38%)	0.461
Duration of mechanical ventilation (days)	Median (range)	0 (0; 6)	0 (2; 12)	0.733
Bronchopulmonary dysplasia⁎	N (%)	2 (8%)	5 (14%)	0.381
Postnatal steroids	N (%)	2 (8%)	0 (0%)	0.171

Chi-square test or Fisher exact test as appropriate.

Student test.

Mann-Whitney test.

Sepsis proven by positive blood cultures; NEC, necrotizing enterocolitis.

Oxygen requirement at corrected age of 36 weeks.

Quality of neonatal DTI data

In the EPO group, corrupted gradient volumes were excluded from the analysis in two subjects (1 and 2 volumes removed). 1–3 volumes (1.57 ± 0.8) were removed for 7 subjects in the placebo group. Network level analysis revealed that the number of DTI frames removed did not correlate with the FAmean or SC values (highest significance, NBS, p = .422 and p = .131, respectively).

Global graph theory characteristics

The network cost of the non-thresholded connectivity matrices was 0.30 ± 0.022 (range: 0.22–0.35). We therefore used 0.22 as the maximum value during the cost thresholding procedure.

Global efficiency, average local efficiency, average clustering coefficient and average nodal strength were depicted as a function of network density separately for the EPO and placebo groups in Fig. 2/a,c,e,g. EPO treated newborns showed 7.75% higher average clustering coefficient, with maximum effect observed at a network density of 0.21 (Fig. 2/e, T = 16.71, p < .0001, estimated marginal means (EMM) EPO: 0.139 (CI, 95%: 0.135–0.142), Placebo: 0.129 (CI, 95%: 0.127–0.132), appearing at CGA = 41.2 weeks). Average clustering coefficient increased with CGA (Fig. 2/f), and the steepness of this was moderately larger after the 42. week in the EPO group (Fig. 2/f). Global network efficiency was 6.41% higher in the EPO group, with maximum effect observed at a network density of 0.05 (Fig. 2/a, T = 7.92, p = .0067, EMM, EPO: 0.083 (CI, 95%: 0.08–0.086), Placebo: 0.078 (CI, 95%: 0.075–0.08), appearing at CGA = 41.2 weeks). Average nodal strength was not affected by EPO treatment in the density range that we used to define statistical significance (Fig. 2/g,h). Average local efficiency was 5.65% higher in the EPO group, maximum effect was observed at a network density of 0.13 (Fig. 2/c, T = 18.16, p < .0001, EMM, EPO: 0.187 (CI, 95%: 0.178–0.196), Placebo: 0.177 (CI, 95%: 0.170–0.185), appearing at CGA = 41.2 weeks), the age dependency showed similar tendency to the clustering coefficient (Fig. 2/d).

Fig. 2

Global graph theory characteristics of the EPO and placebo treated groups. (a) global efficiency, (c) average local efficiency, (e) average clustering coefficient, (g) average strength. Global characteristics are plotted on the left axis as a function of network cost. The right axis shows the T statistic of the multivariate analysis of group differences. Network costs at which significant EPO-placebo differences appear are marked with segmented line and dot (p < .05, FDR-corrected across network nodes). Blue lines: mean values, placebo group, red lines: mean values, EPO group, transparent area reflects the variability of measurements (±1SD range). Network characteristics that were significantly different across group are plotted as a function of CGA (b, d, f, h).

Local graph theory characteristics

Local graph theory characteristics of three nodes in the FAmean network were significantly different between the EPO and placebo groups, the results of the nodal analysis are summarized in Table 2 and in Fig. 3. Graph theoretical measures of two frontal brain regions were significantly different between the EPO-treated and placebo groups. In the FAmean network, the left middle frontal gyrus showed 16.8% higher clustering coefficient in the EPO-treated neonates (p = .014, EMM, EPO: 0.118 (CI, 95%: 0.110–0.125), placebo: 0.101 (CI, 95%: 0.095–0.107), while the right superior frontal gyrus had 13.5% higher clustering coefficient (p < .001, EMM, EPO: 0.118 (CI, 95%: 0.112–0.124), placebo: 0.104 (CI, 95%: 0.098–0.109). At a network density of 0.21, the strength of the right anterior cingulate gyrus was 20.1% lower in the EPO treated group (p < .001, EMM, EPO: 6.35 (CI, 95%: 5.93–6.78), placebo: 5.29 (CI, 95%: 4.93–5.65). In the SC network, the clustering coefficient of the right middle frontal gyrus was 37.7% higher in EPO treated newborns (p < .001, EMM, EPO: 4.86*105 (CI, 95%: 4.31*105–5.41*105)), placebo: 3.53*105 (CI, 95%: 3.07*105–3.99*105).

Table 2

Effect of EPO treatment – local graph theory characteristics.

Dependent variable, AAL region number	Network	Test statistic (T)1	Effect size (B)1	Significance (p)1	Network cost2	Role of node network?
Clustering coefficient, left middle frontal gyrus (7)	FAmean	11.33	0.101	0.014	0.22	Peripheral
Clustering coefficient, right superior frontal gyrus (4)	FAmean	12.56	0.104	<0.001	0.22	Peripheral
Strength, right anterior cingulate gyrus (32)	FAmean	14.49	97.1	<0.001	0.21	Peripheral
Clustering coefficient, right middle frontal gyrus (8)	SC	13.76	3.5*105	<0.001	0.05	Peripheral

Model corrected for CGA. AAL = automated anatomical labeling atlas.

Network density at which maximum effect size is observed.

Fig. 3

Effect of EPO treatment on nodal graph theory characteristics. (a) nodal clustering coefficient of the right superior frontal gyrus, (b) nodal clustering coefficient of the left middle frontal gyrus, (c) nodal strength of the right anterior cingulate gyrus, (d) nodal clustering coefficient of the right middle frontal gyrus, SC-network. In each sub-plot, nodal characteristics are plotted (left axis) as a function of network density, while the right axis shows the T statistic of the multivariate analysis of group differences. Network costs at which significant EPO-placebo differences appear are marked with segmented line and dot (p < .05, FDR-corrected across network nodes). Blue lines: mean values, placebo group, red lines: mean values, EPO group, transparent area reflects the variability of measurements (±1SD range). The right panels in each sub-plot shows the CGA dependency of the main effect.

Effect of gestational age on the structural connectivity networks

Structural connectivity strength increased with the corrected gestational age (CGA) of the neonates. FAmean values were found to be significantly positively correlated with CGA in 455 edges of the non-thresholded structural connectivity network (Fig. S1). SC increased with CGA in 322 network edges of the structural connectivity network (Fig. S1). These results prompted us to correct the following statistical analysis steps for the effect of CGA.

Network based statistics and core structure

We carried out the core analysis and network based statistics at a cost of 0.21. This threshold was selected based on results of the global and local graph theoretical analyses, where FAmean networks thresholded at 0.21 and 0.22 were found to be significantly different between EPO treated group and controls. We found that 18 (2.14% of 841 non-zero network edges of the common structural connectivity network at cost = 0.21) connections in the FAmean network were significantly stronger in the neonates treated with EPO when controlling for the effect of CGA. The biggest difference was measured for intra-lobar medium range, frontal, limbic and temporal connections (Fig. 4/a). The structural connections within this network were characterized by 7.72% higher FAmean in the EPO group than in the placebo-treated neonates (p = .0047, estimated marginal means (EMM) from a linear regression model controlling for the effect of CGA, EPO group: 0.265 (CI, 95%: 0.255–0.275), Placebo group: 0.246 (CI 95%: 0.238–0.255), appearing at CGA = 41.2 weeks). The network nodes with the highest number of significantly different network edges (>3 significant edges) were the right anterior and middle cingulate gyri, left middle frontal gyrus, right superior frontal gyrus and right putamen (Fig. 4/a).

Fig. 4

Effect of EPO treatment: network based statistics and core structure analysis. (a) structural connectivity network (FAmean) with increased connection strength in EPO-treated neonates. The size of nodes correlates with the number of significantly positive edges arising from the node, while edge width correlates with the test statistic of the NBS analysis. Gray edges mean inter-modular, colored edges intra-modular or intra-lobar connections. LH: left hemisphere, RH: right hemisphere. (b) linear regression plot showing the relationship between average FAmean of the significant graph edges, CGA at the time of MRI and group allocation. (c) network core affected by EPO treatment, NBS results, (d) common network core of the study population. In the middle image, core and peripheral nodes are depicted in red and yellow, respectively, while core connections, feeder connections and local connections are depicted as red, yellow and gray lines. For the abbreviations of ROIs, see Table S1.

Stepwise linear regression analysis did not select any clinical parameters besides CGA (see Table 1 for variables and parameters considered for inclusion in the model) that may have acted as mediators or confounders of this effect.

Next, we reconstructed the core connectivity network in the study population. The core consisted of the 8–8 hemispherically symmetric nodes (Fig. 4/c,d): middle cingulate gyri and posterior cingulate gyrus on the right side, olfactory lobe, caudate, putamen, thalamus, insula, hippocampus and precuneus, which overlaps with the previously reported rich-club structural core (van den Heuvel and Sporns, 2011). We found that the treatment had a predilection towards nodes that were part of the network periphery in the core network analysis based on our preterm population. While the superior frontal gyrus and cingulate has previously been shown to be part of the core structural connectivity network, our core analysis did not classify these as such. Out of the 18 edges affected by EPO treatment, 4 were core connections (22.2% of network, connections of the putamen, caudate, middle cingulate gyrus and precuneus), 5 were feeders (27.8%) while 9 (50%) were peripheral connections. The coreness statistic was not significantly different between EPO and placebo-treated newborns (T = 2.403, p = .1269). The affected core network structure is illustrated in Fig. 4/c.

At a density of 0.21, NBS revealed no differences between the EPO and placebo groups when using the SC connectivity network instead of the FAmean (minimum p = .067).

As the edge strength correlated with CGA across the whole structural connectivity network (Supplementary Fig. S1) and the steepness of the regression line appeared to differ between placebo and EPO treated groups (Fig. 4/b), we tested whether group difference varied as a function of CGA (Group-CGA interaction). A two-group NBS analysis with continuous covariate interaction was performed, which revealed that interaction is not significant in any of the network edges in the FAmean connectivity network (minimum p = .1124). Next, we repeated the statistical test to capture the effects of EPO treatment on the structural connectivity networks, this model was corrected for the Group-CGA interaction term. The effects of EPO on the FAmean connectivity network were seen in an almost identical set of edges to the previously performed tests with covariate interaction (23 edges, no covariate interaction: 18 edges). With covariate interaction, 78.05% of the edges were overlapping (Dice coefficient of significant edges) with the tests not accounting for covariate interaction in the FAmean network. The effects of EPO with and without Group-CGA interaction are demonstrated in Supplementary Fig. S2.

Supplementary Fig. S1

Age related changes of structural connectivity in the EPO-treated and placebo infant groups. (a) Structural connectivity network where corrected gestational age at MRI (CGA) significantly positively correlated with the FAmean of graph edges, (b) structural connectivity network where CGA was significantly correlated with the SC. For the abbreviations of ROIs, see Table S1.

Supplementary Fig. S2

Analysis of the effect of EPO in the presence of group-CGA interaction. The effect of EPO treatment on the FAmean network is demonstrated (a) with Group-CGA interaction term and (b) without covariate interaction term in the model.

Discussion

In this study, we applied connectomic analysis to identify the effect of early high dose EPO administration on the structural brain connectivity network in preterm infants. While we have previously shown that early high dose EPO administration reduces brain injury and improves white matter development in the major white matter tracts (O'Gorman et al., 2015), little is known about the effects of EPO on connectivity networks in the preterm brain.

Our findings reflect an early, weak effect of the EPO treatment on the overall structural brain connectivity of preterm infants. This effect was not limited to any of the large anatomical subdivisions, such as lobes of the brain, but the connectome of EPO treated infants was rather characterized by an increase of local structural connectivity strengths. Global network strength – the sum of connectivity weights of all connections of a node – was not affected by the EPO treatment. The increased average local and global efficiency after EPO treatment reflects the better ability of the connectivity network to transmit information at the local level and a shift towards a more regular network. We examined four local topological properties of the structural brain networks: network integration was measured by nodal strength and efficiency, while segregation was indicated by the clustering coefficient and betweenness centrality (Onnela et al., 2005).

It has been shown that brain hubs form a so-called “rich club” or core structure (van den Heuvel and Sporns, 2011), which is characterized by a tendency for high-degree nodes to be more densely interconnected. This structure is likely present as early as the 30th gestational week (Ball et al., 2014), and during neurodevelopment after very preterm birth, this core is most likely prioritized over local connectivity (Fischi-Gomez et al., 2016; Karolis et al., 2016). In a recent study, a stronger rich-club architecture was found in preterms than in normally developing infants, despite a relative deficiency of white matter resources (Karolis et al., 2016). Our network based statistical analysis revealed that EPO facilitates both the development of peripheral and frontal, temporal, subcortical and limbic core connections. The EPO-induced development of peripheral, i.e., local connectivity, may provide valuable means to counteract the structural connectivity impairments that are specific for preterm birth (Karolis et al., 2016).

Our results reflect a weak, but widespread effect in the structural connectivity network after EPO treatment, which may give rise to more concentrated changes in the role of frontal and limbic regions in the entire network. These regions are presumed to be densely interconnected areas – hubs – early on. Clustering and the formation of hubs is an emergent feature of the human connectome. Clustering coefficient is a measure of brain network segregation and is used to quantitatively describe the brain's ability for specialized processing within such interconnected groups of brain regions (Sporns, 2013). We found increased clustering coefficient in the right superior and middle frontal gyri and left middle frontal gyrus. The strength of the anterior cingulate gyrus, which is part of the preterm core connectivity network, was increased after EPO treatment. An increase in this metric is the direct consequence of a higher prevalence of clustered connectivity around a brain region (Watts and Strogatz, 1998), which may potentially reflect increased anatomical and functional segregation of this area (Rubinov and Sporns, 2010). Results from neuroimaging studies imply that the majority of large-scale pathways are already present by birth (Doria et al., 2010; van den Heuvel et al., 2015), and the newborn connectome features an adult-like small-world modular architecture that could be sensitive to disruptions, as shown, for example, in children after hypoxic ischemic encephalopathy (Tymofiyeva et al., 2012). van den Heuvel el al. found that clustering coefficient averaged over the entire structural connectivity network increased from the 30th gestational week until term equivalent age in preterm infants (van den Heuvel et al., 2015), which is in line with our observations of gradually increasing clustering coefficients. Increased anatomical segregation may support more efficient consolidation of information from different functional networks through long-range connections, and it is safe to assume that during neurotypical development, clustering coefficient of the most important connectivity hubs gradually increase (van den Heuvel and Sporns, 2013). The developmental predilection for white matter injury to occur during prematurity appears to relate to both timing of appearance and regional distribution of susceptible pre-oligoendrocytes (Back et al., 2007; Riddle et al., 2006). The deep white matter myelinates in a caudal-to-rostral, proximal to distal, central to peripheral direction, with sensory pathways myelinating before motor pathways and with deep occipital white matter maturing first and frontal white matter maturing last (Kinney et al., 1988). Hence, the increased clustering coefficients in the frontal white matter in EPO treated infants might also reflect a regional and temporal variability of the EPO effect. Postnatal age has been associated with increasing connectivity in the frontal areas (Pandit et al., 2014). Our results suggest that EPO might enhance the normal development of connections with increasing age in predominantly frontal areas. This is an important finding as previous imaging studies have shown that FA in the frontal white matter lobe was reduced by prematurity (Rose et al., 2008). Such white matter alterations seem to persist into adolescence and adulthood (Vollmer et al., 2017; Yung et al., 2007), and such alterations were correlated with measures of executive function and general cognitive abilities in adolescence and young adulthood (Vollmer et al., 2017).

Based on previous works (Dosenbach et al., 2007), we propose that the newborn and infant brain gradually builds up the ability of its connectivity network to integrate information between its subparts. During this process, the presumed neuroprotective effect of EPO could steer the development into a trajectory that is characterized by a steeper increase of clustering coefficients. The most prominent feature of the developing infant brain in dMRI is the rapid increase of the diffusion anisotropy. We found that the corrected gestational age at MRI correlates with mean FA in many network edges after correcting for multiple comparisons, which is most likely the consequence of the progression of myelination (Berman et al., 2005; Huppi et al., 1998). Fractional anisotropy is known to increase gradually during the last trimester of gestation, and this trajectory is not significantly different in the postnatal development of prematurely born infants (Bonifacio et al., 2010). We found that the connectivity network most affected by the EPO treatment showed a steeper increase compared to the placebo group, however, no significant group-CGA interaction was found. A similar tendency was observed for the clustering coefficient, efficiency and strength during the nodal graph theory analysis. The increasing FA is considered to reflect an increase in white matter structure and thus an increase in the efficacy of long-distance neuronal signal transmission in the neonatal brain (van den Heuvel et al., 2015). The described steeper increase in FA and in clustering coefficient might reflect the long-term effects of EPO on neurogenesis (Juul, 2010). Based on these findings, we may also speculate that a possible trophic effect of EPO manifests as increasing network integration, especially at more advanced CGA.

A recent meta-analysis including 6163 children born very preterm and 5471 term-born children showed that children born very preterm have lower IQs, lower scores on measures of executive functioning and lower scores on measures of processing speed (Brydges et al., 2018). Recent neuroimaging and lesion studies suggest that executive functions depend on distributed networks encompassing both frontal and posterior (mainly parietal) associative cortices (Collette et al., 2006; Jurado and Rosselli, 2007). This is consistent with an imaging study of children born very preterm, in which frontal network alterations have been associated with higher-cognitive function deficits (Fischi-Gomez et al., 2015; Fischi-Gomez et al., 2016). A recent meta-analysis showed that prophylactic EPO improved cognitive outcome at 18 months of infants born very preterm (Fischer et al., 2017). Whether the improved structural connectivity in the EPO treated children seen in this study is reflected by better higher cognitive function such as executive function has to be shown.

We identified the following technical limitations in our study protocol. The low number of diffusion-weighting directions (n = 21) allows estimating only one tensor per image voxel reliably; this inherently affects the result of probabilistic tracking of connections. Newborns underwent MRI in natural sleep, and longer scan time during DTI will greatly increase the ratio of corrupted image frames due to head motion, and this limits the number of applicable diffusion-weighting directions. During pre-processing, we used linear registration to correct eddy-current induced distortions, during which residual misalignment can adversely impact model fitting to the data. Improved correction would be achieved by techniques that model the diffusion signal and correct for eddy-currents and head motion in one resampling step (Andersson and Sotiropoulos, 2016). More accurate modeling of the diffusion propagator function and the application of more complex tractography techniques to the newborn brain would result in more reproducible description of the structural connectome in the future (Kuhn et al., 2016; Lebel et al., 2012; Wakana et al., 2007). We found only a weak effect of the EPO treatment on the overall structural brain connectivity of preterm infants. This might be explained by the timing and short duration of the EPO intervention.

Early high dose administration of EPO in very preterm infants shows a neuroprotective and trophic effect on the global brain connectivity, which effect was found to be marginally larger at a more advanced infant age. The connectome of EPO treated infants was characterized by increased average local efficiency, and increased structural connectivity strengths in edges that contribute to the peripheral connectivity in the presence of a preserved connectivity core. The increased integrity and local efficiency of the structural brain connectome might well be reflected by improved executive functions and general cognitive abilities in later childhood.

Funding

This work was supported by the OPO Foundation, the Forschungszentrum fur das Kind (FZK) Grant, the “Stiftung für wissenschaftliche Forschung an der UZH” and the EMDO Foundation, R.T. was supported by the EMDO Foundation. The study was supported by a grant received from the Swiss National Science Foundation (SNF 3200B0-108176).

The following are the supplementary data related to this article.

Supplementary Table 1

Region of interest abbreviations and locations.

Diffusion tensor tractography script used in the manuscript.

Neonatal fractional anisotropy (FA) template used in the manuscript.

Diffusion tensor image processing script used in the manuscript.