Preterm birth affects the developmental synergy between cortical folding and cortical connectivity observed on multimodal MRI

The survival rates of infants born prematurely have improved as a result of advances in neonatal care, although there remains an increased risk of subsequent disability. Accurate measurement of the shape and appearance of the very preterm brain at term-equivalent age may guide the development of predictive biomarkers of neurological outcome. We demonstrate in 92 preterm infants (born at an average gestational age of 27.0±2.7weeks) scanned at term equivalent age (scanned at 40.4±1.74weeks) that the cortical sulcation ratio varies spatially over the cortical surface at term equivalent age and correlates significantly with gestational age at birth (r=0.49,p<0.0001). In the underlying white matter, fractional anisotropy of local white matter regions correlated significantly with gestational age at birth at term equivalent age (for the genu of the corpus callosum r=0.26,p=0.02 and for the splenium r=0.52,p<0.001) and in addition the fractional anisotropy in these local regions varies according to location. Finally, we demonstrate that connectivity measurements from tractography correlate significantly and specifically with the sulcation ratio of the overlying cortical surface at term equivalent age in a subgroup of 20 infants (r={0.67,0.61,0.86}, p={0.004,0.01,0.00002}) for tract systems emanating from the left and right corticospinal tracts and the corpus callosum respectively). Combined, these results suggest a close relationship between the cortical surface phenotype and underlying white matter structure assessed by diffusion weighted MRI. The spatial surface pattern may allow inference on the connectivity and developmental trajectory of the underlying white matter complementary to diffusion imaging and this result may guide the development of biomarkers of functional outcome.


