Action of Diffeomorphisms on Discrete Samplings

In the following sections, we use \(C^k\) to represent the space of continuous functions with k continuous derivatives, where the domain and range can be inferred by the context. We model a neuronal branch (dendrite or axon) as a regular 3D curve \(c:[0,L] \rightarrow \mathbb ^3\), \(c \in C^k\), and \(\vert \dot \vert > 0\), where L, without loss of generality, is the arc length of the curve (Younes, 2010). When a neuronal curve is traced, it is typically stored as a sequence of points \(\\}_^n\), where the independent variables \(t_i\) can be taken to be the indices of the points. When there is a diffeomorphism between coordinate systems \(\phi : \mathbb ^3\rightarrow \mathbb ^3\), these traces are mapped via the group action:

$$\begin \phi \cdot \_^n&= \_^n \end$$

We want to extend the space of traces, and the associated action, to include derivatives of the underlying curve denoted \(\partial _t c\). This can be done using the jet space \(J^k\). In our setting, \(J^k=[0,L] \times X^\), where an element of \(X^\) is a \(k+1\)-tuple \((x^0,x^1,...,x^k) \in (\mathbb R^3)^\) representing a position and first k derivatives of a curve in \(\mathbb R^3\). A \(C^k\) curve \(c: [0,L] \rightarrow \mathbb R^3\) can be extended to a curve \(\hat: [0,L] \rightarrow X^\) simply by adding derivatives, with \(\hat(t) = (c(t), \partial _t c(t), \ldots , \partial _t^k c(t)) \in X^\) (Olver, 1995).

The \(C^k\) diffeomorphisms have a natural group action on the jet space \(J^k\), ensuring the commutation between the standard action of diffeomorphisms on curves, \((\phi , c) \mapsto \phi \circ c\) and their extensions, such that the identity \(\phi \cdot \hat(t) = \widehat(t)\) holds for all curves c and times t, defining the left-hand side. For example, for \(k=2\), this provides

$$\phi \cdot (t, x^, x^, x^) = (t, \phi (x^), D\phi (x^) x^, D\phi (x^) x^ + D^2\phi (x^) (x^, x^))$$

Neuron traces, as mentioned before, involve a sequence of samples with time-stamps \(\_i)\}_^n\), identified as elements of \((J^k)^n\), the n-fold Cartesian product of \(J^k\). Our diffeomorphisms will act on such a sequence as follows:

Statement 1

For a sequence of time-stamped elements on the jet space, \(T = \)\}_^n\) in \((J^k)^n\), we define the action of diffeomorphisms

$$\begin \phi \cdot T = \)\}_^n \end$$

(1)

The fact that this operation provides an action is an established result (Olver, 1995), and the proof is provided in the Supplement. We will define k’th order discrete mapping to be the action in Eq. 1 of a diffeomorphism on a curve sampling that includes k derivatives. The axioms that define group actions are important to verify because they ensure that applying the identity transformation does not change the object, and that applying a composition of transformations is equivalent to applying the individual transformations successively. Further, group actions can exchange mathematical structure between the acting group and the set being acted upon, and they are at the core of several important theorems (Suksumran, 2016).

The k’th order discrete mapping method allows us to compute the first k derivatives of the transformed curve. We will interpolate the transformed curve using splines of order \(2k+1\) that satisfy the derivative values. For example, zeroth order mapping will produce a first order spline and first order mapping will produce a cubic Hermite spline (Spitzbart, 1960).

Error Analysis of Zeroth and First Order Mapping

Now we will examine the error introduced by zeroth order mapping, which is used by existing neuron mapping methods. First, note that under affine transformations, zeroth order mapping of piecewise linear curves introduce no error, so these results are only useful under non-affine transformations. The following results require that the curve c be parameterized by arc length. However, we note that all continuously differentiable regular curves can be reparameterized by arc length (Smale, 1958). We use \(\vert \cdot \vert\) to denote the Euclidean norm for elements of \(\mathbb ^d\), and the spectral norm for matrices.

Proposition 1 [Zeroth Order Mapping Error Bound]

Say \(\phi : \mathbb ^3 \rightarrow \mathbb ^3\) is a \(C^1\) diffeomorphism and \(c: [0,L] \rightarrow \mathbb ^3\) is a continuous, piecewise linear curve parameterized by arc length with knots \(\ < t_i\}_^n\). For the transformed curve \(f=\phi \, \circ \, c\), the zeroth order mapping defines a first order spline g which satisfies:

$$\begin \max _ \vert f(t)-g(t)\vert&\le \max _, t \in [t_, t_i]} \frac \left( \vert D\phi \, \circ \, c(t) - I \vert \vert t_i-t_ \vert + \vert \epsilon _i - \epsilon _\vert \right) \end$$

(2)

where \(\epsilon _i\triangleq c(t_i)-\phi (c(t_i))\) and \(D\phi \, \circ \, c(t)\) is the Jacobian of \(\phi\) evaluated at c(t).

