Moment-based metrics for molecules computable from cryogenic electron microscopy images

Abstract

Single-particle cryogenic electron microscopy (cryo-EM) is an imaging technique capable of recovering the high-resolution three-dimensional (3D) structure of biological macromolecules from many noisy and randomly oriented projection images. One notable approach to 3D reconstruction, known as Kam’s method, relies on the moments of the two-dimensional (2D) images. Inspired by Kam’s method, we introduce a rotationally invariant metric between two molecular structures, which does not require 3D alignment. Further, we introduce a metric between a stack of projection images and a molecular structure, which is invariant to rotations and reflections and does not require performing 3D reconstruction. Additionally, the latter metric does not assume a uniform distribution of viewing angles. We demonstrate the uses of the new metrics on synthetic and experimental datasets, highlighting their ability to measure structural similarity.

Impact Statement

Single-particle cryogenic electron microscopy (cryo-EM) is a popular method to obtain three-dimensional (3D) reconstructions of biological molecules from noisy two-dimensional (2D) tomographic projection images. Many iterative techniques for this reconstruction require initializations sufficiently close to the unknown structure to obtain high-quality reconstructions. To help select an initialization from a database of known structures, this paper introduces a metric to compare the similarity of known 3D structures with a stack of noisy 2D tomographic projection images of an unknown structure. We show that this metric distinguishes differing structures and present an efficient method to compute it, notably without performing 3D reconstruction.

1. Introduction

Single-particle cryogenic electron microscopy (cryo-EM) enables high-resolution reconstruction of three-dimensional (3D) biological macromolecules from a large collection of noisy and randomly oriented projection images. Kam’s method⁽¹⁾ is one of the earliest methods proposed for homogeneous reconstruction in cryo-EM. It is a statistical method-of-moments approach applied to the cryo-EM reconstruction problem, where the computation of low-order statistics of two-dimensional (2D) images allows for the estimation of 3D structure through solving a polynomial system. Kam’s method has helped push the theoretical understanding of the reconstruction process – under certain conditions, it is a provable algorithm and provides bounds for the estimated structure’s quality in terms of the noise level and the number of images.^{( 2 – 8 )} Kam’s method also enjoys other remarkable properties:

1. 1. It bypasses the need for angular assignment, typically a large computational burden in competing methods.

2. 2. It is a streaming algorithm and is thus theoretically much faster than iterative methods.

3. 3. It can – in theory – break the detection limit of the minimal size of proteins that can be observed in cryo-EM.^{( 9 )}

While theoretically attractive, Kam’s method has not been able to yield high-resolution reconstructions as yet. One direction that is currently being explored to resolve this issue is the development of better priors, for instance, based on the sparsity of the solution as in the study by Bendory et al.^{( 7 )} Separately, there has been considerable, continued interest in utilizing the vast amount of solved structures stored in the Protein Data Bank (PDB)^{( 10 )} to improve cryo-EM reconstructions.

Leveraging the PDB as a prior, we propose a method to match either projection images or molecular volumes to a database of previously solved structures (Section 3). We use this procedure as a rotationally and reflectionally invariant metric that can be directly computed from raw image datasets without needing a 3D reconstruction process. Importantly, our metric neither relies on prior knowledge of rotations nor assumes a uniform rotational distribution, making it applicable to typical datasets.

To demonstrate the efficacy of our metric, we compare it to existing methods and show empirically that it achieves similar performance to alignment-based metrics. As an application, we use our metric to compute a low-dimensional embedding of a subset of the PDB into the Euclidean plane, visually showcasing how structurally similar proteins are embedded near each other (Section 4.2). Further, we apply the version of the metric that can be directly computed from stacks of 2D images and show that it gives an efficient methodology to identify the nearest neighbors in a database to a given set of experimental moments on synthetic and real datasets (Sections 4.3 and 4.4).

2. Background

This section presents the mathematical preliminaries needed to define our metric. Let $ \Phi :{\mathrm{\mathbb{R}}}^3\to \mathrm{\mathbb{R}} $ be the electrostatic potential of a molecule and $ \hat{\Phi}:{\mathrm{\mathbb{R}}}^3\to \mathrm{\mathbb{C}} $ be its Fourier transform, which we define by

$$ \hat{\Phi}\left(\xi \right)={\int}_{{\mathrm{\mathbb{R}}}^3}\Phi (x){e}^{- i\xi \cdot x} dx. $$

A single projection image is given by

$$ I=h\hskip0.35em \ast \hskip0.35em PR\Phi +\varepsilon, $$

and its Fourier transform is

$$ \hat{I}=H\cdot SR\hat{\Phi}+\hat{\varepsilon}, $$

where $ P $ is the projection operator, $ S $ is the slicing operator, $ h $ is a point spread function, $ H $ is the contrast transfer function (CTF), $ \varepsilon $ is noise, and $ R\in \mathrm{SO}(3) $ is a rotation operator. We assume that the Fourier transform $ \hat{\Phi} $ can be expanded in a spherical harmonic expansion:

(1) $$ \hat{\Phi}\left(r,\theta, \varphi \right)=\sum \limits_{\mathrm{\ell}=0}^L\sum \limits_{m=-\mathrm{\ell}}^{\mathrm{\ell}}{A}_{\mathrm{\ell},m}(r){Y}_{\mathrm{\ell}}^m\left(\theta, \varphi \right), $$

