Long-term reorganization of structural brain networks in a rabbit model of intrauterine growth restriction

Characterization of brain changes produced by intrauterine growth restriction (IUGR) is among the main challenges of modern fetal medicine and pediatrics. This condition affects 5-10% of all pregnancies and is associated with a wide range of neurodevelopmental disorders. Better understanding of the brain reorganization produced by IUGR open a window of opportunity to find potential imaging biomarkers in order to identify the infants with a high risk of having neurodevelopmental problems and apply therapies to improve their outcomes. Structural brain networks obtained from diffusion magnetic resonance imaging (MRI) is a promising tool to study brain reorganization and to be used as a biomarker of neurodevelopmental alterations. In the present study this technique is applied to a rabbit animal model of IUGR, which presents some advantages including a controlled environment and the possibility to obtain high quality MRI with long acquisition times. Using a Q-Ball diffusion model, and a previously published rabbit brain MRI atlas, structural brain networks of 15 IUGR and 14 control rabbits at 70 days of age (equivalent to pre-adolescence human age) were obtained. The analysis of graph theory features showed a decreased network infrastructure (degree and binary global efficiency) associated with IUGR condition and a set of generalized fractional anisotropy (GFA) weighted measures associated to abnormal neurobehavior. Interestingly, when assessing the brain network organization independently of network infrastructure by means of normalized networks, IUGR showed increased global and local efficiency. We hypothesize that this effect could reflect a compensatory response to reduced infrastructure in IUGR. These results present new evidence on the long-term persistence of the brain reorganization produced by IUGR that could underlie behavioral and developmental alterations previously described. The described changes in network organization have the potential to be used as biomarkers to monitoring brain changes produced by experimental therapies in IUGR animal model. exploratory analysis of the regional features altered by IUGR performed. Both global and regional network features associated with neurobehavioral tests results. The contribute the on long-term brain changes associated to neurobehavioral dysfunctions in IUGR, showing the feasibility of using brain networks features from diffusion MRI as biomarkers to assess and potentially monitor treatment in IUGR using experimental models.

potential imaging biomarkers in order to identify the infants with a high risk of having neurodevelopmental problems and apply therapies to improve their outcomes. Structural brain networks obtained from diffusion magnetic resonance imaging (MRI) is a promising tool to study brain reorganization and to be used as a biomarker of neurodevelopmental alterations. In the present study this technique is applied to a rabbit animal model of IUGR, which presents some advantages including a controlled environment and the possibility to obtain high quality MRI with long acquisition times.
Using a Q-Ball diffusion model, and a previously published rabbit brain MRI atlas, structural brain networks of 15 IUGR and 14 control rabbits at 70 days of age (equivalent to pre-adolescence human age) were obtained. The analysis of graph theory features showed a decreased network infrastructure (degree and binary global efficiency) associated with IUGR condition and a set of generalized fractional anisotropy (GFA) weighted measures associated to abnormal neurobehavior. Interestingly, when assessing the brain network organization independently of network infrastructure by means of normalized networks, IUGR showed increased global and local efficiency. We hypothesize that this effect could reflect a compensatory response to reduced infrastructure in IUGR. These results present new evidence on the long-term persistence of the brain reorganization produced by IUGR that could underlie behavioral and developmental alterations previously described. The described changes in network organization have the potential to be used as biomarkers to monitoring brain changes produced by experimental therapies in IUGR animal model.