This shows how the error introduced by the state of the art mapping method is related to the displacement magnitude, \(\epsilon\), and the extent to which the Jacobian of the transformation, \(D\phi\), differs from the identity matrix. Note that the bound in Eq. 2 goes to zero as \(\phi\) approaches the identity map (in which case zeroth order mapping has zero error for piecewise linear curves). It depends on the arc lengths of the original curve segments and the spectral norm of \(D\phi\), which is related to the finite time Lyapunov exponent (\(\log \vert D\phi \vert\)), a well-known quantity in field dynamics which characterizes the amount of stretching in a differentiable flow. Also, the bound applies to \(\max _ \vert f(t)-g(t)\vert\), which is not parameterization invariant, and therefore not a strictly geometric quantity. However we note that this quantity is an upper bound of the Frechet distance, which is parameterization invariant.

In this paper we demonstrate first order mapping in an effort to mitigate this mapping error. Such a method has the advantage of having superior error convergence at the knots as a consequence of Taylor’s theorem. Further, we present a set of error bounds that helps clarify the advantage of first order mapping.

Proposition 2 [Comparable Bounds for Zeroth and First Order Mapping]

Say \(\phi : \mathbb ^3 \rightarrow \mathbb ^3\) is a \(C^4\) diffeomorphism and \(c: [a,b] \rightarrow \mathbb ^3\) is a continuous, piecewise \(C^4\) curve parameterized with knots \(\ < t_i\}_^n\). For the transformed curve \(f=\phi \circ c\) defined by coordinate functions \(f=(f^0,f^1,f^2)^T\), the zeroth order mapping defines a first order spline \(g_0\) which satisfies:

$$\begin \begin \max _ \vert f(t)-g_0(t)\vert&\le \frac} \max _} \vert \partial ^_t f^j(t)\vert \left( \frac\right) ^4 + \\&\frac} \left( \frac\right) ^2 \max _, j \in \} \vert \partial ^_t f^j(t_i)\vert \left( \frac\right) + \\&\frac} \left( \frac\right) ^2 \max _, j \in \} \vert \partial ^_t f^j(t_i)\vert \end \end$$

(3)

where \(\delta \triangleq \max _ \vert t_i-t_\vert\) and \(\partial ^_t f^j(t)\) is the k’th derivative of \(f^j\) evaluated at t. Also, the first order mapping defines a third order spline \(g_1\), which satisfies

$$\begin \max _ \vert f(t)-g_1(t)\vert&\le \frac} \max _} \vert \partial ^_t f^j(t)\vert \left( \frac\right) ^4 \end$$

(4)

and we note that the bound in (4) is tighter than the bound in (3). Further, there exists a transformed curve f and a set of knots \(\_^n\) that achieves both bounds exactly.

Thus, we have made a connection between the state of the art (zeroth order mapping) and a higher order method (first order mapping) via worst-case bounds on mapping error. The error bound for first order mapping is smaller than that for zeroth order mapping, though for any given curve, either method may produce smaller error than the other. Proofs for the propositions are in the supplement.

Software Implementation

We implemented both zeroth and first order discrete mapping in our our open-source Python package brainlit. For first order mapping, we compute one-sided derivatives at the knots of the curve from first order splines in accordance with original SWC formulation (Stockley et al., 1993; Cannon et al., 1998). Then, once the knot positions and derivatives are transformed, we generate a continuous curve in the new space using Hermite interpolation. Further details of our implementation can be found in the Methods.

First order mapping involves more computations than zeroth order mapping. First, if the Jacobian of the transformation, \(D\phi\), is not immediately available, it needs to be approximated. Our software uses a finite difference method which, to approximate \(D\phi (x)\), involves calling \(\phi\) four times, three vector addition operations, and three vector scaling operations. Next, the derivatives at the knots need to be approximated and transformed by \(D\phi\). For each line segment, our method computes one vector addition and one vector scaling to approximate one-sided derivatives, and computes two matrix-vector products to transform the derivatives (Algorithm 3). Lastly, generating and evaluating cubic splines involves more computations than first-order splines (Kincaid & Cheney, 2002). Despite these differences, first order mapping still scales linearly with the number of trace nodes, and the number of computations differs from zeroth order mapping by a constant factor of computations. In practice, both zeroth order and first order mapping take on the order of seconds for traces with thousands of nodes, which should not be a bottleneck in neuron morphology studies.

Figure 2 shows examples of our method on simulated data, compared to the zeroth order method, and the “ground truth” where we map a dense sampling of points along the first order spline of the original curve.

Fig. 2

Preserving derivative information can mitigate errors when transforming discretized curves. a-b Applying a nonlinear deformation field to a single line segment (a) using zeroth and first order mapping (b). c-d Applying a nonlinear deformation field to a piecewise linear curve (c) using zeroth and first order mapping (d). Zeroth and first order discrete mapping methods are shown relative to ground truth considered to be the application of the vector field to a dense sampling of the original curves

Application to Real Neurons