where $ \left(r,\theta, \phi \right) $ are spherical coordinates, and $ {Y}_{\mathrm{\ell}}^m $ denotes the complex spherical harmonic:

$$ {Y}_{\mathrm{\ell}}^m\left(\theta, \varphi \right):= {\left(\frac{\left(\mathrm{\ell}-m\right)!\left(2\mathrm{\ell}+1\right)}{4\pi \left(\mathrm{\ell}+m\right)!}\right)}^{1/2}{e}^{im\varphi}{P}_{\mathrm{\ell}}^m\left(\cos \theta \right), $$

where $ {P}_{\mathrm{\ell}}^m $ are the associated Legendre polynomials, $ {A}_{\mathrm{\ell},m}(r) $ are $ r $ -dependent coefficients, and $ L $ is a bandlimit parameter. See Eq. 14.30.1 in^{( 11 )} for the definitions of $ {Y}_{\mathrm{\ell}}^m $ and $ {P}_{\mathrm{\ell}}^m $ .

Let $ \rho :\mathrm{SO}(3)\to \mathrm{\mathbb{R}} $ be the probability density function of the rotational distribution, which without loss of generality is invariant to in-plane rotations and reflections. (Note that by augmenting the dataset with in-plane rotations and reflections of all 2D images, one can always reduce to the case of such an invariant distribution $ \rho $ ; for example, see the study by Ponce and Singer^{( 12 )}). More formally, $ \rho $ is a function on $ SO(3)/O(2)\simeq {\unicode{x1D54A}}^2/\left\{\pm 1\right\} $ _, which is formed by identifying antipodal points on the sphere $ {\unicode{x1D54A}}^2 $ _.^{( 13 )} Thus, we model the density as a function $ \rho :\mathrm{SO}(3)\to \mathrm{\mathbb{R}} $ with an even-degree spherical harmonic expansion:

(2) $$ \rho (R)=\sum \limits_{\mathrm{\ell}=0}^P\sum \limits_{m=-2\mathrm{\ell}}^{2\mathrm{\ell}}{B}_{2\mathrm{\ell},m}\overline{Y_{2\mathrm{\ell}}^m}\left(\theta (R),\varphi (R)\right), $$

where $ \left(\theta (R),\varphi (R)\right) $ represents the third column of the rotation matrix given by $ R $ in spherical coordinates, and $ P\in {\mathrm{\mathbb{Z}}}_{\ge 0} $ is a bandlimit parameter (see Section 4.2 in the study by Sharon et al.^{( 8 )}). The analytical first and second moments $ {\mathbf{m}}_1 $ and $ {\mathbf{m}}_2 $ of the Fourier-transformed projection images with respect to $ \rho $ are

(3) $$ {\displaystyle \begin{array}{c}{\mathbf{m}}_1\left(x,y\right)=\int \left(R\cdot \hat{\Phi}\right)\left(x,y,0\right)\rho (R) dR,\hskip2em \mathrm{and}\\ {}{\mathbf{m}}_2\left(\left(x,y\right),\Big({x}^{\prime },{y}^{\prime}\Big)\right)=\int \left(R\cdot \hat{\Phi}\right)\left(x,y,0\right)\left(R\cdot \hat{\Phi}\right)\left({x}^{\prime },{y}^{\prime },0\right)\rho (R) dR,\end{array}} $$

where $ dR $ denotes integration with respect to the Haar measure on $ SO(3) $ . It will be convenient to write $ \left(x,y\right) $ and $ \left({x}^{\prime },{y}^{\prime}\right) $ in terms of polar coordinates $ \left(r,\phi \right) $ and $ \left({r}^{\prime },{\phi}^{\prime}\right) $ , respectively. In Appendix A.1, we show in (12) and (13) that the first moment only depends on $ r $ , that is, $ {\mathbf{m}}_1={\mathbf{m}}_1(r):{\mathrm{\mathbb{R}}}_{\ge 0}\to \mathrm{\mathbb{C}} $ , and that the second moment only depends on $ r,{r}^{\prime } $ and $ \Delta \phi =\phi -{\phi}^{\prime } $ , that is, $ {\mathbf{m}}_2={\mathbf{m}}_2\left(r,{r}^{\prime },\Delta \phi \right):{\mathrm{\mathbb{R}}}_{\ge 0}\times {\mathrm{\mathbb{R}}}_{\ge 0}\times \left[-2\pi, 2\pi \right]\to \mathrm{\mathbb{C}} $ . We write $ {\mathbf{m}}_1={\mathbf{m}}_1\left(\hat{\Phi},\rho \right) $ and $ {\mathbf{m}}_2={\mathbf{m}}_2\left(\hat{\Phi},\rho \right) $ to denote the first and second moments defined by $ \hat{\Phi} $ and $ \rho $ when discussing multiple structures. The basis of Kam’s method is that the moments in (3) can be estimated from images and related to expansion coefficients for the potential $ \hat{\Phi} $ ; see Appendix A for explicit formulas.

3. Definition of Kam’s metrics

We now use metrics between the moments in (3) to define similarity between proteins and stacks of images. A first function is used to measure the similarity of two known structures by the moments of their potential as defined in (3). The second is used to measure the similarity between a known structure and the unknown structure observed in a dataset of images.