Introduction
During the last trimester of pregnancy there are substantial changes in the appearance and connectivity of the fetal brain. During this short period, the cortex develops from a lisencephalic state with rudimentary associative connections between regions and dramatically increases in volume and surface area (Kapellou et al., 2006) to one that resembles the adult cortex in both appearance and structural connectivity in all but scale (Dubois et al., 2008;Partridge et al., 2004). In babies born very prematurely much of this structural development occurs entirely in the postnatal period, during exposure to the relatively harsh extrauterine environment of the neonatal intensive care unit. Although many babies born very preterm will go on to develop normally, a significant number will have neurological developmental difficulties that persist to childhood and beyond (Doyle and Saigal, 2009;Marlow et al., 2005) and despite advances in neonatal intensive care over the past decade the rates of disability and neurodevelopmental sequelae for infants that survive extreme preterm birth have remained disappointingly unchanged (Moore et al., 2012). As a consequence there is much interest in the relationships between very preterm birth, gross brain injury and focal cellular injury to oligodendroglial precursors and their influence on later neurological outcome (Volpe, 2009). A significant challenge researchers face is to understand how the brains of babies born prematurely differ at term equivalent age with those of term-born infants.
In the second half of gestation, cortical growth proceeds extremely rapidly and with a well-defined organisation. Measurements from histology have been used to assess how the cortical folding pattern and internal pattern of myelination change in late gestation and infancy. Unmyelinated projective and commissural fibres are in place early in brain development (Kolasinski et al., 2013;Kostovic and Jovanov-Milosevic, 2006), as are the positions of major sulci (Chi et al., 1977). Secondary folds and connectivity appear from around 24 weeks followed by the appearance of associative fibres and tertiary folding in the late preterm period. By term equivalent age tertiary folding is advanced and associative fibres are present, resembling the adult morphology (Zilles et al., 1988).
The relationship of cortical folding to underlying white matter connectivity is theorised to be one of improving neuro-communication efficiency (Cherniak et al., 2004) whilst establishing a functional parcellation (Fischl et al., 2008), but the exact mechanism has generated some peculiar controversy over the past two decades. Fundamental to this controversy is the difficulty of establishing the actual mechanisms that give rise to the canonical human brain phenotype, since this requires simultaneous in vivo measurement of both cortical folding and white matter composition and shape. (Van Essen, 1997) proposed that cortical folding may be explained by a process of axonal tension between sulci and that the strength of the axonal connections would thus predict sulcal depth. There are a number of variations on this theory, adjusting the relative weights of genetics, cellular proliferation and the axonal contribution within the confined boundary of the skull. Conflicting evidence for the role of mechanical tension has been produced by a number of authors based on both theoretical and diverse experimental results. Some specific examples are gyrogenesis purely as a result of differential cellular proliferation (Kriegstein et al., 2006;Toro and Burnod, 2005); gyrification occurring, but not as predicted by a pattern of axonal tension in the developing ferret brain (Xu et al., 2010); whilst (Nie, 2012) proposed a model of axonal pushing in the formation of gyri rather than sulci, reinforced with a rather tangential validation in Drosophila. Unsurprisingly, the role of genetics has also been proposed as a major determinant of appearance; (Im et al., 2011) found that identical twins have more similar cortical folding patterns than non-identical twins and specific genes related to variants of polymicrogyria and to brain morphology differences in schizophrenia have been found (Gregorio et al., 2009;Rakic, 2004). However, not all disorders of brain development are related to genetic factors (Francis et al., 2006). In addition, indirect factors such as gender (Gilmore et al., 2007b), diseases such as congenital heart disease (Awate et al., 2010) and preterm birth may also play a role in the eventual cortical appearance and function.
The outward appearance of the brain on Magnetic Resonance Imaging (MRI) has been used to reveal the gyrification pattern of the cortex. Techniques for extracting the cortical surface have previously been proposed: histology based methods frequently use a gyrification index (GI) defined on slices measuring the ratio of the length of the grey matter to that of the shortest route around the coronal surface. (Zilles et al., 1988) used this method to profile the adult cortex anterior to posterior and generate a sectional profile of the human cortical folding pattern. Automatic methods are now popular: using a semi-automatic segmentation, (Dubois et al., 2008) delineated the cortical surface by a surface mesh and analysed the appearance and formation of major sulci over 26-36 weeks gestational age making use of a sulcation ratio in relation to total cortical area. Their data showed similar relations to the allometric relationships described in (Xue et al., 2007) and (Kapellou et al., 2006) with a rapid increase in cortical folding over this period of development. (Xue et al., 2007) proposed an implicit surface formulation for cortical surface extraction, that has the advantage that the local curvature does not require an additional triangulation procedure. From these surfaces it is possible to derive a number of summary descriptions of the surface features, although none have achieved the status of a de facto standard (Batchelor et al., 2002;Yotter et al., 2011). The implicit surface formulation was further adopted by (Awate et al., 2010) who developed summary statistics of the cortical surface by collecting complementary values of curvature into two-dimensional histograms. The representation of the surface in this technique is such that their histograms have two visible peaks distinguishing separately the properties of gyri and sulci; thus they were able to immediately assess their work using a statistical mapping approach. The authors found significant differences in folding pattern between hemispheres and major lobular regions in a group of infants with congenital heart disease. In additional, measurements of underlying sub-cortical grey matter volumes have shown significantly reduced volumes in preterm at term equivalent age infants relative to control term age infants (Boardman et al., 2010;Cardoso et al., 2013).
Further to the controversy on the links between cortical folding and white matter connectivity, the development and presence of myelin has itself been the source of debate. Development of axonal myelin sheaths is fundamental to increased speed and efficiency of electrical activity in the brain, and myelin related disorders are associated with a wide range of conditions. In the developing brain, processes leading to myelination progress in an established spatial pattern, ascending into the corticospinal tracts from as early as 30 weeks gestation and progressing from this region along the central sulcus and anterior and posterior over the first few months of life, with the posterior myelination occurring more rapidly (Brody et al., 1987;Kinney et al., 1988). On MRI, the absence of myelin is a major contributor to the inverted contrast, relative to the adult brain, seen on T 1 and T 2 weighted images and thus there is interest in using MR measurement of myelin as a biomarker of later neurological outcome. Some developmental changes can be observed on Diffusion Weighted Imaging (DWI); in grey matter the increasing cortical connectivity between 30 and 40 weeks gestational age reduces the observed diffusion anisotropy in a characteristic pattern which is thought to represent loss of early-development radial glial scaffolding and replacement by cortical neuronal proliferation and synapse formation (McKinstry et al., 2002). Debate about myelination arises since changes in the apparent connectivity and T 2 shortening can be inferred from DWI prior to histologically visible myelin; which manifests as increased anisotropy and decreased free diffusion. Thus there have been conflicting reports on whether any form of DWI can be used to specifically measure myelin content or change. Diffusion changes are likely to be related to increased structural organisation, axonal development and possibly ensheathment of axons by pre-oligodendrocytes within the developing brain (Volpe, 2009) and thus will have a non-specific relationship to myelin content (Hagmann et al., 2009;Huppi and Dubois, 2006). Nonetheless, changes to the cellular content of white matter have been used to infer the presence of myelin: during the preterm period; at term equivalent age and to follow myelination through infancy (Deoni et al., 2011;Hagmann et al., 2009;Partridge et al., 2004). If myelin content can be accurately measured; delays in the myelination process might be detectable in preterm infants at term equivalent age.
Despite the relative lack of specificity to myelin content, DWI is widely used to assess axonal disease and injury (Ulug et al., 1999) and changes have been observed to be dependent on gestational age at birth in both the corticospinal tracts and corpus callosum with reduced connectivity in both regions Thompson et al., 2011). Furthermore, DWI measures of connectivity have shown a reduction in thalamo-cortical connectivity assessed by diffusion MRI (Ball et al., 2012). Although correlations between indices of gyrification and the observed diffusivity of corresponding grey matter have been observed (Adams et al., 2010;Deipolyi et al., 2005), relationships between the surface pattern and the underlying white matter diffusivity profile have not been established. More recently, work has shown some of the links between white matter structure and functional performance: (Quairiaux et al., 2011) demonstrated how cortical maturation corresponds with whisker control in a murine model whilst (Northam et al., 2012) showed relationships specifically between inter-hemispheric temporal lobe connectivity and language performance in ex-preterm adolescents. Relationships of this type would themselves guide theories of gyrification, which would be crucial for teasing apart the correlative synergy between mechanical, genetic and environmental factors in the relationship between internal connectivity and outward shape and appearance. Fundamentally these relationships must come from preterm and neonatal cohorts in which increases in cortical folding and white matter connectivity are ongoing. Here we hypothesise that the number of white matter computed connections to the cortex by diffusion weighted tractography correlates with the cortical folding pattern of the cortex and we combine an adaptive neonate specific segmentation routine with a level-set based surface analysis technique to investigate correlations between the regional white matter connectivity and the overlying cortical folding pattern.

