NeuroPath2Path: Classification and elastic morphing between neuronal arbors using path-wise similarity

TL;DR

NeuroPath2Path introduces a novel graph-theoretic method for neuronal morphology analysis, enabling classification and morphing of neurons based on path-wise similarity, improving over existing techniques.

Contribution

The paper presents a new approach using rooted path decomposition and elastic deformation for neuronal shape comparison and classification, with enhanced accuracy and visualization.

Findings

Effective neuron classification using path-wise similarity.

Superior performance over existing methods.

Provides a visualization tool for morphological morphing.

Abstract

The shape and connectivity of a neuron determine its function. Modern imaging methods have proven successful at extracting such information. However, in order to analyze this type of data, neuronal morphology needs to be encoded in a graph-theoretic method. This encoding enables the use of high throughput informatic methods to extract and infer brain function. The application of graph-theoretic methods to neuronal morphological representation comes with certain difficulties. Here we report a novel, effective method to accomplish this task. The morphology of a neuron, which consists of its overall size, global shape, local branch patterns, and cell-specific biophysical properties, can vary significantly with the cell's identity, location, as well as developmental and physiological state. Various algorithms have been developed to customize shape based statistical and graph related features for quantitative analysis of neuromorphology, followed by the classification of neuron cell types using the features. Unlike the classical feature extraction based methods from imaged or 3D reconstructed neurons, we propose a model based on the rooted path decomposition from the soma to the dendrites of a neuron and extract morphological features on each path. We hypothesize that measuring the distance between two neurons can be realized by minimizing the cost of continuously morphing the set of all rooted paths of one neuron to another. To validate this claim, we first establish the correspondence of paths between two neurons using a modified Munkres algorithm. Next, an elastic deformation framework that employs the square root velocity function is established to perform the continuous morphing, which, in addition, provides an effective visualization tool. We experimentally show the efficacy of NeuroPath2Path, NeuroP2P, over the state of the art.

Tables

Table 1. Table 1: Importance weight δ 𝛿 \delta values for dataset-1. For space constraint, we provide feature-specific importance weight for classification in case of pairwise classes and all classes separately.

	Tortuo$\left(\kappa \right)$	Bifur-angle$\left(b\right)$	Part-aym$\left(\alpha \right)$	Concur$\left(C\right)$	Seg-len$\left(\beta \right)$	Diverg$\left(\lambda \right)$
Motor1-Motor2	0.0729	0.1072	0.0876	0.2644	0.3794	0.0955
Motor1-Prelimbic	0.0820	0.0737	0.0631	0.2262	0.4483	0.1067
Motor1-Somato	0.0706	0.0929	0.0646	0.2328	0.4175	0.1218
Motor1-Visual1	0.0689	0.1230	0.0659	0.3317	0.2878	0.1227
Motor1-Visual2	0.0499	0.1212	0.0812	0.2612	0.3952	0.0912
Motor2-Prelimbic	0.2719	0.0737	0.0703	0.0899	0.5367	0.0484
Motor2-Somato	0.0452	0.0256	0.0328	0.4008	0.3730	0.1266
Motor2-Visual1	0.0167	0.0773	0.1041	0.2996	0.2969	0.2055
Motor2-Visual2	0.0545	0.0411	0.0564	0.2235	0.5489	0.0756
Prelimbic-Somato	0.1314	0.1168	0.1121	0.0568	0.4621	0.1187
Prelimbic-Visual1	0.2091	0.0821	0.0617	0.1012	0.4615	0.0795
Prelimbic-Visual2	0.2311	0.0551	0.0453	0.0790	0.3556	0.2339
Somato-Visual1	0.0693	0.0598	0.0584	0.1319	0.5816	0.0941
Somato-Visual2	0.0144	0.0310	0.2043	0.6963	0.0159	0.0382
Visual1-Visual2	0.0075	0.0880	0.1482	0.3551	0.3558	0.0473
Overall	0.0785	0.0807	0.1736	0.1956	0.2600	0.2076

Table 2. Table 2: Distance and correspondence between paths. The correspondences between the paths of neuron-1 and neuron-2 are enlisted in the first two columns. The numbers in yellow indicate that the correspondence obtained by NeuroP2P matches with the candidates of best correspondence in the sense of minimum distance. The pairs in blue are subjected for further verification because of large differences in the hierarchy values (Routine 4 4 4 in Algorithm 2 2 2 ).

Index (Neuron-1)	Index (Neuron-2)	Distance	ElasticP2P batabyal2018elasticpath2path	Best match
1	25	1.0595	1	9
2	2	0.8031	2	5
3	27	0.4530	3	5
4	12	0.7910	4	4
5	3	0.5571	5	5
6	26	0.6165	6	5
7	28	0.5690	7	5
8	5	0.5606	8	8
9	23	0.6271	9	9
10	13	0.7904	10	4
11	24	0.6096	11	5
1	7	1.4997	NA	8
2	1	0.9902	NA	9
3	15	1.1052	NA	5
4	11	0.9223	NA	9
5	6	0.5932	NA	5
6	17	1.7108	NA	9
7	18	1.2869	NA	9
8	20	0.8704	NA	9
9	4	0.7708	NA	5
10	14	0.9186	NA	5
11	10	1.1796	NA	9
9	8	0.8504	NA	5
4	9	1.2279	NA	5
10	16	1.1810	NA	5
8	19	1.1658	NA	9
2	21	1.1522	NA	5
5	22	0.8931	NA	5

Table 3. Table 3: Importance weight δ 𝛿 \delta values for dataset-2

	Tortuo$\left(\kappa \right)$	Bifur-angle$\left(b\right)$	Part-aym$\left(\alpha \right)$	Concur$\left(C\right)$	Seg-len$\left(\beta \right)$	Diverg$\left(\lambda \right)$
Ganglion-Granule	0.0278	0.0138	0.0241	0.2102	0.6727	0.0533
Ganglion-Motor	0.0464	0.0797	0.0918	0.1865	0.4730	0.1226
Ganglion-Purkinje	0.27	0.032	0.1437	0.1905	0.2270	0.1368
Ganglion-Pyramidal	0.0433	0.0001	0.0553	0.5121	0.3156	0.0809
Granule-Motor	0.0356	0.0491	0.0372	0.1449	0.6697	0.0636
Granule-Purkinje	0.0147	0.1218	0.1875	0.1073	0.5285	0.0402
Granule-Pyramidal	0.0453	0.0305	0.0372	0.2231	0.5390	0.1248
Motor-Purkinje	0.0094	0.4046	0.0556	0.0007	0.5167	0.0130
Motor-Pyramidal	0.0223	0.0737	0.0465	0.1275	0.6650	0.0650
Purkinje-Pyramidal	0.0219	0.2201	0.1088	0.1528	0.4223	0.0740
Overall	0.1304	0.0372	0.0417	0.1786	0.4489	0.1633

Equations

$${\upsilon }^f$$

$${\upsilon }^g$$