Crucially, the metrics can be computed without performing 3D alignment of the structures, reducing their computational costs compared to other approaches. Moreover, one of the metrics can be directly computed from noisy and CTF-affected projection images. This enables a nearest neighbor search among known structures to determine an initialization for the 3D reconstruction pipeline, especially in the expectation–maximization procedure.^{( 14 , 15 )}

3.1. Kam’s volume metric $ {d}_{\mathrm{vKam}} $

Here, we introduce the first of Kam’s metrics, which measures the similarity of two 3D structures. We use this to perform dimensionality reduction to visualize the relationship between structures from a subset of the PDB.

In detail, given two 3D structures $ {\Phi}_1 $ and $ {\Phi}_2 $ , we define the distance between them through their first and second moments $ {\mathbf{m}}_1 $ and $ {\mathbf{m}}_2 $ under a uniform distribution of viewing directions, which we denote by $ \rho ={\rho}_u $ . We will derive the explicit equations for the uniform case in (17) and (18). We then measure the resulting weighted deviation of the first and second moments by

(4) $$ \parallel {\mathbf{m}}_2\left({\hat{\Phi}}_1,{\rho}_u\right)-{\mathbf{m}}_2\left({\hat{\Phi}}_2,{\rho}_u\right){\parallel}_{L^2\left({\mathrm{\mathbb{R}}}^2\times {\mathrm{\mathbb{R}}}^2\right)}^2+\lambda \parallel {\mathbf{m}}_1\left({\hat{\Phi}}_1,{\rho}_u\right)-{\mathbf{m}}_1\left({\hat{\Phi}}_2,{\rho}_u\right){\parallel}_{L^2\left({\mathrm{\mathbb{R}}}^2\right)}^2, $$

where $ \lambda \ge 0 $ is a hyperparameter that we set to 1 for all experiments. The moments will be represented on a discretized voxel grid, and we therefore replace the continuous norms with discrete norms. More specifically, we will represent the second moment using a grid $ r,{r}^{\prime}\in \left\{{r}_1,\dots, {r}_{\left\lfloor N/2\right\rfloor}\right\} $ and , where $ N $ is the number of pixels of one side of the discretized volume. We define the grid points $ {r}_k=\delta k/N $ , $ \Delta {\phi}_j=2\pi j/N-2\pi $ , for $ k=1,\dots, \left\lfloor N/2\right\rfloor $ , and $ j=1,\dots, 2N $ , where $ \delta $ is the side length of the volume grid in angstroms. We then use the following two approximations to the continuous norms above

(5) $$ \parallel \mathbf{M}{\parallel}_{w_2}^2:= \sum \limits_{\begin{array}{c}j=1,\dots, 2N\\ {}{k}_1=1,\dots, \left\lfloor N/2\right\rfloor \\ {}{k}_2=1,\dots, \left\lfloor N/2\right\rfloor \end{array}}{\left|\mathbf{M}\left(\Delta {\phi}_j,{r}_{k_1},{r}_{k_2}\right)\right|}^2{r}_{k_1}{r}_{k_2},\hskip1em \parallel \mathbf{N}{\parallel}_{w_1}^2:= \sum \limits_{k=1,\dots, \left\lfloor N/2\right\rfloor }{\left|\mathbf{N}\left({r}_k\right)\right|}^2{r}_k. $$

With these norms, we define the metric comparing two sets of moments of two 3D structures by

(6) $$ {d}_{\mathrm{vKam}}\left({\Phi}_1,{\Phi}_2\right)\hskip0.35em := {\left(\parallel {\mathbf{m}}_2\left({\hat{\Phi}}_1,{\rho}_u\right)-{\mathbf{m}}_2\Big({\hat{\Phi}}_2,{\rho}_u\left){\parallel}_{w_2}^2+\lambda \parallel {\mathbf{m}}_1\right({\hat{\Phi}}_1,{\rho}_u\left)-{\mathbf{m}}_1\right({\hat{\Phi}}_2,{\rho}_u\Big){\parallel}_{w_1}^2\right)}^{1/2}, $$

This distance is rotationally invariant since for any rotation $ R $ , we have $ \hat{R\cdot \Phi}=R\cdot \hat{\Phi} $ and the moments $ {\mathbf{m}}_1 $ and $ {\mathbf{m}}_2 $ in (2) satisfy

(7) $$ {\mathbf{m}}_i(R\cdot \hat{\Phi},\rho )={\mathbf{m}}_i(\hat{\Phi},{R}^T\cdot \rho ), $$

as can be seen through a change of variables in (3). When $ \rho ={\rho}_u $ is uniform, clearly $ {R}^T\cdot \rho =\rho $ , which therefore shows rotational invariance of the cost function in (4), up to the discretization of the volume grid. Note that this bypasses the need for an alignment step. We detail the procedure for computing $ {\mathbf{m}}_1 $ , $ {\mathbf{m}}_2 $ and therefore $ {d}_{\mathrm{vKam}} $ in Appendix A.1. Under certain conditions, it has been demonstrated that the second moment of the image collection identifies the 3D structure uniquely^{( 2 – 4 , 6 , 7 )} or up to a finite list of candidate structures.^{( 8 )} In Section 4.2, we show that our metric is alike other similarity scores but remarkably does not rely on alignment.

3.2. Kam’s image metric $ {d}_{\mathrm{iKam}} $