Participants
Ninety-two infants are included in this study, details of the cohort are found in Table 1. The local ethics committee granted permission for this study, and informed parental consent was obtained for each infant. The ethics committee provided approval to administer sedation during the MRI scan.

Magnetic resonance image acquisition
Infants were sedated with an oral dose of chloral hydrate (Rosemont Pharmaceuticals, Leeds, UK) as agreed by the local ethics committee and imaged within a transparent MR-compatible pod. T1-weighted data were acquired on a 1.5 T Siemens's Avanto using TR = 17 ms, TE = 6 ms and flip angle of 21°and resolution of 0.39 × 0.39 × 1 mm.
Eighty-two of the above 92 subjects had additional diffusionweighted data with either 6 (62 infants) or 30 (20 infants) direction diffusion weightings as a result of a change in the diffusion imaging protocol (Table 1). The 6 directions were acquired at b = 600 s. mm − 2 with 6 repeats in addition to six b = 0 s. mm − 2 images at 0.9 × 0.9 × 3 mm resolution. 2x30 directions at b = 600 s. mm −2 were acquired in addition to two b = 0 s. mm −2 images at 0.9 × 0.9 × 3 mm resolution, although one diffusion dataset was removed from the analysis due to severe motion artefacts (see Table 1)). There were no significant fractional anisotropy (FA) differences between the two subgroups of 6 directional and 30 directional acquisitions (For the genu of the corpus callosum; FA 6 = 0.22 ± 0.04, FA 30 = 0.21 ± 0.05 (p = 0.076) and posterior white matter; FA 6 = 0.11 ± 0.034, FA 30 = 0.10 ± 0.039 (p = 0.088)).

White matter injury scoring
Conventional MRI was assessed for white matter (WM) abnormality using an established qualitative scoring system based on a grading system of 5, 3-point scales detailed in (Woodward et al., 2006). White matter abnormality is assessed by: the type and distribution of whitematter signal abnormality; loss of volume of (periventricular) white matter; the presence and location of any cystic abnormalities; ventricular dilation; or thinning of the corpus callosum. A summary white matter score is then generated using from these scales to described the white matter abnormality: none (a score of 5 to 6), mild (7 to 9), moderate (10 to 12), and severe (13 to 15). Images were scored by an experienced neuroradiologist blinded to clinical history and neurodevelopmental outcome (RG). WM abnormality was classified as normal/mild or moderate/severe. Fourteen babies were classified with moderate/severe white matter injury. The remaining 78 babies were classified as normal/mild. Neither gestational age at birth, birthweight, nor imaging subgroup was found to be significantly different between normal/mild or moderate/severe WM appearance (see Table 1).