INTRODUCTION
Intrauterine growth restriction (IUGR) affects 5-10% of all pregnancies and is a major public health issue, being a prevalent condition that has been associated with a wide range of short-and long-term neurodevelopmental and cognitive dysfunctions (Arcangeli et al., 2012;Baschat, 2013), even in adulthood (Løhaugen et al., 2013). With the significant advance of magnetic resonance imaging (MRI) in the recent years, the brain alterations and reorganization underlying these neurofunctional alterations are starting to be elucidated. It has been suggested that brain reorganization starts in utero, where different patterns of cortical development (Egaña-Ugrinovic et al., 2013) and altered quantitative MRI texture predictive of altered neurodevelopment  have been shown in IUGR. At neonatal period IUGR has been reported to have decreased volume in gray matter (GM) (Tolsa et al., 2004) and hippocampus (Lodygensky et al., 2008) and discordant patterns of gyrification (Dubois et al., 2008). At one year of age, persistence of structural changes has been demonstrated, including reduced volumes of GM (Padilla et al., 2011) and decreased fractal dimension in both GM and white matter (WM) that correlate with abnormal neurodevelopment (Esteban et al., 2010). Studies on IUGR at later ages have reported changes in regional brain volumes and cortical thickness in 4 to 7-year-old children (De Bie et al., 2011), reduced volumes for thalamus and cerebellar white matter (Martinussen et al., 2009), and thinning of corpus callosum and general WM reduction  in adolescents. There is a need to better characterize the brain reorganization underlying neurodevelopmental and cognitive dysfunctions in IUGR.
Likewise, the development of imaging biomarkers is an urgent clinical and experimental need (Ment et al., 2009).
The study of brain connectivity holds great promise for the development of pathophysiological insights and biomarkers of human disease characterized with subtle brain changes that are not reflected in conventional MRI techniques (Gratacos, 2012). Indeed, one of the major recent advances in the application of new MRI modalities has been the emerging technique of "connectomics" (Hagmann, 2005), opening the possibility to extract macroscopic circuitry of the connections of the brain, in what has been called "the connectome" (Sporns et al., 2005). In particular, the use of graph theory analyses on brain networks has been demonstrated to be a useful tool to characterize brain organization by a few comprehensible parameters (Bassett and Bullmore, 2009). Different sets of data, including functional MRI and diffusion MRI, have been used to extract macroscopic brain networks and analyze network features in healthy adults, adolescents and infants (Gong et al., 2009a;Hagmann et al., 2008;Hagmann et al., 2010;Iturria-Medina et al., 2008;Yap et al., 2011) and to report altered group connectivity parameters in a wide range of neurological, neurobehavioral and neurodegenerative diseases (Alexander-Bloch et al., 2010;Liu et al., 2008;Lo et al., 2010;Shu et al., 2009;Shu et al., 2011;Wang et al., 2009;Wu et al., 2009). Importantly, connectomics and graph theory features have been shown to be potential tools to develop biomarkers to predict neurological outcomes in adult Li et al., 2009;Wee et al., 2010;Wen et al., 2011) and perinatal diseases Batalle et al., 2013;Tymofiyeva et al., 2012). Particularly, brain networks of one-year-old infants obtained from diffusion MRI have been reported to have reduced level of weighted organization and a pattern of altered regional network features that is associated with latter neurodevelopmental problems , showing their potential to develop imaging biomarkers to detect infants at high risk of having neurodevelopmental problems one year later. Nonetheless, whether persistent brain reorganization produced by IUGR persists at long-term (adolescence and adult period) and whether connectomic analysis could be a suitable tool to characterize the patterns induced by this conditions is still unknown.
Assessing long-term effects of IUGR in the human brain is a challenging task, limited by the influence of uncontrolled environmental factors (Hall and Perona, 2012) and the difficulty of obtaining sufficiently large sample sizes. The induction of IUGR in rabbit models has been proven to reproduce major features of human IUGR (Bassan et al., 2000;Eixarch et al., 2009;Eixarch et al., 2011).
Furthermore, white matter maturation process in rabbit is closer to humans than other species, since it starts in intrauterine period (Derrick et al., 2007). Hence, albeit their obvious limitations, rabbit model may be a useful tool to analyze long-term brain remodeling in IUGR. They could play a key role in the definition of image biomarkers for early diagnosis that are critical to demonstrate changes after the application of experimental therapies, especially when those should be tested in fetuses or neonates.
Besides highly reproducible experimental conditions, high quality MRI with long acquisition times can be performed in isolated whole brain preparations. Using this model, regional brain changes in fractional anisotropy, correlated with poorer outcome in neurobehavioral tests has been reported in newborns , some of them persisting in preadolescent period , where changes in the connectivity of anxiety, attention and memory networks has been shown. Due to the recent development of an MRI rabbit brain atlas , the possibility to obtain whole brain structural networks based on diffusion MRI arises. This opens the opportunity to assess long-term network reorganization associated with functional impairments without a priori hypothesis, taking advantage of the huge potential of graph theory measures to characterize brain functioning and organization (Bassett and Bullmore, 2009) that have been previously used to characterize one-year-old infants with IUGR Batalle et al., 2013).
In the present study, graph theory features from diffusion MRI brain networks were calculated in 15 rabbits with surgically induced IUGR and 14 controls at equivalent preadolescent age, in order to assess the long-term impact of IUGR in brain organization that could underlie behavioral and developmental alterations. The results showed a specific pattern of global network features altered in IUGR, characterized by an impaired network infrastructure, but an increase in the relative terms of organizational efficiency that we hypothesize to be associated with a compensatory effect in IUGR.
An exploratory analysis of the regional features altered by IUGR condition was also performed. Both global and regional network features were associated with neurobehavioral tests results. The results here presented contribute to the knowledge on long-term brain changes associated to neurobehavioral dysfunctions in IUGR, showing the feasibility of using brain networks features from diffusion MRI as biomarkers to assess and potentially monitor treatment in IUGR using experimental models.

METHODS
The design of the study and each of the steps of the procedures are shown in Figure 1. A detailed description of the methodology used is included in this section.

Animals, study protocol and surgical model
Animal experimentation of this study was approved by the Animal Experimental Ethics Committee of the University of , and all efforts were made to avoid or minimize suffering.
Part of the animals used in this study has been previously used in a recent study .
From 13 New Zealand pregnant rabbits provided by a certified breeder, we selected two cases and two controls of each dam at birth. At 70 th postnatal day, all surviving cases and one control for each case were included resulting in a total population of 30 rabbits (15 with induced IUGR and 15 sham controls). Dams were housed for 1 week before surgery in separate cages on a reversed 12/12h light cycle, with free access to water and standard chow. At 25 days of gestation (term at 31 days), a ligation of 40-50% of uteroplacental vessels was performed following a previously described protocol (Eixarch et al., 2009). Briefly, after midline abdominal laparotomy, the gestational sacs of both horns were counted and numbered. Afterwards, only one uterine horn was kept outside the abdomen and the induction of IUGR was performed by ligating 40-50% of the uteroplacental vessels of all the gestational sacs from this horn. After the procedure, the abdomen was closed in two layers and postoperative analgesia (meloxicam) was administered for 48 hours. After surgery, the animals were allowed free access to water and standard chow for 5 days until delivery. Cesarean section was performed at 30 days of gestation and living pups were obtained. All living newborns were weighed and identified by a subcutaneous microchip inserted in their back (Microchip MUSICC, Avid Microchip S.L., Barcelona, Spain). Cases were considered those pups delivered from the ligated horn, whereas controls were those delivered from the contralateral horn (non-ligated). Both cases and controls were housed with a wet nurse rabbit with part of the offspring until the 30th postnatal day when they were weaned. Thereafter both groups of rabbits were housed in groups of three with a reversed 12/12h light cycle with free access to water and standard chow. On the 70th postnatal day, which is considered to be equivalent to preadolescence period in humans in terms of sexual maturity (Moorman et al., 2000), functional tests were applied and the rabbits were anaesthetized and sacrificed thereafter. Left and right common carotid arteries were cannulated and the brains were perfused with phosphate-buffered saline (PBS) followed by 4% paraformaldehyde PBS. Then, brains were excised from the skull and immersed in 4% paraformaldehyde PBS at 4°C for 48h.