We now introduce a metric between the empirical moments computed from a set of experimental projection images and the moments computed from the atomic coordinates of a known structure that compares images to the known structure. We detail the procedure for computing these moments in Appendix A.1.

If the distribution of poses in the experimental dataset would be known to be uniform, the empirical moments could directly be substituted for $ {\mathbf{m}}_1 $ and $ {\mathbf{m}}_2 $ in (6) and the metric could be defined as the deviation between the moments of the two structures. In practice, however, the distribution of angles is not uniform and is unknown. Since the moments are functions of this distribution, it must therefore be inferred.

We will show in (12) and (13) that $ {\mathbf{m}}_1 $ and $ {\mathbf{m}}_2 $ depend linearly on the expansion coefficients $ {B}_{p,u} $ of the distribution of viewing directions. The optimization problem minimizing the discrepancy between the moments of the two structures is, therefore, a linear least-squares problem in $ {B}_{p,u} $ . It follows from Table 3 of^{( 8 )} that this linear least-squares is Zariski-generically full-rank (although not necessarily well-conditioned) for various small bandlimits $ L $ and $ P. $ Solving this optimization problem efficiently eliminates the unknown rotational distribution. We then define the metric between the moments of the structure $ \Phi $ and the experimental moments $ {\tilde{\mathbf{m}}}_1,{\tilde{\mathbf{m}}}_2 $ by

(8) $$ \underset{\rho \in \mathcal{P}}{\min}\parallel {\tilde{\mathbf{m}}}_2-{\mathbf{m}}_2\left(\hat{\Phi},\rho \right){\parallel}_{L^2\left({\mathrm{\mathbb{R}}}^2\times {\mathrm{\mathbb{R}}}^2\right)}^2+\lambda \parallel {\tilde{\mathbf{m}}}_1-{\mathbf{m}}_1\left(\hat{\Phi},\rho \right){\parallel}_{L^2\left({\mathrm{\mathbb{R}}}^2\right)}^2, $$

where $ \lambda \ge 0 $ is a hyperparameter which we set to 1 for all experiments and

(9) $$ \mathcal{P}=\left\{\rho (R):\mathrm{SO}(3)\to \mathrm{\mathbb{R}}\;\mathrm{s}.\mathrm{t}.\rho (R)=\sum \limits_{\mathrm{\ell}=0}^P\sum \limits_{m=-2\mathrm{\ell}}^{2\mathrm{\ell}}{B}_{2\mathrm{\ell},m}\overline{Y_{2\mathrm{\ell}}^m}\Big(\theta (R),\varphi (R)\Big),\mathrm{and}\;{\int}_{\mathrm{SO}(3)}\rho (R) dR=1\right\}, $$

is the set of admissible distributions of viewing directions that are invariant to global reflections and in-plane rotations, where $ \left(\theta (R),\varphi (R)\right) $ is as in (2). To simplify the optimization problem and lead to faster algorithms, note that we do not impose positivity of the distributions $ \rho \in \mathcal{P} $ , though this could be enforced, for instance, by imposing linear constraints $ \rho \left({R}_i\right)\ge 0 $ for a suitable choice of $ {R}_i\in \mathrm{SO}(3) $ . Moreover, the constraint $ {\int}_{\mathrm{SO}(3)}\rho (R) dR=1 $ is equivalent to imposing $ {B}_{0,0}=1 $ _,^{( 8 )} which can be achieved by removing $ {B}_{0,0} $ from the set of optimization variables and fixing its value to $ 1 $ . The values of the bandlimit parameters $ L,P $ and the hyperparameter $ \lambda $ used in our numerical experiments are given in Appendix B.1.

Just as in the previous section, we replace the continuous norms in (8) by discrete norms to define the metric between empirical moments and the moments from a 3D structure as

(10) $$ {d}_{\mathrm{iKam}}\left(\left({\tilde{\mathbf{m}}}_1,{\tilde{\mathbf{m}}}_2\right),\Phi \right):= {\left(\underset{\rho \in \mathcal{P}}{\min}\parallel {\tilde{\mathbf{m}}}_2-{\mathbf{m}}_2\left(\hat{\Phi},\rho \right){\parallel}_{w_2}^2+\lambda \parallel {\tilde{\mathbf{m}}}_1-{\mathbf{m}}_1\Big(\hat{\Phi},\rho \Big){\parallel}_{w_1}^2\right)}^{1/2}. $$

The cost function in (8) is rotationally invariant, in that it does not depend on the orientation of $ \Phi $ , since (7) implies that