We applied our method to 20 reconstructed neurons in SWC format from a whole mouse brain image from the Janelia MouseLight project (Winnubst et al., 2019). We selected the first 20 SWC files that successfully downloaded from MouseLight’s NeuronBrowser repository and did not have repeat trace nodes. Neurons have a tree-like structure, and we split them into non-branching curves in order to apply our mapping methods. We follow a method introduced previously (Athey et al., 2021) where the root to leaf path with the longest arc length is recursively removed until the tree is reduced to non-bifurcating “branches”. The transformed branches are then reconnected to maintain the topology of the original trace, and therefore our methods generalize naturally from 3D curves, to the branching structure of neuron traces.

We generate random transformations using the Large Deformation Diffeomorphic Metric Mapping (LDDMM) framework described in Miller et al. and applied in Tward and Miller (Miller et al., 2006; Tward & Miller, 2017). We generate an initial momentum field by sampling Gaussian noise with zero mean and varying standard deviation, \(\sigma\). The momentum is smoothed to construct a velocity field, and integrated in time according to the conservation laws established in Miller et al. to generate a diffeomorphic transformation (Miller et al., 2006). We generated four diffeomorphisms with \(\sigma\) levels of \(80\), \(160\), \(320\) and \(640\) \(\mu m/\text \). The position and tangent displacement profiles of these four diffeomorphisms are shown in Figure 3a. We centered the neuron traces at the origin then applied the random diffeomorphisms to compare zeroth and first order mapping to ground truth (Fig. 3b-g). Ground truth was generated by upsampling the original traces to a maximum node spacing of \(2 \mu m\) followed by zeroth order mapping.

Fig. 3

Application of zeroth and first order mapping of neuron traces under diffeomorphisms derived from random Gaussian initial momenta. a Different values of \(\sigma\) produced diffeomorphisms with different position and tangent vector displacement profiles. The positions and tangent vectors sampled in the histogram were distributed as a uniform grid with a spacing of \(500 \mu m\). b-g Two examples of the diffeomorphism with \(\sigma =640\) applied to neuron traces to produce zeroth and first order mappings, along with ground truth. Both examples show the original trace and the transformation (b, e), the results of the different transformation methods (c, f), and a zoomed in view of the region outlined by the dotted line to show discrepancies between the methods (d, g). Plot axes are in units of microns

For each neuron trace, we computed the discrete Frechet error from ground truth (Fig. 4a). We also wanted to measure which mapping method better matched the ground truth with respect to a neuron’s distribution of common morphometric quantities, such as path angle, branch angle, tortuosity, and segment length. We used the Kolmogorov-Smirnov test statistic to measure how much the distribution of these quantities differed from ground truth (Fig. 4b). We performed two-sided Wilcoxon signed-rank tests for each comparison and used a Bonferroni correction across the different \(\sigma\) values (Fig. 4b). Lastly, we compared the discrete Frechet errors to the average sampling period of the trace i.e. the average distance between trace nodes (Fig. 4c).

Fig. 4

Comparison of zeroth and first order mapping of neuron traces under random diffeomorphisms. a Discrete Frechet error was computed between the different order mappings, and ground truth. b Distributions of common morphometric quantities were compared to that of ground truth using the Kolmogorov-Smirnov test statistic. Differences between zeroth and first order methods were tested using Wilcoxon signed-rank test with Bonferroni correction across different values of \(\sigma\) (\(*: p \le 0.05\), \(**: p \le 0.01\), \(***: p \le 0.001\), \(****\, p\le 0.0001\)). Box plots show median, upper and lower quartiles and whiskers have a maximum length of 1.5x the interquartile range with other outlier data marked with points. c Relationship between discrete Frechet error and average sampling period (distance between trace points) under the random diffeomorphisms

To explore the effect of downsampling neuron traces on mapped morphologies, we identified non-branching nodes in straight portions of the trace, and measured the impact of removing those nodes from the trace. Specifically, we performed both zeroth and first order mapping on the segment with the node removed, and compared it to the ground truth mapping of the original segment. We determined which fraction of nodes maintained a discrete Frechet error less than one micron, serving as an estimate for the fraction of nodes which are not necessary to maintain the mapped morphology (Fig. 5c). We performed two-sided Wilcoxon signed-rank tests to compare zeroth and first order methods and used a Bonferroni correction across the different values of \(\sigma\).

Fig. 5

Counting how many nodes in MouseLight neuron traces can be removed without affecting the mapped morphology. a For each non-branching node with path angle above 170 degrees, we generated a line segment with that node removed. b We performed first order mapping on the downsampled line segment and compared the result with the ground truth mapping of the original curve. c For each mapped neuron trace, we determined the fraction of nodes where the discrete Frechet error is less than or equal to one micron under the four random diffeomorphisms. Box plots show median, upper and lower quartiles and whiskers have a maximum length of 1.5x the interquartile range with other outlier data marked with points. We performed Wilcoxon signed-rank test between the paired neuron traces at \(\alpha =0.05\) with Bonferroni correction across different values of \(\sigma\) (\(*: p \le 0.05\), \(***: p \le 10^\))