$$d({\upsilon }^{fg})=\sqrt{\frac{1}{{\varphi }_2^b}\sum _{k=1}^{{\varphi }_2^b}{w}_k({\upsilon }_k^f-{\upsilon }_k^g{)}^2}$$

$${\mu }^{fg}$$

$${\chi }_{G_1{G}_2}=\sum _{k=1}^{n_2}{\mu }^{{\gamma }_{1k}{\gamma }_{2k}}$$

$${\delta }^{opt}$$

$$L(\delta )$$

$$\frac{dL}{d\delta }$$

Full text

The morphology of a neuron, which consists of its overall size, global shape, local branch patterns, and cell-specific biophysical properties, can vary significantly with the cell’s identity, location, as well as developmental and physiological state. Various algorithms have been developed to customize shape based statistical and graph related features for quantitative analysis of neuromorphology, followed by the classification of neuron cell types using the features. Unlike the classical feature extraction based methods from imaged or 3D reconstructed neurons, we propose a model based on the rooted path decomposition from the soma to the dendrites of a neuron and extract morphological features on each path. We hypothesize that measuring the distance between two neurons can be realized by minimizing the cost of continuously morphing the set of all rooted paths of one neuron to another. To validate this claim, we first establish the correspondence of paths between two neurons using a modified Munkres algorithm. Next, an elastic deformation framework that employs the square root velocity function is established to perform the continuous morphing, which, in addition, provides an effective visualization tool. We experimentally show the efficacy of NeuroPath2Path, NeuroP2P, over the state of the art.

1 Introduction

Neurons process information by transmitting electrical signals via complex circuitry. The functionality of each neuron depends on a set of intrinsic factors, such as morphology, ionic channel density, gene expression, including the extrinsic ones, such as connectivity to other neurons brown2008quantifying ; ascoli2008petilla . In 1899, Cajal y1972histologie , considered the founder of modern neuroscience, put forward his pioneering work on neuroanatomy with detailed, accurate, and meticulous illustrations, and posited that the shape of a neuron determines its functionality. Experimental results strongly support this idea. Inspired by this fundamental work, the study of neuromorphology primarily aims at analyzing and quantifying the complex shape and physiology of neurons in specific functional regions to identify relationships.

Neurons vary significantly in size, shape, and length. A major obstacle towards understanding the brain is the development of efficient ways to encode these shapes. The anatomical and geometrical features of neurons of any cell-type, for example, pyramidal cells differ based on the regions in which the cells reside bielza2014branching ; brown2008quantifying . Fig. 1 shows regional variation in the structure and geometry of dendritic arbors of pyramidal cells romand2011morphological . It is observed that the number of branches, length, surface area, and volume of apical dendrites is $4-9$ times larger for hippocampal than for cortical regions, whereas in terms of the same features of basal dendritic arbors, it is approximately $3$ times brown2008quantifying . Another source of variation stems from technical imprecision in measurements obtained while performing 3D reconstruction from image stacks using software tracing tools, such as Neurolucida glaser1990neuron . Noise due to technical imprecision includes wide variations in the number of manually or semi-automatically traced 3D locations (approximately between $60$ to $70,000$), the number of ramified branches and bifurcations by different tracers, and deletion of dendritic spines adversely affect the registration of neurons, and thereby induce error in morphological feature quantification. The skeletons of dendritic and axonal branches form a tree topology with a number of bifurcations. The bifurcations at successive stages help in a series of effective and unambiguous signal processing modules, such as active and passive signal propagation, integration, filter, attenuation, oscillation, and backpropagation ascoli2008petilla ; london2005dendritic .

From the soma to the dendritic terminals of a neuron, the diameter of the dendritic shaft tapers jan2010branching ; wen2008cost . The increased diameter of a dendritic shaft near the soma is tailored to faster signal propagation to the soma compared to the dendritic tuft , which helps generate action potential in the soma. Several research works consider the branches in the proximity of soma are more important compared to the distant dendritic tuft and spines in the analysis of neuromorphology kanari2018topological ; cervantes2018morphological ; lopez2011models . The length of dendritic branch segments shows similar behavior when propagating away from soma. For instance, the terminal segments are longer than the intermediate branch segments for basal dendrites in cortical pyramidal cells bielza2014branching . These observations support the Bayesian philosophy which is geared towards the analysis of morphogenesis of neurons lopez2011models . Functions such as synaptic boosting migliore2002emerging , coadaptive local spiking gasparini2004initiation , and global spike amplification williams2004spatial suggest the use of other morphometrics to describe the structural aspects on the functions. For example, packing density of ramified branches and bifurcations of neuron potentially trigger intermittently co-adpative spiking.

The tree-type arbors of neurons and the availability of the inventory of digitally-traced 3D reconstructed neurons, Neuromorpho ascoli2007neuromorpho , provided significant momentum in the last decade for the quantitative and qualitative assessment of neuroanatomy via graph-based morphometrics. In Neuromorpho, the sequentially-aligned slices of microscopic images are registered and traced using software meijering2010neuron , such as Neurolucida glaser1990neuron and Neuromantic myatt2012neuromantic , and the reconstructed images can then be processed through software, such as L-measure scorcioni2008measure to extract an extensive list of morphological metrics. On one hand, there are several research works dedicated to analyze the neuromorphology of specific cell types, such as basal dendrites of cortical pyramidal cells lopez2011models ; bielza2014branching , GABAergic interneuron cells ascoli2008petilla and others. These works account for region-specific variations in the physiology and anatomy of a neuron cell to establish the effect of certain functions on the structure. On the other hand, research efforts, such as blastneuron wan2015blastneuron , neurosol batabyal2017neurosol ,and TMD kanari2018topological , attempt to extract model based features, which are catered to the need for automated classification of different neuronal cells. The motivation behind this avenue of research is that it is impossible to identify and categorize one trillion neuronal cells by adopting manual or even semi-automatic methods.

State-of-the-art methods wan2015blastneuron ; batabyal2018neurobfd for the classification of neurons can be broadly divided in two categories. Research in the $first$ category, which is supervised in nature, employs different feature extraction algorithms followed by suitable classifiers to obtain classification accuracy in percentage. The validation of the methods are performed by adopting a series of statistical tests. However, the significant variation in the neuron skeletons precludes the selection of the optimal set of morphometrics as features. Adoption of feature transformations, such as principal component analysis (PCA) or kernel transformations, may improve the classification accuracy. Nevertheless, these transformations obscures the identity of discriminating features as the transformed space is formed by linear or nonlinear composition of extracted features. In addition, the classification accuracy of categorization does not quite explain the physiological and structural differences between two neurons.