(11) $$ {\displaystyle \begin{array}{l}\underset{\rho \in \mathcal{P}}{\min}\hskip0.3em \parallel {\tilde{\mathbf{m}}}_2-{\mathbf{m}}_2\left(\hat{R\cdot \Phi},\rho \right){\parallel}_{L^2\left({\mathrm{\mathbb{R}}}^2\times {\mathrm{\mathbb{R}}}^2\right)}^2+\lambda \parallel {\tilde{\mathbf{m}}}_1-{\mathbf{m}}_1\left(\hat{R\cdot \Phi},\rho \right){\parallel}_{L^2\left({\mathrm{\mathbb{R}}}^2\right)}^2\\ {}=\underset{\rho \in \mathcal{P}}{\min}\hskip0.3em \parallel {\tilde{\mathbf{m}}}_2-{\mathbf{m}}_2\left(\hat{\Phi},{R}^T\cdot \rho \right){\parallel}_{L^2\left({\mathrm{\mathbb{R}}}^2\times {\mathrm{\mathbb{R}}}^2\right)}^2+\lambda \parallel {\tilde{\mathbf{m}}}_1-{\mathbf{m}}_1\left(\hat{\Phi},{R}^T\cdot \rho \right){\parallel}_{L^2\left({\mathrm{\mathbb{R}}}^2\right)}^2\\ {}=\underset{\rho \in \mathcal{P}}{\min}\hskip0.3em \parallel {\tilde{\mathbf{m}}}_2-{\mathbf{m}}_2\left(\hat{\Phi},\rho \right){\parallel}_{L^2\left({\mathrm{\mathbb{R}}}^2\times {\mathrm{\mathbb{R}}}^2\right)}^2+\lambda \parallel {\tilde{\mathbf{m}}}_1-{\mathbf{m}}_1\left(\hat{\Phi},\rho \right){\parallel}_{L^2\left({\mathrm{\mathbb{R}}}^2\right)}^2.\end{array}} $$

where the last equality follows because $ {R}^T\cdot \rho $ lies in $ \mathcal{P} $ , since rotating a viewing angle distribution over $ \mathrm{SO}(3) $ results in another viewing angle distribution over $ \mathrm{SO}(3) $ .

At the cost of solving the small linear system detailed in Appendix A.3, our method allows for the comparison between a stack of images and a resolved structure, without performing a 3D reconstruction. Furthermore, we precompute the least-squares matrices necessary for optimization, after which the distance function can be calculated in real time. With sufficient storage and precomputation, the procedure is scalable to the entirety of the PDB.

In particular, $ {d}_{\mathrm{iKam}} $ can be used in an efficient scheme to match a stack of synthetic images to the potentials of nearby PDB structures. By selecting a subset of the PDB database, one can efficiently compute $ {d}_{\mathrm{iKam}}\left(\left({\tilde{\mathbf{m}}}_1,{\tilde{\mathbf{m}}}_2\right),\Phi \right) $ for each $ \Phi $ in the subset and find the nearest neighbors. The method for processing image moments in practice is detailed in Appendix A.4, and the computational complexity of the metric is derived in Appendix A.5.

4. Results

4.1. Existing measures of structural similarity

There are several existing methods for reporting structural similarity between two known volumes. We list two approaches based on computing alignment and Zernike moments. We compare both $ {d}_{\mathrm{iKam}} $ and $ {d}_{\mathrm{vKam}} $ to these approaches in the experiments in the following subsections. Note that the following existing metrics are limited to measuring similarity between two structures and cannot compare images to structures, whereas $ {d}_{\mathrm{iKam}} $ can.

1. 1. Euclidean alignment: A classical approach for comparing the similarity of two structures is to sample the volumes on a 3D grid and calculate the Euclidean distance between pairs over rotations and translations. However, this method is expensive to compute since optimization over $ \mathrm{SO}(3) $ is required to align the structures. Accelerated methods for computing these alignments by maximizing the correlation between two volume maps over rotations and translations have been implemented in various programs, for example, via gradient ascent in Chimera.^{( 16 )} Further acceleration can be achieved by calculating volumetric correlations by expanding the volumes in a well-chosen basis and applying dimensionality reduction^{( 17 )} or by maximizing the correlation between common lines in projection images generated from the volumes.^{( 18 )} Similar alignment methods, such as those described in the study by Bartesaghi et al. and Xu et al.^{( 19 , 20 )}, are also used in electron tomography for sub-volume similarity. In this paper, we use a Bayesian optimization algorithm to minimize an Euclidean loss function, as described in the study by Singer and Yang,^{( 21 )} to compute the alignment and minimum distance between two volumes.

2. 2. Zernike moments: Another metric for structural similarity is to expand the molecule’s structure in Zernike polynomials and compute a metric from the Zernike expansion coefficients, as described in the study by Guzenko et al.^{( 22 )}, which is used by the PDB for structural similarity search.

4.2. Applying $ {d}_{\mathrm{vKam}} $ to a PDB subset

To test the ability of $ {d}_{\mathrm{vKam}} $ to discern the similarity between 3D structures, we first generate a database using 1420 structures downloaded from the PDB.^{( 10 )} The subset chosen here was selected by filtering for human proteins with an experimental structure at resolution between 1 and 3 Å and a molecular weight between 150 and 250 kDa. We use this subset because it encompasses a diverse range of shapes and symmetries as well as many homologous structures. Additionally, the weight range reflects a smaller and more challenging protein size for a typical cryo-EM experiment.^{( 23 )} In the future, a larger database containing the entire PDB can also be generated.

Using our database, we first generate a discretized potential for each structure as described in Appendix A.2. The first and second moments of each structure can then be computed using (3). We then compute $ {d}_{\mathrm{vKam}} $ in (6) pairwise for all structures in the database.

To compare the performance of $ {d}_{\mathrm{vKam}} $ against existing metrics, we calculate pairwise scores using $ {d}_{\mathrm{vKam}} $ , Euclidean alignment, and the Zernike metric. We then plot the returned scores against each other and calculate a ranking similarity using normalized discounted cumulative gain^{( 24 )} (NDCG). We use this metric since it is a popular method to quantify the similarity between sets of rankings; its calculation is given in Appendix A.6.