Image segmentation
We employed a segmentation algorithm specifically developed for very preterm neonates (Cardoso et al., 2013). 1 Briefly, the method makes use of a maximum a priori expectation maximisation routine with both iteratively relaxed spatial priors and a prior term on the initial tissue intensities. Robustness to low signal and contrast to noise ratio can be enhanced via a spatial smoothness term by means of a Markov random field (Van Leemput et al., 1999) and inclusion of MRI intensity non-uniformity correction. In addition, due to anatomical variability we assumed that the relative weighting of the tissue priors to the data were not known a priori and we employed an iterative scheme which spatially relaxes the tissue priors towards the data, facilitating the segmentation of pathological cases similar to the technique in (Shiee et al., 2011). As in (Cardoso et al., 2013), age-specific tissue priors are obtained from a publicly available atlas of preterm neonates as described in (Kuklisova-Murgasova et al., 2011). The tissue priors were brought into correspondence using an affine registration (Ourselin et al., 2000) followed by a free-form non-rigid registration algorithm (Modat et al., 2009). Neonate specific partial volume correction was carried out as a final step in order to automatically reduce the a priori probability for partial volume containing voxels to belong to white matter when the probability of being grey and white matter is also high. The segmentation protocol proceeded using the following pipeline: 1. Tissue intensity priors were registered to the individual subject using affine and non-rigid free-form deformations (Modat et al., 2009), the registration result was then used to define the subject-specific brain mask and the spatial priors. 2. Segmentation of six major tissue classes was carried out using a maximum a priori expectation-maximisation strategy (white matter, cortical grey matter, subcortical grey matter (which also included the small amount of myelinated white matter present at term), total cerebrospinal fluid (CSF) space, cerebellum and brainstem) using the manual segmentation priors. This algorithm included intensity non-uniformity correction, spatial smoothness introduced via a Markov Random Field and estimation of an outlier class enabling the rejection of intensity clusters that have a large Mahalanobis distance from the estimated model, thus reducing their influence in the parameter estimation. 3. Explicit white-matter partial volume correction was carried out and the white matter prior adjusted as in (Cardoso et al., 2013). Briefly, the method downweights the prior probability of white matter (WM prior ) in the presence of high grey matter (GM) or CSF probability (p GM and p CSF ) estimated from the segmentation step in 2 using the update: 4. Segmentation of the six major classes in step 2 was repeated including the adjusted white matter prior of step 3 yielding a final probabilistic segmentation.

Cortical surface analysis
The segmentation method allowed an inspection of the shape of the boundary between grey and white matter surfaces. In our study we found the boundary between grey and white matter using a level set (Eq. (1)); solving the Eikonal equation via a fast marching technique for the function ϕ(x,t) by weighting a signed distance function from the boundary with λ d and λ c weighting the relative strengths of the speed of normal front propagation (d(x,t)) and the constraint on the maximum mean curvature (c(x,t)) as in (Xue et al., 2007). The level set was initialised along the white matter boundary where p(WM) = 0.5 and resulted in a signed distance function from the WM/GM boundary. 2 Although topology preserving strategies can be included in the level-set formulation we do not do so here; assuming that both the segmentation results (Cardoso et al., 2013) are good as a result of the MRF and partial volume correction strategy and that selecting the WM/GM boundary reduces the chances of a broken topology.
We defined measures of curvature on this implicit surface, similar to the method described and used in (Awate et al., 2010). The local Hessian matrix of second order derivatives can be found at each point of this implicit surface relative to the orientation of this surface. The resulting eigensystem of eigenvalues κ n and eigenvectors ν n may be solved to yield at most two non-zero eigenvalues (Eq. (2)). Practically, cortical curvatures are calculated as finite difference estimates at all points in the volume and since the level set is smoothly varying and the finite difference kernel relatively large, this should give a good representation of the local curvature at all points. Representation of the surface curvature in the interpolated image domain uses a zero-crossing image defined according to whether there is a sign-change between a voxel and its adjacent neighbours. Since this interpolation is of smaller scale than the finite difference curvature calculation, the representation of the surface curvature on these discrete points should be of justifiable accuracy.
Hence, at every point on the surface we can use the eigenvalues of the local Hessian matrix, κ 1 and κ 2 (with κ 1 N κ 2 ) to summarise the local shape. Explicitly we define the shape index, S (Eq. (3)), describing how cup-like or saddle-like the surface is, and a curvature, C (Eq. (4)), describing the distortion of the surface relative to that of a flat sheet (Koenderink and van Doorn, 1992).
Intuitively C and S are complementary in the sense that they are a polar representation of the eigenvalue sum (κ 1 + κ 2 ) and difference (κ 2 − κ 1 ). The value of C is corrected for brain volume using the strategy described in (Awate et al., 2010), where each value of C is divided by the cube root of the ratio of the subject brain volume to the mean dataset brain volume. The use of S and C is motivated by the separation of gyral and sulcal regions into sign-separated clusters as in (Awate et al., 2010) when the values are summarised in a 2D histogram, one that is filled by an adaptive kernel density estimation procedure (Botev et al., 2010). The sulcation ratio (SR) is defined as the ratio of the number of surface points with positive shape index to the total number of surface points (as illustrated in Fig. 1).
Furthermore we may define the outer boundary of the grey matter by initialising the level set at the grey and white matter interface and evolving this initial surface outward to the p(WM) + p(GM) = 0.5 boundary, maintaining correspondence of the inner white matter boundary prior to parcellation. Thus we have defined two surfaces that can be used as an approximation to the histological gyrification index (Zilles et al., 1988) by dividing one by the other, although alternative methods are also available (Batchelor et al., 2002).