Neurobehavioral tests and cognitive evaluation
In order to assess neurobehavioral changes related to emotion and cognition, two standard test used extensively in rodents and recently adapted to rabbits  were used: Open Field Behavioral Test (OFBT) and Object Recognition Task (ORT). Neurobehavioral tests were performed in a subset of rabbits. Particularly, of the total sample of 29 subjects, 20 rabbits were selected to perform OFBT and ORT.
OFBT evaluates locomotion and exploratory activities competing against fear, anxiety and attention (Bouet et al., 2003;Walsh and Cummins, 1975), based on the procedures originally described by Hall (Hall, 1934). Briefly, the testing area was divided into 36 squares of 23 x 23 cm, the four central squares were considered as the internal area and the remaining squares were defined as the external (peripheral) area. The rabbits were taken out of their cage wrapped with a cloth and placed close to one of the lateral walls (starting point) and behavior was assessed during 10 minutes.
In order to minimize interference due to human contact, each session was videotaped and later evaluated by two blinded observers. Rabbits that escaped from the testing area were excluded from the analysis. The main parameters were recorded including latency of leaving the starting point (seconds) and number of squares explored (total, internal and external). In order to obtain a single dichotomic value summarizing a normal or an abnormal performance of a given subject, normality of the control population was established at the 25 th centile for each of the four main parameters recorded. Having two or more of these items below this 25 th centile was considered to be an abnormal exploration. OFBT was performed successfully in 10 controls and 9 IUGR.
ORT assesses short-term memory, especially recognition (Olton and Feustle, 1981) and attention capacity (Cowan et al., 1999), based on the tendency of rodents to explore new stimuli for a longer time compared to familiar stimuli (Dere et al., 2006;Ennaceur and Aggleton, 1997;Mumby, 2001).
The test was adapted from the original description (Ennaceur and Delacour, 1988), including some modifications in the stimulus used (instead of a visual stimulus, an odor-based stimulus was used).
Briefly, after a familiarization phase where the rabbit was presented with two boxes containing the same odor-based stimuli (apple), the animal was returned to its cage for a 30-minutes retention interval and returned to the testing area for a 5 minutes testing phase, where one of the boxes was removed and replaced by a novel stimulus (orange). Cumulative time exploring each object was recorded and discrimination index (DI) was calculated. DI represents the ability to discriminate the novel from the familiar object and was calculated as follows:

= Novel Object Exploration Time − Familiar Object Exploration Time Novel Object Exploration Time + Familiar Object Exploration Time
Learning criteria was considered when DI was above 0, and dichotomized version of the score was generated accordingly. ORT was performed in all animals with a successful OFBT test. Animals that did not explore the familiar object at least once in the testing phase or did not explore any of the objects in the familiarization phase were excluded from the analysis, as previously suggested (De Bie et al., 2011;Illa et al., 2013). ORT was performed successfully in 8 controls and 6 IUGR. Further details on the specifics of both tests can be consulted at Illa et al. (2013).

MRI acquisition
MRI was performed on fixed brains using a 7T animal MRI scanner (BrukerBioSpin MRI GmbH).
High-resolution three-dimensional T1 weighted images were obtained in the brain samples by a Modified Driven Equilibrium Fourier Transform (MDEFT) sequence with the following parameters: Time of Echo (TE) = 3.5 ms, Time of Repetition (TR) = 4000 ms, 0.7-mm slice thickness with no interslice gap, 70 coronal slices, in-plane acquisition matrix of 188 x 188 and Field of View (FoV) of 28 x 28 mm 2 , resulting in a voxel dimension of 0.15 x 0.15 x 0.7 mm 3 . Diffusion-weighted images (DWI) were acquired using a standard diffusion sequence covering 126 gradient directions with a bvalue of 3000 s/mm 2 together with a reference (b=0) image. Other experimental parameters were: TE = 26 ms, TR = 250 ms, slice thickness = 0.7 mm with no interslice gap, 70 coronal slices, in-plane acquisition matrix of 40 x 40 , FoV of 28 x 28 mm 2 , resulting in a voxel dimension of 0.7 x 0.7 x 0.7 mm 3 . The total scan time for both acquisitions was 13h56m40s. Any potential tissue alteration, mainly significant tissue loss, was considered as exclusion criteria, being one control of the total population excluded for this reason.