The $second$ category follows mostly unsupervised approaches and attempts to compute pairwise distances between neurons. Authors in sarkar2013shape used Fourier based shape descriptors to encode the global shape of a neuron, which however ignored the local features of the neuron arbor. Gillette gillette2015topological ; gillette2015topological1 performed a sequence alignment based algorithm for categorization by decomposing a neuron into a sequence of branches. The approach failed to consider geometric features. Blastneuron wan2015blastneuron adopted a mixed strategy. Using a supervised approach, the method first extracted $21$ global morphological features and $13$ moment invariant features to retrieve a set of targets that closely matches to each test neuron in terms of the anatomical structures. Each target is then RANSAC sampled fischler1981random and aligned optimally to the test neuron, which outputs a distance value. This unsupervised routine decides the output category of the test neuron based on the minimum distance criteria. The method involved initial pruning of branches and resampling of each neuron, which collectively alters the morphology statistics. Moreover, the retrieval accuracy of $233$ projected neurons (PN) of Drosophila drops significantly to $39$% as the number of potential candidates that are to be compared with the target increases. NeuroSoL batabyal2017neurosol offered a graph-theoretic method which is free from registration and resampling. In spite of its appeal of using graph theory, the matrix alignment routine is NP-hard in nature, thereby producing suboptimal results. The problem of comparing a pair of neuron topologies can also be regarded as a graph kernel based similarity measure problem vishwanathan2010graph . However, the rationale behind conventional graph kernels, such as the random walk kernel may be inconsistent with the morphological understanding of a neuron.

Instead of modeling a neuron as a generic graph, the neuron can be modeled as a specialized graph that contains a collection of rooted paths, where each path starts from the soma, called the root node, and ends up in a dendritic terminal. It is important to note that each path acts like an atomic neuron, as it contains the soma and a dendritic end to complete a circuit. Most of the synapses along a path will be nearer the soma than at the end of the path. It is convenient to think problems, such as synaptic plasticity as the evolution of a set of synapses over time along all the paths. During this evolution, there are birth, death and rearrangement of paths. Following the same logic, quantifying the problem of distinguishing two neurons can be equivalently mapped as finding the cost of evolving a set of circuits optimally from one neuron to the other.

Another relevant fact is that path based models basu2011path2path ; batabyal2018elasticpath2path ; kanari2018topological integrate both global (overall shape based approach) and local (vertex or sampled location based approach) features of neuron topology. Topological morphology descriptor (TMD) kanari2018topological aimed at solving the categorization problem, encoded the birth and death of path segments over time in a persistence diagram used as a barcode. The authors showed that TMD exhibits robustness to erroneous 3D sampling and ambiguous branching when the neuron is reconstructed using two different tracing tools. However, the conversion of a discrete 3D reconstructed neuron to the persistence image space is irreversible and many-to-one. Based on the distance used to mark and quantify the birth and death of a branch or component of the neuronal tree, a single persistence image may correspond to multiple neurons. In addition, an appropriate distance measure between persistence diagrams is still unavailable. The work in Path2Path basu2011path2path shows potential to address the neuron cell categorization problem and can be extended to several other related problems, such as synaptogenesis, degeneration of neurons due to neurological diseases, and synaptic plasticity which can be studied by inspecting the path statistics. The work described in this article is motivated by the framework of Path2Path.

1.1 What is Path2Path and its variants?

The principle of Path2Path is based on finding the optimal correspondence between the paths of one neuron to that of the other using a proposed metric. It is an intuitive circuit-based approach that appeals to its electrical engineering inventors. In Path2Path, each sampled location on a path is endowed with 3D coordinate values and two features, concurrence and hierarchy. The concurrence value at each location denotes the number of paths from the soma to dendritic ends that visit that node. The hierarchy value at a location indicates the depth of the location from the soma in terms of the number of bifurcations between the point and the soma. The hierarchy value of a location counts the number of bifurcations one has to cross while traversing from the soma to that location. Using the 3D coordinates, concurrence, and hierarchy values of each location on a path, authors in Path2Path proposed an empirical metric that outputs a distance value between two paths. A path from a neuron corresponds to a path from another neuron if the distance between the paths is minimum over all the paths of the latter neuron.

This approach has several drawbacks. The Path2Path algorithm is dependent on the number of sampled locations of each path and the registration. The selection of the metric is arbitrary in a sense that the metric is null when two paths have the same set of concurrence values but different locations and hierarchy values. Therefore, it does not qualify the axioms of a metric. In addition, the proposed distance measure uses the Euclidean distance between two paths as a part of the distance computation routine, which favors the pair if they are aligned in proximity after registration. The algorithm of finding the correspondence is not one-to-one and it often leads to the degenerate case where all paths from one neuron are matched to only one path in the second neuron. The problem exacerbates when the number of samples in the two paths are unequal. One potential solution is to resample each path using a constant step wan2015blastneuron , but may, unfortunately, eliminate the importance of the locations, such as curvature of a rooted path prior to resampling.

ElasticPath2Path batabyal2018elasticpath2path attempted to address the previously mentioned problems. It introduced a mid-point based resampling routine as opposed to constant-length resampling. To ensure one-to-one correspondence between a pair of paths from two different neurons, the Munkres algorithm munkres1957algorithms is employed. Most importantly, elasticPath2Path envisaged the problem of distinguishing two neurons as a continuous deformation between the corresponding paths of the neurons. Such homeomorphism is computed by applying the square root velocity function (SRVF) srivastava2011shape to the Euclidean coordinates of each sample on a path. The visual deformation of the corresponding paths has an enormous impact in the validation of the path based on customized features and the proposed distance measure. On the flip side, elasticPath2Path failed to address the problem where there is a significant difference in the number of paths between two neurons. As the correspondence is one-to-one, it asserts that elasticPath2Path performs subgraph matching. Both Path2Path and elasticPath2Path did not consider important anatomical morphometrics, such as bifurcation angles and partition asymmetry.

1.2 Key aspects of NeuroPath2Path

$\boldsymbol{1)}$ The inception of NeuroPath2Path comes from the realization of neuromorphogenesis and the self-similar phenotype of neuronal arbor. Since its birth from the soma, a path of a neuron has an exploratory attribute to collect external resources by the minimal-length-maximal-routing sporns2005human strategy. Due to the parsimonious exploitation of intrinsic resources (ion density, ATP and other electrophysiological items), the path, which fails to procure external resources, retracts. The exploratory attribute of a path can be expressed by the concurrence values at each sample point of a path. More paths imply more exploration. As the path matures, it has a competitive attribute miina2002application ; genet2014incorporating ; lopez2011models with respect to the other paths in its neighborhood in order to form a synapse. To account for competition, we count the number of paths in the proximity of each sampled location on a path and assign the count to that location.

$\boldsymbol{2)}$ The fractal dimension puvskavs2015fractal ; brown2008quantifying of a neuronal arbor is considered one of the key morphometrics because the fan-out branches of a neuron bear self-similarity. In Path2Path and elasticPath2Path, the notion of matching the paths ignores this important feature. We extend the use of Munkres algorithm to perform one-to-one matching in a sequential fashion, which replicates the self-similar behavior.

$\boldsymbol{3)}$ As path features, we consider the bifurcation angle, partition asymmetry, and fragmentation score to each 3D location on a path. It is shown that the distribution of bifurcation angles in the basal dendrites of cortical pyramidal cells follows a Von Mises distribution bielza2014branching . An experimentally observed fact is that the mean bifurcation angle of branches ordered in a reversed fashion is discriminative for pyramidal cells in different cortical regions. However, the mean bifurcation angle of branches in standard order remains similar for the pyramidal cells. We take the standard ordering of branches, instead of the reverse order, to discriminate different neuronal cell types. Partition asymmetry brown2008quantifying ; polavaram2014statistical ; samsonovich2006morphological is another visually-significant morphometric. We use the caulescence measure as defined in brown2008quantifying to account for the tree asymmetry.