Definition of lobe-based regions of interest
To analyse the local cortical curvature properties we propagated a map of anatomical labels to each segmentation (Oishi et al., 2011) and found the intersection of this with the GM tissue class. Although gyrification is well into the tertiary phase at term equivalent age, detailed anatomical information may be poorly defined due to the distance of the initial single-subject term neonate atlas described in (Oishi et al., 2011) from the preterm infant population, so here we restricted the parcellation to the grey matter surfaces of major lobes (the frontal, parietal, temporal and occipital lobes) and analysed the surface joint histograms in these regions on both left and right hemispheres. To define the lobe-based parcellation we take the existing 122 labels defined in the atlas and group these into regions corresponding to their associated lobe. In order to align the atlas, we propagated the label map to each individual segmentation using a deformation field determined in combination with the segmentation step (Cardoso et al., 2013). We take the overlapping region of the label map and the underlying grey matter segmentation to define the lobe extent when calculating the SR and GI (Figs. 2 and 1).

Linking DTI and local cortical properties
For those infants with diffusion data we calculated local values of the fractional anisotropy and mean diffusivity in regions of interest in the left and right, anterior and posterior white matter. Two-dimensional regions were manually positioned on mean diffusivity maps with an 11 × 11 pixel circular kernel in the anterior and posterior white matter (see Fig. 4a top and bottom). This kernel size was chosen to be small enough to fit entirely within the supra-ventricular white matter whilst allowing tissue averaged parameters to be found. In addition, regions were selected in both the genu and splenium of the corpus callosum (CC) using the same circular kernel. Average values of the mean diffusivity (MD) and FA were calculated for each region for comparison with cortical folding indices derived from the adjacent lobar region. For the DWI subset with higher angular resolution (30 diffusion directions), we correlated diffusion properties with cortical folding properties by registration of the diffusion space images to the structural T1-weighted image via the DWI MD image using asymmetric affine transformation 3 We made use of this deformation to further propagate subsequent tractography volumes to the T1 space. Tractography volumes were formed by defining three seed regions within the posterior limbs of the internal capsule and CC of each of these 19 infants (see Fig. 4). These seed regions were subsequently used as starting regions for the probabilistic tractography algorithm (Behrens et al., 2003). In addition, we constrained the tractography with termination masks for each corticospinal tract (CST) on the superior longitudinal fasciculi anterior and posterior of the central sulcus. Although a direct biological interpretation remains difficult (Jones et al., 2013), the probabilistic tractography algorithm nonetheless obtains a connected pattern representing the ease at which the cortex can be reached from each of these seed regions, either by progressing parallel to coherent white matter regions, or via a stochastic diffusive process in regions of lower FA. We use the term connectivity to describe the ease that areas outside the seed region can be reached as defined by the streamline count in (Behrens et al., 2007), but do not correct for distance from the seed point in order to reduce unexpected bias in the result interpretation.
The intersection of the tractography volumes with the cortex in the space of the T1 volumes was subsequently used to define the regions of interest for analysis, including the number of surface voxels in the region of interest, the number of tract counts reaching the region of interest and calculation of the regional SR. Thus each infant had four regions of interest combining both diffusion and cortical folding information  (see Fig. 3). Furthermore we investigated the effect of thresholding the connectivity to reduce the impact of noise from weakly connected regions (Section 3.2.2).

Spatial trends in the cortical folding pattern
Population trends in the gyrification index Fig. 5 shows trends in the overall GI from posterior to anterior. The one-dimensional GI data are resampled onto a standard onedimensional space between the anterior and posterior poles (defined as those regions for which the GI falls below 1). The labelled limits of each lobe are defined by the average position of the anterior-most and posterior-most labelled voxels in each region of interest. High GI suggests increased complexity in the posterior occipital lobe, falling off towards the operculum and rising again in the frontal lobe. Results along the medial-lateral axis are highly symmetric from left-side though to right-side.