Pre-processing and tractography
Brain tissue was segmented from the background in the T1 volume by means of an Otsu threshold method (Otsu, 1975). In the case of DWI, brain tissue was segmented from the background by means of an in-house algorithm previously described  based on the high SNR of the brain tissue on the average diffusion volume. The orientation diffusion function (ODF) of each voxel was reconstructed by means of a Q-Ball approach (Tuch, 2004) and used to reconstruct fiber tracts by means of a deterministic tractography algorithm with an angle threshold of 30 degrees implemented in Diffusion Toolkit (http://trackvis.org/dtk/) (Wang et al., 2007). Extracted brain volume size and number of fibers reconstructed for each subject were assessed.

Brain parcellation
Automatic brain parcellation of the subjects' brain was performed taking advantage of a recently  (Warfield et al., 2002). The elastic transformation obtained matching the template atlas and each subject's T1 was applied to the ROI labels, obtaining a parcellation of each brain in 60 ROIs. In order to align the labels obtained for each subject in the T1 volume to its corresponding DWI, affine registration (Studholme et al., 1999) between T1 and the baseline image was performed with IRTK (www.doc.ic.ac.uk/~dr/software/). Discrete values of the labels were preserved by nearest neighbor interpolation in both transformations. ROIs comprising only WM tissue were discarded, leaving a total of 44 regions for each subject (Table 1), each of which considered as a node in the brain networks obtained.

Network extraction
Brain network of each subject was extracted by means of an in-house algorithm combining the fiber tracts obtained by Q-Ball tractography and the ROIs resulting from the automatic brain parcellation previously performed. Two nodes (regions) i and j were considered to be connected by an edge eij when there existed at least one fiber bundle f with end-points in i and j regions, with selfloops excluded. In addition to the binary network produced with this approach, weights were also assigned to each edge eij. There is still not a gold standard in the literature to assign connectivity weight between two regions, hence, two of the most used approaches in the literature were followed to weight the connectivity between each pair of regions, w GFA , average generalized fractional anisotropy (GFA) along all the fibers connecting a pair of regions; and w FD , fiber density (FD) defined as Hagmann et al. (2008) with slight modifications: is the FD-weight given to the connection between nodes i and j, V i is the volume of the ROI i, F ij is the set of all fibers connecting ROI i and j, and l(f) is the length of fiber f along its trajectory, introduced in the denominator to eliminate the linear bias towards longer fibers introduced by the tractography algorithm. In addition, the resulting weighted networks were normalized by the total weight of all the connections in the network, to assess the brain organization independently of the network average strength.
In summary, five networks were extracted for each subject: binary, GFA-weighted, FD-Weighted, normalized GFA-weighted and normalized FD-weighted.

Network analysis
Network analysis by means of graph theory allows obtaining a set of features that summarize the organization of a brain network, represented as an adjacency matrix (binary or weighted). In the present study, we applied graph theory features over each of the five brain network classes computed for each subject. All the features were computed using Brain Connectivity Toolbox (Rubinov and Sporns, 2009).
Global functioning of each network was assessed by its infrastructure (average degree for binary networks and average strength for weighted networks), integration (binary and weighted global efficiency) and segregation (binary and weighted local efficiency). In addition, binary and weighted characteristic path length and average clustering coefficient were also assessed.
Regional characteristics were evaluated by means of nodal degree and strength, features that assess the total connectivity of a node in a network in its binary and weighted versions respectively.
In addition, nodal efficiency was also assessed in its binary and weighted version. Nodal efficiency measures the efficiency of the sub-network associated to a given node. Nodes with a high nodal efficiency indicate a high fault tolerance of the network to the elimination of the given node, which is given by a high clustering of the neighborhood of this node (Achard and Bullmore, 2007).
In order to assess the small-world behavior of the computed networks (Watts and Strogatz, 1998), binary characteristic path length (L p ) and binary clustering coefficient (C p ) were compared with the average characteristic path length and clustering coefficient of one hundred equivalent random networks with the same size and degree distribution of each subject's network. The ratio of the values obtained for the original subject network and its equivalent random version allowed to obtain the , and σ = γ λ ⁄ . Small-worldness value (σ) is above 1 in small-world regimes (Humphries and Gurney, 2008).
Formulation of the graph theory features used to assess each network is based on the definitions compilated by Rubinov and Sporns (2009).

Levels of analysis: raw, cost-corrected, cost-integrated and normalized strength
We analyzed the obtained networks with different approaches in order to catch different aspects of their organization. The first straightforward approach is to directly apply binary and weighted graph theory measures to the raw binary, GFA-and FD-weighted networks obtained for each subject.
However, some of these measures could be closely associated to the measured infrastructure of the network, i.e., their average network strength and degree (Ginestet et al., 2011;van Wijk et al., 2010).
In order to assess organization independently from network density (i.e., average network degree), a cost-thresholding approach was followed (Achard and Bullmore, 2007) obtaining binary networks for a set of costs in the rank from 0 to 0.39 at 0.01 steps, using GFA-and FD-weighted networks to select the "strongest" connections that yield to a network with the desired network cost for each subject.
0.39 is the smallest network density value of any subject included in the study, hence being the largest network density value that it is possible to study for all subjects fairly. This approach allows comparing network features independently of differences in the raw number of connections (cost) between subjects. However, it increases the number of comparisons and the complexity of the results interpretation. In order to disentangle network density from topology in a single value for each graph theory feature, a cost-integrated version of each topological version was calculated as the mean value across the assessed rank. In addition, this cost-correction approach has an additional shortcoming in the fact that it leads to different connections sets across subjects. In order to assess the behavior of the networks using the same set of connections for all the subjects, an experimental approach is used in Supplementary Materials inspired in the methodology used by Gong et al. (Gong et al., 2009b).
With the goal to assess the pure weighted organizational characteristic of the network topology, normalized weighted networks were also calculated for each subject in their GFA and FD versions.
This normalization was performed by means of dividing each connection weight by the total weight on the network. Hence, normalized weights were calculated as w X "X" is GFA or FD. By this means, each normalized weight w GFA N (i, j) and w FD N (i, j) corresponds to the percentage of connection weight used in the link between i and j relative to the total amount of weights in the network. This way the measures performed on normalized GFA-and FD-weighted networks were independent to the total amount of connectivity strength each subject had, and then assessing only its distribution of weights.

Statistics
Comparisons among groups were performed by general linear models (GLM) with gender as cofactor. Interaction of group and gender was first included into the model, but as it did not show any significant effect was excluded of the final model. Significance was declared at p<0.05 (uncorrected).
Normality was assessed with Shapiro-Wilk Test and homoscedasticity with Levene's Test, and when the null hypothesis rejected, log-transformation was performed prior to GLM analysis. Descriptives of the variables were expressed as mean (standard deviation) for normal distributions and median (interquartile range) for non-normal distributions. Differences between cases and controls of dichotomized values of OFBT and ORT were assessed with a binary logistic regression with group and gender as co-factors. Association of network features with OFBT and ORT was performed by means of a partial correlation controlling the effect of gender. Association with dichotomized OFBT and ORT was performed by a GLM with gender as co-factor. The most discriminative network features of dichotomized OFBT and ORT were assessed by means of a conditional step forward binary logistic regression with an entry probability of 5%. Note that the population available with OFBT and ORT was 19 subjects (10 controls and 9 IUGR) and 14 subjects (8 controls and 6 IUGR) respectively. Regional features were corrected for multiple comparisons with a False Discovery Rate (FDR) approach (Benjamini et al., 2006), controlling alpha error to 10% and 5% for each feature calculated among all 44 regions. Regional alterations were shown in the rabbit brain atlas template using BrainNet viewer (www.nitrc.org/projects/bnv/) (Xia et al., 2013). The software package SPSS 18.0 (SPSS, Chicago, IL) was used for the statistical analyses. Computational algorithms were implemented using MATLAB (2009b, The MathWorks Inc., Natick, MA).

Population
After applying exclusion criteria, the final population under study consisted in 14 controls and 15 subjects with induced IUGR. As expected, birth weight in IUGR subjects was significantly decreased vs. 66.7%, p=0.059). In summary, these results showed that IUGR rabbits have a significant degree of anxiety, attention and memory problems when compared to controls.

