Nagesh Adluru^{1}, Hyunwoo J Kim^{2}, Richard J Davidson^{3}, Andrew L Alexander^{4}, Sterling C Johnson^{5}, and Vikas Singh^{6}

This work presents novel statistical image analysis methods to characterize complex morphological brain changes using MRI data. Specifically, our procedure utilizes the fundamental representations of "longitudinal change" -- voxel-wise Jacobian matrices obtained from image registration. Currently their univariate summaries (for example determinants) are ubiquitously used in neuroimaging studies. Operating directly with representations of Jacobians namely Cauchy deformation tensors, which are elements of an abstract mathematical manifold of symmetric positive definite matrices, yields promising improvements in statistical power in detecting subtle but statistically significant effects. The key technical contributions are computational algorithms for estimating multivariate general linear models with manifold-valued response variables.

Let us consider a standard voxel-wise analysis setup. A general linear model (GLM) identifies associations between covariates and the measurements at each voxel. For each voxel $$$i$$$, we solve $$$y_i=\alpha+\mathbf{\beta}^T\mathbf{x}_i+\varepsilon_i$$$, for $$$\alpha$$$, the $$$y$$$-intercept and $$$\mathbf{\beta}$$$, a $$$p$$$-dimensional coefficient vector for the covariates $$$\mathbf{x}_i=\{x_i^{1},x_i^{2},\ldots,x_i^{p}\}$$$. $$$\varepsilon_i$$$ is the error in the voxel measurement $$$y_i$$$. Using data from $$$N\geq (p+1)$$$ subjects, the coefficients $$$\alpha$$$ and $$$\beta$$$ can be estimated by minimizing the residual $$$\displaystyle\sum_{j=1}^N\left(y_i^j-\alpha-\mathbf{\beta}^T\mathbf{x}_i^j\right)^2$$$. When the response is a CDT (SPD matrix) we cannot directly perform an element-by-element subtraction of the matrices for computing residuals^{1,6,7}. To see this, consider a manifold visualized in Fig. 2 in the context of the subtraction operation. The residual (or distance) between the measurement $$$y$$$ and the model prediction $$$\hat{y}$$$ should in principle be the shortest geodesic curve (black) connecting $$$y$$$ and $$$\hat{y}$$$ that lies on the manifold. A direct element-by-element subtraction estimates the length of the straight line (blue) that lies in the ambient space and not entirely on the manifold. This can be a poor approximation when the curvature of the manifold is high. Such approximations may result in reduction of power for detecting statistical associations between $$$\mathbf{x}$$$ and $$$y$$$. Hence the distance $$$(y-\hat{y})$$$ should be computed by an operation that measures the length of the shortest curve connecting these two matrices *along *the manifold.For simplicity, we present the main ideas in our algorithms for computing residuals using a single covariate $$$(p=1)$$$ MGLM. There are two main pieces in minimizing the residual. (1) $$$\alpha$$$ which corresponds to a distinct unknown point on the manifold serves as an “offset”. This offset will be used to define a tangent space – a plane locally tangent to the manifold (Fig. 3). (2) The coefficient $$$\beta$$$, in our formulation, corresponds to a vector $$$v$$$ in the tangent space. It turns out that the tangent vector $$$v$$$ precisely corresponds to a “geodesic curve” on the manifold (shown in black in Fig. 3). Thus we can parameterize geodesic curves via tangent vectors. Then using concepts from differential geometry^{8-12}, we can compute the primary quantity of interest (residual), directly on the manifold. The GLM estimation procedure then reduces to searching through possible values of the offset and the tangent vector such that residual is as small as possible. The technical details when $$$p>1$$$ (Fig. 3, right) are more complicated but conceptually similar to the $$$p=1$$$ case. Since statistics using manifold-valued data do not satisfy many of the standard parametric-distributional assumptions, $$$p$$$-values can be obtained using permutation testing.

1. Fletcher P. and Joshi S. Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors. Computer Vision and Mathematical Methods in Medical and Biomedical Image Analysis, ed: Springer, 2004, pp. 87-98.

