Accompanying code and materials for Griffiths et al. The Spherical Harmonic Structure of the Human Connectome

Code

coming imminently...

Paper Summary

Abstract

The most salient gross anatomical feature of the human brain is the cerebral hemispheres, a pair of large spheroidal structures that sit atop the evolutionarily older subcortical systems. Here we demonstrate that this spherical topology of the cerebral cortex shapes and constrains brain dynamics, using a combination of cortical shape analysis, diffusion-weighted MRI tractography, and source-localized magnetoencephalography. We find that the eigenvectors of anatomical and functional connectivity closely resemble spherical harmonics, signifying the tight mathematical constraints on brain organization imposed by geometry. Through analytic theory and numerical simulations, we show that these phenomena emerge naturally from the spherical embedding, spatial statistics, and linear dynamics of cortico-cortical connectivity. Our findings offer a deeper theoretical foundation for understanding the geometric principles of brain organization and function.

Figures

Figure 1: Spherical embedding of cerebral connectivity. Each cerebral hemisphere is a complex network with a distinctive geometric and topological structure. Nodes of the network lie on a (2D) manifold in 3D space - the approximately spheroidal surface defined by the cortical grey matter ribbon. The left side of the central image - a coronal slice of a T1-weighted MR volume with grey matter, white matter, and cerebrospinal fluid coloured dark grey, light grey, and black respectively - shows this clearly, with the surface boundary between grey and white matter additionally highlighted in bright white. The cortical surface is well-known to be topologically equivalent to a sphere, as indicated by the cortical inflation diagrams on the left. Unlike its surface-embedded nodes, the edges of the cortical network occupy and traverse the full 3D space on the interior of the grey matter surface - the cortical white matter. Although the white matter appears as an undifferentiated volume on the T1 slice on the left, its complement on the right shows that this tissue in fact comprises a rich and complex set of anatomical fibre pathways, whose 3D trajectories are indicated by the diffusion-weighted MRI (dw-MRI) tractography streamlines in purple. The 3D volume comprising the network edges is fully enclosed by the 2D manifold (the cortical surface) comprising the network nodes. The phenomena studied in this paper are a direct consequences of this special anatomo-geometric organization.

Figure 2: Shared spherical harmonic patterns in cortical geometry, structural, and functional connectivity eigenmodes. A) Spatial topographies of the first 10 eigenvectors of $\mathbf{\hat{\nabla}}_{Sph}^2$, $\mathbf{\hat{\nabla}}_{Pia}^2$, $\mathbf{C}_a$, and $\mathbf{C}_f$. Colours are the normalized eigenvector values at each point on the hemispheric surface, and can be understood as the basis functions or component modes for global neural activity patterns. The prototypical SH-like sequence of topographies is exemplified in $\hat{\nabla}_{Sph}^2$ and $\hat{\nabla}_{Pia}^2$, and followed remarkably closely by $\mathbf{C}_a$ and $\mathbf{C}_f$. B) Corresponding order-summed coefficients ($P_{\ell}$s) from SH decompositions of the $\mathbf{\hat{\nabla}}_{Pia}^2$, $\mathbf{\hat{\nabla}}_{Pia}^2$, $\mathbf{C}_a$, and $\mathbf{C}_f$ modes. For each mode, the five degree levels are indicated by bar colours (blue/orange/green/red/purple for $P_{\ell}$=0/1/2/3/4 respectively), and coefficient values are indicated by bar height. The main bar plots show patterns for the first 10 modes, and the line plot insets on the top right of each row show the patterns extending out to the 30th mode. As with the topographies in panel A, here the $P_{\ell} \mathbf{\hat{\nabla}^2}_{Sph}$ coefficient patterns in the first row define the canonical SH-like pattern. This pattern is a sequence of dominant coefficients at progressively higher degree levels, where the number of orthogonal SH basis functions at each level (1 blue, 3 orange, 5 green, 7 red, 9 purple) is also indicated by the spherical SH patterns at the top right. As expected, $P_{\ell}\mathbf{\hat{\nabla}^2}_{Sph}$ follows the canonical pattern perfectly and $P_{\ell} \mathbf{\hat{\nabla}^2}_{Pia}$ follows it very closely; with our key empirical result being that $P_{\ell}\mathbf{C}_a$ and $P_{\ell} \mathbf{C}_f$ also follow it well for the first 5-10 modes, and follow it approximately for the subsequent 20. C) Diagrammatic summary of orthogonal plane orientations defined by the zero crossings of modes 2-4 in axial/sagittal/coronal views (left) and cortical flat map representation (right). The node line overlays on the flat map are computed directly from the zero values of the $\mathbf{\hat{\nabla}}_{Pia}^2$ mode maps. D) Box-and-whisker representations of $M_C$ null distributions generated from 1000 random permutations of the data matrices for $\mathbf{\hat{\nabla}}_{Pia}^2$, $\mathbf{\hat{\nabla}}_{Pia}^2$, $\mathbf{C}_a$, $\mathbf{C}_f$, $\mathbf{D}^{Inv}_{Pia}$, and $\mathbf{D}^{Inv}_{Sph}$. On each row, the black x indicates the observed value of $M_C$ in the empirical data, which lie substantially and significantly beyond the upper confidence intervals of the surrogate distributions, and are thus highly unlikely to have occurred by chance.