Global network features
As shown in Figure 2, the analysis of the basic infrastructure of the network showed a decreased average degree (p=0.039) in IUGR cases. The analysis of the average strength of the weighted networks did not show any significant effect, although a tendency towards significance is found in FDweighted average strength (p=0.084).
Analysis of small world characteristics (Figure 3) showed that both cases and controls have a small-worldness σ value higher than one for all subjects, a typical characteristic of small-world networks, however, no significant differences were found between cases and controls.
The assessment of brain network organization by means of the analysis of global and local efficiency ( Figure 4) demonstrated several differences between controls and IUGR. Having a significantly reduced average degree it is expectable to find that IUGR also had a reduced binary global efficiency (p=0.041). IUGR group did not show significant differences in any weighted measure of efficiency, however, the analysis of normalized brain networks, independent to differences in the average strength of each subject, showed that IUGR group had a significantly increased normalized GFA-weighted global efficiency (p=0.018), normalized GFA-weighted local efficiency (p=0.039) and normalized FD-weighted local efficiency (p=0.010). Note that while FD indirectly measures the axonal density connecting two regions, GFA measures have been related with the integrity of these axonal connections. These results were also confirmed when looking directly into normalized characteristic path length and average clustering coefficient (Supplementary Material). Finally, when forcing all the subjects' networks to have the same density (cost-corrected approach), differences between cases and controls were found for several cost values in GFA-weighted global efficiency and FD-weighted local efficiency ( Figure 5) and GFA cost-integrated measure of global efficiency was found to have a tendency of being increased in IUGR (p=0.070). Overall, global network features obtained suggested an altered brain network organization characterized by a reduced level of network connectivity and a compensatory effect after assessing pure organizational features by means of different normalization approaches. Although a compensatory effect in local efficiency was observed with both GFA and FD weighting approaches, statistically increased global efficiency was only observed for GFA-weights.

Regional network features
The analysis of regional network features showed a pattern of alterations in IUGR group ( Figure   6, Supplementary Table 1). Particularly, IUGR is characterized by a significantly lower nodal degree in several brain regions, although vermis is the only region withstanding a 5% FDR correction.
Interestingly, right cingulate cortex, lenticular nucleus and vermis also showed a significantly decreased GFA-weighted strength, both in absolute and relative terms (normalized), vermis also resisting a 5% FDR correction. In line with the results obtained for global network features, regions with a higher degree or strength in IUGR group were only found in relative terms.
Concerning binary nodal efficiency, IUGR showed a significant decrease in right parietal cortex and right frontal cortex, but an increase in vermis. FD-weighted network efficiency also showed a significantly increase in left cerebellar hemisphere and vermis. In line with the significant increase found in global normalized weighted features, IUGR showed a wide pattern of regions with a significantly increased normalized efficiency, including basal forebrain, piriform cortex, cerebellar hemisphere and amygdala in normalized GFA-weighted networks (right piriform cortex and cerebellar hemisphere withstanding a FDR correction at 10%). Efficiency of mesencephalon, pons, right parietal and cingulate cortex, hippocampus, cerebellar hemisphere and vermis were also increased in normalized FD-weighted networks; with left cerebellar hemisphere resisting a 10% FDR correction and vermis resisting a 5% FDR correction.
In summary, IUGR mainly showed a pattern of decreased features in raw networks, especially in binary nodal efficiency and degree and GFA-weighted strength, and a pattern of increased features in normalized GFA-and FD-weighted networks, especially in nodal efficiency.

Association of network characteristics with altered neurobehavior
The analysis of dichotomized OFBT normal/abnormal status with GLM showed significantly decreased GFA-weighted average strength (p=0.039), GFA-weighted global efficiency (p=0.035) and GFA-weighted local efficiency (p=0.017) in the abnormal group. A tendency to have increased normalized FD-weighted global efficiency (p=0.058) and cost-integrated GFA global efficiency (p=0.083) was also found in the group with an abnormal exploration. Analysis of ORT normal/abnormal exploration showed only significantly increased cost-integrated GFA global efficiency (p=0.007) in the abnormal group.
The association with normal/abnormal exploration in OFBT and ORT was confirmed by partial correlation of OFBT continuous values and global features (Supplementary Table 2). Briefly, OFBT scores showed significant correlations with several GFA-weighted measures (average strength, characteristic path length, average clustering, global efficiency and local efficiency) and with normalized FD-weighted global efficiency. ORT, however, only showed some tendencies towards significance with GFA-weighted clustering coefficient and GFA cost-integrated local efficiency.
Interestingly, the results showed that for all significant correlations, and for those showing a tendency towards significance, absolute GFA-weighted measures had a significant correlation with neurobehavioral scores in the opposite direction of cost-integrated and normalized measures.
Regional features associated with abnormal neurobehavior were assessed by means of a GLM with dichotomized OFBT / ORT and gender as co-factor (Figure 7). The results showed a significant association of several regional features with neurobehavioral tests. As it could be expected from the strong association with global GFA-weighted network features, abnormal OFBT exploration was associated with a wide rank of GFA-weighted nodal measures, particularly with nodal efficiency features. The association of regional features with abnormal ORT was not so spread, but some consistent patterns of alterations across different kinds of networks were also observed, as it is the case of right thalamus and right caudate nodal strength and efficiency features.
Further analysis to find the network features with a higher potential as an image biomarker of altered neurodevelopment was performed by means of a step forward binary logistic regression. All the assessed global network features were entered in a first analysis, selecting GFA-weighted local efficiency as the more discriminative in OFBT (90% sensitivity, 66.7% specificity, Nagelkerke R 2 =0.335, Chi 2 =5.494, p=0.019). The most discriminative global feature in ORT was cost-integrated global efficiency (88.9% sensitivity, 80% specificity, Nagelkerke R2=0.623, Chi 2 =8.467, p=0.004), selecting also normalized FD-weighted local efficiency which allowed a full separation of normal/abnormal ORT (Chi 2 =18.249, p<0.001). In a second analysis, all the assessed regional network features were introduced. Regional features more discriminative of OFBT were GFAweighted right insular cortex strength and normalized FD-weighted septal nuclei nodal efficiency (90% sensitivity, 89.5% specificity, Nagelkerke R 2 =0.939, Chi 2 =23.101, p<0.001). In the case of ORT, the most discriminant regional network feature was normalized FD-weighted nodal strength of right thalamus, which also fully separated by itself normal/abnormal ORT (Chi 2 =18.249, p<0.001).