In Figure 1, we report the NDCG scores between pairs of metrics. All NDCG scores are close to 1, indicating strong agreement among the three different metrics on which structures are most similar. However, the alignment metric and $ \log \left({d}_{\mathrm{vKam}}\right) $ share the highest average NDCG score. To verify the statistical significance of this agreement, we report a t-test by selecting $ 10 $ different subsets, showing that the NDCG score between $ {d}_{\mathrm{vKam}} $ and the alignment metric is statistically significantly higher (with a p-value $ p\approx 8\times {10}^{-9} $ ) than the NDCG score between the alignment metric and the Zernike metric. We thus conclude that $ {d}_{\mathrm{vKam}} $ provides a fast and accurate alternative for the alignment metric.

Figure 1. Comparison between $ {d}_{\mathrm{vKam}} $ , the Zernike metric, and Euclidean alignment. a-c) A random size 100 subset of the database is selected. Then, pairwise similarity metrics are calculated and plotted, where each point represents a pair of structures. The NDCG score is calculated using the metric on the $ y $ -axis as the predicted metric, and the metric on the $ x $ -axis as the true metric. d) The procedure is repeated with 10 randomly selected size 100 subsets, and the mean ( $ \mu $ ) and standard deviation ( $ \sigma $ ) of the NDCG scores are calculated. The error bars and points visualize $ \mu \pm \sigma $ .

Although it is the most interpretable metric, Euclidean alignment is computationally expensive to execute for all pairs of structures in a database. To achieve a manageable runtime for alignment, we calculate pairwise Euclidean alignment distances for a subset of the database of size 100. Pairwise alignment on this subset took 8 hours on a 2.6 GHz Intel Skylake Central Processing Unit (CPU) running AVX-512 using 16 physical cores and 80 GB random-access memory (RAM). To do pairwise alignment via Bayesian optimization for the entire database of 1420 structures would require 46 days of computation, whereas using $ {d}_{\mathrm{vKam}} $ (including precomputation) to calculate pairwise distances between all 1420 structures in the database requires 3 minutes on the same hardware. Despite containing an alignment component, the Zernike metric is also fast, taking 3 minutes to compute pairwise distances for the entire database.

After observing high agreement between $ {d}_{\mathrm{vKam}} $ and the other metrics, we compute a 2D embedding of the similarity between structures in our database using t-distributed stochastic neighbor embedding (t-SNE)^{( 25 )} (see Figure 2). Analogous t-SNE plots for the alignment metric and Zernike metric are reported in Appendix B.2. We find that $ {d}_{\mathrm{vKam}} $ provides interpretable results in identifying similar molecules from their moments without the need for alignment. In particular, we observe that both homologous (i.e., structures with similar sequences) and similar-shaped structures are shown to be clustered together.

Figure 2. 2D embedding of protein structures based on their similarity using $ {d}_{\mathrm{vKam}} $ . The analytical moments of 1420 proteins were computed and compared using (6), and t-SNE was applied for visualization. Each node represents a single structure and is colored by the number of atoms. Distinct clusters containing homologous or similarly shaped structures suggest that $ {d}_{\mathrm{vKam}} $ provides interpretable results.

4.3. Database search using $ {d}_{\mathrm{iKam}} $ with synthetic cryo-EM data

We next demonstrate the ability of $ {d}_{\mathrm{iKam}} $ to accurately find a match for the moments computed from projection images to a database of analytical moments computed from the atomic coordinates of known structures. To test our metric, we use the same dataset as the previous section, selecting the protein structure of a Mas-related G-protein-coupled receptor (available as entry PDB-7VV3^{( 26 )}) from our database described in Section 4.2. We use this entry because our database includes several similarly shaped yet nonidentical structures, on which we examine our metric’s performance.

We generate a synthetic cryo-EM dataset as illustrated in Figure 3: We take $ 25000 $ clean projection images from a nonuniform distribution over $ \mathrm{SO}(3) $ at viewing angles given by a mixture of three von Mises–Fisher distributions.^{( 27 )} To simulate cryo-EM data, the images are then corrupted with one of 100 unique radial CTFs, after which we add white noise with a signal-to-noise ratio (SNR) of $ 0.1 $ . We define the SNR by taking the signal as the average squared intensity over each pixel in all the clean images and setting the noise variance to the appropriate ratio of the signal. These simulated images are generated using the ASPIRE software package^{( 28 )} and have parameters consistent with many experimental datasets.

Figure 3. Visualization of the generation of simulated images. (a) Protein structure of PDB-7VV3. (b) Clean projection images from PDB-7VV3 generated with a nonuniform viewing angle distribution. (c) Projection images corrupted with a CTF and white noise with $ SNR=0.1 $ . (d) Distribution of nonuniform viewing angles.

We then compute the moments of the simulated images as will be shown in (12) and (13) and compare them to the database of moments using the image-to-volume metric described in (8). We also report the effect of varying the number of images on the metric’s performance in Appendix B.3. Using our metric, we can rank the similarity of the image’s moments to our database as shown in Figure 4. We show that the most similar score (i.e., the smallest value in image Kam’s metric) corresponds to the ground truth structure used to generate the images. Furthermore, based on our results, the next top 116 structures correspond to structures with similar volumes and sequences. These results demonstrate that we are able to compare directly between noisy, CTF-corrupted images and known structures. This approach could be especially valuable if there is no known model for initialization in 3D reconstruction or if the molecule generating the images is unknown.^{( 29 )}

