Automated Stitching of Microtubule Centerlines across Serial Electron Tomograms

Tracing microtubule centerlines in serial section electron tomography requires microtubules to be stitched across sections, that is lines from different sections need to be aligned, endpoints need to be matched at section boundaries to establish a correspondence between neighboring sections, and corresponding lines need to be connected across multiple sections. We present computational methods for these tasks: 1) An initial alignment is computed using a distance compatibility graph. 2) A fine alignment is then computed with a probabilistic variant of the iterative closest points algorithm, which we extended to handle the orientation of lines by introducing a periodic random variable to the probabilistic formulation. 3) Endpoint correspondence is established by formulating a matching problem in terms of a Markov random field and computing the best matching with belief propagation. Belief propagation is not generally guaranteed to converge to a minimum. We show how convergence can be achieved, nonetheless, with minimal manual input. In addition to stitching microtubule centerlines, the correspondence is also applied to transform and merge the electron tomograms. We applied the proposed methods to samples from the mitotic spindle in C. elegans, the meiotic spindle in X. laevis, and sub-pellicular microtubule arrays in T. brucei. The methods were able to stitch microtubules across section boundaries in good agreement with experts' opinions for the spindle samples. Results, however, were not satisfactory for the microtubule arrays. For certain experiments, such as an analysis of the spindle, the proposed methods can replace manual expert tracing and thus enable the analysis of microtubules over long distances with reasonable manual effort.

The operator diag(v) creates a diagonal matrix from a vector v ∈ R D by diag(v) = We denote a vector containing only ones as 1.
Further prerequisites To account for outliers, Myronenko and Song introduce an additional term P (x, M + 1) = 1/N to the joint distribution. This outlier term can be weighted with 0 ≤ w ≤ 1. The marginal probability of x then takes the form: We use the same outlier term in all three settings. For convenience, we will often write sums as matrices using the following equalities: P (m|x n )x n y m R = X P YR (20) Algorithm for linear alignment from orientation E-step The posterior P ( m| x) is computed by To compute P ( m, x), we add an outlier term P ( x|M + 1) = 1/N and weight the priors with 0 ≤ w ≤ 1 so that P (m) = (1 − w)/M and P (M + 1) = w/N . P ( m, x) can be computed from the priors and P ( x| m) which is given by Equation 6 for m ≤ M . The joint distribution becomes: P ( x, M + 1) = w/N , for m = M + 1.
P ( x) can be computed from P ( m, x) by marginalizing out m: The posterior becomes M-step In the M-step, we must write out the expectation of the negative log-likelihood function Q (Equation 3) and minimize it with respect to the transformation parameters R and the distribution parameter κ. Q becomes: We split the term tr(( X P Y) R) into terms that are dependent on and independent ofR:

S3
where X 2d is a matrix containing the first two columns of X and X z contains the last column of X.
The optimal rotation can be obtained with a Singular Value Decomposition from the term X 2d P Y 2d (equivalently to Myronenko and Song [28], see Umeyama [42] and Myronenko and Song [43] for details on the method). κ cannot be obtained in closed form, but it can be computed with Newton's method. We iterate Because we only compute a rotation, the endpoints are not aligned. To align them, we first note that even though we did not explicitly take the positions into account, the likelihood that two orientations match also applies to the corresponding endpoint positions. To compute an alignment of the positions, we compute the squared distances of the positions, which we weight by the posterior computed from the obtained optimal parameters Θ * , P ( m| x, Θ * ): To minimize the distances, we first eliminate the brackets and take the derivative with respect to t: Solving for t yields: We summarize the algorithm in Figure S2.
Algorithm for linear alignment from position and orientation E-step The joint distribution of x, x and m, m is (Equation 2 and 6): To compute the posterior P (m, m|x, x), we must again apply Bayes' rule. We add an outlier term as in the previous algorithm and weight the priors as before with 0 ≤ w ≤ 1. The posterior then is: S4 M-step The expectation of the complete negative log-likelihood function is given by To solve for each of the parameters during the M-Step, we first compute t by taking the partial derivative of Q with respect to t, as we did in the last section, to obtain (equivalently to Myronenko and Song [28]) We then plug t back into Equation 32. The log-likelihood becomes wherex n = x n −µ x andŷ m = y m −µ y . We then write out ||x n −sRŷ m || 2 as a product (x n −sRŷ m ) (x n − sRŷ m ), eliminate the brackets and write Q in matrix form: Q(s,R, σ 2 , κ) = + 1 2σ 2 (tr(X diag(P 1)X) − 2s · tr((X P Ŷ ) R ) + s 2 tr(Ŷ diag(P1)Ŷ))− −κ tr(( X P Y) R) + N P log σ 2 − N P log( κ e κ − e −κ ) + const.
We could not find a closed form solution for any of the parameters. Instead we search the optimum for the parameters numerically, except for t, which we already eliminated. In general, using numerical optimization in this case is problematic due to the constraintsR R = I and det(R) = 1. In our case, however, the rotation matrix is only two-dimensional and can therefore be parameterized by a single angle ρ that can be treated as one additional parameter. To find a numerical solution for this problem, usually first and second derivatives of the objective function are required. These are not hard to find from Equation 35, but it is tedious work. For the convenience of the reader, we give details on the derivative with respect to ρ below (Section Q derivatives) and list all first and second order derivatives in Figure S15. Figure S3 summarizes the algorithm.
The numerical optimization in the M-step is an unfortunate consequence of the scaling parameter s. Without the scaling, the objective function Q could probably be minimized by formulating the transformation as dual quaternions and solving the transformation parameters as described by Walker [44]. While this still would not provide a final optimization in closed form, it might be possible to calculate closed form update equations for a generalized expectation maximization. Note that we cannot use other approaches that use scaled vectors such as Wells [26], because our orientations are unit vectors.