DISCUSSION
This study describes the use of diffusion MRI brain networks to analyze the brain reorganization in a rabbit model of IUGR at 70 days of age, equivalent to preadolescence in humans. The results show that IUGR brain reorganization persists into pre-adolescence equivalent age with significant differences at various levels in network infrastructure. GFA-weighted features and regional network features were associated with poor neurobehavioral performance. Furthermore, the increased efficiency of normalized GFA-weighted networks associated with IUGR suggests the existence of reorganizational compensatory processes in post-natal life.

Global network features
Comparison of the brain networks obtained with random equivalent networks showed small-world characteristics. This result is coherent with previously found in humans, animals and a wide range of networks (Bullmore and Sporns, 2009). Small-world networks are suggested to be a balance between random and lattice networks, having a high level of organization but also allowing to communicate any given pair of nodes with relatively few intermediate steps.
No differences in small-world features were found in IUGR, but the analysis of other global network features showed a pattern of brain reorganization persisting at preadolescent equivalent age. Particularly, reduced average degree is found in IUGR, supporting the idea that IUGR have an impaired network infrastructure. As expected, having a reduced network infrastructure may be producing a dysfunctional connectivity, expressed by a reduced binary global efficiency. Supporting this concept, we observed also a trend of IUGR to have reduced FD-weighted average strength. Note that strength can be understood as the total power of brain connections between regions. However, although having a less connected brain network, GFA-and FD-weighted efficiencies were not found significantly decreased in IUGR. Broadly, global efficiency assess how easy is on average to connect a pair of regions, while local efficiency has been related with the level of clusterization of the network. The results obtained in raw GFA-weighted efficiencies are in contraposition with previous findings suggesting significantly reduced fractional anisotropy (FA) weighted global and local efficiency found at one year of age in human IUGR . One possible explanation could be related with changes in postnatal brain maturation after IUGR. Previous studies have reported that there is a compensation in the long-term myelination and WM volume in a guinea pig model of IUGR (Tolcos et al., 2011). In this line, previous studies assessing WM alterations in a rabbit animal model of IUGR using voxel based analysis have reported a less prominent pattern of alterations in pre-adolescence age  than in neonatal age . On the other hand, in this study we could not reproduce previous observations in one-year-old infants' networks showing reduced average degree in IUGR subjects . We hypothesize that this apparent discrepancy could be explained by evolutionary changes of brain network reorganization associated with IUGR at different ages. Alternatively, we cannot exclude that this effect could be specific of IUGR in this rabbit animal model.
An interesting finding of this study was the observation of increased efficiencies in the normalized networks of IUGR compared with controls. Consistent with this finding, the evolution of graph theory measures in a cost-thresholding approach showed a discordant pattern in IUGR, particularly an increased FD local efficiency and an increased GFA global efficiency for several network costs. We hypothesize that this effect might reflect a compensatory mechanism in IUGR subjects, whereby available network resources seem to be optimized and reorganized more efficiently from a network perspective. This compensatory effect in brain connectivity could also explain some other intuitively non-expected findings in IUGR in the long-term, such as the reported, regional cortical thickening in small for gestational age children (De Bie et al., 2011;Martinussen et al., 2005).