Figure 4. Histogram ranking of dissimilarities computed using $ {d}_{\mathrm{iKam}} $ on simulated noisy projection images generated from PDB-7VV3.

We report alignment scores between molecules in our database to PDB entry 7VV3, compare these to our metric’s scores, and plot the results in Figure 5. Most notably, when the protein structure becomes less similar to the ground truth (7VV3), the alignment metric begins to lose discriminative power. Figure 5 shows structures with varying degrees of dissimilarity as having the same score ( $ \sim $ 100). In contrast, our metric retains discriminative power, ranking structures with similar sequences/functions before structures with similar shapes.

Figure 5. Comparison between the rankings given by $ {d}_{\mathrm{iKam}} $ (computed from simulated images) and the minimum Euclidean distance after alignment (computed from volumes). The structures shown are superimposed with the ground truth after alignment in panels (a)–(d). The points on the graph that correspond to these structures are colored and labeled. The ground truth corresponds to the green cross in the lower left.

Note that the Euclidean alignment metric shows stagnation whereas Kam’s metric does not.

Alignment via Bayesian optimization between one structure and the 1420 structures in the database took 95 minutes using the hardware described in Section 4.2. Aside from the computational cost, the interpretation of the optimal rotation returned by alignment becomes unclear when comparing two structures that are not volumetrically similar. On the other hand, our metric does not return an alignment between two structures, which could render it less useful when an explicit alignment must be computed. Without this alignment, it may become harder to visually compare their volumes.

It is computationally costly to generate and perform moment estimation on synthetic images for every molecule in the database. As such, to compare the performance of our metric against the Zernike metric, we select from our database a random subset of 100 structures. For each structure, we repeat the process we perform on PDB-7VV3: First, we generate a nonuniform distribution over $ {\unicode{x1D54A}}^2 $ as a mixture of three von Mises–Fisher distributions with random means, weights, and covariance matrices. We then generate 25000 images, corrupt with SNR = 0.1 and radial CTFs, compute the moments, and search across the database.

For every structure, we recover the ground truth as one of the first six lowest-scoring molecules. Moreover, 88 of the 100 tests recovered the ground truth as the lowest-scoring molecule. To evaluate how well the metrics agree on structural similarity, we compute the size of the intersection between the top ten structures returned by our metric and those returned by the Zernike metric. As shown in Figure 6, we find that the metrics agree on two to three structures, and a large number of structures agree only on the ground truth structure. When they occur, disagreements between the metrics are likely due to the presence of near-identical molecules in the database.

Figure 6. Comparison between $ {d}_{\mathrm{iKam}} $ _, computed from simulated images, and the Zernike metric, computed from volumes. Here, we repeat simulated experiments 100 times. Then, the size of the intersection of the top ten structures returned by $ {d}_{\mathrm{iKam}} $ and the Zernike metric is plotted as a histogram.

4.4. Toward matching experimental datasets by $ {d}_{\mathrm{iKam}} $

While our simulated result shows success in matching a synthetic cryo-EM dataset to PDB structures, many experimental cryo-EM datasets are corrupted by a large number of unmodeled effects that we have not considered. Among the real-data effects are scattering potential’s corruption by a solvent effect,^{( 30 )} the B-factor,^{( 31 )} a global scaling ambiguity, imperfect centering, junk particles, non-radial CTF, and imperfect noise model. Our simulation falls short on these counts.

In a first step toward applying $ {d}_{\mathrm{iKam}} $ to real experimental datasets, we compare the moments of a stack of images deposited in the Electron Microscopy Public Image ARchive (EMPIAR)^{( 32 )} to the moments of its preprocessed 3D reconstructions given by the program CryoSPARC^{( 33 )}. We select the dataset EMPIAR-10076,^{( 34 )} a heterogeneous dataset containing five major structures. The dataset is well characterized, and each image in the dataset has been classified into one of the five major states^{( 34 )} or “junk” particles, which we discard. We use the classification to generate five separate datasets, allowing us to compute five different moments, one for each of the major states. This test case allows us to examine our metric matching on a real dataset, while bypassing some of the issues associated with comparing datasets and volumes obtained in different experimental conditions.

We downsample the image stack to $ 64\times 64 $ , center using the deposited shift, and mask the images with a circular binary mask of radius $ 0.8 $ times half the side length of the image. We then estimate the moments for each structure and compare them to moments computed analytically from preprocessed volume reconstructions of the five major structures, as well as two other distinctly shaped ribosomes from the Electron Microscopy Data Bank^{( 35 )} (EMDB), EMD-8457 and EMD-2660, used as a baseline. Scaling issues between the moment computed from the images and the moment computed from the volume are resolved by examining the diagonal entries of the second moments. Specifically, we find a multiplicative scaling factor that best matches the diagonal of the image-computed second moment and those of the volume-computed second moment under a uniform distribution with respect to the $ {l}^2 $ norm.

As shown in Figure 7, it is observed that Kam’s metric recovers the ground truth structure at the lowest distance for the experimental images corresponding to structure 001. We note that the scores for molecules 001 and 002, as well as molecules 003 and 004, are almost identical in value. Also, we find that the analytical moments are closer to each other than to the experimentally determined moments. Finally, the metric reports the baseline structures, which are very different in shape and size, at the largest distances.