$\boldsymbol{4)}$ We provide visualization of the continuous deformation between a pair of neurons and enumerate path similarity statistics to justify the correspondences between the paths. In contrast, conventional methods perform feature customization and extraction, and the classification, in supervised or unsupervised settings, depends on the abstract feature space and the strength of the classifier. In those methods, the mapping between the space of 3D reconstructed neurons and the feature space is irreversible and abstract. Therefore, apart from the statistical quantification and analysis, it is ambiguous whether improved accuracy of the categorization stems from the trained classifier or the discriminating strength of the extracted features or both. In NeuroPath2Path, the classification problem is modeled as a variant of the transport problem. First, the correspondence of paths between a pair of neurons are decided in the feature space. Next, the correspondence is utilized to deform one neuron to the other. The distances computed between the paths and the deformation together justify the validity of the correspondence.

$\boldsymbol{5)}$ With suitable feature selection, NeuroP2P framework can be applied to perform morphological analysis of any cell type with ramified branching arbors, such as microglia and astrocytes. The continuum that is present in the evolution from one cell type to the other can be utilized in the analysis of cell differentiation. As an example, under certain constraints, the strategy of continuous morphing with branch splitting (explained later) can retrieve the intermediate states of a neuron cell while it evolves from a neural progenitor cell to its fully developed state. In short, NeuroP2P can serve as an effective tool for cell-specific informatics, which is not restricted to classification only.

2 Path modeling of a neuron

As mentioned in the introduction, a digitally-traced 3D sampled neuron can be modeled as a graph. Let the graph be represented by \mathcal{G}=\big{(}\mathcal{V},\mathcal{E},\mathcal{W}\big{)}, where $\mathcal{V}$ is the number of 3D locations as vertices and $\mathcal{E}$ is the set of edges connecting the vertices with the corresponding weights $\mathcal{W}$ harary1969graph . $\mathcal{G}$ is said to be simple if it does not contain multiple edges between any two vertices. A graph is called undirected if there is no preferred direction associated to an edge. A sequence of contiguous edges is called a path if no vertex and edge are repeated in that sequence. A path of length $k$ has $k$ number of edges or equivalently $(k+1)$ vertices. A sequence of contiguous edges is called a trail if no edge is repeated. If all the vertices except the start and the end of a trail are distinct, it is called a loop. A simple graph without a loop is termed as a tree. If the degree of each vertex is fixed, tree has the fastest growth by volume, hence smallest curvature lin2010ricci . A graph is said to be single-connected if there exists at least one path between a pair of vertices. In case of a neuron, $\mathcal{G}$ is a simple, undirected, weighted, and single-connected tree.

A path can be considered as an open curve, $f_{i}(t),~{}t\in[0,1]$, as defined in differential geometry. The cardinality of the set of vertices, or ,equivalently, the total number of 3D locations, is given by $|\mathcal{V}|=N$. Here, there are $n$ dendritic terminals, which implies that the total number of paths rooted at the soma is $n$. Let $\Gamma$ be the set containing all the paths $f_{i},~{}i\in\{1,2,..,n\}$, which is a linear subspace of the classical Wiener space. Each path is sampled with the number of samples as $\phi$ with the sampled path denoted by $\tilde{f}_{i}$. We extract $K$ features for each sample on $\tilde{f}_{i}$, which can be compactly given by the feature matrix for $f_{i}$ as $\Theta_{i}\in\mathbb{R}^{\phi\times K}$. Let $\Theta$ be the ordered set containing the feature matrix for all the paths, $\Theta=\{\Theta_{1},\Theta_{2},...,\Theta_{n}\}$, where $\Theta_{i}$ corresponds to the $i^{th}$ path $\tilde{f}_{i}$. The path model of a neuron $\mathcal{G}$ can be mathematically represented as $H=\{\Gamma,\Theta,\mu\}$, where $\mu$ is a measure that we define in the next section. Note that we use the set of paths, $\Gamma$, as an ordered set which has a one-to-one correspondence with the elements in $\Theta$. The standard branch order of a path, $f_{i}$, is defined as the order in which the locations of bifurcation on a path are visited from the root to the end of the path. Similarly, the reverse branch order is defined when the direction of traversal is reversed. For interclass comparison of neurons, we use the standard order. Whereas, for intraclass comparison, we follow the reverse order.

3 Proposed methodology

Our proposed method, which is sequential, scalable and modular, consists of four key stages as depicted in Fig. 4. In the first stage, centrally curated files of 3D-traced neurons in SWC format (or equivalent formats) are read and then preprocessed to extract only the dendritic arbors, including the soma. Several preprocessing modules, such as range-wise calibration, bifurcation location determination, and synaptic tip identification are employed to aid in preparing an assembly of rooted (soma) paths.

In the second stage, a set of features are extracted from each path, forming a feature descriptor of the path, $\Theta_{i}\forall i$. An exhaustive list of the features that are used in our method is provided in Section 3.1, and the systematic quantification of the features is provided in the Appendix. Notice that each path descriptor can be populated with additional structural and geometric features in order to perform fine-grained analysis.

The central aspect of the third stage is finding an appropriate cost function, as illustrated in Section 3.2. The cost function assimilates several anatomical features (such as segment length and bifurcation angle) and physiologically relevant factors (such as the competitive behavior, decaying anatomical importance of a path from the soma to synapse). A rigorous optimization framework is also formulated to find the relative contributions of such factors. In short, this stage delivers a distance measure between a pair of paths to the last stage.

With a distance measure between neuron paths in hand, this measure is augmented in the final stage as elucidated in section 3.3. We theoretically establish the correspondence of paths of a pair of neurons by repeatedly applying the Munkres algorithm. In contrast to the conventional approaches where the distance between neurons inherently accounts for sub-graph matching, we propose a full-tree matching algorithm. The repeated application of the Munkres algorithm reveals the fractal or self-similar nature of a pair of neurons. Equivalently, the following question may be posed: how many identical copies, taken together, of the first neuron can match with the second neuron, assuming the second neuron is much larger than the first one? Once the correspondence is found, neurons are diffeomorphically transformed to each other by morphing corresponding paths. This visual representation aids in justifying the correspondence of paths.

3.1 Feature extraction on a path

We extract a set of discriminating features from each path $f_{i}\in\Gamma$ of $H$, which are bifurcation angle ($b_{i}$), concurrence ($C_{i}$), hierarchy ($\xi_{i}$), divergence ($\lambda_{i}$), segment length ($\beta_{i}$), tortuosity ($\kappa_{i}$), and partition asymmetry ($\alpha_{i}$). Therefore, $\Theta_{i}=[b_{i},C_{i},\xi_{i},\lambda_{i},\beta_{i},\kappa_{i},\alpha_{i}]\in\mathbb{R}^{\phi\times 7}$. Each feature encodes a specific structural property of a neuronal arbor, as described in the Appendix. A schematic of different features along with the systematic quantification is shown in Fig.3.

