The study was approved by the appropriate Institutional Review Boards and Ethics Boards of the University of British Columbia (UBC) and the University of North Carolina (UNC). All structural data were obtained as part of fMRI studies whose results are reported elsewhere (e.g., [34]).

### MR Imaging at the University of British Columbia

All subjects gave written informed consent prior to participating. Nine volunteers with clinically diagnosed PD participated in the study (5 men, 4 women, mean age 68.1 ± 6.8 years, 7 right-handed, 2 left-handed). All subjects had mild to moderate PD (Hoehn and Yahr stage 2–3) [35] with mean symptom duration of 3.6 ± 2.6 years. We recruited ten healthy, age-matched control subjects without active neurological disorders (3 men, 7 women, mean age 55 ± 12.4 years, 9 right-handed, 1 left-handed). Exclusion criteria included atypical Parkinsonism, presence of other neurological or psychiatric conditions and use of antidepressants, sleeping tablets, or dopamine blocking agents.

MRI was conducted on a Philips Achieva 3.0 T scanner (Philips, Best, The Netherlands) equipped with a 6 channel Sense head-coil. A high resolution, three dimensional (3D) SPGR image data set of the whole brain consisting of 170 axial slices at a FOV of 256 × 200 mm^{2} was acquired for WM/GM segmentation purposes and as anatomical reference (inversion prepared 3D T1TFE, TR = 7.746 ms, TE = 3.55 ms, inversion delay = 880 ms, flip angle = 8.00°, voxel dimensions 1.0 × 1.0 × 1.0 mm^{3}).

The thalami were one of eighteen specific regions of interest (ROIs) that were manually drawn on each unwarped, aligned structural scan using the Amira software (Mercury Computer Systems, San Diego, USA). Although the thalami were manually segmented on the axial slices, they were carefully examined in the coronal and saggital planes to ensure accuracy. The trained technician performing the segmentation was blinded to the disease state.

### MR Imaging at the University of North Carolina

All subjects gave written informed consent prior to participating. Nine volunteers with clinically diagnosed mild to moderate PD (Hoehn and Yahr stage 2–3 – mean symptom duration of 2.1 ± 2.0 years) participated in the study (5 men, 4 women, mean age 58 ± 12 yrs, all right-handed). We also recruited eight healthy, age-matched control subjects without active neurological disorders (5 men, 3 women, mean age 49 ± 14 yrs, 8 right-handed). Images were acquired on a 3.0 Tesla Siemens scanner (Siemens, Erlangen, Germany) with a birdcage-type standard quadrature head coil and an advanced nuclear magnetic resonance echoplanar system. The head was positioned along the canthomeatal line. Foam padding was used to limit head motion. High-resolution T1 weighted anatomical images were acquired (3D SPGR, TR = 14 ms, TE = 7.7 ms, flip angle = 25°, voxel dimensions 1.0 × 1.0 × 1.0 mm, 176 × 256 voxels, 160 slices).

ROIs (including thalami) were drawn manually by the same trained research associate with assistance from multiple on-line and published atlases (e.g. [36]).

### Thalami Shape Analysis

As described in the technical appendix, the analysis of each shape results in a unique feature vector, of length *n* = 1440. The left and right thalami were analyzed separately.

For comparison, we examined for any differences in volume. The volume of each thalamus was estimated as the number of voxels that each ROI contained multiplied by the volume of a single voxel.

To assess the significance of group differences between feature vectors, we used a permutation test to generate a null distribution of Euclidean distances between feature vectors. The permutation test does not require a priori assumptions about the data distribution, and is thus preferred over T-test and F-test [37]. We assessed the differences in left vs. right thalami in controls, left vs. right in PD subjects, PD vs. controls for the left thalamus, and PD vs. controls for the right thalamus. Although the boundaries of the thalami were determined by visual inspection, in prior work we compared feature vectors derived from thalami segmented from structural scans obtained before and after giving L-dopa medication (as part of another fMRI study) [38]. As expected, no significant differences could be detected in the two groups, suggesting that independent manual segmentation did not incur significant systematic errors.

### Shape Analysis – technical aspects

