Principal Curvatures Estimation with Applications to Single Cell Data

Our method starts with Local PCA as described in [13]. Given a point cloud $x_{1},\cdots,x_{m}$, we select a neighborhood $\mathcal{N}_{x_{i},\epsilon_{\text{PCA}}}:=\{x_{j}:0<\|x_{j}-x_{i}\|<\epsilon_{\text{PCA}}\}$ around each point $x_{i}$ for a hyperparameter $\epsilon_{\text{PCA}}>0$ that has to be determined. Each data matrix containing the neighbors of a point $x_{i}$ is shifted to be centered around $x_{i}$ to get a matrix $X_{i}=\left[x_{i_{1}}-x_{i},\ldots,x_{i_{N_{i}}}-x_{i}\right]$ where $N_{i}:=|\mathcal{N}_{x_{i},\epsilon_{\text{PCA}}}|$. Then, the columns of $X_{i}$ are rescaled to $B_{i}=X_{i}D_{i}$ by applying a diagonal weighting matrix $D_{i}$ to emphasize the importance of local data. Finally, SVD decomposition yields a numerical approximation of the tangent plane. Since, we are interested mainly in surfaces, the first two eigenvectors are selected as a basis for the local tangent space and the third one serves as a normal vector to the surface.

AdaL-PCA uses the explained variance ratio for the first two singular values given by

to select a suitable parameter $\epsilon_{\text{PCA}}$. This ratio describes the fraction of data variance captured by the tangent plane approximated by the span of the first two singular vectors. We set a threshold $\gamma$ for the ratio $\rho(r)$ and compute the largest $r$-neighborhood that explains a fraction $\gamma$ of the data variance. That is,

We use a similar method to select a radius $\tau_{i}$ for estimating the curvature around each data point $x_{i}$. In this case, we need a neighborhood large enough to capture the “bending” of the surface. This is done by computing the lowest value reached by the graph of the explained variance ratio $\rho(r)$,

As illustrated in Fig.1 at a point $p$ on a torus, this approach is motivated by the fact that as the $\tau$-neighborhood increases past a certain threshold, the variance in data can no longer be explained by the selected tangent plane. We refer the reader to Algorithm 1 for a summary of AdaL-PCA’s key steps and emphasize that both $\epsilon_{\text{PCA}}$ and $\tau$ are adjusted at each data point to capture the local geometry.

The directional curvatures $\kappa_{i}(T)$ at a point $x_{i}$ in a direction $T$ are approximated (see for instance [14]) by

Here $T$ is replaced by the entries of $X_{i}$ in the proper $\tau_{i}$-neighborhood and $N$ is the orthonormal vector to the frame obtained by local PCA. The principal curvatures $\kappa_{1}$ and $\kappa_{2}$ correspond, respectively, to the highest and lowest values of the directional curvatures.¹¹1Note that this is sometimes taken as a definition of principal curvatures. In practice, we select a percentage (20%) of the highest (respectively, lowest) curvatures and average them (using a Gaussian kernel) to approximate $\kappa_{1}$ (respectively, $\kappa_{2}$). By selecting the directional vectors $T$ corresponding to the highest curvatures $\kappa_{T}$, this averaging yields principal directions, while we obtain Gaussian curvature by the product $\kappa_{1}\kappa_{2}$. The time complexity of our current implementation is $O(n_{\tau}m(m^{2}+\text{log }n_{\tau}))$, where $n_{\tau}$ is an upper bound on the cardinality of the $\tau$-neighborhoods (in general $n_{\tau}\ll m$). This can be improved significantly in practice with fast PCA algorithms for scalability [15], [16].

We compare AdaL-PCA’s estimates of Gaussian curvature against two contemporary methods, Hickok & Blumberg [10] and Diffusion Curvature [9]. We also quantify AdaL-PCA’s recovery of ground-truth mean curvature as a validation of the fidelity of its principal curvatures.

To assess the ability of various models to recover the Gaussian and mean curvatures, we generate datasets from three canonical 2-dimensional manifolds: the torus, ellipsoid, and the hyperbolic paraboloid (saddle). Tables I and II were generated from 5000 points sampled uniformly from these surfaces. To study the robustness of each method to noise, we corrupt each point:

where each $\epsilon_{1},\ldots\epsilon_{N}\overset{\mathrm{iid}}{\sim}\mathcal{N}(0,\sigma)$ and $\sigma$ between 0.0 and 0.5.

Note that both Hickok & Blumberg’s method and Diffusion Curvature require manually specified parameters, which must be tuned for each dataset. By contrast, AdaL-PCA’s heuristics adapt the method to each dataset. This results in an improved performance observed in Fig. I and Fig. II. Diffusion Curvature is an unsigned measure of local curvature for point clouds sampled from a manifold. Although it differs from Gaussian curvature, numerical experiments detailed in [9] suggest a correlation. Therefore, we report only the Pearson correlation for diffusion curvature.

We apply our model to single-cell data for cell state differentiation direction discovery. We use RNA sequencing data for human embryonic stem cells available at [17], collected over 27 days during which cells start as embryonic stem cells and then progressively differentiate into different cellular lineages. Low-dimensional manifold visualization of this data using PHATE (Fig. 4, A) shows that embryonic cells (days 0-3, displayed in red) branch into two lineages: endoderm (upper split) and ectoderm (lower split) around day 6. Further differentiation occurs during days 12-27. This is reflected in Fig. 4B with relatively constant zero curvature values at days 0-3 and a transition into a region of high variations in curvature. We observe starting from day 3 a transition into very low negative values of curvature and then a rapid progression into higher values close to zero as we approach day 27. This is consistent with the fact that the region 0-3 days corresponds to the stem state and the region 12-27 to the differentiated state. In addition to the signed curvature that provides a better appreciation of the cellular differentiation into several lineages (cell types), the principal directions in Fig. 4C, D obtained from projecting the three-dimensional principal directions using the PHATE embedding allow us to track the state towards which the cells differentiate, adding directional information.

We estimated the curvature of a publicly available single-cell induced pluripotent stem cell (iPSC) reprogramming. In this dataset, mass cytometry is used to quantitatively measure 33 protein biomarkers in 2005 mouse fibroblast cells induced to undergo reprogramming into stem cell state. Low-dimensional PHATE visualization of this data shows fibroblasts progressing to a point of divergence where two lineages emerge, one that successfully undergoes reprogramming and another that undergoes apoptosis (cell death). Our model correctly identifies the initial branching point as having negative values of Gaussian curvature indicating saddle-like divergent paths out of the branching point (Fig. 5). Moreover, the principal directions on the diverging branch correctly identify the directions in which the cell lineages diverge.