Algorithm for elastic alignment
For the elastic transformation model, y is transformed via a displacement vector ν(y) (Myronenko and Song [28]) Here, the task is to find a function ν that gives enough flexibility to capture the deformation of the data but still maintains the overall neighborhood structure of the points. Myronenko and Song derive the displacement vectors by solving the regularized expectation of the negative log-likelihood function where Φ penalizes high frequencies of ν. See ) and the rotational part of the transformation is not explicitly accessible. In practice, we assume that the elastic transformation contains only a minimal rotation that can be ignored, because we have applied a linear alignment before. With this assumption we can rewrite Equation 37 with the orientation as p old (y m , y m |x n , x n )||x n − (y m − ν(y m ))|| 2 + N P log σ 2 + const.
Because Q d is independent of the displacement vectors µ, the optimization can be performed exactly as described in Myronenko and Song [28]. The only difference here is that the prior P (m, m|x, x) depends on κ, which must also be updated at each M-step. κ can be updated in each iteration again with Newton's method. With this strategy, the orientation only influences the result via the posterior P old (m, m|x, x).

E-step
The posterior again is computed by Bayes' rule The algorithm is summarized in Figure S4.

Implementation details
To optimize Q for the linear alignment from position and orientation (Algorithm S3), we use IPOpt [48], a free optimization library. As the starting point for the optimization, we use the parameters from the previous step, except for σ 2 . In our experiments, σ 2 is not optimized correctly if the initial value is placed on the wrong side of the peak of the derivative of Q (see Figure S14). To avoid this, we initialize σ 2 as σ 2 = 1 2N P (tr(X diag(P 1)X) − 2s tr((X P Ŷ ) R ) + s 2 tr(Ŷ diag(P1)Ŷ)),

S6
which would be the minimum of σ 2 if the other parameters were fixed. In our experiments this was sufficient to guarantee that the initial σ 2 is placed on the right side. Furthermore, the smaller σ becomes the sharper the peak that is the minimum of the objective function becomes (see Figure S14). This sharp peak might make it hard to estimate a reasonable step size when optimization is done numerically. We believe that this is the reason why IPOpt does not always find an optimum in later expectation maximization iterations in our implementation. However, for our test data, we found that about 4000 iterations in IPOpt were sufficient to find a new parameter set that decreased the log-likelihood.
The running time of the algorithm for elastic alignment is O(M 3 ). It can can reduced to linear by computing P with a fast Gauss transform and G by a low rank approximation as described by Myronenko and Song [28]. We did not implement this technique, which is why our implementation has cubic running time.
To compute the partial derivative with respect toR, we note thatR is parametrized via one rotation angle ρ asR = cos ρ − sin ρ sin ρ cos ρ in the 2d case. Writing the trace as a sum and taking partial derivatives with respect to this angle we find that and ∂R ∂ρ =R = − sin ρ − cos ρ cos ρ − sin ρ The second derivative can be found equivalently.
We give the partial derivatives of Equation 43 with respect to s, ρ, σ 2 and κ in Figure S15.