Gestational age dependent features of the cortical surface
Overall SR correlates significantly with gestational age at birth on image data acquired at term equivalent age (r = 0.49, p b 0.0001, normal/mild group only (r = 0.38, p = 0.0001, for all infants); Fig. 6a), correcting for age at scan. For the normal/mild group, left hemisphere SR values are higher than the right (difference in means = 0.031 (0.026 − 0.036 95 % confidence interval (ci))). The regional correlation of SR with gestational age at birth also varies consistently with gestational age at birth with the exception of the left occipital lobe (Fig. 6b), symmetric along the medial-lateral axis (Table 2).
After correction for gestation, occipital lobar SR values were higher than all other lobes and temporal SR values were lower than all other lobes. Values for frontal and parietal regions were intermediate and similar. Regional SR values were similarly ordered between left and right sides. In detail, after adjustment for gestational age at birth, the occipital lobes are found to have higher SR than all other lobes (p b 0.001 for all lobes). Conversely, the SRs of the temporal lobes are found to be lower than all other lobes (p b 0.001 for all lobes). Frontal and parietal lobes are found to be non-significantly different (F R − P R falls between − 0.02 − 0.011 (95 % ci), p = 0.15) for left frontal to left parietal and borderline significant for right frontal to right parietal (F L − P L falls between 0.0003 − 0.016 (95 % ci), p = 0.039).

Correlations with white matter DTI
In 73 of the 82 infants with both diffusion data and classified in the normal/mild group, we inspect variations in the white matter FA and SR. FA is found to be specifically lower in the anterior (0.11 ± 0.04) and posterior (0.10 ± 0.03) white matter relative to both the genu (0.22 ± 0.04) and splenium (0.25 ± 0.06) of the CC, correcting for gestational age at scan. A positive correlation with gestational age at birth was found for genu (r = 0.26, p = 0.02) and for splenium (r = 0.52, p b 0.001) in contrast to a lack of relationship in anterior and posterior WM (r = − 0.02, p = 0.83 and r = 0.01, p = 0.87 respectively). Conversely mean diffusivity in genu and splenium was negatively correlated with gestation (r = − 0.33, p = 0.002 and r = − 0.27, p = 0.013, respectively). SR of individual lobar regions was not significantly correlated with regional FA or MD controlling for gestational age at birth (results not shown).

Correlations with white matter tractography
In the 20 babies with 30-direction diffusion data, we investigated three tract systems: left and right CST emanating from the posterior limb of the internal capsule and the tract system originating from the CC (Fig. 4b). Tract connectivities are all corrected by the number of seed points used in their generation. All results are shown corrected for gestational age at birth, gestational age at scan, birthweight and brain volume defined as the sum of the grey matter and white matter volumes determined by segmentation (in addition, the number of connections reaching the cortex for each system correlate with gestational age at birth as follows: CC: r = 0.35, p = 0.16, Left CST: r = 0.11, p = 0.67, right CST: r = 0.48, p = 0.04 and for birthweight: CC: r = 0.40, p = 0.10, left CST: r = 0.13, p = 0.59, right CST: r = 0.35, p = 0.15). Lastly, the tract-average FA was not found to correlate significantly with gestational age at birth or birthweight.
Figs. 8a,b,c plot the relationships between the number of connections to the cortex found between each tract system (thresholded with at a connectivity value of 50), the SR found on the connected region of cortical surface and the average FA of this tract system respectively. Figs. 8d,e,f plots the inter-relationships between the regioncortex tractography and the overlying surface properties: correlating, for each system, the SR of the connected surface (SR), the number of unique points on the surface (Cn) and the number of times the surface is reached by tractography (Cx) against one another. Table 3 column one shows the correlation coefficients found for the corresponding graphs in Fig. 8 between the number of connections to the cortex found for each tract system thresholded with every other at a connectivity value of 50 (see Fig. 9). Only the correlation between the CC and the left CST reaches significance (Cx vs. Cx). Table 3 column two shows the correlations between the SR calculated from the intersection of each tract system with the cortical surface: left and right CST cortical SRs show a significant correlation. The SR ratio found between the CC and right and left CST also show significant correlation. Table 3 column three compares the average tract FA values revealing a high similarity between left and right CST but not between CST and the CC. FA is found to be significantly higher in the CST than in the CC ((average FA = 0.18 in the CC and FA = 0.24 in the CST, thresholding at a connectivity index of 1000) see Fig. 7). Table 3 column four shows correlations coefficients between the number of connections to the cortex and the associated SR of that connected surface region. Tract connectivity values are shown to correlate specifically with the SR ratio of the connected surface; correlations between tract systems and unconnected surfaces do not correlate significantly. Table 3 columns five  and six show correlations between the number of unique surface voxels reached by tractography and the SR; and correlations between the number of unique surface voxels reached against the tract connectivity. To inspect the relationship of the correlation coefficients to one-another, they may be modified using the Fisher transformation z = tanh −1 (r) to make the distribution of correlation coefficients more normal. Carrying out this procedure does not reveal any significant differences between the results in columns 1 and 2 and of Table 3 for an equivalent p b 0.05. For Column 3, the FA value between left and right CST is significantly higher than those of the left and right CST to CC values. For Column 4 only the CC result remains significantly different from the mixed tract/surface results whilst it is not itself different from the relationships in the CSTs. Comparable relationships hold for Columns 5 and 6 emphasising differences found in the CC-to-surface system.
In order to evaluate the sensitivity of the analysis to the connectivity threshold we repeated the above analysis with a variable connectivity threshold at values of [1,2,5,10,50,100,200,1000] since each seed region is constrained to produce at most 5000 streamline counts. Fig. 9 illustrates the relative insensitivity of the derived cortical-diffusion correlation to the connectivity threshold. Correlations between cortical folding and diffusion measurements remain stable until beyond a threshold of 100 after which the correlation coefficient is seen to decline, motivating our threshold of 50. For the number of datasets used in this section, correlations above about r = 0.49 are significant at the p = 0.05 level.

