Cell shape characterization, alignment, and comparison using FlowShape

Abstract Motivation The shape of a cell is tightly controlled, and reflects important processes including actomyosin activity, adhesion properties, cell differentiation, and polarization. Hence, it is informative to link cell shape to genetic and other perturbations. However, most currently used cell shape descriptors capture only simple geometric features such as volume and sphericity. We propose FlowShape, a new framework to study cell shapes in a complete and generic way. Results In our framework a cell shape is represented by measuring the curvature of the shape and mapping it onto a sphere in a conformal manner. This single function on the sphere is next approximated by a series expansion: the spherical harmonics decomposition. The decomposition facilitates many analyses, including shape alignment and statistical cell shape comparison. The new tool is applied to perform a complete, generic analysis of cell shapes, using the early Caenorhabditis elegans embryo as a model case. We distinguish and characterize the cells at the seven-cell stage. Next, a filter is designed to identify protrusions on the cell shape to highlight lamellipodia in cells. Further, the framework is used to identify any shape changes following a gene knockdown of the Wnt pathway. Cells are first optimally aligned using the fast Fourier transform, followed by calculating an average shape. Shape differences between conditions are next quantified and compared to an empirical distribution. Finally, we put forward a highly performant implementation of the core algorithm, as well as routines to characterize, align and compare cell shapes, through the open-source software package FlowShape. Availability and implementation The data and code needed to recreate the results are freely available at https://doi.org/10.5281/zenodo.7778752. The most recent version of the software is maintained at https://bitbucket.org/pgmsembryogenesis/flowshape/.


Supplementary figures
Comparison of the average shape for the ABpl cell. Left: Average after alignment with maximal crosscorrelation (our method). Right: Average after alignment with the first order ellipsoid. Averaging the cell based on the ellipsoid loses characteristic features of the cell that are maintained by our alternative approach. Appendix: Technical details

Spherical harmonics
As we represent the shape as a scalar function on the sphere using ρ, we can apply concepts from spectral analysis to parameterize the shape. The spherical harmonics (SH) decomposition is the analogue of the Fourier transform, applied to a spherical domain. It decomposes a function into an infinite sum of basis functions. This decomposition can be found if ρ is square-integrable The real SH form an orthonormal basis of the space of square-integrable functions L 2 R (S 2 ). This means that every square-integrable function ρ ∈ L 2 R (S 2 ) can be decomposed as a linear combination of SH. Usually, the SH are given as complex functions, but they can be transformed into a real-valued form. For a real function, using the real SH has the added benefit that the coefficients are also guaranteed to be real.
Since ρ is a function on the sphere, we can parameterize it with spherical angular coordinates (θ, φ), with θ the colatitude (0 ≤ θ ≤ π) and φ the longitude (0 ≤ φ < 2π). The SH decomposition is where ρ m ℓ are real coefficients and Y m ℓ are the SH.

Definitions
The SH are a set of functions defined on the surface of a sphere. The complex SH Y m ℓ : S 2 → C of degree ℓ and order m are defined as N ℓm is the normalization factor P m ℓ are the associated Legendre polynomials, which can be defined in terms of the ordinary Legendre polynomials P ℓ (x) The ordinary Legendre polynomials can be expressed as The real SH Y m ℓ : S 2 → R can then be derived as follows For every degree ℓ there are 2ℓ + 1 such functions, with order m = −ℓ, . . . , +ℓ.
The orthonormality of SH means that the inner product satisfies