3.2 Path alignment and path distance measure, $\mu$

Given an unequal number of samples in a pair of paths, finding the appropriate distance between two paths or open curves is challenging. Due to the resampling bias imposed by a given tracer, in general, a path contains erroneous sampled locations which could alter the path statistics. For example, adding an extra leaf vertex changes the concurrence values of all the locations on a path. Unlike conventional approaches that used different resampling procedures, such as mid-point based resampling, RANSAC sampling, and spectral sampling, we use the help of the branch order as mentioned in section 2 for suboptimal alignment.

Consider two neurons, $G_{1}$ and $G_{2}$, with the corresponding path models given as $H_{1}$ and $H_{2}$, respectively. Let $f$ and $g$ be the two paths that are arbitrarily selected from $\Gamma_{1}$ and $\Gamma_{2}$, respectively. Without loss of generality, let us assume that $f$ and $g$ contain $\phi_{1}^{b}$ and $\phi_{2}^{b}$, the number of locations from which the current paths bifurcate. In the case, where $\phi_{1}^{b}<\phi_{2}^{b}$, we append $(\phi_{2}^{b}-\phi_{1}^{b})$ zeros at the end (standard branch order) or at the front (reverse branch order) of a feature vector on $f$.

Experimental evidence  bielza2014branching suggests that the importance of a bifurcation location on a path decays as one travels the path from the soma to the synaptic end. We utilize this relative importance by way of hierarchy values of the bifurcation locations on a path. Let the sequential order of hierarchy values from the root to the terminal on $f$ be $\xi_{f}=[\xi_{1},\xi_{2},...,\xi_{\phi^{b}_{1}}]$. Using $\xi_{f}$, the $k^{th}$ importance weight is given by $w_{k}=\frac{1}{\xi_{k}+\epsilon}/\sum_{j=1}^{\phi^{b}_{1}}\frac{1}{\xi_{j}+\epsilon}$. $\epsilon$ is introduced to avoid the indeterminate case. According to the hierarchy, it is obvious that $\xi_{1}<\xi_{2}<...<\xi_{\phi^{b}_{1}}$. Thus, $w_{1}>w_{2}>...>w_{\phi^{b}_{1}}$.

Let us consider a feature $\upsilon\in\{b,C,\lambda,\kappa,\beta,\alpha\}$. The values of the feature on the paths, $f$ and $g$, are defined by

[TABLE]

The distance between $\upsilon^{f}$ and $\upsilon^{g}$, weighted by the importance factor, is given by

[TABLE]

This distance is computed for each $\upsilon\in\{b,C,\lambda,\beta,\kappa,\alpha\}$. The overall distance between the paths $f$ and $g$ can be expressed as a weighted average of individual distances.

[TABLE]

For simplicity, we take $\delta_{i}=\frac{1}{6}\forall i$ and consider the final distance as the intrinsic distance between the neurons. For classification, we determine $\delta$ through optimization using $maximizing-interclass-minimizing-intraclass$ distance strategy (See algorithm 2 in the Appendix). We term $\delta$ as the relative importance of features.

3.3 Path assignment and self-similarity

Let the number of paths in $H_{1}$ be $|\Gamma_{1}|=n_{1}$. Similarly, for $H_{2}$, this value is $|\Gamma_{2}|=n_{2}$. Without loss of generality, let us assume $n_{1}\leq n_{2}$. Using eq. 3, the cost matrix of paths between $G_{1}$ and $G_{2}$ becomes $\mathcal{D}$ ($\mathcal{D}_{ij}=\mu^{ij},~{}i\in\Gamma_{1},j\in\Gamma_{2}$). By applying an analogy for the path assignment as a job assignment problem with $n_{1}$ workers and $n_{2}$ jobs, we adopt the Munkres algorithm to find the optimal assignment of jobs to the workers from $\mathcal{D}$. In most cases, including inter- and intra-cellular neurons, the job assignment problem is an unbalanced $n_{1}<n_{2}$. We append $(n_{2}-n_{1})$ zero rows to $\mathcal{D}$ to serve as dummy workers. ElasticPath2Path batabyal2018elasticpath2path employed this technique and resulted in an output of $n_{1}$ optimally matched paths between $G_{1}$ and $G_{2}$. However, this is essentially subgraph matching, which may lead to misclassification while dealing with two structurally similar, but different, cell types. For example, hippocampal CA3 pyramidal and cerebellar Purkinje cells have similar dendritic branch patterns, but significantly different number of paths. To resolve this problem we devise an algorithm, given in the Appendix, by applying Munkres algorithm repeatedly to obtain a full-tree matching. To meet such criterion, the algorithm gives $n_{2}$ pair of paths. Let the pair be $(\gamma_{11},\gamma_{21}),...,(\gamma_{1n_{2}},\gamma_{2n_{2}})$, where $\gamma_{1i}\in\Gamma_{1}$ and $\gamma_{2j}\in\Gamma_{2}$. Recall that $n_{1}<n_{2}$, which implies that some of the $\gamma_{1i}$ are repeated while forming the pair. Finally, the distance between $G_{1}$ an $G_{2}$ is given by

[TABLE]

Let $\lfloor\frac{n2}{n1}\rfloor=T$. Then, this procedure to find the correspondence is termed as $T-$regular matching, which in turn can be thought of $T$ nearly self-similar structures akin to a fractal system. The detailed algorithm is provided in the Appendix.

There are four modules that are sequentially executed in the algorithm. The first module mathematically deciphers the relatively self-similar anatomy of a larger neuron compared to a smaller one, yielding the number of copies of the smaller one needed to stitch together to approximately obtain the larger one. The routine runs for $\lfloor\frac{n_{2}}{n_{1}}\rfloor$ times, which indicates that each path in neuron 1 (containing $n_{1}$ paths), is matched with $\lfloor\frac{n_{2}}{n_{1}}\rfloor$ paths of neuron 2 (containing $n_{2}$ paths). Here $n_{2}>n_{1}$.

The second module runs for the remaining unpaired paths of neuron 2. The assigned correspondence is added to the list of paired paths from the first module. However, not all the pairs are anatomically consistent. This is dictated by an internal constraint of Munkres algorithm, in which the assignment is carried out without replacement. In the Munkres algorithm, if one ‘worker’(a path from neuron 1) is assigned a ‘job’ (a path from neuron 2), then the ’job’ is not available for further assignment. Therefore, if the distance between two paths is significantly large, it demands further inspection whether the pair of paths is morphologically different to each other or the algorithmic constraint induces the large distance value. This motivates us to introduce the third module.