Regional network features
Concerning regional differences found in IUGR, only vermis features withstood a FDR correction at 5%, while only features of cerebellar hemispheres and piriform cortex withstood a FDR correction at 10% (Supplementary Table 1). Particularly, degree and normalized GFA-weighted strength of vermis were significantly reduced in IUGR, and nodal efficiency of both raw and normalized FDweighted vermis were significantly increased in IUGR. We could interpret that IUGR subjects have a vermis with reduced neighborhood but with a stronger connectivity (higher clustering). Again, this would support the hypothesis of a compensatory effect and goes in line with the results obtained on global network features. Interestingly, alterations in cerebellum and vermis connectivity in infants with IUGR have previously been reported by network analysis Batalle et al., 2013), as well as reduced WM volumes in cerebellar areas in short- (Padilla et al., 2011) and long-term follow-up studies (De Bie et al., 2011;Martinussen et al., 2009). Interestingly, cerebellar areas are implicated in motor learning, memory, cognition and behavior (Baillieux et al., 2008), and have been reported to increase their volume in infants with IUGR after intervention by means of individualized intensive care nursery associated with an improvement in executive function (McAnulty et al., 2013).
Considering other animal models of IUGR, in guinea pig it also has been reported a reduction in cerebellar WM volume (Tolcos et al., 2011), and a reduction in Purkinje neuronal population (Mallard et al., 2000). However, we must note that in previous voxel-based analysis studies in IUGR rabbit model Illa et al., 2013) we only observed mild differences in cerebellar areas.
We hypothesize that this could be partially explained by stronger alterations in the sub-network associated to cerebellum (described by network features), than in the structure of cerebellum by itself.
Another factor that could explain this discrepancy can be the technical characteristics of voxel-based analysis technique, especially in cerebellar areas, where the limited resolution is especially critical for cerebellar WM.
Beside those associated with cerebellar areas, other regional features were also significantly different in IUGR, although only normalized FD-weighted efficiency of right piriform cortex withstood a 10% FDR correction. It seems interesting that regional features of normalized networks were found significantly increased in IUGR mainly in specific basal and deep gray matter areas, such as amygdala, hippocampus, thalamus, medulla oblongata, and some of the cortical areas controlling essential functions for the rabbit, as is the case of piriform cortex, which has been strongly associated to the smelling processing (Kadohisa and Wilson, 2006). One possible interpretation of these results could be that a compensatory reorganization in IUGR brain might occur preferentially in brain areas regulating function with critical importance in survival, such as memory, attention or smelling processing. Considering other features that did not withstood FDR correction, we must note that the regional features altered in IUGR were coherent with the global results obtained, mainly showing reduced values in binary and raw GFA-weighted measures (especially in degree/strength) and increased values in normalized measures (especially in efficiency). Interestingly, changes in regional features of raw networks were manly identified in brain regions located in the right hemisphere of the brain, suggesting a certain asymmetry in the reorganization associated with IUGR. Reports on asymmetrical brain alterations produced by IUGR are scarce in the literature, and further long-term studies in humans would be required to confirm this effect and discard that is exclusive to animal models. Another interesting finding was left cingulate cortex and right lenticular nucleus decreased nodal degree and relative and absolute GFA-weighted nodal strength. Cingulate cortex features were also correlated with OFBT performance, which is in line with previous results in a similar sample of IUGR rabbit model, where FA changes in the cingulate cortex were found to be correlated to a series of neurobehavioral domains, especially those of OFBT . Note that cingulate cortex has been strongly associated with anxiety processes (Kim and Whalen, 2009), which is thought to have an important role on the behavioral response in OFBT. In addition, cingulate cortex volumetric differences in humans (Emond et al., 2009;Spampinato et al., 2009) as well as histological changes in rodents (Miller et al., 2012) have been associated with traits of attention deficit and anxiety.
Concerning the alterations in lenticular nucleus, FA changes have been previously reported in the putamen of neonatal  and long-term rabbit model of IUGR , as well as disrupted regional network features associated with putamen and globus pallidum in 1-yearold infants . Moreover, features associated to other subcortical GM nuclei including caudate and thalamus were also found affected in IUGR rabbits in this study. These findings are in line with previous evidence suggesting the important role of striatal injury as a risk factor of behavioral disturbances associated to IUGR (Toft, 1999), and alterations in cortico-striato-thalamic network has been associated to cognitive disorders such as attention deficit hyperactivity disorder (Castellanos et al., 1994;Faraone and Biederman, 1998), which in turn is more prevalent in IUGR (Heinonen et al., 2010). Concerning the specific associations with OFBT and ORT by means of a step forward binary logistic regression, features of insular cortex and septal nuclei were found strongly associated with OFBT and features of thalamus with ORT. Septal nuclei participation in the regulation of anxiety and depression in experimental models has been extensively documented (Estrada-Camarena et al., 2002), as well as insular cortex, which has been proposed to have a key role in anxiety proneness (Paulus and Stein, 2006). Importantly, thalamus has been previously associated to memory impairments in animal models (Mitchell and Dalrymple-Alford, 2005) and with poor performance in object recognition tasks in patients with schizophrenia (Heckers et al., 2000), disorder that has also been associated to low birth weight (Nilsson et al., 2005).
In summary, a pattern of regional network features altered in IUGR were identified, some of them being specifically associated with abnormal neurobehavior. Noteworthy, alterations in network features of cerebellum stand out as the more strongly associated with IUGR condition, and septal nuclei, insular cortex and thalamus with poor performance in neurobehavioral tests. These results, and its association with previous literature, reinforce the potential of network features based on these regions to be used as biomarkers of altered neurodevelopment.

