A Continuous Model of Cortical Connectivity

A Con tin uous Mo del of Cortical Connectivit y Daniel Mo yer, Boris A. Gutman, Josh ua F ask owitz, Neda Jahanshad, and P aul M. Thompson Imaging Genetics Cen ter, Universit y of Southern California moyerd@usc.edu Abstract. W e present a con tinuous model for structural brain connec- tivit y based on the Poisson point process. The mo del treats each stream- line curve in a tractography as an observed even t in connectome space, here a pro duct space of cortical white matter b oundaries. W e approxi- mate the mo del parameter v ia kernel density estimation. T o deal with the heavy computational burden, we dev elop a fast parameter estima- tion method b y pre-computing associated Legendre pro ducts of the data, lev eraging prop erties of the spherical heat kernel. W e sho w how our ap- proac h can b e used to assess the quality of cortical parcellations with resp ect to connectivt y . W e further present empirical results that suggest the “discrete” connectomes derived from our mo del ha ve substantially higher test-retest reliabilit y compared to standard metho ds. Keyw ords: Human Connectome, Diffusion MRI, Non-Parametric Esti- mation 1 In tro duction In recen t years the study of structural and functional brain connectivit y has expanded rapidly . F ollo wing the rise of diffusion and functional MRI connec- tomics has unlo c ked a wealth of knowledge to b e explored. Almost synonymous with the connectome is the net work-theory based representation of the brain. In m uch of the recen t literature the quantitativ e analysis of connectomes has fo cused on region-to-region connectivit y . This paradigm equates physical brain regions with no des in a graph, and uses observed structural measurements or functional correlations as a pro xy for edge strengths b etw een no des. Critical to this representation of connectivit y is the delineation of brain re- gions, the parcellation. Multiple studies hav e shown that the c hoice of parcel- lation influences the graph statistics of b oth structural and functional netw orks [15,17,18]. It remains an open question whic h of the prop osed parcellations is the optimal represen tation, or even if such a parcellation exists [14]. It is th us useful to construct a more general framew ork for cortical connectiv- it y , one in which an y particular parcellation of the cortex may be expressed and its connectivity matrix deriv ed, and one in which the v ariability of connectivity Earlier forms of this pap er included erroneous results due to a softw are bug. These errors, whic h artificially raised ICC scores, hav e b een corrected in this version. measures can b e mo deled and assessed statistically . It is also imp ortan t that this framework allow comparisons b et ween parcellations, and representations in this framework must b e b oth analytically and computationally tractable. Since sev eral brain parcellations at the macroscopic scale are possible, a representation of connectivit y that is indep enden t of parcellation is particularly app ealing. In this pap er, we develop such a general framework for a parcellation inde- p enden t connectivity representation, building on the work of [8]. W e describ e a con tinuous p oin t pro cess mo del for the generation of observed tract 1 (streamline) in tersections with the cortical surface, from which w e may recov er a distribution of edge strengths for an y pair of cortical regions, as measured by the in ter-region tract count. Our model is an intensit y function o ver the pro duct space of the cortical surface with itself, assigning to every pair of p oin ts on the surface a p oint connectivity . W e describe an efficient metho d to estimate the parameter of the mo del, as well as a metho d to recov er the region-to-region edge strength. W e then demonstrate the estimation of the mo del on a T est-Retest dataset. W e pro vide repro ducibilit y estimates for our metho d and the standard direct count metho ds [10] for comparison. W e also compare the representational p ow er of common cortical parcellations with resp ect to a v ariet y of measures. 2 Con tin uous Connectivit y Mo del The key theoretical comp onen t of our work is the use of p oint pro cess theory to describ e estimated cortical tract pro jections. A p oin t pro cess is a random pro cess where an y realization consists of a collection of discrete p oin ts on a measurable space. The most basic of these pro cesses is the Poisson pro cess, in whic h even ts occur indep enden tly at a specific asymptotic intensit y (rate) λ ov er the chosen domain [12]. λ completely characterizes each particular pro cess, and is often defined as a function λ : Domain → R + , which allows the pro cess to v ary in intensit y b y lo cation. The exp ected count of any sub-region (subset) of the domain is its total in tensity , the integral of λ ov er the sub-region. In this pap er, our domain is the connectivity space of the cortex, the set of all pairs of p oin ts on the surface, and the even ts are estimated tract intersections with the cortical surface. 2.1 Mo del Definition Let Ω b e union of tw o disjoint subspaces eac h diffeomorphic to the 2-sphere represen ting the white matter b oundaries in each hemisphere. F urther consider the space Ω × Ω , which here represents all p ossible end p oint pairs for tracts that reac h the white matter b oundary . W e treat the observ ation of suc h tracts as ev ents generated by an inhomogeneous (symmetric) Poisson pro cess on Ω × Ω ; in our case, for ev ery even t ( x, y ) we hav e a symmetric even t ( y , x ). 1 It is critical to distinguish b etw een white matter fib ers (fascicles) and observ ed “tracts.” Here, “tracts” denotes the 3d-curv es reco vered from Diffusion W eigh ted Imaging via tractograph y algorithms. Assuming that each even t is indep enden t of all other even ts except for its symmetric even t (i.e., eac h tract is recov ered indep endently), we mo del con- nectivit y as a intensit y function λ : Ω × Ω → R + , such that for an y regions E 1 , E 2 ⊂ Ω , the n umber of even ts is Poisson distributed with parameter C ( E 1 , E 2 ) = Z Z E 1 ,E 2 λ ( x, y ) dxdy. (1) Due to prop erties of the Poisson distribution, the exp ected num b er of tracts is exactly C ( E 1 , E 2 ). F or any collection of regions { E i } N i =1 = P , we can compute a w eighted graph G ( P , λ ) by computing each C ( E i , E j ) for pairs ( E i , E j ) ∈ P × P . Eac h no de in this graph represents a region, and each weigh ted edge represents the pairwise connectivity of the pair of no des (regions) it connects. W e call P a parcellation of Ω if S i E i = Ω and T i E i has measure 0 ( { E i } is almost disjoin t). 2.2 Reco very of the Intensit y F unction A sufficien t statistic for Poisson pro cess mo dels is the intensit y function λ ( x, y ). Estimation of the function is non-trivial, and has been the sub ject of muc h study in the spatial statistics communit y [3]. W e choose to use a non-parametric Kernel Density Estimation (KDE) approac h due to an efficient closed form for bandwidth estimation describ ed b elo w. This pro cess is self-tuning up to a choice of desiderata for the bandwidth parameter. W e first inflate each surface to a sphere and register them using a spherical registration (See section 3.1); how ever this entire method can b e undertaken without group registration. W e treat each hemisphere as disjoint from the other, allo wing us to treat Ω × Ω as the pro duct of spheres ( S 1 ∪ S 2 ) × ( S 1 ∪ S 2 ). Throughout the rest of the pap er D denotes a dataset containing observed tract endp oin ts ( x, y ) i , and ˆ λ denotes our estimation of λ . The unit normalized spherical heat k ernel is a natural choice of k ernel for S 2 . W e use its truncated spherical harmonic representation [1], defined as follo ws for an y tw o unit vectors p and q on the 2-sphere: K σ ( p, q ) = H X h 2 h + 1 4 π exp {− h ( h + 1) σ } P 0 h ( p · q ) Here, P 0 h is the h th degree asso ciated Legendre p olynomial of order 0. Note that the non-zero order p olynomials hav e coefficient zero due to the radial symmetry of the spherical heat k ernel [1]. How ever, since we are estimating a function on Ω × Ω , we use the product of t wo heat k ernels as our KDE k ernel κ . F or any t w o p oin ts p and q , the kernel v alue asso ciated to a end p oin t pair ( x, y ) is κ (( p, q ) | ( x, y )) = K σ ( x, p ) K σ ( y , q ). It is easy to show that R Ω × Ω K σ ( x, p ) K σ ( y , q ) dpdq = 1. The spherical heat kernel has a single shap e parameter σ which corresp onds to its bandwidth. While in general tuning this parameter requires the re-estimation of ˆ λ at every iteration, by rewriting our kernel we can memoize part of the com- putation so that we only need to store the sum of the outer pro ducts of the harmonics. W riting out κ (( p, q ) | D ) = P ( x i ,y i ) ∈ D K σ ( x i , p ) K σ ( y i , q ), we ha ve the follo wing: κ (( p, q ) | D ) = H X h H X k   2 h + 1 4 π   2 k + 1 4 π  exp {− σ ( h 2 + h + k 2 + k ) } | {z } Independent of D , ev aluated every iteration × X ( x i ,y i ) ∈ D P 0 h ( x i · p ) P 0 k ( y i · q ) | {z } Independent of σ, ev aluated once  Th us, ev aluations of the k ernel at any p oin t ( p, q ) can b e done quic kly for se- quences of v alues of σ . W e then are left with the choice of loss function. Denoting the true intensit y function λ , the estimated intensit y ˆ λ , and the leav e-one-out estimate ˆ λ i (lea ving out observ ation i ), In tegrated Squared Error (ISE) is de- fined: I S E ( σ | D ) = Z Ω × Ω ( ˆ λ ( x, y | σ ) − λ ( x, y )) 2 dxdy ≈ Z ˆ λ ( x, y | σ ) 2 dxdy − 2 | D | X ( x i ,y i ) ∈ D ˆ λ i ( x i , y i ) + C onstant. Hall and Marron [9] suggest tuning bandwidth parameters using ISE. In practice, w e find that replacing each lea v e-one-out estimate with its logarithm log ˆ λ i ( x i , y i ) yields more consisten t and stable results. 2.3 Selecting a parcellation Giv en an estimated intensit y ˆ λ and tw o or more parcellations P 1 , P 2 , . . . , w e w ould lik e to know which parcellation and asso ciated graph G ( P , ˆ λ ) b est repre- sen ts the underlying connectivit y function. There are at least t wo p erspectives to consider. Appro ximation Error : Because eac h P i co vers Ω (and P i × P i = Ω × Ω ), eac h G ( P 1 , ˆ λ ) can b e view ed as a piece-wise function g : Ω × Ω → R + , where g ( x, y ) = 1 | E i || E j | C ( E i , E j ) such that x ∈ E i and y ∈ E j . In other words, g is the constant approximation to λ ov er every pair of regions. A natural measure of error is another form of In tegrated Squared Error: E r r ( ˆ λ, G ( P 1 , ˆ λ )) = Z Z Ω × Ω ( g ( x, y ) − λ ( x, y )) 2 dxdy . (2) This is analogous to squared loss ( ` 2 -loss). Lik eliho o d : An alternative viewp oin t leverages the p oin t pro cess mo del to measure a lik eliho o d: log L ( P ) = X E i ,E j ∈ P log Poisson( |{ ( x, y ) i ∈ D : x ∈ E i , y ∈ E j }| ; C ( E i , E j )) . (3) Here, the independence assumption plays a critical role, allowing pairs of regions to b e ev aluated separately . Unfortunately this is biased to ward parcellations with more, smaller regions, as the Poisson distribution has tied v ariance and mean in one parameter. A p opular likelihoo d-based option that somewhat counterbal- ances this is Ak aike’s Information Criterion (AIC), AI C ( P ) = − 2 log L ( P ) + | P | ( | P | − 1) . (4) AIC balances accuracy with parsimon y , p enalizing ov erly parameterized mo dels - in our case, parcellations with to o many regions. 3 Application to CoRR T est-Retest Data W e demonstrate the use of our framework on a test-retest dataset. W e measure connectome reproducibility using In traclass Correlation (ICC ) [13], and compare three parcellations using m ultiple criteria (See Equations 2,3, and 4). 3.1 Pro cedure, Connectome Generation, and Ev aluation Our data are comprised of 29 sub jects from the Institute of Psychology , Chinese Academ y of Sciences sub-dataset of the larger Consortium for Reliability and Repro ducibilit y (CoRR) dataset [19]. T1-w eighted (T1w) and diffusion w eigh ted (D WI) images were obtained on 3T Siemens T rioTim using an 8-channel head coil and 60 directions. Each sub ject w as scanned t wice, roughly t wo w eeks apart. T1w images were pro cessed with F reesufer’s [4] recon-all pip eline to obtain a triangle mesh of the grey-white matter boundary registered to a shared spherical space [5], as w ell as corresp onding v ertex lab els p er sub ject for three atlas-based cortical parcellations, the Destrieux atlas [6], the Desik an-Killiany (DK) atlas [2], and the Desik an-Killiany-T ourville (DKT31) atlas [11]. Probabilistic streamline tractograph y was conducted using the D WI in 2mm isotropic MNI 152 space, using Dipy’s [7] implemen tation of constrained spherical decon volution (CSD) [16] with a harmonic order of 6. As p er Dipy’s ACT, we retained only tracts longer than 5mm with endp oin ts in lik ely grey matter. W e pro vide the mean ICC score computed b oth with and without en tries that are zero for all sub jects. When estimating ˆ λ the kernels are divided b y the n umber of tracks, and we use a sphere with unit surface area instead of unit radius for ease of computation. W e threshold each of the kernel integrated connectomes at 10 − 5 , which is approximately one half of one unit trac k density . W e then compute three measures of parcellation representation accuracy , namely ISE, Negativ e Log Likelihoo d, and AIC scores. 3.2 Results & Discussion T able 1 shows a surprisingly lo w mean ICC scores for regular coun t matrices. This ma y b e b ecause ICC normalizes each measure by its s 2 statistic, meaning that A tlas DK Destrieux DKT31 Num b er of Regions 68 148 62 Coun t ICC 0.2093 0.1722 0.2266 In tensity ICC (F ull) 0.4868 0.4535 0.4388 In tensity ICC (w/Threshold) 0.5613 0.6481 0.4645 T able 1. This table shows mean ICC scores for each connectome generation metho d. The count method - the standard approac h - defines edge strength b y the fib er endpoint coun t. The integrated in tensity metho d is our prop osed metho d; in general it returns a dense matrix. Ho wev er, man y of the v alues are extremely lo w, and so we include results b oth with and without thresholding. Highest ICC scores for each atlas are b olded. T yp e DK Destrieux DKT31 ISE 1 . 8526 × 10 − 5 2 . 1005 × 10 − 5 2 . 1258 × 10 − 5 Negativ e LogLik 85062.5 355769.4 88444.5 AIC Score 174680.95 733294.8 185253.5 Retest ISE 1 . 0517 × 10 − 5 1 . 0257 × 10 − 5 1 . 1262 × 10 − 5 Retest Negativ e LogLik 85256.0 357292.9 88434.9 Retest AIC Score 175068.1 736341.9 185234.3 T able 2. This table sho ws the means o ver all sub jects of three measures of parcellation “go odness”. The retest versions are the mean of the measure using the parcellation’s regional connectivity matrix (or the count matrix) from one scan, and the estimated in tensity function from the other scan. en tries in the adjacency matrices that should b e zero but that are sub ject to a small amoun t of noise – a few erroneous tracks – hav e v ery lo w ICC. Our method in effect smo oths tracts endp oin ts into a density; end p oin ts near the region b oundaries are in effect shared with the adjacen t regions. Thus, ev en without thresholding we damp en noise effects as measured by ICC. With thresholding, our metho d’s p erformance is further impro ved, handily b eating the coun ting metho d with resp ect to ICC score. It is imp ortan t to note that for many graph statistics, changing graph top ology can greatly affect the measured v alue [18]. While it is imp ortan t to ha ve consistent non-zero measurements, the difference b et w een zero and small but non-zero in the graph con text is also non-trivial. The consistency of zero-v alued measuremen ts is th us very imp ortan t in connectomics. T able 2 suggests that all three measures, while clearly differen t, are consisten t in their selection at least with respect to these three parcellations. It is somewhat surprising that the Destrieux atlas has quite low lik eliho o d criteria, but this ma y b e due to the (quadratically) larger n umber of region pairs. Both lik eliho od based retest statistics also choose the DK parcellation, while ISE chooses the Destrieux parcellation by a small margin. It should b e noted that these results m ust be conditioned on the use of a probabilistic CSD tractograph y model. Differen t mo dels may lead to differen t in tensity functions and resulting matrices. The Fig. 1. A visualization of the ICC scores for connectivity to Bro dmann Area 45 (De- strieux region 14) for the Count connectomes ( left ) and the prop osed Integrated In- tensit y connectomes ( right ). Blue denotes a higher score. Fig. 2. A visualization of the marginal connectivity M ( x ) = R E i ˆ λ ( x, y ) dy for the Left P ost-central Gyrus region of the DK atlas (Region 57). The region is shown in blue on the inset. Red denotes higher connectivity regions with the blue region. biases and merits the differen t mo dels and metho ds (e.g. gray matter dilation for fib er counting vs streamline pro jection) remain imp ortan t op en questions. 4 Conclusion W e hav e presented a general framework for structural brain connectivity . This framew ork pro vides a representation for cortical connectivity that is indep enden t of the c hoice of regions, and thus may be used to compare the accuracy of a giv en set of regions’ connectivity matrix. W e provide one p ossible estimation metho d for this represen tation, leveraging spherical harmonics for fast parameter estimation. W e ha ve demonstrated this framew ork’s viabilit y , as w ell as provided a preliminary comparison of regions using sev eral measures of accuracy . The results presented here lead us to conjecture that our connectome esti- mates are more reliable compared to standard fib er counting, though we stress that a muc h larger study is required for strong conclusions to b e made. F ur- ther adaptations of our metho d are p ossible, suc h as using F A-weigh ted fib er coun ting. Our future work will explore these options, conduct tests on larger datasets, and inv estigate the relative differences b et ween tracking metho ds and parcellations more rigorously . Ac kno wledgments This w ork was supported b y NIH Grant U54 EB020403, as w ell as the Rose Hills F ellowship at the Univ ersity of Southern California. The authors would like to thank the reviewers as well as Greg V er Steeg for multiple helpful conv ersations. References 1. Ch ung, M.K.: Heat k ernel smoothing on unit sphere. In: Biomedical Imaging: Nano to Macro, 2006. 3rd IEEE International Symposium on. pp. 992–995. IEEE (2006) 2. Desik an, R.S., et al.: An automated labeling system for sub dividing the human cerebral cortex on MRI scans into gyral based regions of interest. NeuroImage 31(3), 968–980 (2006) 3. Diggle, P .: A kernel metho d for smo othing p oin t pro cess data. Applied Statistics pp. 138–147 (1985) 4. Fisc hl, B.: F reesurfer. NeuroImage 2(62), 774–781 (2012) 5. Fisc hl, B., et al.: High-resolution intersub ject av eraging and a co ordinate system for the cortical surface. Human Brain Mapping 8(4), 272–284 (1999) 6. Fisc hl, B., et al.: Automatically parcellating the human cerebral cortex. Cerebral Cortex 14(1), 11–22 (2004) 7. Garyfallidis, E., et al.: Dip y , a library for the analysis of diffusion M RI data. F ront. Neuroinform 8(8) (2014) 8. Gutman, B., et al.: Registering cortical surfaces based on whole-brain structural connectivit y and contin uous connectivity analysis. In: MICCAI 2014, pp. 161–168. Springer (2014) 9. Hall, P ., Marron, J.S.: Exten t to whic h least-squares cross-v alidation minimises in tegrated square error in nonparametric density estimation. Probability Theory and Related Fields 74(4), 567–581 (1987) 10. Jahanshad, N., et al.: Alzheimer’s Disease Neuroimaging I (2013) genome-wide scan of healthy human connectome discov ers SPON1 gene v ariant influencing dementia sev erity . Pro c Natl Acad Sci USA 110(12), 4768–4773 11. Klein, A., T ourville, J., et al.: 101 lab eled brain images and a consistent human cortical lab eling proto col. F ront. Neurosci 6(171), 10–3389 (2012) 12. Moller, J., W aagep etersen, R.P .: Statistical inference and sim ulation for spatial p oin t pro cesses. CRC Press (2003) 13. P ortney , L.G., W atkins, M.P .: Statistical measures of reliability . F oundations of Clinical Researc h: Applications to Practice 2, 557–586 (2000) 14. de Reus, M.A., V an den Heuvel, M.P .: The parcellation-based connectome: limita- tions and extensions. NeuroImage 80, 397–404 (2013) 15. Satterth waite, T.D., Dav atzik os, C.: T ow ards an individualized delineation of func- tional neuroanatom y . Neuron 87(3), 471–473 (2015) 16. T ournier, J.D., Y eh, C.H., Calamante, F., Cho, K.H., Connelly , A., Lin, C.P .: Re- solving crossing fibres using constrained spherical deconv olution: v alidation using diffusion-w eighted imaging phantom data. NeuroImage 42(2), 617–625 (2008) 17. W ang, J., et al.: Parcellation-dependent small-world brain functional netw orks: A resting-state fMRI study . Human Brain Mapping 30(5), 1511–1523 (2009) 18. Zalesky , A., et al.: Whole-brain anatomical net works: do es the choice of no des matter? NeuroImage 50(3), 970–983 (2010) 19. Zuo, X.N., et al.: An op en science resource for establishing reliability and repro- ducibilit y in functional connectomics. Scientific Data 1 (2014)