Implementation
Since ρ is a function on the faces, we will take the spherical coordinates of the triangle barycenters to be the sampling points The barycenter of a triangle is calculated by taking the mean of the Cartesian coordinates of the three vertices. Since this mean does not lie on the sphere, we project it back to the surface. We use the generalized iterative residual fitting (IRF) procedure, originally proposed by Chung et al. (2008) and later generalized by Elahi et al. (2017). In this method, the problem of finding the SH coefficient is efficiently solved by partitioning the problem into subspaces. Here, we use one subspace per degree ℓ. Let c l be the column vector of estimates with length 2ℓ + 1, and Y ℓ a matrix with size n × (2ℓ + 1) Set the initial residual r 0 = ρ. The first subproblem is then to findĉ 0 that minimizesĉ The estimate is Y 0ĉ0 so that the new residual is This procedure continues iteratively. At step j we findĉ j that minimizeŝ Updating the residual with At each iteration j, eq. (1) is solved using weighted least squareŝ using the diagonal mass matrix M as weights to account for unequal triangle sizes. The estimate forĉ is then just the concatenation of all theĉ j . This whole procedure is then iterated again, until the norm of the residual can not be decreased further.

Fourier methods and alignment
The convolution theorem says that we can find the convolution of two functions g and h by simply multiplying them in the Fourier domain. We will use this principle in two ways: first, to efficiently evaluate filters on spherical functions, and second, to calculate cross-correlations.

Filters
Given kernel h and a function f , the convolution is (Driscoll and Healy, 1994) Note that the kernel only has coefficients where m = 0, which is equivalent to saying it is axially symmetric. For example, a smoothing filter can be expressed by the heat kernel G k This heat kernel is the spherical analogue of the well-known Gaussian distribution (Bigot et al., 2008). The Laplacian of Gaussian is often used to filter for 2D images (Han and Uyyanonvara, 2016). It is straightforward to derive an analogous filter for shapes By doing the filtering on the spherical domain instead of the original mesh, the conformal scale factor is ignored. The result is that the filter kernel is locally scaled, but because the map is conformal, it remains isotropic.

Möbius centering
One issue with conformal maps is that they are not unique. There are many conformal maps from the unit sphere to itself, known as the Möbius transformations. General Möbius transformations on the sphere can be found as compositions of inversions and rotations. Inversions can be understood as fixing the poles that lie on some axis through the origin and then 'pushing' the geometry towards one of the poles. To obtain a unique, canonical mapping, we use the Möbius balancing algorithm described by Baden Baden et al. (2018). This algorithm finds the inversion that optimally distributes the area distortion over the sphere. This is beneficial as it facilitates using a minimal number of SH. Further, when trying to find a correspondence between two shapes, we only have to search for the optimal rotation. This is an easier problem than also having to consider the best Möbius transformation to match the shapes.

Rotations
The inner product for two functions f, g ∈ L 2 (S 2 ) is Because of the orthonormality of SH, this can also be expressed as the inner product of their SH coefficients f m ℓ and g m Given two shapes, we want to find the rotation that best aligns them. To do this we map them both to the sphere and calculate their curvature functions f and g. We then try to maximize their cross-correlation This is similar to a convolution except that the argument R is a rotation, so C is a function from the space of rotations SO(3) to R. We want to find a rotation R that maximizes this cross-correlation. Evaluating C(R) on all possible rotations would be extremely costly. Fortunately, this can be efficiently computed thanks to a Fourier transform (Kostelec and Rockmore, 2008;Baden et al., 2018;Healy et al., 2003). Here the forward transform is the SH decomposition. The inverse transform is done using the generalized Fourier transform over SO(3). The element in the Fourier domain of SO(3) is formed by taking the outer product of each subspace ℓ of the SH, forming a sequence of (2ℓ + 1) × (2ℓ + 1) matrices. Rotations are parameterized using ZYZ Euler angles (α, β, γ), with α, γ ∈ [0, 2π) and β ∈ [0, π). A rotation matrix is thus expressed as R = Rz(γ)Ry(β)Rz(α) The output of generalized Fourier transform is a 2B × 2B × 2B cube of Euler angles (α, β, γ). B = ℓmax, the band-limit and each point has a value C(R α,β,γ ).
After finding the grid point that maximizes the correlation, we do a local quadratic fit to estimate the location of the peak at sub-grid precision. At a bandwidth of B = 32, the grid is sufficiently densely sampled so that the average error is about 1°.
To evaluate rotations efficiently in the space of SH, rotations matrices are calculated using the method described by Pinchon and Hoggan (2007), as implemented by Cohen and Welling (2016).