Figure 3: Inverse and exponential distance-connectivity relationships yield SH eigenmode patterns. Considering the idealized case (A) of a network embedded on a spherical cortical surface: When connectivity weights are given by inverse distance (B), then when taking the matrix eigenvectors and plotting on the spherical surface, or summarizing with order-summed SH coefficients (C), SH patterns will always be seen, irrespective of the spatial distances and or other linear scaling. This is because (B; Eqs. 3-4) the spherical surface inverse-distance matrix is a numerical evaluation of the integral operator with inverse-distance Green function kernel, and this operator has the same eigenvectors as its inverse, the spherical Laplacian operator. When connectivity weights are instead given by an exponential function of distance (D)), the Laplacian structure still holds approximately for this regime (c.f. Eq. 5), but SH eigenvectors are only observed numerically over a restricted range of distances and decay exponents. Panel E shows $\mathbf{M}_C$ for 100 combinations of the exponential decay constant $b$ and brain diameter $d$, half of which show strong SH structure ($\mathbf{M}_C$ > 0.8) and half show weak SH structure ($\mathbf{M}_C$ < 0.6). Example cases from these weak and strong SH regimes are marked respectively with green star and red cross, and shown in full (distance/weight functions, connectivity matrices, mode topographies, $P_l$ profiles) through panels D and F. The degradation in SH structure is manifest in this case as an intrusion of higher-frequency (3 node line) mode patterns at the 5th and 6th modes (dotted boxes).

Figure 4: Spherical harmonic eigenmodes `propagate up' from structure to activity in oscillatory neural mass model simulations. A) Schematic of how the model is constructed and used to generate simulated neural activity. Again, we focus on intra-hemispheric cortical-only connectivity, to allow direct comparison with cortical geometry. Each region in the parcellation is represented as a node in a dynamic network model, where the weights in the network are given by same group-mean normalized tractography streamline counts anatomical connectivity matrix $\mathbf{C}_a$ discussed in previous sections. The dynamics at each node are given by a nonlinear second order differential equation (Fizhugh-Nagumo model; Equations 7-8). B) By driving all nodes with a stochastic white noise process, we generate simulated resting state regional neural activity time series. From these, we computed band-limited power AEC matrices, analogous to the empirical MEG data studied above. These simulated AEC matrices ($\mathbf{C}^{Sim}_f$) are then entered into the analysis pipeline outlined in Fig. S1, yielding spatial eigenvector patterns, SH coefficient matrices, $P_{\ell}$ s, and the matching metric $M_C$. C) Left and right columns show AEC calculations for two different bandpass filters - one centred on the peak oscillatory frequency (40Hz), and one centred far away from the dominant frequency (8-12Hz). As is clear, for the 35-45Hz AEC FC, the same SH structure described in Figure 2 is also observed in these synthetic AEC eigenvectors. This is not the case, however, for AEC at 8-12Hz, which appears largely random in both the connectivity matrix and the $P_\ell$ coefficients. This shows that the SH spatial patterns of interest are not merely due to correlations between background stochastic noise processes, but reflect frequency-specific oscillatory coupling.

Suppementary Figure 1: Analysis workflow. A) The principal data used were i) spherical surface and ii) cortical (pial) surface geometry (both as represented by the numerically evaluated discrete Laplace-Beltrami Operator, $\mathbf{\hat{\nabla}}_{Sph}^2$ and $\mathbf{\hat{\nabla}}_{Pia}^2$, iii) normalized DW-MRI tractography streamline counts, an estimate of anatomical connectivity ($AC$) strength, iv) empirical and v) simulated MEG functional connectivity ($FC$). Each of these was processed using standard tools (see Methods) to yield either a parcel-by-parcel or vertex-by-vertex adjacency matrix $\mathbf{Y}$. B) Eigendecomposition of $\mathbf{Y}$ yields a spatial pattern of loadings across brain regions for each eigenvector $\mathbf{u}_N$. C) These spatial patterns are mapped on to cortical surface vertices using the surface parcellation. The correspondence of Freesurfer pial and spherical surfaces allows the re-representation of each spatial pattern $\mathbf{u}_N$ in a 2D spherical coordinate system with coordinates $\theta,\phi$. This in turn allows us to compute SH coefficients for each eigenvector pattern. SH decomposition up to degree $\ell$ and order $m$ yields an ($\ell$+1) x $m$ matrix of coefficients with 2 $\ell$ +1 nonzero values on each row. Summing each row gives the five ($\ell$+1) order-summed coefficients $P_{\ell}$. D) The pattern of these coefficients for a given set of eigenvectors is the principal quantity of interest in the present study. Presence of SH-like patterns in these coefficients is then quantified using $M_C$ and $M_{Dice}$, and permutation methods used to compare these numbers against synthetic null distributions.