Figure 7. $ {d}_{\mathrm{iKam}} $ visualization and ranking results for experimental data corresponding to structure 001 (a) Experimental images from EMPIAR-10076 corresponding to structure 001 downsampled to $ 64\times 64 $ pixels, centered, and with binary mask applied. (b) Comparison between diagonal entries of the second moment computed from the reconstructed volumes 001 and 002 and the moment estimated from experimental images corresponding to structure 001. (c) Comparison between diagonal entries of the second moment computed from the reconstructed volumes 003 and 004 and the moment estimated from experimental images corresponding to structure 001. (d) The five reconstructions (000–004) and two baseline structures (EMD-8457 and EMD-2600) ranked using $ {d}_{\mathrm{iKam}} $ , ordered from left to right.

In Figure 8, we plot the distances between the five reconstructions (or in the case of $ {d}_{\mathrm{iKam}} $ , their experimental images) and the seven candidate structures given by both of our metrics. The exact values for $ {d}_{\mathrm{iKam}} $ are given in Appendix B.4. There is also scaling ambiguity in $ {d}_{\mathrm{vKam}} $ since our reconstructions are preprocessed; hence, we use the same approach as above: We scale each candidate structure’s moment by a multiplicative scaling factor that best matches the candidate structure’s diagonal entries of the second moment with those of the ground truth structure. Analyzing the trends in each row, we observe that the metrics seem to agree on the general ranking of the molecules. While the structures 001, 002 and 003, 004 are very similar, $ {d}_{\mathrm{vKam}} $ shows that the metric distinguishes between them given accurate moment estimation, whereas $ {d}_{\mathrm{iKam}} $ loses some discriminative power. However, when it comes to distinct molecules such as EMD-8457 and EMD-2660, both metrics agree on their rankings.

Figure 8. Visualization of $ \log \left({d}_{\mathrm{vKam}}\right) $ and $ \log \left({d}_{\mathrm{iKam}}\right) $ values on the seven candidate structures. Here, EMD-8457 and EMD-2660 are listed as 8457 and 2660 for brevity. Note that there are five ground truth structures but seven candidate structures since EMD-2660 and EMD-8457 are baseline structures for which there are no images in the experimental dataset.

5. Limitations and future work

$ {d}_{\mathrm{iKam}} $ currently falls short of being directly applicable to experimental datasets. As stated in Section 4.4, there are several unmodeled effects not considered in this work that could lead to unexpected results for real data. The net effect of ignoring these experimental considerations is to bias our moment estimator, which may explain the inability of $ {d}_{\mathrm{iKam}} $ to detect the smaller differences between structures 001 and 002, as well as 003 and 004. Developing an estimator that is robust to outliers (such as junk particles) could help alleviate this.

While we address a few of these parameters, we do so with prior knowledge. For example, the shifts used to center images are a byproduct of the reconstruction process. In future work, we aim to develop methods to correct these effects directly from the raw images. Likewise, here we have controlled for experiment-specific artifacts by using images and structures resolved from the same experiment, whereas in the future we wish to compare across all structures. Furthermore, in the future we seek to compare moments computed from real data directly to the PDB, by appropriately correcting for the discrepancies between PDB and reconstructed structures.

Even with our current mitigations, issues such as the B-factor and inaccuracies in the noise model remain completely unmodeled. Further studies will be required to investigate which of these omissions is important and which can safely be made. Then, our method could be modified to account for the important effects.

6. Discussion

We introduced structural similarity metrics for proteins based on moments, inspired by the moment computation in Kam’s method. $ {d}_{\mathrm{vKam}} $ compares known 3D structures according to the difference between the moments of their potentials. We showed that the metric accurately captures similarity according to the rotationally aligned Euclidean metric, an interpretable but expensive-to-compute molecular similarity metric. Therefore, $ {d}_{\mathrm{vKam}} $ allows for the efficient comparison of large number of known structures. A potential application is to improve the similarity search presently in the PDB, which uses the Zernike metric – a fast but less principled metric that involves learning weights and which our results suggest performs worse than ours.

A second metric, termed $ {d}_{\mathrm{iKam}} $ , allows for the computation of a similarity score between an unknown structure present in a large cryo-EM dataset and a solved structure. The computation of this metric does not require a 3D reconstruction process for the image stack and therefore is very efficient. We showed on simulated projection images that our method could discriminate between even very similar proteins with reasonably sized datasets. If it were to work on experimental datasets, $ {d}_{\mathrm{iKam}} $ could become a versatile tool for 3D reconstruction. Typical reconstruction algorithms used in practice are only locally optimal and thus require good initialization, which $ {d}_{\mathrm{iKam}} $ could provide by returning the homologous structures present in the PDB. By extending the database to the entirety of the PDB and including structure predictions, both solved and predicted structures could be quickly compared against.

Beyond its application to experiments, $ {d}_{\mathrm{iKam}} $ demonstrates that Kam’s method is a feasible strategy for high-resolution reconstruction. Recent works have improved the viability of Kam’s method by using sparsity^{( 7 )} or neural network^{( 36 )} priors; likewise, the search over the PDB using Kam’s metric can be interpreted as simply running Kam’s method under a very strong prior, where only a finite number of structures appear with nonzero probability. Our results suggest that, if one could formulate a tractable prior over the manifold of proteins, Kam’s method could yield high-resolution reconstructions.