Methodological issues and future work
Although the study of brain networks in animal models has been seminal to the field (Sporns et al., 2005), brain networks obtained from diffusion MRI in small animals are scarce. To the best of our knowledge, this is the second study of this kind in small mammals (Iturria-Medina et al., 2011), and the first showing brain networks in rabbits. This animal model has some specific methodological difficulties that had to be overcome. The recent successful development of an MRI atlas for the rabbit brain , allowed to automatically segment and analyze brain networks from this model. It is common in human studies to confine deterministic tractography only to WM tissue. However, due to the difficulty of properly defining a WM mask in rabbits, the deterministic tractography was computed for the whole brain. This could create some redundant fibers, and increase tractography noise, but it was partially overcome by the use of a high angular resolution diffusion MRI acquisition, and the use of graph theory, particularly weighted measures, to analyze the results obtained. This limitation also ramifies in the calculation of FD-weighted networks. Instead of the formulation proposed by Hagman et al. (2008) to calculate FD between two regions (using the surface of the WM-GM interface of the regions involved in a given link in the denominator), we modified this formula introducing the volume of each region instead of the surface of the WM-GM interface. In addition to specific limitations of the rabbit model, some others common to the study of brain networks must be discussed. Using weighted measures of the brain networks obtained has a series of advantages, including a better management of weak and potentially non-significant links (Rubinov and Sporns, 2009); however, how the connectivity between regions must be quantified and its correlation with the underlying anatomical substrate is still an open question that has been approached in very different ways (Griffa et al., 2013). To this regard, instead of selecting a single approach, we decided to use two very different ways of quantifying the weight of brain connectivity: average GFA along the fibers, and fiber density between regions; characterizing different aspects of connectivity that we expected to be complementary. Furthermore, we also computed normalized versions of the networks, allowing disentangling the measures obtained from the variability of absolute weighted strength of each individual network. Notwithstanding this, connectivity could still be dependent to some extent to the variability in network density (average degree). To obtain measures independent of network density, analyzing the evolution of a brain network across a set of network costs is a common approach, although part of the intrinsic connectivity patterns of a network are lost by forcing all the subjects to have the same network density. This issue is especially important when isolated nodes in the network of subjects start appearing at some levels of cost. Although the effect of isolated nodes is partially mitigated by the use of efficiency network features instead of characteristic path length and clustering coefficient measures (Rubinov and Sporns, 2009), brain networks with isolated nodes seem to be biologically implausible. Note that to address this issue some authors decide to restrict the range of network costs assessed in order to ensure fully connected networks without isolated nodes (Bassett et al., 2008). However, we decided to report the whole range of comparable network costs, given that even if this may lead to implausible brain networks at certain costs, the evolution of efficiencies in restrictive network costs is still giving relevant useful information about the weighted topology of the networks under study, directly associated to the backbone structure of each subject brain network. It is also important to note that although different normalization approaches might be capturing different aspects of brain network organization, compensatory effects found in normalized weighted networks (increased global and local efficiencies in IUGR group) were confirmed by the cost-corrected approach, hence minimizing the chance of the reorganizational patterns being dependent of the variability in network density. In addition, we acknowledge that in order to support some of the results, especially the compensatory effect found after normalization of brain networks, it would be very advisable to report the evolution of the network features at previous ages, being of special interest the neonatal period. If the rabbit MRI atlas is successfully adapted to neonatal age, research on the brain networks of neonatal population of IUGR rabbit model will be the likely focus of future studies. Finally, we acknowledge that the relatively reduced sample size used might be reducing the statistical power of the comparisons between cases and controls, preventing to generalize some of the findings obtained, especially concerning regional characteristics that did not withstand a FDR correction.
Finally, we would like to stress the opportunity that represents the study of different neurological conditions with brain networks from animal models. In the specific case of IUGR, animal models are crucial not only to the better characterization of the pathophysiology of the condition, which is intrinsically difficult to assess, but also to the development of reliable image biomarkers of altered neurodevelopment. Hence, the application of brain networks in animal models might allow assessing and monitoring changes after possible interventions consisting on experimental drugs or enriched environment, in line with the recent reports showing that IUGR long-term prognosis significantly improves after intensive care nursing (McAnulty et al., 2013). Further work will involve the deeper characterization of changes in brain networks organization by means of the analysis in differences in individual links between regions (Meskaldji et al., 2011;Zalesky et al., 2010), the study of the predictive power of network features using machine learning, and the assessment of the results in long-term effects produced by IUGR in humans.

CONCLUSIONS
The evidences presented here support the hypothesis that previously described neurodevelopmental changes produced by IUGR in the long-term could be associated with underlying brain reorganization. This reorganization was characterized by an impaired network infrastructure, which was accompanied by an increase on the relative organization of GFA-and FD-weighted networks. In addition, a pattern of altered regional features was identified, among which changes in cerebellar areas stand out. Furthermore, global and regional network features were associated to behavioral and cognitive tests especially designed for a rabbit model. We hypothesize that IUGR at

ACKNOWLEDGEMENTS
We would like to acknowledge Guadalupe Soria, Raul Tudela and 7T MR animal platform of IDIBAPS for the diligent and careful performance of MRI acquisitions. The Image Registration Toolkit was used under License from Ixico Ltd.

CONFLICT OF INTEREST
The authors do not have conflicts of interest to declare. Note that the blue object in the left panel corresponds to the blanket used to bring the subjects to the open field, serving also as the starting familiar point for the animal. (c) Subjects are sacrificed and brain samples are obtained. MRI acquisition consisting of 126 directions DWI and anatomical T1 volumes is performed. Brain parcellation is performed in T1 volume and Q-Ball reconstruction and tractography is performed on DWI volumes. Parcellation is registered to DWI space by means of an affine registration and brain networks are extracted. (d). Binary, GFA-weighted and FDweighted networks are obtained, as well as normalized GFA-weighted and normalized FD-weighted versions. Global and regional graph theory features are used to characterize the five brain networks obtained for each subject.

Characteristic path length and average clustering coefficient
Characteristic path length, which is inversely associated to global efficiency, was significantly increased in IUGR binary raw network (p=0.048), but significantly decreased in its normalized GFA version (p=0.017). Average clustering coefficient, which is closely related to local efficiency, was also significantly increased in IUGR normalized networks in its FD version (p=0.016) and had a tendency towards significance in its normalized GFA-weighted version (p=0.072).

Figure S1. Characteristic path length and average clustering of IUGR and controls
This approach lead to a binary connected network with a network density of 22%. Multiplying this binary network for each subject GFA-or FD-weighted network we obtained individual networks with the same network cost and the same position of connections, only changing the individual weight of each connection.
Assessing the differences between cases and controls, we observed a tendency of IUGR to have a decreased FD-weighted average strength and global efficiency (p=0.079 and p=0.096 respectively). Note that with this approach we are assessing the differences in the average connection weight of the network most plausible connections (i.e. the network "backbone"). This result is coherent with previously reported tendency of IUGR to have decreased FD-weighted average strength in their raw networks, supporting the hypothesis that there is a reduced infrastructure in the "backbone" network of these subjects independently of their network density. Varying the threshold to generate this backbone binary network at different levels of network density we observe similar results than obtained for the maximum set of networks at 22% of network density (see Figure S2). Figure S2. Network cost, average strength, global efficiency and local efficiency of the network "backbone" compared between cases (red) and controls (blue) as a function of GFA-and FDweighted connectivity. + p<0.1. a 10 controls and 9 IUGR. b 8 controls and 6 IUGR. Highlighted in bold those correlations statistically significant (p<0.05), in italics those with a tendency towards significance (p<0.1).