UBC Faculty Research and Publications

Connectivity-Based Parcellation of Functional SubROIs in Putamen using A Sparse Spatially Regularized… Zhang, Yiming; Liu, Aiping; Tan, Sun Nee; McKeown, Martin J.; Wang, Z. Jane (Zheng Jane) May 31, 2016

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


52383-Zhang_Y_et_al_Connectivity_based_parcellation_functional.pdf [ 1.52MB ]
JSON: 52383-1.0319906.json
JSON-LD: 52383-1.0319906-ld.json
RDF/XML (Pretty): 52383-1.0319906-rdf.xml
RDF/JSON: 52383-1.0319906-rdf.json
Turtle: 52383-1.0319906-turtle.txt
N-Triples: 52383-1.0319906-rdf-ntriples.txt
Original Record: 52383-1.0319906-source.json
Full Text

Full Text

Connectivity-Based Parcellation of Functional SubROIsin Putamen using A Sparse Spatially RegularizedRegression ModelYiming Zhanga, Aiping Liuc,∗, Sun Nee Tanb, Martin J. McKeownb, Z. JaneWangaaDepartment of Electrical and Computer Engineering,University of British Columbia, V6T 1Z4, Vancouver, CanadabDepartment of Neurology and Pacific Parkinson’s Research Centre,University of British Columbia, V6T 1Z4, Vancouver, CanadacDepartment of Biomedical Engineering,Hefei University of Technology, 193 Tunxi Road, 230009 Hefei,ChinaAbstractIn this paper, we present a novel framework for parcellation of a brain region intofunctional subROIs (Sub-Region-of-Interest) based on their connectivity pat-terns with other brain regions. By utilising previously established neuroanatomyinformation, the proposed method aims at finding spatially-continuous, func-tionally consistent subROIs in a given brain region. The proposed frameworkrelies on 1) a sparse spatially-regularized fused lasso regression model for en-couraging spatially and functionally adjacent voxels to share similar regressioncoefficients; 2) an iterative merging and adaptive parameter tuning process; 3) aGraph-Cut optimization algorithm for assigning overlapped voxels into separatesubROIs. Our simulation results demonstrate that the proposed method couldreliably yield spatially continuous and functionally consistent subROIs. Weapplied the method to resting-state fMRI data obtained from normal subjectsand explored connectivity to the putamen. Two distinct functional subROIscould be parcellated out in the putamen region in all subjects. This approachprovides a way to extract functional subROIs that can then be investigatedfor alterations in connectivity in diseases of the basal ganglia, for example in∗Corresponding authorEmail address: aipingl@ece.ubc.ca (Aiping Liu )Preprint submitted to Biomedical Signal Processing and Control October 18, 2016Parkinson’s disease.Keywords: fMRI, brain connectivity, brain parcellation, spatial regularization1. IntroductionFunctional magnetic resonance imaging (fMRI) is a functional neuroimagingtechnique that indirectly measures brain activity by detecting associated alter-ations in blood oxygenation (BOLD signal). In the past, most fMRI studiesfocused on detection of localized neural activities by modeling the relationshipbetween fMRI signals and experiment stimulus, i.e., activity studies [1, 2, 3].However, the human brain relies on efficient networks of interacting brain re-gions [4]. Hence, interests in studying the associations between brain regionshave grown, i.e., connectivity studies [5, 6, 7]. Connectivity studies can be ex-plored using task-related as well as resting state fMRI data, with the latterlooking at spontaneous interactions between different brain regions without re-quiring active engagement from the subject. Resting state fMRI may thereforebe more suitable for studies involving aging and diseased populations [8] whooftentimes suffer from sensory, motor and/or cognitive impairments renderingthem incapable of performing challenging tasks.Connectivity studies can be conducted at either the voxel or ROI (regions-of-interest) level. Voxel-based approaches usually involve a large number ofvariables, are computationally-inefficient, and must account for the massiveamount of multiple comparisons. Such voxel-based approaches are typicallydone by spatially transforming all brain volumes to the same anatomical tem-plate. However, subtle misregistration can make the assumption that, afterregistration, a given voxel will represent the same functional region across allsubjects tenuous. ROI-based connectivity analysis may reduce the number ofmultiple comparisons, and does not necessarily require spatial transformation,but requires careful consideration as to the definition of an ROI. AnatomicalROIs may be used to infer functional ROIs [9], but a single anatomical ROI, suchas the putamen or amygdala, may in fact encompass distinct functional subROIs2Figure 1: Illustration of functional subROIs in the putamen brain region. The figure is madeaccording to the graph shown in [18].[10]. A number of attempts have been made to utilize data-driven approachesto subdivide a given ROIs into subROIs based on functional connectivity. Onebroad approach is based on cluster analysis [11] which first retrieves connectivityfeatures of each voxel within an ROI using general linear regression models orPearson’s pairwise correlation coefficients, and then applies clustering accordingto the functional distances defined by the extracted features. Various clusteringmethods have been adopted in previous studies, such as the fuzzy clusteringmethod [11], k-means clustering method [12], a self-organized mapping method[10] and maximum margin clustering method [13]. However, in order to acquirespatially continuous results, most clustering methods require meticulous de-noising preprocessing methods as they are very sensitive to outliers in the data.Another popular category of data-driven approaches is based on graph theory,where each voxel represents one node in the graph, and methods such as thenormalized cut approach [14] and modularity detection method [15] are used toseparate the graph (and hence ROI) into distinct subROIs. Similar to clustering3methods, graph theory methods do not take into account spatial information,and thus it is often difficult to obtain spatially continuous subROIs using graphtheory methods. Furthermore, most graph theory methods are only concernedwith the connectivity map between voxels within an ROI without incorporatingconnectivity from other ROIs, which might limit their usefulness.In this paper, we propose a novel framework which defines subROIs in oneROI based on their functional connectivity to other ROIs. The proposed methodemploys a fused lasso regression model [16] with a spatial regularization penaltyincorporated. The fused lasso approach encourages sparsity of the regressioncoefficients as well as sparsity of their successive differences between coefficients.We further introduce the normal lasso penalty for all voxels and the fused lassopenalty on spatially adjacent pairs of voxels.Our framework differentiates from other approaches in the literature [11] [12][13][14][15]in two main aspects. Firstly, in our proposed framework, we utilizethe functional connectivity between voxels in the task ROI and the average timeseries of other related ROIs (where the task ROI and the related ROIs belongto one neural control loop). In the literature [11][14][15], people only considerthe functional connectivity between voxels within the task ROI. Secondly, weincorporate spatial information that is often ignored in the literature [11][14][15]into our problem formulation. In addition, in order to limit the amount ofbias that spatial constraint could introduce, we proposed a novel algorithmfor adaptive, data-dependent parameter selection which allows us to only addspatial constraint when it is ‘necessary’.Functionally, in the basal ganglia-cortical loops in animals, the putamen isconnected to several cortices with a clear topography [18]. We have chosenthree reference brain regions namely the orbitofrontal (OF) cortex, cingulategyrus (CG) and sensorimotor cortex(SMA) to assess connectivity to the ROI ofinterest, the putamen. Specifically, the DLS is more strongly connected to theSMA while the DMS has more reciprocal connections with the OF and CG [19][18]. However, the spatial boundary between the DLS and the DMS is blurry andtheir exact locations are unknown due to overlapped connections; Some voxels4of the DMS also have weak connections with the SMA while the DLS, too, mayreceive weak connections from the CG and OF. This scenario is illustrated inFigure 1. In addition, as we observed in real fMRI data, due to the spatial signalnoise, head movement and other possible artifacts, the data could be spatiallycorrupted with some outlier voxels. As a result, many current parcellationmethods are not able to deal with such corrupted voxels to obtain a spatiallycontinuous parcellation. Therefore, we plan to design an algorithm to parcellatethe putamen region not only according to functional connectivity features fromprior knowledge but also through integration of spatial information.This pilot study aimed to investigate a novel technique to parcellate theputamen into two subregions with distinct functional and structural connectionsto the cortices of the brain. The putamen and caudate are two structurallydistinct brain regions in the brainstem with the former lying more laterallyand inferior to the latter. Due to the proximity of these two brain regions toeach other and their shared neuronal connections, together, they are knownas the striatum. Within the striatum (caudate and putamen), it has spatiallysegregated functional topography. The dorsolateral striatum (DLS) consists ofthe dorsal and lateral aspects of the caudate and putamen, and is associatedwith control of habitual, automatic movements. On the other hand, more medialareas within the caudate and putamen are known as the dorsomedial striatum(DMS) and this region is functionally related to learning and execution of goal-oriented movements. As an exploratory first step, in this paper, we chose toparcellate the putamen into the DLS and DMS subregions.In the remainder of the paper, we will present the proposed method in Sec-tion 2. In Section 3.1, we investigated on a synthetic dataset to compare theresults of the proposed method with those of clustering and graph theory meth-ods. We tested the proposed method on real resting state fMRI data set in100Section 3.2. In Section 4, we presented a summary of our results with a conclu-sion.52. MethodIn this section, we will describe the proposed framework to separate a givenROI into functionally consistent and spatially continuous subROIs. We definethe ROI to be separated as the task ROI and other ROIs used to estimate theconnections with the task ROI as reference ROIs. We intended to find subsets ofadjacent voxels (i.e., subROIs) in the task ROI that share similar connectivitypatterns to other reference ROIs and iteratively merge them into groups. Inthe proposed algorithm, according to prior neuroanatomical knowledge, thereare several reference ROIs that share similar connectivity patterns with onefunctional subROI in the task ROI (putamen in this paper) so that we canobtain the functional boundary between these subROIs.The proposed framework can be summarized in Table 1. We will now elab-orate on the individual components of the proposed framework in the followingsubsections. First, we start by describing the spatially regularized fused lassomodel.2.1. Spatially Regularized Fused Lasso MethodLet X = [x1, x2, · · · , xn] be a (T ×n)-dimensional data matrix with n denot-ing the number of voxels and T denoting the length of time points. X representsthe fMRI signals in the task ROI and xi, i = 1, 2, · · · , n represents fMRI timecourse of voxel i. Y is a (T × 1) vector representing the signals of one referenceROI which is acquired by averaging time courses across all voxels contained inthat reference ROI. Let β be a (n × 1) vector where each element in β repre-sents the connectivity coefficient between one voxel and the reference ROI. Astandard linear model with one reference ROI Y is:Y = Xβ +  (1)with error  assumed to be zero mean and constant variance.The linear regression model tries to obtain the weight for each voxel in thetask ROI in order to get the best linear fit for the reference signal, i.e., smallest6Table 1: Connectivity-based Brain Parcellation Framework.Input : X: T × n matrix represents task ROI signal.Y1, Y2, . . . , Ym: Each Yi is a T × 1 vector representing timecourse for reference ROI iStep 1 : According to each reference ROI Yi, iterate until finalseparation being obtained.1) Perform spatially regularized fused lasso algorithm betweenvoxels in the task ROI and one reference ROI.2) Merge adjacent voxels with similar connectivity weightacquired from 1) into one group.3) Treat groups from 2) as new voxels, go back to 1)Step 2 : With more than one reference ROI, combine allparcellation results from all reference ROIs by performingGraph-Cut algorithm.Output : Two spatially continuous functional subROIs from Xmean square error. So the problem could be formed as,βˆ = arg minβ‖Xβ − Y ‖2 (2)For fMRI data, it is common that the number of voxels n in task ROI is largerthan the length of time points T . In this case, the problem of multicollinearitywill arise, and an optimal solution to above problem will not be unique. In orderto eliminate the effect of multicollinearity and make the solution unique, weintroduce the normal Lasso (Least Absolute Shrinkage and Selection Operator)penalty∑nj=1 |βj | ≤ s1 where s1 is a tuning parameter [20]. This penaltywill encourage sparsity in the coefficient vector β, forcing small elements tobe exactly zero. In addition, it is well known in fMRI analysis that similarconnectivity patterns tend to exist in clusters of spatially adjacent voxels ratherthan small isolated groups of voxels [21]. However, with only the Lasso penalty,we cannot incorporate spatial information. Therefore, we further introduce a7spatially regularized fused Lasso penalty∑nj=1∑j<i Cji|βj − βi| < s2, wheres2 is a tuning parameter andCji = 1, if voxel j and i are adjacent0, otherwise (3)This is a spatially regularized fused lasso penalty that comes from the normalfused lasso penalty [16], and it is also a special case of a graph-guided fused lassomodel [22]. The spatially regularized fused lasso penalty will encourage sparsityin the difference between spatially adjacent voxels, forcing adjacent voxels tobelong to the same cluster.With the Lagrange multipliers’ method, we could form the optimizationproblem of the spatially regularized fused lasso model as,arg minβ{‖Xβ − Y ‖2 + λ n∑j=1|βj |+ γn∑j=1∑j<iCji|βj − βi|}(4)where the second term encourages sparsity among elements in β, and the thirdterm encourages sparsity in differences between spatially successive elements.In fMRI, usually n > T , i.e., the number of voxels within a ROI, n, is largerthan the number of fMRI time points T . However, regarding to different sizesof ROIs and lengths of fMRI experiments, it is possible that in some ROIs wecan have n < T . For instance, in the fMRI data we have been working on, thelength of fMRI time course T (T = 240) is larger than the number of voxels n(n is around 200) in ROI Caudate. It is commonly recognized that if n < T andif time courses of different voxels are independent to each other, the introduc-tion of Lasso penalty to the regression model is redundant. However, accordingto realistic assumptions of fMRI data, the fMRI time series of different voxelswithin one ROI are not independent and indeed neighbour voxels are generallyexpected to have high correlations. Thus, many predictors suffer from collinear-ity which would allow the coefficients of predictors to be sparse and coefficientdifferences between adjacent voxels to be sparse. Therefore, the introductionof Lasso penalty and fused-lasso penalty is not unnecessary and the proposedalgorithm is still feasible when n < T .8Equation (4) leads to a quadratic programming problem with non-smoothconstraints. In order to solve this problem efficiently, we implemented a Smooth-ing Proximal Gradient (SPG) method recently proposed [23].2.2. Merge Voxels into GroupsWith the lasso penalty and spatially regularized fused lasso penalty, theinferred connectivity weights β and their differences should be sparse. In thesecond step, we merged adjacent voxels that share similar connectivity weightsinto groups. To better illustrate our method, we utilize some definitions fromgraph theory in our description. We mapped the task ROI as a connected graph,where each voxel is taken as one node in the graph, and the edge between twonodes exists only when these two nodes are adjacent. In this network, adjacencyis defined by any pairs of nodes that are horizontally, vertically, or diagonallyadjacent in a 3-dimensional space, i.e., the Euclidean distance between the twonodes is not greater than√3. We then merged those adjacent nodes correspond-ing to the same connectivity weights into the same group.After acquiring new groups, we averaged the signals contained in each groupand treated the averaged signals as new predictor variables. Then, we recal-culated the connectivity weights of the updated predictor variables using thespatially regularized fused lasso model. After merging of the voxels (groups),edges in the updated graph exist only when the minimal distance between vox-els in the two groups is not greater than√3. This process is repeated until thenumber of groups reaches the expected number of subROIs, which is 2 in ourcase. However, with real fMRI data, noise cannot be ignored and tuning param-eters would also affect the results such that the algorithm may not converge.Therefore, if the number of groups remains unchanged after one iteration, wewill increase the penalty parameters to enforce the regularization and relax themerging precision. Merging precision is controlled by a rounding parameter δwhich is defined as the precision of judging the equality of two connectivityweights (For instance, merging precision δ = 104 means that the connectivityweights will be compared by 4 digits after the decimal point).9The detailed algorithm of parcellation with one reference ROI is describedin Algorithm 1.Algorithm 1 Spatially Regularized Fused Regression Based Algorithm forBrain Region Parcellation Algorithm.Require:task ROI voxel level signal X,the mean time course Y of the reference ROI,the lasso penalty parameter λ, the fused lasso penalty parameter γ and themerging precision parameter δ.Initialization:(1)Run spatially regularized fused regression model on all voxels in taskROI, get connectivity weights β for all voxels.(2)Round β by δ, merge adjacent voxels with same connectivity weights intoone group, get m groups. Before iteration, save original merging precisionparameter, δinitial ← δ.Ensure:1: while m 6= 2 do2: Take groups as nodes, use mean value to represent groups signal. Runspatial regularized fused lasso regression on all nodes;3: Round β by δ, merge adjacent nodes who have same connectivity weightsinto one group, get m groups.4: if m = 1 then5: γ ← γ/2, λ← λ/2, δ ← δ × 10, redo 2 and 3;6: else7: if m remain unchanged then8: γ ← γ + 1, λ← λ+ 1,δ ← δ/10, redo 2 and 3;9: end if10: else11: δ ← δinitial, continue12: end if13: end while102.3. Assign Voxels Using the Graph-Cut OptimizationIn our study, three reference ROIs were used for parcellation of the taskROI, i.e., the putamen. After obtaining parcellation based on each reference ROIrespectively, there were some overlapped voxels with different group assignmentsaround the boundaries of all parcellation results. In order to further define theboundary between the DLS and DMS, we combined all information from thereference ROIs, and applied the Graph-Cut algorithm to assign those voxelsinto two groups.The Graph-Cut algorithm is mostly applied in the field of computer vision,where it could be used for image smoothing, segmentation, and many othercomputer vision problems which could be categorised under energy minimiza-tion problems. It is basically a combinatorial optimization technique basedon network flow theory (i.e., max-flow min-cut theory [24]). Once an energyminimization problem is formulated in terms of a network flow problem, manyenergy functions can be solved through the Graph-Cut algorithm.In our framework, there were two determined subROIs and several voxelsaround the boundary to be assigned, so our goal is to assign those voxels into two200groups. Firstly, we formalized the voxel assignment problem as the minimizationof the following energy function,minfE(f) =∑vDv(fv) +∑{v,q}∈NV (fv, fq) (5)where f represents labels for all voxels, fv is the label of voxel v,∑vDv(fv) isdata energy which measures the cost of labeling voxel v as fv,∑{v,q}∈N V (fv, fq)is the smooth energy that measures the extent to which the labeling f is notpiecewise smooth and N is the set of all pairs of adjacent voxels. In summary,the first term in equation (5) tries to ensure that each voxel has been labeled cor-rectly while the second term maintains the continuity of parcellation. Secondly,in order to solve the above energy minimization problem, we have to formulatethe problem as a network flow problem and define Dv(fv), V (fv, fq). Here weconstruct the voxel assignment problem as network flow problem based on [25]11Figure 2: Illustration of graph example in the voxel assignment problem, where each smallcube represents one voxel in an overlapped part of the task ROI. α and β are labels of twodetermined subROIs in the task ROI. Link ei,j between voxels exist if voxel i and voxel jare neighbours. tαi and tβi represent links between voxel i and determined subROI α and βrespectively.and solve it using the αβ-swap algorithm described in [26]. In our applicationof Graph-Cut, we define,Dv(fv) = 1− ρ(v, fv) (6)V (α, β) = T (α 6= β) (7)where ρ(v, fv) is the Pearson’s correlation between voxel v and region fv, regionsignal is calculated as mean time course of the obtained subROI fv, and T (.) isPotts model.In the αβ-swap algorithm, initial label f must be given. Here we initializethe label fv for voxel v as α if the number of voxel v labeled as α (nα) islarger than the number of that labeled as β (nβ), vice versa. If nα = nβ , weset the initial label randomly. Then we try to find a labeling fˆ that minimizeequation (5) over all labeling by swapping the labels for all voxels one by one.As described in Figure 2, let Pαβ be the set of voxels to be assigned. Thereare two links connecting each voxel with terminal α and β, noted as tαp and tβpwhere p indicates the voxel p. There is another possible link connecting adjacent12voxels p and q together, noted as e{p,q}. According to [25], the weights assignedto the links are,tαp = Dp(α) +∑q∈Nv, q /∈PαβV (α, fq) for v ∈ Pαβ (8)tβp = Dp(β) +∑q∈Nv, q /∈PαβV (β, fq) for v ∈ Pαβ (9)e{p,q} = V (α, β) for {v, q} ∈ N and v, q ∈ Pαβ (10)where Nv is the set of adjacent voxels of v, N is the set of all pairs of adjacentvoxel. Therefore, the energy function (5) can be solved through maximum-flowalgorithm on this graph [27].2.4. Extension to Regions with Three or More SubROIsAlthough the aforementioned algorithm was applied to our fMRI datasetwhere the task ROI is assumed to have two subROIs, the method can be easilyextended to other brain regions that may contain three or more subROIs, withsome small modifications: First, for the parcellation of the task ROI according toone reference ROI, we only need to modify the stop criteria. The algorithm willstop merging groups when the number of groups reaches the expected numberof subROIs. Second, the Graph-Cut approach is still applicable in this case. Inthe case with two subROIs, the αβ swap algorithm only applies on one pair oflabels, while in any case with three or more subROIs, the αβ swap algorithmneeds to run on all pairs of labels. For each pair, the process is exactly the sameas described earlier.3. Results3.1. Synthetic Data SetIn this section, we first tested the proposed method on synthetic datasetsby comparing its performance with those of k-means clustering as well as themodularity detection approach.13Figure 3: Illustration of the synthetic data generation. 1,2,3 indicate three connected referenceROIs. Y1, Y2, Y3 represent signals in reference ROIs 1,2,3 respectively. XA represents signal inregion A. XB represents signal in region B. Y and Z are source signals generated from normaldistribution with zero mean and unit variance. In order to comply with real fMRI signal, wefurther perform temporal Gaussian smoother on the source signal. Here we set β = 1 andα = 0.2.In order to ensure our simulations are compliant with real fMRI conditions,we generated three reference ROIs which are connected to each other and onecubical task ROI in which one spatially continuous subROI had strong con-nectivity to all three reference ROIs and the other part had weak connectionswith all reference ROIs. The boundary was set to separate these two subROIs.The process of generating the synthetic dataset was shown in Figure 3. The10×10×10 cubic task ROI contains 1000 voxels, and the sample length of eachvoxel was 300. Similarly, each reference ROI contained 300 voxels with datalength of 300.In the simulation, two different scenarios were considered. In the first case,the dataset was generated with the same signal to noise ratio (SNR) for allvoxels (syn-data 1). Here SNR = σ2signal/σ2noise where σ2 was the variance ofthe data. In the second simulation, we randomly chose 75 voxels in each subROIof the task ROI in synthetic data set 1. Then, we assigned those voxels withsmaller SNR, i.e., partially corrupted by outliers (syn-data 2). In order to ensure14the boundary is clear in synthetic dataset 2, we chose the outlier voxels that areaway from the boundary. Results from the two datasets were compared withthose of the k-means clustering method and the modularity detection method.Each procedure was repeated fifty times and the averaged performance of eachalgorithm was evaluated by looking at its error percentage, which is defined as,Error Percentage =# of all negative assigned voxels# of all voxels in task ROI(11)To implement the k-means algorithm according to [12], we first computedthe Pearson’s correlation coefficients between the task ROI and each referenceROI at the voxel level. Then, we transformed the Pearson’s correlation ρ intoFisher’s z statistics,z = 0.5× log 1 + ρ1− ρ (12)Then, each voxel of the task ROI will have a corresponding feature vector whereeach element represents the z score between current voxel and one voxel in areference ROI. k was set to 2 and Euclidean distance measure was chosen inthe simulation. For modularity detection method in [15], we constructed theundirected weighted graph by considering all voxels as nodes and their pairwisePearson’s correlation coefficients as edge weights. Since all algorithms performwell on syn-data 1, in this paper, we only demonstrated one representative resultfrom syn-data 2 in Figure 4.In addition, we performed k-means clustering with spatial regularization tofurther compare with our spatially regularized algorithm. In this paper, ac-cording to [28], we perform the spatially regularized k-means clustering methodby extending the feature vector of each voxel in task ROI with three spatialcoordinates, multiplied with a parameter p indicating the weights given to thespatial information. Figure 5 shows the performance of algorithm on syn-data2with different values of parameter p, form which we could infer that with spa-tial information incorporated result could be improved. However, the choice ofparameter p is crucial and a larger value would cause loss of robustness.From the results, we noted that, without the incorporation of spatial infor-15Table 2: Synthetic data setsSNR Configurationsyn-data1 0.5 normal spatial noisesyn-data2 0.5(full)+0.1(part) data with few outliersTable 3: Best error percentage from synthetics data setssyn-data1 syn-data2Spatially Regularized Fused Lasso 0 1%K-means Clustering 0 10.90%Modularity detection 0 11.30%Table 4: Average Error Percentage from Synthetics Data Setssyn-data1 syn-data2Spatially Regularized Fused Lasso 0.10% 2.9%K-means Clustering 0.11% 11.83%Modularity detection 0.09% 12.47 %mation, it was difficult to deal with noise that was not normally distributed andthe outliers can not be eliminated. As a result, a spatially continuous parcella-tion could not be guaranteed. Nevertheless, the spatially regularized fused lassomodel with the spatial regularization term can provide reliable parcellation withspatially continuous and functionally consistent subROIs.3.2. Real FMRI Data Set3.2.1. Data Description9 healthy subjects were recruited from the Pacific Parkinson’s Research Cen-tre (PPRC)/Movement Disorders Clinic at the University of British Columbia(UBC). All experiments were approved by the Ethics Board at UBC, and allsubjects provided informed consent prior to experiment participation. Duringthe experiment, all subjects were required to lie on their back in the scanner16(a) (b)(c) (d)Figure 4: One example of parcellation results from syn-data2. (a) ground truth. (b) results ofthe proposed spatially regularized fused lasso algorithm. (c) results of the k-means clustering.(d) results of the modularity detection.with their eyes closed and keep awake.MRI data were collected from a Philips Achieva 3.0 T scanner (Philips, Best,Netherlands) equipped with a headcoil. A whole brain three-dimensional T1-weighted image consisting of 170 axial slices with high resolution were acquired300to facilitate the anatomical localization for each individual. Blood oxygenationlevel-dependent (BOLD) contrast echo-planar (EPI) T2*-weighted images weretaken with the following specifications: repetition time 1985 ms, echo time 37ms, flip angle 90◦, field of view (FOV) 240.00mm, matrix size 128×128, andpixel size 1.9 mm×1.9 mm. The duration of each functional run was 8 minduring which we obtained 240 time points of 36 axial slices with 3 mm thicknessand 1 mm gap thickness. The FOV was set to include the cerebellum ventrallyand the dorsal surface of the brain.17Figure 5: Simulation results from spatially regularized k-means methods with different param-eter p chosen. The results are shown from average error rate perspective and also minimumerror rate perspective. The upper line (blue) shows the average error rate changing withspatial regularization parameter p, while the lower line (red) shows the minimum error ratetrending.Figure 6: Initial parcellation results for 9 normal subjects according to 3 reference ROIs. Eachcolumn shows the results from the same subject according to 3 different reference ROIs, whileeach row represents results from the same reference ROI across different subjects. In eachscatter plot, voxels are shown as small circles with different colors representing groups, andthe red color demonstrates significantly stronger correlations. It is worth noting that in ourfMRI data, as a result of not using the common spatial template, the spatial coordination ofthe putamen region changes from subject to subject.18After acquiring the raw fMRI data, we utilized our own optimized prepro-cessing pipeline of the fMRI data[29], which includes SPM-based slice timing cor-rection, isotropic reslicing (based on algorithms generated by our team), SPM-based motion correction and FSL-derived registration. 54 regions-of-interest(ROIs) were derived from the open source Freesurfer program (http://freesurfer.net/)for time course extraction. Rather than warping brain images into a commontemplate for connectivity analysis, all data analysis were done on the unwarpedimages (in the native space) on a subject-by-subject basis.3.3. Region of Interest(ROI) ExtractionWe employed the open source Freesurfer program and the available stan-dard atlases for cortical parcellation and subcortical segmentation. Freesurferuses a high-dimensional registration process by utilizing probabilistic atlases forcortical and subcortical labelling. The labels are then propagated back to thesubjects native space by relying on the template label and the subjects trans-formed voxel value obtained from T1-weighted structural scans. Based on ourexperience, this ROI-based segmentation method is superior to manually-drawnROIs and voxel-based method as it minimises registration error particularly insubcortical brain regions such as the putamen.For highly accurate structural analyses of subcortical structures, we havein the past used Large Deformation Diffeomorphic Metric Mapping (LDDMM)[28]. However, we note that we are using the Freesurfer segmentation of thehigh-resolution structural scans (1x1x1mm) to provide ROI-based binary masksto the relatively low-resolution fMRI scans (3x3x3 mm). We have found theFreesufer segmentations well-suited for this purpose. The extracted ROIs arevisually checked by experienced neurologists if needed.Table 5 lists all four of the regions-of-interest (ROIs) extracted from bothsides of the brain. For simplicity, the right column contains only Freesurfer andHuman Motor Area Template (HMAT) labels from the left side of the brain.HMAT labels are in bold. In this study, the orbitofrontal gyrus (OF) is definedby both the medial and lateral orbital frontal cortices, the cingulate gyrus in-19Table 5: List of ROIs Used in This PaperRegions of Interest (ROI) Freesurfer and HMAT labelsSensorimotor Area (SMA) ctx-lh-postcentral, L S1 L M1Cingulate gyrus (CG) ctx-lh-caudalanteriorcingulate,ctx-lh-posteriorcingulateOrbitofrontal gyrus (OF) ctx-lh-medialorbitofrontal,ctx-lh-lateralorbitofrontalPutamen Left-Putamenclude the anterior and posterior cingulate gyri, and finally the somatosensoryarea consists of the primary motor cortex as well as the somatosensory cortex.3.3.1. Algorithm ApplicationIn the real fMRI data application, in this pilot study, we only applied ouralgorithm on the left hemisphere. We choose OF, CG and SMA as referenceROIs as discussed before, to separate the putamen region into two functionalstriatal subROIs - the DLS and DMS. We assumed the motor control loops ofthe DLS and DMS were independent. Before performing any analysis, we firstremoved the global mean artifact in the putamen region and all three referenceROIs through linear regression. Then, in order to exploit the partial corre-lation between each reference ROI and the putamen region, we also removedartifacts from the SMA in the OF and CG and artifacts from the CG and OFin the SMA using the same method. The DLS is assumed to be functionallyconnected to the SMA. Therefore, we take SMA as a reference ROI to definethe DLS in the putamen region. In order to determine DLS, we calculatedthe correlation coefficients between each voxel contained in the putamen andmean time course of SMA. After transferring correlation coefficients to Fisher’sz-statistics, a one-sided Welch’s t-test is applied on magnitude of z-statistics totest whether one region has significantly stronger (under 5% significance level)correlation with the SMA and the region with stronger correlation is definedas the DLS. Similarly, to define the DMS in putamen region, OF and CG are20selected as the reference ROIs according to prior knowledge. The Graph-Cutalgorithm is essential in this case to combine two parcellation results based onconnections to the CG and OF. Figure 6 shows initial parcellation results from 9normal subjects respectively with SMA, CG and OF as the reference ROI. Eachpoint represents one voxel in 3-dimensional coordinates in Figure 6. Red andgreen colors denote two functional subROIs respectively with red representingstronger correlation. Figure 7 demonstrates the final definition of the DLS andDMS from 9 normal subjects with red representing the DMS and green repre-senting the DLS. In Figure 8, we showed our parcellation result of one subject(N003) with registration on the structural image.3.3.2. Performance and Group Consistency EvaluationSince we have no access to the ground truth of parcellation of real fMRIdata, we have to indirectly validate our results anatomically or functionally.Firstly, under the assumption that functional subROIs exist in the putamenregion, it is reasonable to assume that the anatomical boundaries across normalsubjects are stable. However, since our data do not use a common spatialtemplate for all subjects, spatial coordination of the putamen region changessubtly from subject to subject. As a result, it is not easy to evaluate theconsistency of anatomical segmentation. The final parcellation results in Figure7 show that the DMS region (red) is consistently located on the right side ofthe putamen region while all DLS points (green) are on the left side. Ourresults show that clear boundaries are consistently observed in normal subjects,however the specific positions of the boundaries vary from subject to subjectprobably due to the inter-subject variability. In addition, since our proposedalgorithm already incorporates prior knowledge (as shown in Figure 1) into theoptimization process, the definition of the DLS and DMS is strongly consistentwith prior knowledge. The results demonstrate that we could obtain spatiallycontinuous and functionally consistent subROIs and suggests that the proposedmethod could provide a promising way for brain parcellation.The adaptive changes of parameters of the proposed algorithm is shown in21Figure 7: Final parcellation results from 9 normal subjects, where red points represent DMSregion defined by combining initial results from CG and OF, green dots represent DLS regiondefined by initial results from SMA. Since DLS and DMS are defined by different parcellationregion, the two groups could be totally spatially separated or overlapped. For reference, thefirst graph(upper left) shows the full putamen region in subject N003 without any parcellation.Note that although the original accurate anatomical definition of putamen region may changesubtly from subject to subject, their geometrical shapes remain almost the same.(a) Coronal (b) Saggital (c) AxialFigure 8: One example of parcellation results in structural image from N003. Here the redregion represents the DLS while the yellow region shows the DMS.Figure 9 with the red line representing the number of groups, green and bluelines representing tuning parameters. Recall that in our algorithm, if the numberof groups could not be reduced, we will increase the penalty parameters. It’sworth noting that in the last few iterations, γ and λ increase dramatically sincethe algorithm is trying to combine those isolated outliers into the groups which22Figure 9: The iteration process in one parcellation (parcellating putamen according to SMAof one normal subject). The red line shows the change of number of groups in each iteration,the green and blue lines represent the change of the spatially regularized lasso parameter γand the normal lasso parameter λ respectively.demands large tuning parameters in the last several iterations,.4. Discussion and ConclusionIn this paper, we presented a novel framework for parcellating one brainROI into distinct functional subROIs using fMRI based functional connectivitypatterns with other reference ROIs. The proposed approach applies the spatiallyregularized fused lasso model to obtain functional connectivity between voxels400in the task ROI with other reference ROIs. With the learned connectivitypatterns, we grouped adjacent voxels with similar connectivity weights together.Then, groups containing similar connectivity adjacent voxels will be treated asnew voxels and the process repeated until the expected number of subROIs areobtained as described in Table 1. Lastly, the Graph-Cut algorithm is applied toassign voxels by combining the parcellation results from several reference ROIs.Several parameters are included in this proposed method. First, the param-eter λ controls the sparsity of elements in the connectivity coefficient vector, β.23This parameter can be set to 0 in order to consider only the sparsity difference.Assigning a larger λ to the algorithm does mean a higher tolerance of noise,however, it can also increase the possibility of merging all voxels into one group.In our simulation on synthetic data, we set it to 1 while on real dataset, it wasset to 10. Second, the parameter γ controls the sparsity of difference betweensuccessive elements in the connectivity vector β. It is important to note that γcan not be set as 0, as otherwise, the algorithm cannot converge. If the initial γvalue is set too high, too many voxels will be merged in the first step, resultingin possibly biased final results. In our experiment, we set γ to 10 on syntheticdataset whereas on real fMRI dataset, due to the existence of different spatialnoises in different subjects, we adjusted the parameter on a subject by subjectbasis. An appropriate choice of this parameter is recognized by well-balancedgrouping results; there should not be a single group that contains less than 10voxels or be located in the middle of other groups. The third parameter is themerging precision parameter, δ, which as its name suggests, controls the preci-sion of merging two adjacent voxels with similar connectivity weights. If δ isset too small, it risks merging of all voxels into only one group whereas a largeδ would make it difficult for the algorithm to converge. In our applications onsynthetic and real data, we set the initial value of δ to 104. It is worth notingthat even though we have adaptively adjusted all parameters, we still have tochoose the initial value carefully. Inappropriate choices of initial parametersmay lead to ill-balanced results.We have tested the proposed algorithm on synthetic datasets and comparedthe result with two other state-of-the-art methods. It is shown that the proposedalgorithm is comparable with these algorithms in general. In particular, theproposed method is better-suited to deal with cases which do not have normallydistributed noise or datasets corrupted with outliers. In a real fMRI data, weapplied the proposed framework to nine healthy subjects in order to define twofunctional subROIs of the putamen known as the DLS and DMS. The resultsare consistent with prior knowledge.We are motivated to develop this novel method to address the need in neurol-24ogy and neuroscience research, i.e., to parcellate the putamen into two functionalparts (DLS and DMS) based on its connectivity with three other reference ROIs.The proposed algorithm is applicable to a set of neurological problems wherethere exist one task ROI to be parcellated into several functional subROIs anda few reference ROIs that have different connectivity patterns with these func-tional subROIs. The proposed algorithm could be applied to both task-relateddata sets as well as resting fMRI data sets. Furthermore, with suitable parame-ter adjustments, we do not require the sophisticated signal denoising procedure.However, it is worth mentioning that, although our algorithm can better dealwith outliers compared with other two methods, outlier detection and removingis still essential for many connectivity analysis problems.Since the proposed method is based on functional connectivity, people mayargue that the results for subsequent connectivity analysis may be biased. Wecould reduce such risks of bias by separating time signals into two parts (e.g.,the odd-indexed points as one part and the even-indexed points as the otherpart), and use one part to obtain consistent parcellation results and then usethe other part of data for further connectivity analysis. We could also benefitfrom this two-part idea for the validation of the parcellation results: We canparcellate the functional ROI with one part of the fMRI time-series signals andvalidate the functional parcellations with the other part of the fMRI signals.Consistency measures such as the discrete entropy can be employed to evaluatethe consistency of the parcellation results from using different parts of the fMRIsignals.The proposed method can be applied to study geriatric and neuro-diseasedpopulations. The extracted functional subROIs themselves are of great interestto study the influence of aging and neurodegenerative diseases. In future study,we plan to use the proposed method to investigate the changes of DLS subROIsin subjects with Parkinson’s disease, where it may serve as a potential biomarkerof Parkinson’s disease.25AcknowledgementThis work was supported by Natural Sciences and Engineering ResearchCouncil of Canada (NSERC) and Canadian Institutes of Health Research (CIHR)and National Natural Science Foundation of China (Grant No. 81571760 and61501164).References[1] M. J. McKeown, S. Makeig, G. G. Brown, T.-P. Jung, S. S. Kindermann,A. J. Bell, T. J. Sejnowski, Analysis of fmri data by blind separation intoindependent spatial components, Tech. rep., DTIC Document (1997).[2] K. J. Friston, L. Harrison, W. Penny, Dynamic causal modelling, Neuroim-age 19 (4) (2003) 1273–1302.[3] A. Mclntosh, F. Gonzalez-Lima, Structural equation modeling and its ap-plication to network analysis in functional brain imaging, Human BrainMapping 2 (1-2) (1994) 2–22.[4] M. P. Van Den Heuvel, H. E. Hulshoff Pol, Exploring the brain network:a review on resting-state fmri functional connectivity, European Neuropsy-chopharmacology 20 (8) (2010) 519–534.[5] A. Liu, X. Chen, M. J. McKeown, Z. J. Wang, A sticky weighted regressionmodel for time-varying resting state brain connectivity estimation, IEEETransactions on Biomedical Engineering (2014) (in press).[6] A. Liu, X. Chen, Z. J. Wang, Q. Xu, S. Appel-Cresswell, M. J. McKe-own, A genetically-informed, group fmri connectivity modeling approach:application to schizophrenia, IEEE TRANSACTIONS ON BIOMEDICALENGINEERING 61 (3) (2013) 946 – 956.[7] A. Liu, J. Li, Z. J. Wang, M. J. McKeown, A computationally efficient,exploratory approach to brain connectivity incorporating false discovery26rate control, a priori knowledge, and group inference, Computational andmathematical methods in medicine 2012.[8] A. L. Cohen, D. A. Fair, N. U. Dosenbach, F. M. Miezin, D. Dierker, D. C.Van Essen, B. L. Schlaggar, S. E. Petersen, Defining functional areas inindividual human brains using resting functional connectivity mri, Neu-roimage 41 (1) (2008) 45–57.500[9] J. Rademacher, V. Caviness, H. Steinmetz, A. Galaburda, Topographicalvariation of the human primary cortices: implications for neuroimaging,brain mapping, and neurobiology, Cerebral Cortex 3 (4) (1993) 313–329.[10] M. Beckmann, H. Johansen-Berg, M. F. Rushworth, Connectivity-basedparcellation of human cingulate cortex and its relation to functional spe-cialization, The Journal of neuroscience 29 (4) (2009) 1175–1190.[11] R. Baumgartner, G. Scarth, C. Teichtmeister, R. Somorjai, E. Moser, Fuzzyclustering of gradient-echo functional mri in the human visual cortex. parti: Reproducibility, Journal of Magnetic Resonance Imaging 7 (6) (1997)1094–1101.[12] T. Kahnt, L. J. Chang, S. Q. Park, J. Heinzle, J.-D. Haynes, Connectivity-based parcellation of the human orbitofrontal cortex, The Journal of Neu-roscience 32 (18) (2012) 6240–6250.[13] L. Wang, Q. Liu, H. Li, D. Hu, Functional connectivity-based parcella-tion of human medial frontal cortex via maximum margin clustering, in:Intelligent Science and Intelligent Data Engineering, Springer, 2013, pp.306–312.[14] X. Shen, X. Papademetris, R. T. Constable, Graph-theory based parcella-tion of functional subunits in the brain from resting-state fmri data, Neu-roimage 50 (3) (2010) 1027–1035.[15] K. A. Barnes, A. L. Cohen, J. D. Power, S. M. Nelson, Y. B. Dosenbach,F. M. Miezin, S. E. Petersen, B. L. Schlaggar, Identifying basal ganglia di-27visions in individuals using resting-state functional connectivity mri, Fron-tiers in systems neuroscience 4(2010).[16] R. Tibshirani, M. Saunders, S. Rosset, J. Zhu, K. Knight, Sparsity andsmoothness via the fused lasso, Journal of the Royal Statistical Society:Series B (Statistical Methodology) 67 (1) (2005) 91–108.[17] P. Griffiths, R. Perry, A. Crossman, A detailed anatomical analysis of neu-rotransmitter receptors in the putamen and caudate in parkinson’s diseaseand alzheimer’s disease, Neuroscience letters 169 (1) (1994) 68–72.[18] A. J. Gruber, R. J. McDonald, Context, emotion, and the strategic pursuitof goals: interactions among multiple brain systems controlling motivatedbehavior, Frontiers in behavioral neuroscience 6 (2012).[19] P. Redgrave, M. Rodriguez, Y. Smith, M. C. Rodriguez-Oroz, S. Lehericy,H. Bergman, Y. Agid, M. R. DeLong, J. A. Obeso, Goal-directed andhabitual control in the basal ganglia: implications for parkinson’s disease,Nature Reviews Neuroscience 11 (11) (2010) 760–772.[20] R. Tibshirani, Regression shrinkage and selection via the lasso, Journal ofthe Royal Statistical Society. Series B (Methodological) (1996) 267–288.[21] G. Tononi, A. R. McIntosh, D. P. Russell, G. M. Edelman, Functionalclustering: identifying strongly interactive brain regions in neuroimagingdata, Neuroimage 7 (2) (1998) 133–149.[22] S. Kim, K.-A. Sohn, E. P. Xing, A multivariate regression approach toassociation analysis of a quantitative trait network, Bioinformatics 25 (12)(2009) i204–i212.[23] X. Chen, Q. Lin, S. Kim, J. G. Carbonell, E. P. Xing, et al., Smoothingproximal gradient method for general structured sparse regression, TheAnnals of Applied Statistics 6 (2) (2012) 719–752.28[24] C. H. Papadimitriou, K. Steiglitz, Combinatorial optimization: algorithmsand complexity, Courier Dover Publications, 1998.[25] Y. Boykov, O. Veksler, R. Zabih, Fast approximate energy minimization viagraph cuts, Pattern Analysis and Machine Intelligence, IEEE Transactionson 23 (11) (2001) 1222–1239.[26] Y. Boykov, V. Kolmogorov, An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision, Pattern Analysis andMachine Intelligence, IEEE Transactions on 26 (9) (2004) 1124–1137.[27] V. Kolmogorov, R. Zabin, What energy functions can be minimized viagraph cuts?, Pattern Analysis and Machine Intelligence, IEEE Transactionson 26 (2) (2004) 147–159.[28] F. Deleus, M. M. Van Hulle, A connectivity-based method for definingregions-of-interest in fmri data, Image Processing, IEEE Transactions on18 (8) (2009) 1760–1771.[29] Chen, Jingyun, et al. ”Freesurfer-initialized large deformation diffeomor-phic metric mapping with application to Parkinson’s disease.” SPIE Med-ical Imaging. International Society for Optics and Photonics, 2009.29


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items