Spatial registration of neuron morphologies based on maximization of volume overlap

Background Morphological features are widely used in the study of neuronal function and pathology. Invertebrate neurons are often structurally stereotypical, showing little variance in gross spatial features but larger variance in their fine features. Such variability can be quantified using detailed spatial analysis, which however requires the morphologies to be registered to a common frame of reference. Results We outline here new algorithms — Reg-MaxS and Reg-MaxS-N — for co-registering pairs and groups of morphologies, respectively. Reg-MaxS applies a sequence of translation, rotation and scaling transformations, estimating at each step the transformation parameters that maximize spatial overlap between the volumes occupied by the morphologies. We test this algorithm with synthetic morphologies, showing that it can account for a wide range of transformation differences and is robust to noise. Reg-MaxS-N co-registers groups of more than two morphologies by iteratively calculating an average volume and registering all morphologies to this average using Reg-MaxS. We test Reg-MaxS-N using five groups of morphologies from the Droshophila melanogaster brain and identify the cases for which it outperforms existing algorithms and produce morphologies very similar to those obtained from registration to a standard brain atlas. Conclusions We have described and tested algorithms for co-registering pairs and groups of neuron morphologies. We have demonstrated their application to spatial comparison of stereotypic morphologies and calculation of dendritic density profiles, showing how our algorithms for registering neuron morphologies can enable new approaches in comparative morphological analyses and visualization. Electronic supplementary material The online version of this article (10.1186/s12859-018-2136-z) contains supplementary material, which is available to authorized users.


Choosing a discrete set of parameters for exhaustive search
The set of discrete values for a transformation parameter is chosen centering on an existing estimate, when one is available, or about [0, 0, 0] (for translation and rotation) or [1, 1, 1](for scaling) when there is none. This is done for each parameter axis separately, by successively adding or subtracting a certain step size from the center for translation and rotation, and by successively dividing and multiplying the center with a step size for scaling until the bounds are reached. This is illustrated in the following sections.

Scaling:
Given a previous estimate w = [s x , s y , s z ], a step size ∆ (see below), and bounds [l, u] for each parameter, the scaling parameter axes are sampled as follows: where lf x = t x * ∆ f loor (log ∆ l/tx ) , similarly for l y and l z ; and uc x = t x * ∆ ceil(log ∆ u/tx ) , similarly for uc y and uc z .

Choosing step sizes
The step size essentially is the inter-sample-interval for the discrete set of parameters chosen for exhaustive searches. To maximize the chances of an exhaustive search that uses a discrete set of parameter values to approximately find the global minima, the step size must be chosen so that the the difference in spatial dissimilarity corresponding to adjacent parameter values is small. Differences in spatial dissimilarity corresponding to adjacent parameter values depend on the voxel size used to discretize the volumes of the morphologies. For example, in most cases, when translations smaller than the voxel size are applied, only a small fraction of the set of voxels representing the morphology change and hence the spatial dissimilarity between the morphologies also changes by a small amount. Therefore, we choose a step size depending on the voxel size being used.
For a voxel size g, the step size ∆ for discretizing the parameter space is chosen as follows: where maxdi st is the distance of the farthest point of the morphology from its centroid. Similar to translation, the step size for rotation (scaling) is chosen as the angle (scale) which moves a point by a distance of one voxel size. Since this angle (scale) depends on the distance of the point being moved from the centroid of the morphology, we choose the smallest of such angles (scales), the one which corresponds to the farthest point of the morphology from its centroid.

Error in estimation
Using an exhaustive search over a discrete set of transforms to find a transform that minimizes spatial dissimilarity between morphologies inherently introduces an uncertainty since it does not search through all possible transforms, but only a subset of them. However, with exhaustive searches that use equally spaced set of transforms, this uncertainty is bounded above by the spacing between transforms and thus can be made very small by choosing a closely spaced set of transforms. Additionally, the estimate could contain an uncertainty that is due to volume discretization of the morphologies. To illustrate this effect, consider the case where two morphologies have just one point each. Since Reg-MaxS calculates spatial dissimilarity based on voxel overlap, any transform that brings the two points within the same voxel would be a solution. For cases where the points are near voxel boundaries, an exhaustive search over a discrete set of equally spaced transformations would prefer a transform that brings the two points into the same voxel over a transform that doesn't bring them into the same voxel but brings them closer. When morphologies have more points, such cases arise more often and this uncertainty becomes more prominent when working with experimentally acquired morphologies with hundreds or thousands of points. Thus Reg-MaxS uses mainly the uncertainty in an estimate caused by volume discretization to determine the plausible region of transforms for the next exhaustive search.
Further, for rotational and scaling differences, points at different distances from the center will have different estimation uncertainties. However, farther points dominate the estimation procedure as they are more sensitive to transformations. Hence, we use the farthest point when calculating estimation uncertainty.

Translation:
Since we apply translations along orthogonal axes, the errors can be calculated independently. This implies that, for each parameter, we need to consider the effect of discretization only along its corresponding axis. The two extreme cases of maximum positive and negative errors occur when discretization causes the separation between two points to maximally decrease and increase respectively. These are illustrated in Fig.  1(a) for the case where each morphology has just one point each. These two points are separated along an axis (X, Y or Z) by a distance of l and the voxel size used for discretization is g. In this case maximal errors in estimating translation parameters can be calculated as follows:.
Translation difference between discretized morphologies δ t = l Smallest possible difference before discretization δ min = l − g ⇒ Maximum negative error δ − = g Largest possible difference before discretization δ max = l + g ⇒ Maximum positive error δ + = g Thus we can expect an error of plus-or-minus one voxel size when estimating translation parameters.