Cortical folding and preterm birth
The new findings from this study are that white matter connectivity measured using tractography correlate significantly and specifically with the sulcation ratio of the overlying cortical surface at term equivalent age in a cohort of preterm-born infants. This result and related work are likely to be commensurate with one of the menagerie of theories describing the relationship of cortical folding to underlying white matter connectivity, either by processes of axonal tension (Van Essen, 1997), cellular proliferation (Kriegstein et al., 2006), brain shape, or a combination of these, although we should emphasise that findings in a preterm-born cohort are likely to deviate from those found in normal  development. Significant results on the first three rows of Table 3 column four indicate that the method is sensitive to correlations between the calculated tract system and the associated SR whilst the weaker correlations in the final three rows suggests that the method is also specific to the associated tract system or cortical region, independent of gestational age at birth. Importantly, despite both SR and FA both showing strong dependence on gestational age at birth the correlation of the two (WM connectivity and overlying cortical surface) remains significant independent of gestational age at birth, birthweight and brain volume. The trends in sulcation ratio with lobe ( Fig. 6) and GI with slice position (Fig. 5) suggest increased folding in the posterior occipital lobe, falling off towards the operculum and rising again in the frontal lobe which can be very tentatively compared to adult histological studies (Zilles et al., 1988), although again we emphasise that this is a preterm cohort. We demonstrate that the cortical SR and GI vary spatially over the cortical surface at term equivalent age with an anterior posterior trend of occipital SR greater than frontal SR and a symmetric medial lateral trend. In addition, we demonstrate that the cortical SR correlates significantly with gestational age at birth at term equivalent age emphasising that differences in brain development brought about by very preterm birth remain detectable in a dose dependent manner at 40 weeks gestation; these results extend what has been found in previous preterm studies (Dubois et al., 2008;Gilmore et al., 2007a) to the term equivalent age cohort.
In the underlying white matter both the FA and MD of local WM regions correlated significantly with gestational age at birth at term equivalent age. In addition the FA and MD in these local regions varies dependent on location, for instance the FA in the CST is higher than the FA in the CC which itself is higher than the FA in the posterior and anterior white matter although caution should be exercised in the interpretation of the lowest values for splenium FA which may have been contaminated with CSF partial volume. The high similarity in FA values between left and right CST seen in Table 3 column 3 but not between CST and the CC may be indicative of the difference in myelination status, cell density and composition of these two regions at term-equivalent age (Deipolyi et al., 2005;Gilmore et al., 2007a). Although a number of authors have sought to link DWI to myelin-related changes, it should be emphasised that changes to the diffusion properties (such as FA) are unable to distinguish between changes to axonal number, density or myelin content. Nonetheless, it is known that at term-equivalent age the ascending CST have higher myelin content than the CC (Brody et al., 1987) and this is likely to be a contributing factor to the observed Table 2 Values for SR and corresponding standard deviation (SD) by lobe (adjusted for gestation, corrected for age at scan) and linear correlation individual lobar SR with gestational age at birth. Overall mean SR was significantly correlated with gestation (r = 0.51; p b 0.001).
FA differences. When combined, these results suggest a close relationship between the cortical surface phenotype and underlying white matter structure assessed by diffusion weighted MRI in this preterm cohort. Whether the differences observed in this work are related to specific brain insults as a result of prematurity or more subtle longerterm differences related to brain development is unclear and future studies linking the imaging to records of the first few weeks of life would have the potential to guide the establishment of links between imaging and outcome. Furthermore, although differences seen in this work are most likely due to the impact of premature birth, it would be interesting to seek similar relationships in a different cohort, for instance in term-born controls.