In the third module, we inspect the pair of paths having distances more than a threshold. The threshold is selected based on the skewness, median and standard deviation of the distance values. As mentioned earlier, in order to find the distance of a feature on two paths (eq. 2), we append zeros to the path having relatively fewer number of locations than the other. The choice of traversal order dictates to which side the zeros are appended. Notice that more zeros lead to higher distance value between paths, and this happens only when there is significant mismatch in the highest level of hierarchy. This fact can be interpreted from the morphological viewpoint. A path with a large number of bifurcation locations (so, large maximum hierarchy value), called a central path of a neuron, exploits the environment of the neuron extensively when compared to path with fewer number of bifurcations. Unless otherwise required, a path with large hierarchy values should not be compared with a path with much smaller maximum hierarchy value. The highest level of hierarchy values of two paths are given by $h_{1}$ and $h_{2}$ with $h_{1}<h_{2}$. We set a criteria that if $|h1-h2|>\frac{max[h1,h2]}{2}$, we do not consider the distance between the pair, and opt for the best match in terms of minimum distance for each path of the pair separately. This is outlined in the reassignment module. The reassigned pairs are added to the list of paired paths serving as the list of correspondence.

3.4 Path morphing

Once the correspondence of paths between neurons is established, it is imperative to know the structural similarity between the paths - whether a pair of paths are structurally similar to each other, or the pair is structurally incoherent but the algorithm outputs such a pair due to its internal constraints. This is achieved in two ways: with a visual representation by morphing the paths of one neuron to that of the other using an elastic framework, and by extracting path statistics.

A rooted path of a neuron can be considered as an open curve as shown in Fig. 5batabyal2018elasticpath2path ; srivastava2011shape . Each location on the path can be considered as a function of a parameter,$t\in[0,1]$. The square root velocity function (SRVF) that is applied on a location $f(t)$ is defined as $q(t)=\frac{\dot{f}(t)}{\sqrt{||\dot{f}(t)||}}$. For a pair of paths $i$ and $j$, we obtain $q_{i}$ and $q_{j}$, which assists in retrieving the intermediate deformations as linear combinations of $q_{i}$ and $q_{j}$ given by $q_{ij}^{n}=q_{i}(1-n)+nq_{j};~{}n\in[0,1]$. $n$ denotes the intermediate algorithmic time steps. Although the deformations are exhibited using the 3D coordinates of the locations of a path, the deformations can also be computed in the feature domain. An example of the continuous morphing process between two pyramidal neurons from the secondary visual cortex of the mouse is shown in Fig. 6. The $15$ paths of the former neuron merge with $11$ paths of the latter upon termination of the morphing process. This implies that more than one path of the first neuron have the same final destination path of the second neuron. It is noted that our algorithm does not consider the costs that are incurred by the merging or splitting of paths during progression. The assessment of such costs requires biophysical measurements of neurons, such as metabolic cost of merging or splitting of branches. Therefore, the cost between paths in eq. 3 is proportional to the cost of structural disparity instead of biophysical costs.

The prime question is: why do we need to inspect intermediate deformations? Statistical assessment of anatomical similarities between paths is sufficient to validate the correspondence that is obtained from the Munkres algorithm. However, to make the correspondence necessary, the intermediate deformations should comply with key cell-type characteristics srivastava2011shape . So we use the SRVF framework to show the deformations so that any noticeable incoherence can be attributed to the feature selection, distance measurement, or both algorithms even though we might obtain improved classification accuracy in the end.

4 Datasets and results

We validate the approach on two datasets that are collected from a centrally curated on-line repository of 3D reconstructed neurons, Neuromorpho.org ascoli2007neuromorpho . To demonstrate the strength of our approach, one dataset is compiled for intraclass and the other one for interclass analysis and comparison.

4.1 Dataset-1 (Intraclass)

This dataset contains 3D-traced neurons from $6$ distinct regions of the mouse neocortex. The regions with their cortical locations are visual-1 or primary visual (occipital), visual-2 or secondary visual (occipital), prelimbic (prefrontal), somato-1 or primary somatosensory (somatosensory), motor-1 or primary motor (frontal), and motor-2 or secondary motor (frontal).

We experiment with $62$ neurons of motor-1, $68$ neurons of motor-2, $24$ neurons of prelimbic, $204$ neurons of somato-1, $237$ neurons of visual-1, and $30$ visual-2 neurons with $625$ neurons in total. The neurons vary widely in their morphological characteristics, such as the number of paths in each neuron. The histogram of paths corresponding to each category is shown in Fig. 7.

Next, we investigate the relative importance of each feature (mentioned in section 3.1) in terms of $\delta$ when comparing a set of classes. For space constraint, we provide $\delta$ values separately for each pair of classes and all the classes taken together. The relative importance is listed in Fig. 8 by a pool of pie charts. A set of class-specific inferences regarding the relative importance is enlisted in the figure description. Whereas the pie charts present a comprehensive view of feature strength. In practice, however, the values are required to report the distance between a pair of neurons. The values are reported in Table 1.

It is worthwhile to note that although there is significant variance in feature strength when all pairs are considered, the distribution approximates a uniform distribution when all classes are taken. This result endorses the selection of our features for all-class classification tasks. It is also important to mention that this framework can incorporate any set of path-specific features, not restricted to our selected features only.

In order to verify the consistency of the path correspondence obtained from Munkres algorithm (provided in section 3.3), we statistically evaluate each pair of paths in the correspondence list using pyramidal neurons from two different regions (the somatosensory cortex and secondary visual cortex.) The neuron from the somatosensory cortex (neuron-2) contains $28$ rooted paths, while the other (neuron-1) has $11$ rooted paths. Table 2 provides the exhaustive list of path correspondences, distances between corresponding paths, correspondences obtained by a competitive approach named ElasticPath2Path batabyal2018elasticpath2path , and the best correspondences of the paths of neuron-2 with that of neuron-1. Notice that the best correspondence of a path $f$ of neuron-2 is the path $g$ of neuron-1, which yields minimum distance with $f$. Whereas, the Munkres algorithm works on the criteria where the sum of path distances (in our case $11$ paths at a time) is minimized.

In Table 2, the two columns on the left enumerate the pair that consists of the path number of neuron-1 and that of neuron-2. A subset of paths of neuron-1 is repeated because neuron-2 (with $28$) has more paths than neuron-1 (with $11$). So from neuron-1 to neuron-2, the correspondence is a surjective mapping. This is in contrast with ElasticPath2Path, where the mapping is bijective and, as result of that, the algorithm outputs only $11$ pairs in the correspondence list. The rest of the $28-11=17$ paths are left unmatched, yielding a solution of the subgraph matching problem. The unmatched paths are marked with ‘NA’ in the fourth column.

The last column, tagged as the best match, identifies only $\{4,5,8,9\}$ path indices out of $11$ paths of neuron-1. Nevertheless, this best matching algorithm also elicits a potential solution for subgraph matching. There are certain extreme cases where all paths of one neuron are matched with only one path of the other neuron, posing degenerate solutions of the neuron matching problem. We mark the correspondences in yellow, where the results of our algorithm and best match coincide.