Scaling:
Consider subfigure (a) in S4_Figure with similar notations as above. Let the points be at a large distance D from the centroid of the morphology. For a given value of l, we have:

Rotation:
Consider subfigure (b) in S4_Figure, with similar notations as above. Since rotation about an axis causes a point to move in a 2D plane, we need to consider the effect of discretization along two axes. The greatest possible increase and decrease of separation between two points will happen when they are separated along the diagonal of a voxel face as shown in subfigure (b) in S4_Figure.
We calculate error for a point farthest from the centroid of the morphology, which is at a distance D. Assuming l D, since sin −1 x ∼ x for small x and for a voxel size of g, we have, Rotation difference between discretized morphologies δ r = l D Smallest possible difference before discretization δ min = l−g Significance of algorithm parameters and how to set them Bounds for transformation parameters: As part of Reg-MaxS, we estimate rotation, translation or scaling errors between morphologies by estimating corresponding parameters. Since the algorithm uses exhaustive search over the parameter space at successively decreasing voxel sizes, it needs lower and upper bounds to sample the parameter space.
These bounds define ranges in which these parameters are expected to be found. If the actual solution is outside these bounds, the solution given by Reg-Max-S might be a local minimum within the specified parameter bounds or at the bound closer to the actual solution, which could be sub-optimal. Setting these parameter bounds appropriately, either using good contextual estimates or by experimentation, can significantly improve the spatial overlap of the final solution reached.
The runtime of exhaustive search, which is the largest component of Reg-MaxS runtime, increases roughly linearly with the number of sample points of the parameter space searched. The total number of parameter sample points searched is the product of the number of sample points along each axis. For the highest voxel size, these are determined by the input bound parameters and the step size. For exhaustive searches at lower voxel sizes, this is fixed as the bounds are chosen depending on the current estimate error. Thus, these bound parameters affect the runtime of Reg-MaxS directly.
Since we expect hemispheric reflections to be taken care of beforehand, multiple morphologies of a stereotypic neuron obtained using similar imaging and reconstruction methods are not expected to have large translational, rotational or scaling differences.

Tolerance for transformation parameters:
Tolerances need to be set for each of the translational, rotational and scaling parameters. They need to be set to values lower than their corresponding step sizes at the lowest voxel size. After successively refining parameter estimates at several decreasing voxel, a final exhaustive search is run over the parameter space using these tolerances as the step size. Note that these do not set the final accuracy of the estimate. They are set by parameter step sizes calculated at the lowest voxel size.

Sequence of Voxel Sizes:
The algorithm requires a sequence of voxel sizes which are used to discretize morphologies and are the most important set of parameters for Reg-MaxS. The lowest voxel size influences the final accuracy of parameter estimates. Hence the lowest voxel size must be chosen small, enough to preserve the features of the arbor being studied but large enough to average out the inter-specimen variability in dendrites.
For exhaustive searches at all voxel sizes other than the highest, the bounds of each parameter are chosen to be directly proportional to the next highest (previous) voxel size, while the step size is proportional to the current voxel size (See "Choosing step sizes"). Thus the number of sample points of the parameter space to search at each voxel size is directly proportional to the cube of the ratio of successive voxel sizes. Thus this ratio affects the runtime of Reg-MaxS and must be set to a small value. However, if this ratio is too small, parameter refinement at each voxel size is small and the algorithm would require refinements at several voxel sizes to achieve enough accuracy.

Initial Reference:
The choice of the initial reference has a major impact on the outcome of the both the PCA-based method and Reg-MaxS-N. This initial reference has to be the best representative morphology of the neuron in the set being studied for the algorithms to perform the best. It must have the most complete arborization among the morphologies being registered and also the least amount of artifacts. Figure 1: Illustrating errors in estimating best transformation parameters due to volume discretization Squares indicate representative voxels used to discretize the morphologies. Points in blue are discretized points, red and green are the points before discretization for the cases that lead to maximum negative and maximum positive errors, respectively. (a) For the case of estimating translation and scaling parameters. The two discretized points of the morphologies are separated along an axis (X, Y or Z) by a distance of l. (b) For the case of estimating rotational parameters. This is for the case where the two points are at an approximate large distance of D from the centroid of the morphology and are separated diagonally along a voxel by a distance of l.

Algorithms
Algorithm 1: Multi-scale estimation algorithm Data: the two morphologies, a set of n(> 0) monotonically decreasing voxel sizes [g 1 , g 2 , · · · , g n ], bounds for each parameter being estimated and a minimum tolerance(∆ min ) Result: transform parameter estimates 1 i ← 1 2 cur r entBounds ← initial parameter bounds 3 foreach g i i n [g 1 , g 2 , · · · , g n ] do 4 stepSi z e ← step size at g i 5 Estimate transform parameters which maximize spatial overlap with cur r entBounds, stepSi z e and g i Apply this transform to cur r entMor ph to generate tempMor ph 10 cur r entMor ph ← tempMor ph 11 par ams ← parameters of the transform chosen above 12 if tempMor ph has higher spatial dissimilarity compared to cur r entMor ph then