1 Introduction
In the usual tensorbased morphometry (TBM), the spatial derivatives of deformation fields obtained during nonlinear image registration for warping individual magnetic resonance imaging (MRI) data to a template is used in quantifying neuroanatomical shape variations [3, 73, 20]. The Jacobian determinant of a deformation field is most frequently used in quantifying the brain tissue growth or atrophy at a voxel level. [22, 55, 26] used the Jacobian determinant of the 2D deformation field as a measure of local areachange at each pixel in 2D crosssections of the corpus callosum. [74, 20]
applied the Jacobian of 3D deformations as a measure of the regional growth. Subsequently, the statistical parametric maps are obtained by fitting the tensor maps as a response variable in a linear model at each voxel, which results in a massive number of univariate test statistics.
Recently, there have been attempts at explicitly modeling the structural variation of one region to another [64, 11, 79, 80, 51, 39, 38]. This provides additional information that complement existing univariate approaches. In most of these multivariate approaches, anatomical measurements such as mesh coordinates, cortical thickness or Jacobian determinant across different voxels are correlated using models such as canonical correlations [4, 64], crosscorrelations [11, 79, 80, 51, 39, 38], partial correlations, which are equivalent to the inverse of covariances [8, 6, 31, 41, 49]. However, these multivariate techniques suffer the smalln largep problem [32, 69, 76, 49, 17]
. Specifically, when the number of voxels are substantially larger than the number of images, it produces an underdetermined linear model. The estimated covariance matrix is rank deficient and no longer positive definite. In turn, the resulting correlation matrix is not considered as good approximations to the true correlation matrix.
The small large problem can be remedied by using sparse methods, which regularize the underdetermined linear model with additional sparse penalties. There exist various sparse models: sparse correlation [49, 17], sparse partial correlation [8, 41, 49], sparse canonical correlation [4] and L1norm penalized loglikelihood [7, 6, 31, 42, 56, 78]. Sparse model is usually parameterized by a tuning parameter that controls the sparsity of the representation. Increasing the sparse parameter makes the solution more sparse. So far, all previous sparse network approaches use a fixed parameter that may not be optimal. Depending on the choice of the sparse parameter, the final statistical results will be different. Instead of performing statistical inference at one fixed sparse parameter that may not be optimal, we introduce a new framework that performs statistical inferences over the whole parameter space using persistent homology [12, 28, 70, 33, 18, 17, 47, 48].
Persistent homology is a recently popular branch of computational topology with applications in protein structures [66], gene expression [25], brain cortical thickness [18], activity patterns in visual cortex [70], sensor networks [23], complex networks [40] and brain networks [47, 48]. However, as far as we are aware, it is yet to be applied to sparse models in any context. This is the first study that introduce persistent homology in sparse models. The proposed persistent homological framework is similar to the existing multithresholding framework that has been used in modeling connectivity matrices at many different thresholds [1, 38, 72, 48]. However, such approach has not been applied in sparse networks before. In a sparse network, sparsity is controlled by the sparse parameter and the estimated sparse matrix entries. So it is unclear how the existing multithresholding framework can be applicable in this situation. In this paper, we establish that thresholding the sparse parameter is equivalent to thresholding correlations under some conditions. Thus, we resolve the unclarity of applying the existing multithreshold method to the sparse networks.
The main methodological contributions of this paper are as follows. (i) We introduce a new sparse model based on Pearson correlation. Although various sparse models have been proposed for other correlations such as partial correlations [8, 41, 49] and canonical correlations [4], the sparse version of the Pearson correlation was not often studied.
(ii) We introduce persistent homology in the proposed sparse model for the first time. We explicitly show that persistent homological structures can be found in the sparse model. This paper differs substantially from our previous study [48], which studies the persistent homology in graphs and networks. However, sparse models and sparse networks were never considered previously.
(iii) We show that the identification of persistent homological structures can yield greater computational speed and efficiency in solving the proposed sparse Pearson correlation model without any numerical optimization. Note that most sparse models require numerical optimization for minimizing L1norm penalty, which can be a computational bottleneck for solving large scale problems. There are few attempts at speeding up the computation for sparse models. By identifying block diagonal structures in the estimated (inverse) covariance matrix, it is possible to bypass the numerical optimization in the penalized loglikelihood method [56, 78]. LASSO (least absolute shrinkage and selection operator) can be done without numerical optimization if the design matrix is orthogonal [75]. The proposed method substantially differs from [56, 78] in that we do not need to assume the data to follow normality since there is no need to specify the likelihood function. Further the cost functions we are optimizing are different. The proposed method also differs from [75] in that our problem is not orthogonal.
As an application of the proposed method, we applied the techniques to the quantification of interregional white matter abnormality in stressexposed children’s magnetic resonance images (MRI). Early and severe childhood stress, such as experiences of abuse and neglect, have been associated with a range of cognitive deficits [60, 67, 53] and structural abnormalities [43, 37, 36]. However, little is known about the underlying biological mechanisms leading to cognitive problems in these children [61] due to the difficulties in the existing methods that do not have enough discriminating power. However, we demonstrate that the proposed method is very well suited to this problem.
2 Methods
2.1 Sparse Correlations
Correlations.
Consider measurement vector
on node . If we center and rescale the measurement such thatthe sample correlation between nodes and is given by . Since the data is normalized, the sample covariance matrix is reduced to the sample correlation matrix.
Consider the following linear regression between nodes
and :(1) 
We are basically correlating data at node to data at node . In this particular case, is the usual Pearson correlation. The least squares estimation (LSE) of is then given by
(2) 
which is the sample correlation. For the normalized data, regression coefficient estimation is exactly the sample correlation. For the normalized and centered data, the regression coefficient is the correlation. It can be shown that (2) minimizes the sum of least squares over all nodes:
(3) 
Note that we do not really care about correlating to itself since the correlation is then trivially .
Sparse Correlations. Let be the correlation matrix. The sparse penalized version of (3) is given by
(4) 
The sparse correlation is given by minimizing . By increasing , the estimated correlation matrix becomes more sparse. When , the sparse correlation is simply given by the sample correlation, i.e. . As increases, the correlation matrix shrinks to zero and becomes more sparse. This is separable compressed sensing or LASSO type problem. Further, there is no need to numerically optimize (4) using the coordinate descent learning or the activeset algorithm often used in compressed sensing [59, 31]. The minimization of (4) can be done by the proposed softthresholding method analytically by exploiting the topological structure of the problem. This sparse regression is not orthogonal, i.e. , the Dirac delta, so the existing softthresholding method for LASSO [75] is not applicable.
Theorem 1.
For , the solution of the following separable LASSO problem
is given by the softthresholding
(5) 
Proof.
Write (4) as
(6) 
where
Since is nonnegative and convex, is minimum if each component achieves minimum. So we only need to minimize each component . This differentiates our sparse correlation formulation from the standard compressed sensing that cannot be optimized in this component wise fashion. can be rewritten as
We used the fact
For , the minimum of is achieved when , which is the usual LSE. For , Since is quadratic in , the minimum is achieved when
(7) 
The sign of depends on the sign of . Thus, sparse correlation is given by a softthresholding of :
(8) 
∎
Theorem 1 is heuristically explained in the conference paper
[17] without the rigorous proof or statement. This paper extends [17] with clearly spelled out softthresholding rule and the detailed proof. The estimated sparse correlation (8) basically thresholds the sample correlation that is larger or smaller than by the amount . Due to this simple expression, there is no need to optimize (4) numerically as often done in compressed sensing or LASSO [59, 31]. However, Theorem 1 is only applicable to separable cases and for nonseparable cases, numerical optimization is still needed.Since different choices of sparsity parameter will produce different solutions in sparse model , we propose to use the collection of all the sparse solutions for many different values of for the subsequent statistical inference. This avoids the problem of identifying the optimal sparse parameter that may not be optimal in practice. The question is then how to use the collection of in a coherent mathematical fashion. For this, we propose to apply persistent homology [27, 47, 48] and establish Theorem 2.
2.2 Persistent Homology in Graphs
Using persistent homology, topological features such as the connected components and cycles of a graph can be tabulated in terms of the Betti numbers. The zeroth Betti number and the first Betti number , which are topological invariants, respectively counts the the number of connected components and holes in the graph [28]. The network difference is then quantified using the Betti numbers of the graph [47, 48]
. The graph filtration is a new graph simplification technique that iteratively builds a nested subgraphs of the original graph. The algorithm simplifies a complex graph by piecing together the patches of locally connected nearest nodes. The process of graph filtration is mathematically equivalent to the single linkage hierarchical clustering and dendrogram construction
[47, 48].Consider a general setting of a weighted graph with node set and edge weights , where is the weight between nodes and . Weighted graph is formed by the pair of node set and edge weights . The edge weights in many brain imaging applications are usually given by some similarity measures such as correlation or covariance between nodes [47, 52, 57, 58, 71].
Given a weighted network , we induce binary network by thresholding the weighted network at . The adjacency matrix of is defined as
(9) 
Any edge weight less than or equal to is made into zero while edge weight larger than is made into one. The binary network is a simplicial complex consisting of simplices (nodes) and simplices (edges), a special case of the Rips complex [33]. Then it can be easily seen that for in a sense the vertex and edge sets of are contained in those of . Therefore, just as in the case of Rips filtration, which is a collection of nested Rips complexes, we can construct the filtration on the collection of binary networks:
(10) 
for Note that is the complete weighted graph while gives the node set . By increasing the value, we are thresholding at higher correlation so more edges are removed and thin out the connections. Such the nested sequence of the Rips complexes (10) is called a Rips filtration, the main object of interest in persistent homology [27]. The sequence of values are called the filtration values. Since we are dealing with a special case of Rips complexes restricted to graphs, we will call such filtration graph filtration. Figure 1 illustrates an example of a graph filtration with 4 nodes. Sequentially we are deleting edges based on the ordering of the edge weights. Since the graph filtration is a special case of the Rips filtration, it inherits all the topological properties of the Rips filtration. Given a weighted graph, there are infinitely many different filtrations. In Figure 1 example, we have two filtrations and among many other possiblities. So a question naturally arises if there is a unique filtration that can be used in characterizing the graph. Let the level of a filtration be the number of nested unique sublevel sets in the given filtration.
Theorem 2.
For graph with unique edge weights, the maximum level of a filtration on the graph is . Further, the filtration with filtration level is unique.
Proof.
For a graph with nodes, the maximum number of edges is , which is obtained in a complete graph. If we order the edge weights in the increasing order, we have the sorted edge weights:
where . The subscript denotes the order statistic. For all , is the complete graph of . For all , . For all , , the vertex set. Hence, the filtration given by
(11) 
is maximal in a sense that we cannot have any additional level of filtration. ∎
Among many possible filtrations, we will use the maximal filtration (11) in the study since it is uniquely given. The finiteness and uniqueness of the filtration levels over finite graphs are intutively clear by themselves and are already applied in software packages such as javaPlex. [2]. However, we still need a rigorous statement to specify the type of filtration we are using out of many.
2.3 Persistent Homology in Sparse Regression
We introduce a persistent homological structure in sparse correlations now as follows. Let be the adjacency matrix obtained from sparse correlation (8):
Let be the graph defined by the adjacency matrix . Then we have the main result of this paper, which relies on the results of Theorem 1 and Theorem 2.
Theorem 3.
For centered and normalized data , be the order statistic of , i.e. the sorted sequence of in increasing order. Then graph obtained from the sparse regression (4) forms the maximal graph filtration
(12) 
Proof.
The proof follows by simplifying the adjacency matrix into a simpler but equivalent adjacency matrix . From Theorem 1, if and 0 otherwise. Thus, the adjacency matrix is equivalent to the adjacency matrix :
(13) 
Let be the graph defined by adjacency matrix . Graph is formed by thresholding edge weights given by the absolute value of sample correlations . From Theorem 2, such graph must have maximal filtration:
(14) 
Since , graph also must have the identical maximal filtration. This proves the statement. ∎
Theorem 3 is illustrated in Figure 2. In obtaining the topological structure of a graph induced by sparse correlation, it is not necessary to solve the sparse regression by the direct optimization, which can be very time consuming. Identical topological information can be obtained by performing the softthresholding on the sample correlations. Figure 2 illustrates how Theorems 3 is used to construct a sparse correlation network using a 4nodes example. In the application, nodes will be used.
The resulting maximal filtration can be quantified by plotting the change of Betti numbers over increasing filtration values [28, 33, 47]. The first Betti number counts the number of connected components of the given graph at the filtration value [48]. Given graph filtration , we plot the first Betti numbers over filtration values (Figure 1). The number of connected components increase as the filtration value increases. The pattern of increasing number of connected components visually show how the graph structure changes over different parameter values. The overall pattern of Betti (number) plots can be used as a summary measure of quantifying how the graph changes over increasing edge weights. The Betti number plots are related but different from barcodes in litearture. The Betti number is equal to the number of bars in the barcodes at the specific filtration value. It is not necessary to perform filtrations for infinitely many possible values in plotting the Betti numbers. From Theorem 2, the maximum possible number of filtration level for plotting the Betti numbers is , where is the number of unique edge weights. For a tree, which is a graph with no cycle, we can come up with a much stronger statement than this.
Theorem 4.
For a tree with nodes and unique positive edge weights , the plot for the first Betti number () corresponding to the maximal graph filtration is given by the coordinates
Proof.
For a tree with nodes, there are total edges. Then from Theorem 2, we have the maximal filtration
(15) 
Since all the edge weights are above filtration value , all the nodes are connected, i.e., . Since no edge weight is above the threshold , . At each time we threshold an edge, the number of components increases exactly by one in the tree. Thus, we have
∎
For a general graph, it is not possible to analytically determine the coordinates for its Bettiplot. The best we can do is to compute the number of connected components numerically using the single linkage dendrogram method (SLD) [48], the DulmageMendelsohn decomposition [62, 16] or existing simplical complex approach [23, 12, 28]. For our study, we used the SLD method.
2.4 Statistical Inference on Betti number plots
The first Betti number will be used as features for characterizing network differences statistically. We assume there are subjects and nodes in Group 1. For subject , we have measurement at node . Denote data matrix as , where is the measurement for subject at node . We then construct a sparse network and corresponding Betti number using . Thus, is a function of . Consider another Group 2 consists of subjects. For Group 2, data matrix is denoted as , where is the measurement for subject at node . Group 2 will also generate single Betti number plot as a function of
. We are then interested in testing if the shapes of Betti number plots are different between the groups. This can be done by comparing the areas under the Betti plots. So the null hypothesis of interest is
(16) 
while the alternate hypothesis is
This inference avoids the use of multiple comparisons. The null hypothesis (16) is related to the following pointwise null hypothesis:
(17) 
If the hypothesis (17) is true, the hypothesis (16) is also true (but inverse is not true). Thus, testing the area under the curve is related to testing the height of the curve at every point. The advantage of using the area under the curve is that we do not need to worry about multiple comparisons associated with testing (17). The area under the curve seems a reasonable approach to use for Bettiplots. A similar approach has been introduced in [14] in removing the multiple comparisons and produce a single summary test statistic.
There is no prior study on the statistical distribution on the Betti numbers so it is difficult to construct a parametric test procedure. Further, since there is only one Bettiplot per group, it is not even possible to construct a statistic without resampling techniques. So it is necessary to empirically construct the null distribution and determine the pvalue by resampling techniques such as the permutation test and jackknife [17, 15, 29, 48]. For this study, we use the jackknife resampling.
For Group 1 with subjects, one subject is removed at a time and the remaining subjects are used in constructing a network and a Bettiplot. Let be the data matrix, where the th row (subject) is removed from . Then for each th subject removed, we compute , which is a function of and . Repeating this process for each subject, we obtain Bettiplots . For Group 2, the th row (subject) is removed from the original data matrix and obtain data matrix as . For each th subject removed, we compute , which is a function of and . Repeating this process for each subject, we obtain Bettiplots . There are 23 maltreated and 31 control children in the study, so we have 23 and 31 Jackknife resampled Bettiplots. Subsequently we compute the areas under the Bettiplots by discretizing the integral and doing the Riemann sum. The area differences between the groups are then tested using the Wilcoxon ranksum test, which is a nonparametric test on median differences [34].
We did not use the permutation test. For the permutation test to converge for our data set, it requires tens of thousands permutations and it is really time consuming even with the proposed timesaving softthresholding method. The proposed method takes about a minute of computation in a desktop but tenthousands permutations will take about seven days of computation. Hence, we used a much simpler Jackknife resampling technique. The procedure is validated using the simulation with the known ground truth. MATLAB codes for constructing network filtration, barcodes and performing statical inference on are given in http://brainimaging.waisman.wisc.edu/~chung/barcodes with the postprocessed Jacobian determinant and FA data that was used for this study.
Simulations. We performed two simulations. In each simulation, the sample sizes are in Group 1 and in Group 2. There are nodes. In Group 1, data at node for subject is simulated as independent standard normal for the both simulations.
Study1 (no group difference): In Group 2, we simulated data at node for subject as . Tiny noise is added to perturb Group 1 data a little bit. It is expected there is no group difference. Following the proposed framework, we constructed the sparse correlation networks and constructed a Bettiplot. Then performed the Jackknife resampling and constructed 20 Bettiplots in each group. The rank sum test was applied and obtained the value of 0.78. So we could not detect any group difference as expected.
Study 2 (group difference):
We first simulate data as independently for all the nodes. Then for four nodes indexed by , we introduce additional dependency:
We added small noise to perturb the node values further. This dependency gives any connection between 1 to 5 to have high correlation. Figure 3 shows the simulated correlation matrix. Following the proposed framework, we constructed the sparse correlation networks and constructed a Bettiplot. Then performed the Jackknife resampling and and constructed 20 Bettiplots in each group. The rank sum test was applied and obtained the value less than 0.001. This significance corresponds to the horizontal gap between the Bettiplots after the filtration value 0.7 (Figure 3 right).
3 Application
3.1 Imaging Data Set and Preprocessing
The study consists of 23 children who experienced documented maltreatment early in their lives, and 31 agematched normal control (NC) subjects. All the children were recruited and screened at the University of Wisconsin. The maltreated children were raised in institutional settings, where the quality of care has been documented as falling below the standard necessary for healthy human development. For the controls, we selected children without a history of maltreatment from families with similar current socioeconomic statuses. The exclusion criteria include, among many others, abnormal IQ (), congenital abnormalities (e.g., Down syndrome or cerebral palsy) and fetal alcohol syndrome (FAS). The average age for maltreated children was 11.26 1.71 years while that of controls was 11.58 1.61 years. This particular age range is selected since this development period is characterized by major regressive and progressive brain changes [50, 36]. There are 10 boys and 13 girls in the maltreated group and 18 boys and 13 girls in the control group. Groups did not differ on age, pubertal stage, sex, or socioeconomic status [36]. The average amount of time spent in institutional care by children was 2.5 years 1.4 years, with a range from 3 months to 5.4 years. Children were on average 3.2 years 1.9 months when they adopted, with a range of 3 months to 7.7 years. Table 1 summarizes the participant characteristics.
T1weighted MRI were collected using a 3T General Electric SIGNA scanner (Waukesha, WI), with a quadrature birdcage head coil. DTI were also collected in the same scanner using a cardiacgated, diffusionweighted, spinecho, singleshot, EPI pulse sequence. The details on image acquisition parameters are given in [36]. Diffusion tensor encoding was achieved using twelve optimum noncollinear encoding directions with a diffusion weighting of 1114 s/mm and a nonDW T2weighted reference image. Other imaging parameters were TE = 78.2 ms, 3 averages (NEX: magnitude averaging), and an image acquisition matrix of 120 120 over a field of view of 240 240 mm. To minimize field inhomogeneity and image artifacts, high order shimming and fieldmap images were collected using a pair of nonEPI gradient echo images at two echo times: TE1 = 8 ms and TE2 = 11 ms. For MRI, a study specific template was constructed using the diffeomorphic shape and intensity averaging technique through Advanced Normalization Tools (ANTS) [5]
. Image normalization of each individual image to the template was done using symmetric normalization with crosscorrelation as the similarity metric. Then the Jacobian determinants of the inverse deformations from the template to individual subjects were computed at each voxel. The Jacobian determinants measure the amount of voxelwise change from the template to the individual subjects. White matter was also segmented into tissue probability maps using templatebased priors, and registered to the template
[9]. For DTI, images were corrected for eddy current related distortion and head motion via FSL software (http://www.fmrib.ox.ac.uk/fsl) and distortions from field inhomogeneities were corrected using custom software based on the method given in [44] before performing a nonlinear tensor estimation using CAMINO [21]. Subsequently, we have used iterative tensor image registration strategy given in [82] and [45] for spatial normalization. Then Fractional anisotropy (FA) were calculated for diffusion tensor volumes diffeomorphically registered to the study specific template.The proposed methods have been applied to MRI and DTI of stressexposed children in characterizing the white matter structural differences against the normal controls.
Maltreated  Normal controls  

Sample size  23  31 
Sex (males)  10  18 
Age (years)  11.26 1.71  11.58 1.61 
Duration (years)  2.5 1.4 (0.25 to 5.4)  
Time of adoption (years)  3.2 1.9 (0.25 to 7.7) 
3.2 Results: Proposed Sparse Correlation
We threshold the white matter density at 0.7 and obtained the isosurface. The resulting isosurface is not the gray and white matter tissue boundary and it is located inside the white matter. We are interested in the white matter changes along the tissue boundary. The surface mesh has 189536 mesh vertices and the average internodal distance of 0.98mm. Since Jacobian determinant and FA values at neighboring voxels are highly correlated, 0.3 of the total mesh vertices are uniformly sampled to produce nodes. This gives average internodal distance of 15.7mm, which is large enough to avoid spurious high correlation between two adjacent nodes (Figure 4). The isosurface of the white matter template was extracted using the marching cube algorithm [54]. The number of nodes are still larger than most region of interest (ROI) approaches in MRI and DTI, which usually have around 100 regions [81]. This resulted in sample covariances and correlation matrices, which are not full rank. We constructed the sparse correlation based network filtrations from the softthresholding method using Theorem 3 (Figure 5). Subsequently, Theorem 4 is used to plot the corresponding Bettiplots (Figure 6). Since each group produces one Bettiplot, the leaveoneout Jackknife resampling technique was performed to produce 23 and 31 resampled Bettiplots respectively for the two groups. We then computed the areas under the Bettiplots. Using the ranksum test, we detected the statistical significance of the area differences between the groups (value 0.001). The Bettiplots for normal controls show much higher Betti numbers at any given threshold.
Biological Interpretation. In the Bettiplots (Figures 6), we obtain more disconnected components for controls than for children in the early stress group for any specific value. It can only happen if Jacobian determinants show higher correlations in the maltreated children across the white matter compared to the controls. So when thresholded at a specific correlation value, more edges are preserved in the maltreated children resulting in decreased number of disconnected components. Thus, the children exposed to early life stress and maltreatment show more dense network at a given value. This is clearly demonstrated visually in Figure 5. If the variations of Jacobian determents are similar across voxels, we would obtain higher correlations. This suggests more anatomical homogeneity across whole brain white matter in the maltreated children. Our finding is consistent with the previous study on neglected children that shows disrupted white matter organization, which results in more diffuse connections between brain regions [36]. Lower white matter directional organization in white matter may correspond to the increased homogeneity of Jacobian determinants and FAvalues across the brain regions. Similar experiences may cause some areas to be connected to other regions of the brain at a higher degree inducing higher homogeneity in the regions. This type of relational interpretation can be obtained from the traditional univariate TBM at each voxel.
3.3 Comparison Against Sparse Covariance
We compared the performance of the proposed sparse correlation method to the existing sparse (inverse) covariance method via the penalized loglikelihood [7, 6, 31, 42, 56], where the loglikelihood is regularized with a L1norm penalty:
(18) 
is the covariance matrix and is the sample covariance matrix. is the sum of the absolute values of the elements. The penalized loglikelihood is maximized over the space of all possible symmetric positive definite matrices. (18) is a convex problem and it is numerically optimized using the graphicalLASSO (GLASSO) algorithm [7, 6, 31, 42]. The tuning parameter controls the sparsity of the offdiagonal elements of the covariance matrix. By increasing , the estimated covariance matrix becomes more sparse.
We also performed the graph filtration technique to the estimated sparse covariance matrix . Let be the adjacency matrix defined from the estimated sparse covariance:
(19) 
The adjacency matrix induces graph consisting of number of partitioned subgraphs:
(20) 
where and are vertex and edge sets of the subgraph respectively. Unlike the sparse correlation case, we do not have full persistency on the induced graph . The partitioned graphs can be proven to be partially nested in a sense that only the partitioned node sets are persistent [17, 42, 56], i.e.
(21) 
for and all . Subsequently the collection of partitioned vertex set is also persistent. On the other hand, edge sets may not be persistent. It is unclear if there exists a unique maximal filtration on the vertex set.
The maximal filtration can be obtained as follows. Let be another adjacency matrix given by
(22) 
where is the sample covariance matrix. It can be shown that the adjacency matrix similarly induces graph [17, 56]:
(23) 
for some edge set . The subgraphs and have identical vertex set but different edge sets. Then from Theorem 2, we have maximal filtration on the graph , where the edge weights are given by the sample covariances. Theorem 2 requires the edge weights to be all unique, which is satisfied for the study data set. Then similar to Theorem 3, the Bettiplots are determined by ordering the entries of the sample covariance matrices. The resulting barcode is displayed in Figure 6. The sparse covariance was also able to discriminate the groups statistically (value ). The changes in the first Betti number are occurring in a really narrow window but was still able to detect the group difference using the areas under the Betti number plots (Figure 6). However, the sparse correlations exhibit slower changes in the Betti number over the wide window, making it easier to discriminate the groups.
3.4 Comparison Against Fractional Anisotropy in DTI
For children who suffered early stress, white matter microstructures have been reported as more diffusely organized [36]. Therefore we predicted less white matter variability in both the Jacobian determinants and also in fractional anisotropy (FA) values as well. The DTI acquisitions were done in the same 3T GE SIGNA scanner; acquisition parameters can be found in [36]. We applied the proposed persistent homological method in obtaining the filtrations for sparse correlations and covariances in the same 548 nodes on FA values (Figure 4). The resulting filtration patterns show similar patterns of a rapid increase in disconnected components for sparse correlations (Figure 5 and 6). The Jackknife resampling followed by the ranksum test on the area differences shows a significant group difference for sparse correlations (value 0.001). These results are due to a consistent abnormality among the stressexposed children that is observed in both MRI and DTI modalities. The stressexposed children exhibited stronger white matter homogeneity and less spatial variability compared to normal controls in both MRI and DTI measurements. However, the covariance results fail to discriminate the groups at level (value 0.043) indicative of a poor performance compared to the sparse correlation method.
3.5 Robustness on Node Size Changes
Depending on the number of nodes, the parameters of graph vary considerably up to 95 and the resulting statistical results will change substantially [30, 35, 81]. On the other hand, the proposed method is very robust under the change of node size. For the node sizes between 548 and 1856 (0.3 and 1 of original 189536 mesh vertices), the choice of node size did not affect the pattern of graph filtrations, the shape of Bettiplots, or the subsequent statistical results significantly. For example, the graph filtration on 1856 nodes shows a similar pattern of dense connections for the maltreated children (Figure 5). The resulting Bettiplots also show similar pattern of the group separation (Figure 6). The statistical results are also somewhat consistent. For both the Jacobian determinant and FA values, the group differences in Bettiplots obtained from sparse correlations and covariances are all statistically significant (value ) in both 548 and 1856 nodes except one case. For the case of the 548 nodes covariance on FA values, we did not detect any group differences at 0.01 level (value = 0.043). On the other hand, we detected the group difference for the 1856 nodes case at 0.001 level. The proposed framework is very sensitive, so it can detect really narrow but very consistent Bettiplot differences.
3.6 Effect of Image Registration
We checked how much impact image registration has on the robustness of the proposed method. Anatomical measurements across neighboring voxels are highly correlated within white matter so we do not expect image misalignment will have much effect on the final results. To determine the variability associated with the image registration, the displacement vector fields from the template to individual brains were randomly perturbed by adding Gaussian noise to each component. This is sufficiently large noise and causes up to 4mm misalignment for some nodes. Then following the proposed pipeline, the Jacobian determinants are correlated across 548 nodes and Bettiplots are computed. Figure 7 shows five perturbation results. The thick line is without any perturbation. The perturbed Bettiplots are very stable and close to the Bettiplots without any perturbation (thick lines). Th height differences in the perturbed Bettiplots are less than in average, which is negligible in the subsequent analysis. In fact, the resulting values are similar to each other and all the perturbed results detected the group difference (value 0.001). Thus, we conclude that the proposed topological framework is robust under image misalignment.
4 Conclusions and Discussions
By identifying persistent homological structures in sparse Pearson correlation, we were able to exploit them for speeding up computations. A procedure that usually takes 56 hours was completed in few seconds without utilizing additional computational resources. Although we have only shown how to identify persistent homology in the sparse Pearson correlation, the underlying principle can be directly applicable to other sparse models and image filtering techniques. These include the least angle regression (LARS) implementation in more general LASSO [13], heat kernel smoothing [19], and diffusion wavelets [46], which all guarantee the nested subset structure over the sparse parameters and kernel bandwidth. We will leave the identification of persistent homology in other frameworks for future studies.
We found that Bettiplots on correlations can visually discriminate better than Bettiplots on covariances. In Figure 6
, almost all topological changes associated with the covariance occur in really small range between 0 and 0.1. However, unlike covariances, correlations are normalized by the variances so the topological changes are more spread out between 0 and 1. This has the effect of making the Bettiplots shape differences spread out more uniformly and wide in the unit interval. This is most clearly demonstrated in the covariance
vs. correlation on FA (second column). The Bettiplots of covariances are difficult to discriminate visually because the Bettiplots are squeezed into small range between 0 and 0.1 but the Bettiplots of correlations are more discriminative since the Bettiplots are more spread out. The visual discriminative power comes from the normalization associated with the Pearson correlation. The change in the metric affects the filtration process itself since it is based on the sorted edge weights. Subsequently, the shape of Bettiplots and the statistical inference results also change.While massive univariate approaches can detect signal locally at each voxel, the proposed graph approach can detect signal globally over the whole brain region. Even though the information obtained by the two methods are complementary, they are somewhat exclusive. The proposed approach tabulates the changes of the number of connected components in the thresholded networks via Bettiplots, which cannot be done at individual node level. Therefore, there is no easy straightforward way of combining or comparing the results from the two methods. The Bettiplots is a global index that is defined over a whole graph so it cannot be directly applicable to nodelevel analysis. However, just like any global graph theoretic indices such as smallworldness and modularity [65, 10], it can be applied to subgraphs around a given node. Thus, it might be possible to measure logical topological characteristic around the node. This is the beyond the scope of the paper and we left it as future research.
This paper is not concerned with white matter anatomical connectivity. Here, we focus on a different issue, namely the degree of interregional dependency of image measurements such as Jacobian determinant and fractional anisotropy across brain regions. The proposed method is general enough to run on any type of volumetric imaging data that is spatially normalized. As an application of the proposed method, we were able to demonstrate developmental differences in brain development among stressexposed children, who are at known risk for cognitive delays. Our Jacobian determinant results are consistent with DTI.
There are recent fMRI studies showing head motion to introduce systematic biases in functional connectivity [63, 68, 77]. Motion makes it appear as if longrange connections are weaker than they really are, and shortrange connections are stronger than they really are [24]. However, unlike fMRI, structural image volumes are acquired across such a long time frame, we do not expect the head motion to introduce spurious correlations. Further, we are also not aware of any study that establishes a relationship between head movement and maltreated children. We do not consider the head motion is a concern for our anatomical studies.
Acknowledgements
This work was supported by NIH grants MH61285, MH68858, MH84051, NIH Fellowship DA028087, NIHNCATS grant UL1TR000427 and the Vilas Associate Award. The authors like to thank Hyekyoung Lee of Seoul National University and Matthew Arnold of University of Bristol for the valuable discussions on persistent homology and sparse regressions, and SeungGoo Kim of Max Planck Institute and Nagesh Adluru of University of WisconsinMadison for help with image preprocessing.
References
 [1] S. Achard, R. Salvador, B. Whitcher, J. Suckling, and E.D. Bullmore. A resilient, lowfrequency, smallworld human brain functional network with highly connected association cortical hubs. The Journal of Neuroscience, 26:63–72, 2006.
 [2] H. Adams, A. Tausz, and M. VejdemoJohansson. javaPlex: A research software package for persistent (co) homology. In Mathematical Software–ICMS 2014, pages 129–136. Springer, 2014.
 [3] J. Ashburner and K. Friston. Voxelbased morphometry  the methods. NeuroImage, 11:805–821, 2000.
 [4] B.B. Avants, P.A. Cook, L. Ungar, J.C. Gee, and M. Grossman. Dementia induces correlated reductions in white matter integrity and cortical thickness: a multivariate neuroimaging study with sparse canonical correlation analysis. NeuroImage, 50:1004–1016, 2010.
 [5] B.B. Avants, C.L. Epstein, M. Grossman, and J.C. Gee. Symmetric diffeomorphic image registration with crosscorrelation: Evaluating automated labeling of elderly and neurodegenerative brain. Medical Image Analysis, 12:26–41.

[6]
O. Banerjee, L. El Ghaoui, and A. d’Aspremont.
Model selection through sparse maximum likelihood estimation for
multivariate Gaussian or binary data.
The Journal of Machine Learning Research
, 9:485–516, 2008.  [7] O. Banerjee, L.E. Ghaoui, A. d’Aspremont, and G. Natsoulis. Convex optimization techniques for fitting sparse Gaussian graphical models. In Proceedings of the 23rd international conference on Machine learning, page 96, 2006.
 [8] P.J. Bickel and E. Levina. Regularized estimation of large covariance matrices. The Annals of Statistics, 36:199–227, 2008.
 [9] M.F Bonner and M. Grossman. Gray matter density of auditory association cortex relates to knowledge of sound concepts in primary progressive aphasia. The Journal of Neuroscience, 32:7986–7991, 2012.
 [10] E. Bullmore and O. Sporns. Complex brain networks: graph theoretical analysis of structural and functional systems. Nature Review Neuroscience, 10:186–98, 2009.
 [11] J. Cao and K.J. Worsley. The geometry of correlation fields with an application to functional connectivity of the brain. Annals of Applied Probability, 9:1021–1057, 1999.
 [12] G. Carlsson and F. Memoli. Persistent clustering and a theorem of J. Kleinberg. arXiv preprint arXiv:0808.2241, 2008.
 [13] M.K. Carroll, G.A. Cecchi, R. Irina, R. Garg, and A.R. Rao. Prediction and interpretation of distributed neural activity with sparse models. NeuroImage, 44:112–122, 2009.
 [14] M.K. Chung. Statistical Morphometry in Neuroanatomy. Ph.D. Thesis, McGill University, 2001. http://www.stat.wisc.edu/~mchung/papers/thesis.pdf.
 [15] M.K. Chung. Statistical and Computational Methods in Brain Image Analysis. CRC Press, 2013.
 [16] M.K. Chung, N. Adluru, K.M. Dalton, A.L. Alexander, and R.J. Davidson. Scalable brain network construction on white matter fibers. In Proc. of SPIE, volume 7962, page 79624G, 2011.
 [17] M.K. Chung, J.L. Hanson, H. Lee, Nagesh Adluru, Andrew L. Alexander, R.J. Davidson, and S.D. Pollak. Persistent homological sparse network approach to detecting white matter abnormality in maltreated children: MRI and DTI multimodal study. MICCAI, Lecture Notes in Computer Science (LNCS), 8149:300–307, 2013.
 [18] M.K. Chung, V. Singh, P.T. Kim, K.M. Dalton, and R.J. Davidson. Topological Characterization of Signal in Brain Images Using MinMax Diagrams. MICCAI, Lecture Notes in Computer Science (LNCS), 5762:158–166, 2009.
 [19] M.K. Chung, K.J. Worsley, M.N. Brendon, K.M. Dalton, and R.J. Davidson. General Multivariate Linear Modeling of Surface Shapes Using SurfStat. NeuroImage, 53:491–505, 2010.
 [20] M.K. Chung, K.J. Worsley, T. Paus, D.L. Cherif, C. Collins, J. Giedd, J.L. Rapoport, and A.C. Evans. A unified statistical approach to deformationbased morphometry. NeuroImage, 14:595–606, 2001.
 [21] P.A. Cook, Y. Bai, S. NedjatiGilani, K.K. Seunarine, M.G. Hall, G.J. Parker, and D.C. Alexander. Camino: Opensource diffusionMRI reconstruction and processing. In 14th Scientific Meeting of the International Society for Magnetic Resonance in Medicine, 2006.
 [22] C. Davatzikos, M. Vaillant, S.M. Resnick, J.L. Prince, S. Letovsky, and N Bryan. A computerized approach for morphological analysis of the corpus callosum. Journal of Computer Assisted Tomography, 20:88–97, 1996.
 [23] V. de Silva and R. Ghrist. Homological sensor networks. Notic Amer Math Soc, 54:10–17, 2007.
 [24] B. Deen and K. Pelphrey. Perspective: brain scans need a rethink. Nature, 491:S20–S20, 2012.
 [25] M.L. Dequeant, S. Ahnert, H. Edelsbrunner, T M.A. Fink, E. F. Glynn, G. Hattem, A. Kudlicki, Y. Mileyko, J. Morton, A. R. Mushegian, L. Pachter, M. Rowicka, A. Shiu, B. Sturmfels, and O. Pourquie. Comparison of pattern detection methods in microarray time series of the segmentation clock. PLoS One, 3:e2856, 2008.
 [26] A. Dubb, R. Gur, B. Avants, and J. Gee. Characterization of sexual dimorphism in the human corpus callosum. Neuroimage, 20:512–519, 2003.
 [27] H. Edelsbrunner and J. Harer. Persistent homology  a survey. Contemporary Mathematics, 453:257–282, 2008.
 [28] H. Edelsbrunner, D. Letscher, and A. Zomorodian. Topological persistence and simplification. Discrete and Computational Geometry, 28:511–533, 2002.
 [29] B. Efron. The jackknife, the bootstrap and other resampling plans, volume 38. SIAM, 1982.
 [30] A. Fornito, A. Zalesky, and E.T. Bullmore. Network scaling effects in graph analytic studies of human restingstate fMRI data. Frontiers in Systems Neuroscience, 4:1–16, 2010.
 [31] J. Friedman, T. Hastie, and R. Tibshirani. Sparse inverse covariance estimation with the graphical lasso. Biostatistics, 9:432, 2008.
 [32] K.J. Friston, A.P. Holmes, K.J. Worsley, J.P. Poline, C.D. Frith, and R.S.J. Frackowiak. Statistical parametric maps in functional imaging: a general linear approach. Human Brain Mapping, 2:189–210, 1995.
 [33] R. Ghrist. Barcodes: The persistent topology of data. Bulletin of the American Mathematical Society, 45:61–75, 2008.
 [34] J. D. Gibbons and S. Chakraborti. Nonparametric Statistical Inference. Chapman Hall/CRC Press, 2011.
 [35] G. Gong, Y. He, L. Concha, C. Lebel, D.W. Gross, A.C. Evans, and C. Beaulieu. Mapping anatomical connectivity patterns of human cerebral cortex using in vivo diffusion tensor imaging tractography. Cerebral Cortex, 19:524–536, 2009.
 [36] J.L. Hanson, N. Adluru, M.K. Chung, A.L. Alexander, R.J. Davidson, and S.D. Pollak. Early neglect is associated with alterations in white matter integrity and cognitive functioning. Child Development, 84:1566–1578, 2013.
 [37] J.L. Hanson, M.K. Chung, B.B. Avants, K.D. Rudolph, E.A. Shirtcliff, J.C. Gee, R.J. Davidson, and S.D. Pollak. Structural variations in prefrontal cortex mediate the relationship between early childhood stress and spatial working memory. The Journal of Neuroscience, 32:7917–7925, 2012.
 [38] Y. He, Z. Chen, and A. Evans. Structural insights into aberrant topological patterns of largescale cortical networks in Alzheimer’s disease. Journal of Neuroscience, 28:4756, 2008.
 [39] Y. He, Z.J. Chen, and A.C. Evans. Smallworld anatomical networks in the human brain revealed by cortical thickness from MRI. Cerebral Cortex, 17:2407–2419, 2007.
 [40] D. Horak, S. Maletić, and M. Rajković. Persistent homology of complex networks. Journal of Statistical Mechanics: Theory and Experiment, 2009:P03034, 2009.
 [41] S. Huang, J. Li, L. Sun, J. Liu, T. Wu, K. Chen, A. Fleisher, E. Reiman, and J. Ye. Learning brain connectivity of Alzheimer’s disease from neuroimaging data. In Advances in Neural Information Processing Systems, pages 808–816, 2009.
 [42] S. Huang, J. Li, L. Sun, J. Ye, A. Fleisher, T. Wu, K. Chen, and E. Reiman. Learning brain connectivity of Alzheimer’s disease by sparse inverse covariance estimation. NeuroImage, 50:935–949, 2010.
 [43] A.P. Jackowski, C.M. de Araújo, A.L.T. de Lacerda, de Jesus M., and J. Kaufman. Neurostructural imaging findings in children with posttraumatic stress disorder: Brief review. Psychiatry and Clinical Neurosciences, 63:1–8, 2009.
 [44] P. Jezzard and S. Clare. Sources of distortion in functional MRI data. Human Brain Mapping, 8:80–85, 1999.
 [45] S. C. Joshi, B. Davis, M. Jomier, and G. Gerig. Unbiased diffeomorphic atlas construction for computational anatomy. NeuroImage, 23:151–160, 2004.
 [46] W.H. Kim, D. Pachauri, C. Hatt, M.K. Chung, S. Johnson, and V. Singh. Wavelet based multiscale shape features on arbitrary surfaces for cortical thickness discrimination. In Advances in Neural Information Processing Systems, pages 1250–1258, 2012.
 [47] H. Lee, M.K. Chung, H. Kang, B.N. Kim, and D.S. Lee. Computing the shape of brain networks using graph filtration and GromovHausdorff metric. MICCAI, Lecture Notes in Computer Science, 6892:302–309, 2011.
 [48] H. Lee, H. Kang, M.K. Chung, B.N. Kim, and D.S Lee. Persistent brain network homology from the perspective of dendrogram. IEEE Transactions on Medical Imaging, 31:2267–2277, 2012.
 [49] H. Lee, D.S.. Lee, H. Kang, B.N. Kim, and M.K. Chung. Sparse brain network recovery under compressed sensing. IEEE Transactions on Medical Imaging, 30:1154–1165, 2011.
 [50] R.K. Lenroot and J.N. Giedd. Brain development in children and adolescents: insights from anatomical magnetic resonance imaging. Neuroscience & Biobehavioral Reviews, 30:718–729, 2006.
 [51] J.P. Lerch, K. Worsley, W.P. Shaw, D.K. Greenstein, R.K. Lenroot, J. Giedd, and A.C. Evans. Mapping anatomical correlations across cerebral cortex (MACACC) using cortical thickness from MRI. NeuroImage, 31:993–1003, 2006.
 [52] Y. Li, Y. Liu, J. Li, W. Qin, K. Li, C. Yu, and T. Jiang. Brain Anatomical Network and Intelligence. PLoS Computational Biology, 5(5):e1000395, 2009.
 [53] M.M. Loman, A.E. Johnson, A. Westerlund, S.D. Pollak, C.A. Nelson, and M.R. Gunnar. The effect of early deprivation on executive attention in middle childhood. Journal of Child Psychiatry and Psychology, pages 224–236, 2010.
 [54] W.E. Lorensen and H.E. Cline. Marching cubes: A high resolution 3D surface construction algorithm. In Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques, pages 163–169, 1987.
 [55] A.M.C. Machado and J.C. Gee. Atlas warping for brain morphometry. In SPIE Medical Imaging, Image Processing, pages 642–651, 1998.
 [56] R. Mazumder and T. Hastie. Exact covariance thresholding into connected components for largescale graphical LASSO. The Journal of Machine Learning Research, 13:781–794, 2012.
 [57] AR Mclntosh and F. GonzalezLima. Structural equation modeling and its application to network analysis in functional brain imaging. Human Brain Mapping, 2:2–22, 1994.
 [58] M.E.J. Newman and D.J. Watts. Scaling and percolation in the smallworld network model. Physical Review E, 60:7332–7342, 1999.
 [59] J. Peng, P. Wang, N. Zhou, and J. Zhu. Partial correlation estimation by joint sparse regression models. Journal of the American Statistical Association, 104:735–746, 2009.
 [60] S.D. Pollak. Mechanisms linking early experience and the emergence of emotions: Illustrations from the study of maltreated children. Current Directions in Psychological Science, 17:17, 370–375, 2008.
 [61] S.D. Pollak, C.A. Nelson, M.F. Schlaak, B.J. Roeber, S.S. Wewerka, K.L. Wiik, K.A. Frenn, M.M. Loman, and M.R. Gunnar. Neurodevelopmental effects of early deprivation in postinstitutionalized children. Child Development, 81:224–236, 2010.
 [62] A. Pothen and C.J. Fan. Computing the block triangular form of a sparse matrix. ACM Transactions on Mathematical Software (TOMS), 16:324, 1990.
 [63] J.D. Power, K.A. Barnes, A.Z. Snyder, B.L. Schlaggar, and S.E. Petersen. Spurious but systematic correlations in functional connectivity MRI networks arise from subject motion. Neuroimage, 59:2142–2154, 2012.
 [64] A. Rao, P. Aljabar, and D. Rueckert. Hierarchical statistical shape analysis and prediction of subcortical brain structures. Medical Image Analysis, 12:55–68, 2008.
 [65] M. Rubinov and O. Sporns. Complex network measures of brain connectivity: Uses and interpretations. NeuroImage, 52:1059–1069, 2010.
 [66] A. Sacan, O. Ozturk, H. Ferhatosmanoglu, and Y. Wang. Lfmpro: A tool for detecting significant local structural sites in proteins. Bioinformatics, 6:709–716, 2007.
 [67] M.M. Sanchez and S.D. Pollak. Socioemotional development following early abuse and neglect: Challenges and insights from translational research. Handbook of Developmental Social Neuroscience, 17:497–520, 2009.
 [68] T.D. Satterthwaite, D.H. Wolf, J. Loughead, K. Ruparel, M.A. Elliott, H. Hakonarson, R.C. Gur, and R.E. Gur. Impact of inscanner head motion on multiple measures of functional connectivity: relevance for studies of neurodevelopment in youth. Neuroimage, 60:623–632, 2012.
 [69] J. Schäfer and K. Strimmer. A shrinkage approach to largescale covariance matrix estimation and implications for functional genomics. Statistical Applications in Genetics and Molecular Biology, 4:32, 2005.
 [70] G. Singh, F. Memoli, T. Ishkhanov, G. Sapiro, G. Carlsson, and D.L. Ringach. Topological analysis of population activity in visual cortex. Journal of Vision, 8:1–18, 2008.
 [71] C. Song, S. Havlin, and H.A. Makse. Selfsimilarity of complex networks. Nature, 433:392–395, 2005.
 [72] K. Supekar, V. Menon, D. Rubin, M. Musen, and M.D. Greicius. Network analysis of intrinsic functional brain connectivity in Alzheimer’s disease. PLoS Computational Biology, 4(6):e1000100, 2008.
 [73] P.M. Thompson, J.N. Giedd, R.P. Woods, D. MacDonald, A.C. Evans, and A.W Toga. Growth patterns in the developing human brain detected using continuummechanical tensor mapping. Nature, 404:190–193, 2000.
 [74] P.M. Thompson and A.W. Toga. Anatomically driven strategies for highdimensional brain image warping and pathology detection. Brain Warping, pages 311–336, 1999.
 [75] R. Tibshirani. Regression shrinkage and selection via the LASSO. Journal of the Royal Statistical Society. Series B (Methodological), 58:267–288, 1996.
 [76] P.A. ValdésSosa, J.M. SánchezBornot, A. LageCastellanos, M. VegaHernández, J. BoschBayard, L. MelieGarcía, and E. CanalesRodríguez. Estimating brain functional connectivity with sparse multivariate autoregression. Philosophical Transactions of the Royal Society B: Biological Sciences, 360:969–981, 2005.
 [77] K.R.A. Van Dijk, M.R. Sabuncu, and R.L. Buckner. The influence of head motion on intrinsic functional connectivity MRI. Neuroimage, 59:431–438, 2012.
 [78] Daniela M Witten, Jerome H Friedman, and Noah Simon. New insights and faster computations for the graphical lasso. Journal of Computational and Graphical Statistics, 20:892–900, 2011.

[79]
K.J. Worsley, A. Charil, J. Lerch, and A.C. Evans.
Connectivity of anatomical and functional MRI data.
In
Proceedings of IEEE International Joint Conference on Neural Networks (IJCNN)
, volume 3, pages 1534–1541, 2005. 
[80]
K.J. Worsley, J.I. Chen, J. Lerch, and A.C. Evans.
Comparing functional connectivity via thresholding correlations and singular value decomposition.
Philosophical Transactions of the Royal Society B: Biological Sciences, 360:913, 2005.  [81] A. Zalesky, A. Fornito, I.H. Harding, L. Cocchi, M. Yücel, C. Pantelis, and E.T. Bullmore. Wholebrain anatomical networks: Does the choice of nodes matter? NeuroImage, 50:970–983, 2010.
 [82] H. Zhang, O. van Kaick, and R. Dyer. Spectral methods for mesh processing and analysis. In EUROGRAPHICS, pages 1–22, 2007.
Comments
There are no comments yet.