Recall that neuron-1 has $11$ paths and neuron-2 contains $28$ paths. Careful observation of the first column of the table suggests that the set of numbers $\{1,2,...,11\}$ is repeated twice in the serial order followed by $6$ path indices which are $\{9,4,10,8,2,5\}$. Here, Munkres algorithm is applied thrice. Each time Munkres algorithm outputs $11$ pairs of paths for correspondence. Therefore, the first two passes encompass $11*2=22$ pairs leaving $28-22=6$ paths of neuron-2 unassigned. Before applying the third pass, the cost matrix $\mathcal{D}$ is cropped with a dimension $\mathcal{R}^{11\times 6}$. The cropped cost matrix is then transposed ($\mathcal{R}^{6\times 11}$), zero-appended ($\mathcal{R}^{11\times 11}$) and subjected to Munkres. The above observation also indicates that $2$ self-similar copies of neuron-1 approximates neuron-2 in the sense of minimum path to path distance. Therefore, the relative fractal index of neuron-2 with respect to neuron-1 is $2\frac{6}{11}$ or $2.545$.

The question is: can the arithmetic average (which is $0.92$) of the ‘Distance’ column of Table 2 be regarded as the final distance between the neuron-1 and neuron-2? Unfortunately, it is not. The reason is explained in section 3.3 and reiterated briefly in the following sentences. After computing the correspondences (column-1 and column-2), we identify the defective sets of pairs for which there are significant differences in the hierarchy levels. The larger the difference, the larger the number of zeros that are appended to each feature on the path, raising the chances of technical error in the final distance value. In the table, the defective pairs are emboldened with blue color. We delete these pairs and replace the correspondences of path indices $18$ and $9$ (neuron-2) with their best matches from neuron-1. Note that path indices $7$ and $4$ of neuron-1 have already been matched with other paths of neuron-2, which are $7--28$, $4--11$, and $4--12$. Therefore, those paths are not subjected to re-assignment. After inserting the best matches for the path indices $18$ and $9$ (which are $9$ and $5$ from neuron-1 respectively), the corresponding distance values are noted. This is described in the fourth routine, ‘Reassignment’ of algorithm 2. The final distance between the two neurons turns out as $0.90$ (rounded off). The competitive approach, ElasticP2P produces a distance value of $0.67$, which implies that the two neurons are more similar. This is discordant with the fact that the two neurons are sampled from two different regions and have two distinct arbor types. This disagreement can be explained due to subgraph matching nature of ElasticP2P. Neuron-1 with $11$ paths is well-matched with a part of neuron-2. However, the rest $17$ paths of neuron-2 are structurally dissimilar with neuron-1. In this case, our method, NeuroPath2Path performs significantly better in distinguishing two neurons in terms of distance owing to its full-graph matching property.

For classification, we compute the importance values $\delta$ from (3) and show them in Table 1. The importance values are applied to compute the distance between a pair of neurons. Using our distance function, we resort to the K nearest neighborhood classifier. We randomly partition the dataset into our training and test set using a constant ratio and rerun the experiment $5$ times. The ratio that we maintain is $0.1$ and $0.2$ as train and test datasets. As the number of paths that a neuron has is a distinguishable feature for certain classes, we devise a strategy to test each neuron from the testing dataset. For a neuron with number of paths as $n_{P}$, we seek candidate neurons from the training set with the number of paths ranging in $[n_{P}-L,n_{P}+L]$. Overall, NeuroPath2Path contains two hyperparameters, $K$ (number of nearest neighbors) and $L$. THis step is followed by the identical testing procedure while considering the interclass dataset. We fixed $L=50$ for our experiments. As noted before, we adopt the reverse and standard branch orders for dataset-1 and dataset-2, respectively.

With the ratio of train and test as $8:2$, the confusion matrix of NeuroPath2Path for an instance of random partition of data is shown in Fig. 14. It can be seen that, while NeuroPath2Path distinguishes Motor-1, Visual-1, and Visual-2 quite well, the class of Somato-1 is significantly misclassified with Motor-1, Motor-2, and Visual-1, leading to a decline in the classification score.

Next, in Fig. 11, we illustrate the comparative performance of NeuroPath2Path against TMD and NeuroSoL. The train to test ratio is set at $8:2$. NeuroSoL shows an erratic behavior as $K$ increases. TMD offers a consistent margin of classification accuracy per $K$. Here, at a given value of $K$, margin implies the difference between the maximum and minimum scores of $5$ experiments which are independently instantiated by randomly partitioning the dataset with $8:2$ train$:$test ratio. Fig. 11 suggests that NeuroPath2Path achieves peak performance when $K$ is set as $7$, but with a noticeable margin.

To scrutinize the performance of TMD and NeuroPath2Path, we routinely inspect the class-wise retrieval accuracy, a crucial metric which is obscured in Fig. 11 due to the averaging effect. The result is shown in Fig. 12. In a majority of cases, despite comparable overall classification scores, TMD tends to be affected by class imbalance, leading to significantly poor accuracy for few classes.

4.2 Dataset-2 (Interclass)

The second dataset consists of 3D-reconstructed neurons that are traced from five major cell types of the mouse: ganglion, granule, motor, Purkinje, and pyramidal. We experiment with an imbalanced pool of $500$ ganglion cells, $490$ granule cells, $95$ motor cells, $208$ purkinje cells, and $499$ pyramidal cells, where the corresponding SWC files are obtained from the neuromorpho repository. The cell-specific distribution of paths is shown in Fig 13.

For classification, we compute the important weights $\delta$ of each features, and due to space constraints, the $\delta$ values are enumerated in Table 3 for pairwise classes and the case with all the classes taken together. The importance-weighted distance value, $\mu^{fg}$ in (3) is used to compute the distance of a pair of neurons. We empirically find that the nonlinear transformation of $\mu^{fg}$, given by $\frac{1}{1+exp(-\mu^{fg})}$, yields an improved classification performance.

With a train$:$test ratio as $8:2$, one instance of the confusion matrix, obtained by NeuroPath2Path is provided in Fig.14.

We demonstrate the effectiveness of NeuroPath2Path over two state-of-the-art approaches - Topological Morphological Descriptor (TMD) kanari2018topological and NeuroSoL batabyal2017neurosol . For each value of $K$, we randomly partition the dataset $5$ times maintaining a constant $9:1$ ratio between the train and test datasets. In short, for every $K$, we obtain $5$ accuracy scores, which are plotted in Fig. 15.

TMD appears to be very consistent in accuracy and range scores, achieving an accuracy of $85.02\%$ when $K=5$. However, while computing the confusion matrices of the classification scores obtained by TMD, we notice that in the majority of instances, the correct classification of motor cells is abnormally low and approaches $0\%$ in some cases. It is important to notice that Dataset-2 has an imbalance in terms of the number of examples in each cell category, with motor cells containing the lowest ($95$) and ganglion cells containing the highest ($500$) number of examples. This fact is unobserved in Fig. 15 due to the averaging effect. We adopt the metric, class-wise accuracy of retrieval, and present the results in Fig. 16. It is evident that NeuroPath2Path exhibits strong resilience against the class imbalance problem.

Similar to the train and test ratio of Dataset-1, we conduct experiments using the ratio of $0.1$ and $0.2$ separately. The classification scores are given in Fig. 17.

5 Conclusion