2. Lenglet C., Rousson M., Deriche R. and Faugeras O. Statistics on the manifold of multivariate normal distributions: Theory and application to diffusion tensor MRI processing. Journal of Mathematical Imaging and Vision, 2006, vol. 25, pp. 423-444.

3. Friedrichs K. Symmetric positive linear differential equations. Communications on Pure and Applied Mathematics, 1958, vol. 11, pp. 333-418.

4. Bhatia R. and Holbrook J. Riemannian geometry and matrix geometric means. Linear Algebra and its applications, 2006. vol. 413, pp. 594-618.

5. Arsigny V., Fillard P., Pennec X. and N. Ayache. Geometric means in a novel vector space structure on symmetric positive-definite matrices. SIAM Journal on Matrix Analysis and Applications, 2007, vol. 29, pp. 328-347.

6. Jayasumana S., Hartley R., Salzmann M., Li H. and Harandi M. Kernel methods on the Riemannian manifold of symmetric positive definite matrices. Computer Vision and Pattern Recognition (CVPR), 2013, pp. 73-80.

7. Moakher M. A differential geometric approach to the geometric mean of symmetric positive-definite matrices. SIAM Journal on Matrix Analysis and Applications, 2005, vol. 26, pp. 735-747.

8. Helgason S. Differential geometry, Lie groups, and symmetric spaces. Academic press, 1979, vol. 80.

9. Eguchi T., Gilkey P. and Hanson A. Gravitation, gauge theories and differential geometry. Physics reports, 1980, vol. 66, pp. 213-393.

10. Do Carmo M. Differential geometry of curves and surfaces. Prentice-hall Englewood Cliffs, 1976, vol. 2.

11. Kobayashi S. and Nomizu K. Foundations of differential geometry. Vol 1: New York, 1963.

12. Spivak M. A comprehensive introduction to differential geometry. vol. 2, IV[A] Boston, Mass., 1970.

An example panel of data generated in morphometric studies using MRI data. **(a)** The moving brain image. **(b)** Warped grid to move image in **(a)** to that in **(d)**. **(c)** Vector field representing local deformations. **(e)** $$$|J$$$ map derived from **(c)**. **(f)** Zoomed in view of the $$$|J|$$$ map. **(g)** CDT map ($$$\sqrt{J^TJ}$$$) map derived from $$$J$$$. **(h)** Zoomed in view of the CDT map. **Significance:** Among the different features (vector field, $$$J$$$, $$$|J|$$$) of brain morphology that can be analyzed, CDTs are the focus of this work. Full resolution image can be viewed online.

Illustration of geodesic curve (black) on a manifold $$$\mathcal{M}$$$. **Significance:** The blue straight line connecting $$$y$$$ and $$$\hat{y}$$$ does not lie on the manifold and can be a poor approximation for representing residual errors when estimating regression models for statistical analysis. The approximation gets worse as the curvature of the manifold gets higher. Full resolution image can be viewed online.

Demonstration of linear regression on manifolds with $$$N=4$$$ samples. An example of univariate ($$$p=1$$$) and bivariate regression ($$$p=2$$$) is presented on left and right respectively. **Significance:** The black curve denotes a curve on the manifold, parameterized by the blue tangent vector $$$v$$$. The various $$$y$$$s are the measurements and the sum of the brown arrows is the residual error. Full resolution image can be viewed online.

Spatio-temporally unbiased estimation of the global coordinate system for the longitudinally acquired imaging data. In our work we analyzed $$$T_1$$$ weighted MRI data from $$$N=243$$$ individuals from two different visits. **Significance:** Visits are averaged first which are then used to estimate the global average. Full resolution image can be viewed online.

Preliminary analyses comparing results from manifold general linear modeling (MGLM) and standard GLM for detecting associations (at $$$p\leq 0.05$$$) between amyloid burden and the longitudinal morphometric changes in the brain as captured by the CDTs. The amyloid burden was measured using the average distributional volume ratio (DVR) of $$${}^{11}$$$C PiB-PET data from 8 different bilateral gray matter regions. The $$$p$$$-values were obtained using 20,000 permutations. **Significance:** CDTs tend to result in lower $$$p$$$-values than the (log) Jacobian determinants. Full resolution image can be viewed online.