Mean curvature and the Dirac equation
We assume a shape to be analyzed is a differentiable two-dimensional manifold M. The manifold has an immersion in R 3 f : M → R 3 .
The differential df of this map is a linear map, from the tangent space at every point to tangent vectors in R 3 . Analogous to the Jacobian matrix, it gives the best linear approximation to the surface at a point.
We also require that the surface is topologically equivalent to the 2sphere, so it is simply connected, without boundary and of genus 0.
The Bonnet problem is concerned with the existence of surfaces that have the same mean curvature at corresponding points. In general, there exist so-called Bonnet pairs. These pairs of surfaces are related by an isometry but have the same mean curvature at corresponding points . For surfaces with spherical topology, this is not the case. A proof that these can be uniquely described by their mean curvature is given by Lawson and de Azevedo Tribuzy (1981).

Quaternions
The algorithms discussed here rely on some quaternionic calculus. We refer the reader to Vicci (2001) for an overview of quaternions and only review the necessary properties. The quaternions H are an extension of the complex numbers, where three basis quaternions i, j and k take the place of the imaginary unit i. A quaternion can be represented as We identify vectors in R 3 with the imaginary quaternions So that we can write any quaternion as a scalar plus vector part: q = a+v. It is well known that quaternions can be used to represent rotations in three dimensions. To be specific, let q be a unit quaternion (|q| = 1). Then q represents a rotation of the vector u as where q * denotes the conjugate quaternion. If q has a magnitude different from 1, then q u q * represents a rotation and uniform scaling. Writing q as its scalar and vector parts, the uniform scaling factor is |q| 2 , the rotation angle is θ and the rotation axis is the normalized vectorv.

Spin transformations
The immersionf is called a spin transformation of f if there exists a smooth quaternion-valued function λ : M → H such that This relation is conformal by construction, as the tangent spaces are transformed by a local scaling and rotation as in eq. (2). The area is locally scaled by |λ * df λ| 2 |df | 2 = |λ| 4 .
We require that the surface is topologically equivalent to a sphere. This guarantees that any two conformal immersions are spin equivalent Lawson and de Azevedo Tribuzy, 1981). By the uniformization theorem, a conformal map always exists. This also holds in the discrete setting (Springborn, 2019).

Dirac equation
Not every function λ corresponds to a valid spin transformation, as the differential λ * df λ might fail to be integrable. In order for λ to be integrable it has to satisfy the Dirac equation where ρ : M → R is a scalar function and D f is the extrinsic Dirac operator, as originally described by . Solving eq. (4) by prescribing some ρ gives a new immersionf . The mean curvature functions H are related bỹ where |df | is the length element. The intrinsic Dirac operator D is defined as (Ye et al., 2018) D such that solving the intrinsic Dirac equation gives a new immersionf as above. The mean curvature is related bỹ This gives an effective algorithm for recovering a shape from a sphere: save the function ρ as in eq. (6) on the sphere, solve the intrinsic Dirac equation eq. (5), and recoverf by integrating df . This operation is conformal by construction, so it only works if the original spherical map is conformal.

Willmore energy and the power spectrum
The Willmore energy of a surface M is the squared norm of the mean curvature (Crane et al., 2013) Since we are dealing with closed surfaces of genus 1, we have W ≥ 4π. When M is exactly a sphere, W = 4π. For this reason, some authors useW = W − 4π ≥ 0 instead. This quantity can be used to measure how much a surface deviates from a perfect sphere. It is also invariant under Möbius transformations (White, 1973). W also has a physical interpretation, as for a thin elastic sheet (i.e. the cellular cortex), W is proportional to the total bending energy (Müller and Röger, 2014).
If we have the conformal map to the sphere, we want to find an expression for W in terms of the curvature function ρ. Squaring both side of eq. (6), now with |ds| as the length element on the sphere yields The area scaling cancels out, so that we can evaluate the Willmore energy directly on the sphere If we have the SH decomposition of ρ with coefficients ρ m ℓ , we can then apply Perseval's theorem So the Willmore energy is simply the sum of squares of the SH coefficients. This means we can interpret the power spectrum of ρ as the distribution of Willmore energy over each frequency ℓ This power spectrum is rotation-invariant (Kazhdan et al., 2003).