NeuroPath2Path follows a graph-theoretic approach that utilizes path-based modeling of neuron anatomy and provides a visualization tool by way of a geometric model that aids in performing continuous deformation between two neurons. NeuroPath2Path offers several advantages. The decomposition of a neuron into paths can be viewed as an assembly of individual circuits from the terminals to the soma, integrating semi-local features that act as path descriptors. Next, instead of subgraph matching, NeuroPath2Path does not leave a single path unassigned, culminating in a full-graph matching algorithm. The matching algorithm presents the notion of relative fractality and path correspondence, and incorporates physiological factors, such as decaying importance of features along the path and exploratory /competitive behavior for resource exploitation.

NeuroPath2Path also precisely investigates the feasibility of algorithmic constraints (such as on Munkres algorithm) on the structural repertoire of neuronal arbors, and thereby enforcing criteria, such as hierarchy mismatch. During classification, NeuroPath2Path delivers resilience to the class imbalance problem. In the future, in order to explore the full potential of the approach beyond classification, we aim to augment NeuroPath2Path in two major domains - morphological analysis and structural transformation of microglia cells, and in progressive degradation of neuronal paths in neurodegenerative diseases.

Appendix

Description of features

We extract a set of discriminating features on each path $f_{i}\in\Gamma$ of $H$, which are bifurcation angle ($b_{i}$), concurrence ($C_{i}$), hierarchy ($\xi_{i}$), divergence ($\lambda_{i}$), segment length ($\beta_{i}$), tortuosity ($\kappa_{i}$), and partition asymmetry ($\alpha_{i}$). •  Bifurcation angle is a key morphometric that dictates the span and the spatial volume of an arbor. It is hypothesized that the span of an arbor at each level of bifurcation depends on the bifurcation of its previous level lopez2011models ; batabyal2018neurobfd ; bielza2014branching , suggesting the influence of Bayesian philosophy. This organizational principle is utilized in several stochastic generative models lopez2011models for the synthesis of specific neuron cell types. The sequence of bifurcation angles at bifurcation vertices located on a path of a neuron captures local geometry. For example, a sequence of non-increasing bifurcation angles from the root to the dendritic terminal of a path indicates the pyramidal shape geometry of the neuron. For a location with multifurcation, we use the maximum of the bifurcation angles computed using pairwise branches originated from that location towards the dendritic terminals.

•  Concurrence, hierarchy and divergence encode the effect of phenomenological factors, which are exploration (ex. Purkinje fanning out rostrocaudally) and competition (ex. retinal ganglion cells), that contribute in the growth of a neuron. The definition of concurrence and hierarchy are already given in section 1.1. The divergence of a location on a path, $f_{i}$ is proportional to the repulsive force that the location experiences from its neighborhood path segments. Let $C_{f_{i}}$ be the sequence of concurrence values of the path $f_{i}\in\Gamma$ when one visits the locations from the root to the dendritic terminal. As an open curve, each path can be parameterized with the parameter $t\in[0,1]$. $C_{f_{i}}(t_{s})=k;~{}t_{s}\in[0,1]$ indicates that $k(\leq|\Gamma|)$ paths share the location $t_{s}$ on $f_{i}$. The divergence $\lambda$ of a location $f_{i}(t_{s})$ is defined as \lambda(f_{i}(t_{s}))=$$1_{\{f_{j}~{}|~{}|f_{j}(t)-f_{i}(t_{s})|\leq\epsilon,f_{j}\neq f_{i},f_{j}\nsucc f_{i}\}}. Here, $1$ is the indicator function computing the number of such $f_{j}$s which follow the conditions $|f_{j}(t)-f_{i}(t_{s})|\leq\delta,f_{j}\neq f_{i}$ and $f_{j}\nsucc f_{i}$. The first condition implies that a location of $f_{j}$ has to be in the $\epsilon$ neighborhood of $f_{i}$. $f_{j}\nsucc f_{i}$ indicates that the location of bifurcation at which $f_{j}$ deviates from $f_{i}$ does not appear after $f_{i}(t_{s})$ on the path $f_{i}$.

•  Tortuosity and partition asymmetry are two important anatomical features of a neuron. Tortuosity refers to the amount of ‘zig-zag’ or bending of a path. Let us take a segment on a path $f_{i}$ as $f_{i}([t_{1},t_{2}]);~{}0\leq t_{1}<t_{2}\leq 1$. Let there be $m-1$ locations in $[t_{1},t_{2}]$. The tortuosity of the segment is defined as $\kappa=\frac{\sum^{m}_{j=1}||f_{i}(t_{j+1})-f_{i}(t_{j})||_{2}}{||f_{i}(t_{2})-f_{i}(t_{1})||_{2}}$ with $t_{m+1}=t_{2}$. Partition asymmetry accounts for how the size of a neuron tree varies within the neuron. We use a variant of caulescence, proposed in brown2008quantifying , as a measure of tree asymmetry. Caulescence at a bifurcation location is evaluated by way of $\alpha=\frac{|l-r|}{l+r}$, where $l$ is the size of the left tree and $r$ of the right tree of the bifurcation vertex. We define the size of a tree by the number of paths or equivalently the number of dendritic terminals. Note that the quantity $(l+r)+1$ is the concurrence value of the bifurcation vertex.

Weight determination

Let the combined distance vector containing the individual feature distances be $D^{fg}=[d(b^{fg})~{}d(C^{fg})~{}d(\lambda^{fg})~{}$ $d(\kappa^{fg})~{}d(\beta^{fg})~{}d(\alpha^{fg})]^{T}$. The corresponding unknown weight vector is $\delta=[\delta_{1},...\delta_{6}]$. While comparing two neurons of sizes $N$ and $M$ with $N<=M$, the distance computation after applying the Munkres algorithm repeatedly will produce $M$ pairs of paths, indicating $M$ such $D^{fg}$s. The desired characteristic of each component of $\delta$ is positivity. In addition, we enforce $\sum\delta_{i}~{}=~{}1$, implying a probability estimate. $\delta_{i}$ thus indicates the relative importance of the feature $\upsilon_{i}$.

We adopt the constrained maximizing-interclass -minimizing-intraclass distances strategy to find our desired $\delta$. Mathematically,

[TABLE]

The first term in the above equation encompasses all the distances between neurons from pairwise classes. The second term encodes the intraclass distances, implying the distances between neurons for each class. The third term enforces positivity of each weight $\delta_{i}$. This is a logarithmic barrier penalty term that restricts the evolution of $\delta$ at intermediate iterations to the region where $\delta>\bar{0}$. The last term accounts for the probabilistic interpretation of $\delta$. $S$ is the number of classes.

Eq. 5 is solved by using gradient descent. The equation and its derivative can be simply written as,

[TABLE]

We use this derivative term in the following algorithm 1 to obtain optimal $\delta$.

Distance between neurons

The algorithm to find distance between a pair of neurons consists of four stages - finding self-similarity (routine-1), remaining path assignment (routine-2), finding pairs with hierarchy mismatch (routine-3) and reassignment of the defective pairs (routine-4).