Let Ψ(*θ*, *φ*) be a function defined on the unit spherewith *θ* and *φ* as the zenithal and azimuthal angles, respectively. The SPHARM representation for this function is given by (1) where {Y}_{lm}^{\ast}\left(\theta ,\phi \right) is the complex conjugate of the *m*
^{th} order spherical harmonic of degree *l*. *l* ranges from 0 to *L* [16]. Increasing the value of *L*, also called the bandwidth, improves the representation accuracy at the cost of higher computation time. This definition can also be extended to real valued 3D distributions Ψ(*r*, *θ*, *φ*) (2), where *r* is the distance from the origin to a given voxel. *k* is an index introduced to account for possible degeneracy due to the additional dimension [39].

{c}_{l}^{m}={\displaystyle \underset{0}{\overset{2\pi}{\int}}d\phi}{\displaystyle \underset{0}{\overset{\pi}{\int}}{Y}_{lm}^{\ast}(\theta ,\phi )\Psi (\theta ,\phi )\mathrm{sin}\phantom{\rule{0.5em}{0ex}}(\theta )d\theta}

(1)

{c}_{kl}^{m}={\displaystyle \underset{0}{\overset{\infty}{\int}}{r}^{2}dr}{\displaystyle \underset{0}{\overset{2\pi}{\int}}d\phi {\displaystyle \underset{0}{\overset{\pi}{\int}}\sqrt{2}\frac{\mathrm{sin}\phantom{\rule{0.5em}{0ex}}(\pi kr)}{r}{Y}_{lm}^{\ast}(\theta ,\phi )\Psi (r,\theta ,\phi )\mathrm{sin}\phantom{\rule{0.5em}{0ex}}(\theta )d\theta}}

(2)

As explained later in this section, rotationally invariant features can be derived from this spherical harmonic representation. In our application, we also need the features to be invariant to any translation of the entire ROI in 3D space. To achieve this, we move the origin of the function Ψ(*r*, *θ*, *φ*), to the centroid of Ψ_{
s
}(*r*, *θ*, *φ*), where Ψ_{
s
}(*r*, *θ*, *φ*) is given by (3).