Implementation
In the discrete setting, a smooth surface is replaced by a mesh. This mesh consists of the sets V, E and F: the vertices, edges and faces respectively. We will only consider meshes where each face is a triangle, so that each face defines a plane and has a well-defined normal. Faces are oriented such that the normals consistently point outwards. The algorithm closely follows the methods proposed by Crane et al. (2011);Ye et al. (2018Ye et al. ( , 2021.

Curvature
The integrated mean curvature is defined over the edges (Hoffmann and Ye, 2020) where θ ij is the bending angle between two face normals n i and n j , and e ij is the edge shared by those faces (see fig. 7). The integrated mean Fig. 7: Bending angle θ ij between face normals n i and n j .
curvature of a face is simply the sum of its edge curvatures As this is an integrated quantity, it can be converted back to a pointwise quantity by dividing by the face area

Hyperedges and the discrete Dirac operator
We can associate with each edge the so-called hyperedge a quaternion with scalar part equal to its integrated mean curvature and imaginary part equal to its immersion in Im H. Let λ be a function λ : F → H. We then define the face-based intrinsic Dirac operator as Note that this differs from the definition given by Ye et al. (2018) since we drop the cosine factor. We can then solve the intrinsic Dirac equation for a function ρ : where D is an |F |×|F | quaternion matrix and the values of ρ are subtracted on the diagonal. We then have a discrete spin transformation (cf. eq. (3)) The edgesẽ ij can be found as the vector (imaginary) part ofẼ ij . The length element |df | in eq. (6) is discretized over every triangle as where A ′ i is the area of the corresponding triangle on the sphere.

Numerical methods
Since most numerical packages can not work directly with quaternionic matrices, we represent a matrix Q ∈ H m×n as a real block matrix Q ′ ∈ R 4m×4n . Each quaternion is replaced by a 4 × 4 block that has the same properties under ordinary matrix multiplication The conjugate quaternion is equivalent to the transpose of the matrix.
Solving eq. (7) directly has two problems. First, ρ might not correspond to a valid spin transformation at all. In that case the Dirac equation has no solution. For example, setting ρ = 0 everywhere on the sphere is clearly invalid, since it is impossible to flatten a sphere, removing all the curvature. So the constraint is relaxed to an eigenvalue problem (Crane et al., 2011) (D − ρ)λ = γλ , where γ ∈ R is the smallest eigenvalue. The Dirac equation will then be satisfied where ρ is shifted by a small constant (D − (ρ + γ))λ = 0 .
The second problem is that the solutions might not be smooth. To constrain the solution-space, it is solved in the vertices instead of the faces. A triangle mesh usually has |F | ≈ 2|V |, so there are fewer degrees of freedom. If we let A be the face-to-vertex averaging matrix, we then solve the following eigenvalue problem A T DρDρAλ = γMλ , for the smallest eigenvalue γ. M is the diagonal vertex mass matrix. The solution is averaged back onto the faces by calculating A T λ.
We can solve this problem efficiently with the inverse power method Since B is positive definite and sparse, the conjugate gradient method is an effective way to solve this system. Instead of using a random initial guess, we set λ to unity. This scheme converges very quickly, so we typically can stop after only three iterations.
After finding λ, the hyperedges are calculated by eq. (8). Integrating back to vertex positions then becomes a system of equations v i − v j =ẽ ij = ImẼ ij .
While the Dirac equation guarantees that the system is integrable, numerical errors can still be introduced. Therefore, the system might not have an exact solution. To account for this, the system is solved by least squares.