Metrics for comparing neuronal tree shapes based on persistent homology

As more and more neuroanatomical data are made available through efforts such as NeuroMorpho.Org and FlyCircuit.org, the need to develop computational tools to facilitate automatic knowledge discovery from such large datasets becomes more urgent. One fundamental question is how best to compare neuron structures, for instance to organize and classify large collection of neurons. We aim to develop a flexible yet powerful framework to support comparison and classification of large collection of neuron structures efficiently. Specifically we propose to use a topological persistence-based feature vectorization framework. Existing methods to vectorize a neuron (i.e, convert a neuron to a feature vector so as to support efficient comparison and/or searching) typically rely on statistics or summaries of morphometric information, such as the average or maximum local torque angle or partition asymmetry. These simple summaries have limited power in encoding global tree structures. Based on the concept of topological persistence recently developed in the field of computational topology, we vectorize each neuron structure into a simple yet informative summary. In particular, each type of information of interest can be represented as a descriptor function defined on the neuron tree, which is then mapped to a simple persistence-signature. Our framework can encode both local and global tree structure, as well as other information of interest (electrophysiological or dynamical measures), by considering multiple descriptor functions on the neuron. The resulting persistence-based signature is potentially more informative than simple statistical summaries (such as average/mean/max) of morphometric quantities—Indeed, we show that using a certain descriptor function will give a persistence-based signature containing strictly more information than the classical Sholl analysis. At the same time, our framework retains the efficiency associated with treating neurons as points in a simple Euclidean feature space, which would be important for constructing efficient searching or indexing structures over them. We present preliminary experimental results to demonstrate the effectiveness of our persistence-based neuronal feature vectorization framework.

where u − v ∞ = max{|u.x − v.x|, |u.y − v.y|} denotes the L ∞ distance between two points; and γ ranges over all bijections between D 1 and D 2 . It is known that the bottleneck distance d B (D 1 , D 2 ) can be computed in O(m 1.5 log m) time, where m is the total number of points in D 1 ∪ D 2 .
Given two functions f, g : |T | → R, suppose that g is a perturbation of f with bounded distance in L ∞ -norm, that is f − g ∞ := max x∈|T | |f (x) − g(x)| measures the amount of perturbation of g from f . The Stability Theorem [1] states that for a function f and its perturbation g, the bottleneck distance between their persistence diagram summaries is bounded from above by the size of the perturbation; that is, d B (Dgf, Dg g) ≤ f − g ∞ . This result is later generalized to more general persistence modules, and show that the bottleneck distance between two persistent diagrams is bounded by the so-called interleaving distance between the corresponding persistence modules that generate them [2,3].
In our setting, given two neuron trees T 1 and T 2 with descriptor functions f 1 : |T 1 | → R and f 2 : |T 2 | → R, we cannot directly compare these two descriptor functions since they are defined on different domains (T 1 and T 2 , respectively). We instead use the so-called functional distortion distance [4] d F D (f 1 , f 2 ) to measure how different the functions f 1 and f 2 are. Intuitively, d F D considers all pairs of mappings between T 1 and T 2 as a way to align them, say φ : |T 1 | → |T 2 | and ψ : |T 2 | → |T 1 | to align T 1 to T 2 (via φ) as well as align T 2 to T 1 (via ψ). It then compares f 1 and f 2 composited with these maps so that they are then defined on a common domain. Each such pair of maps (alignment) (φ, ψ) will give a cost, measuring how well f 1 and f 2 are aligned under these two maps, and d F D (f 1 , f 2 ) returns the minimum cost under all possible such alignments (pairs of maps). We refer the readers to [4] see the formal definition. It follows from results of [4] that Figure 1: T 1 is a noisy version of T 1 : there could be local combinatorial changes such as the subtrees B, C and D merge at slightly different height, and there could also be spurious noisy branches. However, such changes do not perturb the tree metric much: in this specific example, the metric distortion is bounded by ε. As a result, their corresponding persistence diagram summaries are also close with d B (D, D ) ≤ ε.
This stability result applies to any descriptor functions. For example, suppose we consider the Euclidean distance functions f 1 and f 2 on the two trees in Figure 1. Then their persistence diagram summaries are close (at most ε), despite that there are noisy branches, as well as combinatorial changes in the tree structures from T 1 to T 2 -Indeed, it is not hard to establish maps φ : T 1 → T 2 and ψ : T 2 → T 1 and show that the cost of the Euclidean distance function incurred by them is at most ε, thus upper bounds d F D (T 1 , T 2 ) by ε too.
As another example, if we use the geodesic distance to the root as the descriptor function, then by using results from [5], the bottleneck distance between resulting persistence diagrams is stable w.r.t. changes in the input neuron trees as measured by the Gromov-Hausdorff distance between these trees. The Gromov-Hausdorff distance is popular way to measure the level of near-isometry between two metric spaces [6,7]. The Gromove-Hausdorff distance between the two neuron trees in Figure 1 is at most ε, implying that using the geodesic distance as descriptor functions f : T 1 → IR and f : T 1 → IR, we have d B (Dgf, Dgf ) ≤ ε as well.
In our framework, to improve computational efficiency, we vectorize the persistent diagrams describe in (Step 2), and the natural L p -distance between them are sum-based, instead of maxbased (as in bottleneck distance). One can extend the bottleneck distance to the so-called degree p-Wasserstein distance d W,p (D 1 , D 2 ) between two persistence diagrams D 1 and D 2 which we will introduce shortly in Eqn (2). The stability for the Wasserstein distance of persistence diagrams is not as well understood as in the bottleneck distance case (which is in fact the case when p = ∞), although there are some results for some special cases [8].
Stability of the persistence feature vectors. We now discuss the stability of feature vectorization step. Specifically, suppose we are given two persistence diagrams D 1 and D 2 , with corresponding feature functions ρ 1 = ρ D 1 , ρ 2 = ρ D 2 : IR → IR induced from D 1 and D 2 as described in Section 2.2 of the main text.
First, the degree-p Wasserstein distance between D 1 and D 2 , for 1 ≤ p < ∞, is defined as: where D 1 := D 1 ∪ L and D 2 := D 2 ∪ L are the persistence diagrams augmented with points in the diagonal L as before.
Theorem 0.1 The L 1 -distance between feature functions ρ 1 and ρ 2 is stable w.r.to the 1-Wasserstein distance between the diagrams D 1 and D 2 generating them. Let ∆ denote the largest persistence value of any point in D 1 ; that is, ∆ = max u,u ∈D 1 |u − u|. Specifically, We now prove this theorem. First, we need the following result bounding the distance between two 1-dimensional Gaussians [9].