Methodological considerations
There are a number of limitations to this work. With regards to the cohort, we did not have ethical approval to collect and store information about the health status of the mother and thus some factors that might be related to preterm birth are unknown. Methodologically, the sulcation ratio summarises abruptly the diverse pattern seen in the cortical folding histograms of Fig. 4. Richer descriptions of the folding pattern are in principle obtainable beyond the summary statistic used in this work but nonetheless we demonstrate correlations with this simple measure. The DWI is limited by the 3 mm through-plane resolution; more advanced higher-resolution diffusion sequences will enable more accurate descriptions of the underlying white matter properties. Furthermore, the probabilistic diffusion tensor tractography algorithm used in this work (Behrens et al., 2007) is one of many, and alternatives could yield more successful representations of the white matter connection pattern. Furthermore we do not currently correct for distance from the seed point since it is not clear which would be the most appropriate unbiased strategy and since this would make a number of assumptions about the underlying space and nature of the connection pattern. The interpretation of results obtained by diffusion tractography also remains practically difficult given that the diffusion measurement is an aggregate of many possible tissue compositions and organisations; multiple configurations of white matter orientation and density can give rise to overlapping values of FA and non-representative principal diffusion direction. Nonetheless, within the confines of a strict definition of what a diffusion tractogram can be said to represent and given consensus on unambiguous terminology (Jones et al., 2013), reproducible and thus potentially useful results can be obtained pending more advanced analysis techniques. A pragmatic, mechanistic view is particularly tolerable if, for instance, studies are searching for biomarkers where a strict causative biological pathway is not necessary.

Conclusion
Despite the limitations, we have shown that relationships exist between the white matter connectivity profile and the overlying cortical folding pattern. Thus in principle, subtleties in the combined connectivity- Fig. 8. Relationships between diffusion tractogaphy from region-of-interest to cortex and the overlying surface sulcation ratio. Inter-tract relationships for a) the number of times the cortex is reached, b) the sulcation ratio of each tractography-defined surface and c) the average tract FA. Inter-parameter relationships for d) the number of times the cortex is reached against the sulcation ratio, e) the number of unique surface voxels reached against the sulcation ratio and f) the number of times the cortex is reached against the number of unique surface voxels reached.

Table 3
Correlation coefficients between tract systems seeded in the CC and left and right CST. Significant correlations (p b 0.05) between the number of connections (subscribing to a mechanistic computational interpretation) traced to the cortex (Cx), the SR of the surface that the tract system crosses, as described in Linking DTI and local cortical properties section) and the average, weighted FA found on the tract system and the number of unique surface voxels reached (Cn cortical folding pattern might be linked to clinical diagnosis or pathology, for instance, describing the biological manifestation of the MR visible grading of white matter injury, second, the methodology might allow analysis of the relative importance of projective, commissural and associative fibres on the connected local folding pattern with respect to normal human ontogeny. If the cortical surface pattern is different between groups of infants with different white matter appearance or other manifestations of prematurity, this technique might detect these differences and in turn be predictive of subsequent cortical appearance in later childhood. In particular, those infants with pathology could be further analysed for specific differences in cortical folding related to lesion location, aided by the region specific nature of the SR and association with diffusion derived indices of connectivity. The results in this work provide a framework for further investigation of the relevance, or otherwise, of variegated spatial patterns in both the preterm-atterm cortex and white matter. Further work, in combination with cohorts at later ages, will determine to what extent the combined preterm cortical folding-connectivity pattern is predictive of the emerging adult appearance and whether this phenotype persists in studies of expreterm children or adults and this may have future application to the development of cortical folding pattern based biomarkers of functional neurological outcome.

Funding
This work was supported by UK registered charity SPARKS to Andrew Melbourne (grant held by Neil Marlow). Giles Kendall is a National Institute for Health Research (NIHR) funded Clinical Lecturer.
M. Jorge Cardoso is supported by EPSRC grant EP/H046410/1. Nicola Robertson and Neil Marlow receive part funding from the Department of Health NIHR Biomedical Research Centres funding scheme at University College London Hospital and University College London. Sebastien Ourselin receives funding from the EPSRC (EP/H046410/1, EP/ J020990/1, EP/K005278), the MRC (MR/J01107X/1), the EU-FP7 project VPH-DARE@IT (FP7-ICT-2011-9-601055), the NIHR Biomedical Research Unit (Dementia) at UCL and the National Institute for Health Research University College London Hospitals Biomedical Research Centre (NIHR BRC UCLH/UCL High Impact Initiative). The authors declare no conflicts of interest.