{\Psi}_{s}\left(r,\theta ,\phi \right)=\{\begin{array}{lll}1\hfill & if\hfill & \Psi \left(r,\theta ,\phi \right)\ne 0\hfill \\ 0\hfill & if\hfill & \Psi \left(r,\theta ,\phi \right)=0\hfill \end{array}

(3)

Since direct computation of (2) is highly inefficient [40], we use an alternate approach by representing the data as a set of spherical functions obtained by intersecting the 3D data with spherical shells. Alternatively, for each value of *r*, Ψ(*r*, *θ*, *φ*) can be visualized as a spherical shell comprising the function values at a distance *r* from the origin. *r* can then be incremented in steps of *t* to encompass the entire ROI. If the initial representation of the function is in the form of a cubic grid (regular isotropic voxels in our case), volumetric interpolation is required to resample the ROI in the spherical coordinate space.

When analyzing multiple subjects' ROIs simultaneously, we define the maximum radius, *R*
_{
max
}, as the minimum radial distance in voxel count that encompasses all non-zero values of all subjects' ROIs being analyzed. To represent the values from the cubic grid of all ROIs with sufficient accuracy, *2R*
_{
max
}shells are used. To achieve scale invariance, the shells must be distributed evenly throughout the spatial extent of each ROI. Since the ROI size across subjects is not uniform, shell spacing *t* must be adjusted for each subject separately. This procedure ensures that each shell captures similar features from the 3D ROI irrespective of its scale.

Surface sampling along each of these shells is performed on an equiangular spherical grid of dimensions *2L* × *2L* [40]. The common bandwidth *L* for all shells of all functions is chosen to satisfy the sampling criterion for the largest shell in this set of ROI, namely the one with radius *R*
_{
max
}. The surface area for this shell represents the maximum surface shell area that needs to be sampled by the equiangular grid; hence, any value of *L* satisfying the required equiangular sampling (*2L* × *2L*) at this shell will be sufficient to represent data from smaller radii shells. The minimum value for *L* is obtained by equating the surface area of this largest shell to the equiangular sampling grid (4). Higher values of L are not used, since it increase computation time with no added benefit. Also, this will result in longer feature vectors, complicating the analysis. Furthermore, when the represented object is a discrete array, higher values of *l*, resulting from a larger *L*, may correspond to sampling noise [39]. Recognizing that in applications pertaining to discrimination, high accuracy in the SPHARM representation is not a necessity, we chose to use the minimum value for *L* as that obtained by (4).

4\pi {R}_{\mathrm{max}\phantom{\rule{0.5em}{0ex}}}^{2}=2L\times 2L,L={R}_{\mathrm{max}\phantom{\rule{0.5em}{0ex}}}\sqrt{\pi}

(4)

To obtain the SPHARM representation for all shells, a discrete SPHARM transform is performed at each value of *r* to obtain {c}_{rl}^{m} (5). Features derived from this representation, however, do not provide a unique function representation [41, 42]. For instance, rotating the inner and outer shells by different amount will result in different spatial distribution of function values. However, in this approach, the derived features are insensitive to these rotational transforms, thus resulting in the same feature values for dissimilar spatial distributions.

\begin{array}{c}{c}_{rl}^{m}={\displaystyle \underset{0}{\overset{2\pi}{\int}}d\phi}{\displaystyle \underset{0}{\overset{\pi}{\int}}{Y}_{lm}^{\ast}(\theta ,\phi )\Psi (r,\theta ,\phi )\mathrm{sin}\phantom{\rule{0.5em}{0ex}}(\theta )d\theta}\\ r=[1,2,3,\mathrm{...},2{R}_{\mathrm{max}\phantom{\rule{0.5em}{0ex}}}]\end{array}

(5)

Burel and Henocq's original equation (2) does not have this problem, since a part of the basis function is a function of *r*. However, since (2) is computationally intractable [17], we proposed an efficient approach that uses a radial transform (6), derived from (2), to obtain a unique function representation. The transform (6) retains the relative orientation information of the shells, thus the features derived will be sensitive to independent rotations of the different shells, thereby ensuring that *unique* feature representation is obtained.

\begin{array}{c}{c}_{kl}^{m}={\displaystyle \sum _{r=1}^{2{R}_{\mathrm{max}\phantom{\rule{0.5em}{0ex}}}}{r}^{2}\sqrt{2}\frac{\mathrm{sin}\phantom{\rule{0.5em}{0ex}}(\pi kr)}{r}{c}_{rl}^{m}}\\ k=[1,2,3,\mathrm{...},2{R}_{\mathrm{max}\phantom{\rule{0.5em}{0ex}}}]\end{array}

(6)

The range of *k* could be changed to obtain different lengths of the final feature vector. However, to avoid unnecessarily increasing the feature vector length or losing any important information caused by reducing the range, we choose to keep the range of *k* the same as that of *r*, i.e. 2*R*
_{
max
}.

From the obtained representation (6), we then compute similarity transform invariant features using (7) for each value of *l* and *k* [39] with *p* and *q* are used to index these features. Note we reshape *I* into a single row vector of dimensions *D* = *L* × 2*R*
_{
max
}for later analysis.

I\left(p,q\right)={\displaystyle \sum _{k=l}^{k=2{R}_{\mathrm{max}\phantom{\rule{0.5em}{0ex}}}}{\displaystyle \sum _{m=-l}^{m=l}{c}_{kl}^{m}\left({c}_{kl}^{m}\right)*}},p=\mathrm{1...}L,q=l\mathrm{...2}{R}_{\mathrm{max}\phantom{\rule{0.5em}{0ex}}}

(7)

In order to provide a scalar estimate how different a given thalamus shape was, we calculated the mean of the row vectors, I (Eqn 7) separately for both the left and right thalami. A distance metric, estimating how "abnormal" a given shape was estimated by determining the Euclidean distance between the given feature vector and the mean vector. For example, the distance for right thalamus for the *j*
^{th}subject was estimated as:

Dis{t}_{right}^{j}=\sqrt{{\displaystyle \sum _{k=1}^{D}({I}_{k}^{{j}^{j}}}-{I}_{right}{)}^{2}}

(8)