A common lines approach for ab-initio modeling of cyclically-symmetric molecules
Gabi Pragier, Yoel Shkolnisky
TL;DR
This paper introduces a robust angular reconstitution algorithm for cyclically-symmetric molecules in cryo-electron microscopy, leveraging self common lines in Fourier space to accurately determine 3D structures from noisy 2D projections.
Contribution
It presents a novel method that exploits symmetry and geometrical constraints of self common lines for improved orientation estimation in cryo-EM.
Findings
Effective in high noise conditions
Works for molecules with high-order rotational symmetry
Validated on simulated and experimental data
Abstract
One of the challenges in single particle reconstruction in cryo-electron microscopy is to find a three-dimensional model of a molecule using its two-dimensional noisy projection-images. In this paper, we propose a robust "angular reconstitution" algorithm for molecules with -fold cyclic symmetry, that estimates the orientation parameters of the projections-images. Our suggested method utilizes self common lines which induce identical lines within the Fourier transform of each projection-image. We show that the location of self common lines admits quite a few favorable geometrical constraints, thus allowing to detect them even in a noisy setting. In addition, for molecules with higher order rotational symmetry, our proposed method exploits the fact that there exist numerous common lines between any two Fourier transformed projection-images of such molecules, thus allowing to determine…
| 1/2 | 1/4 | 1/8 | 1/16 | |
|---|---|---|---|---|
| 500 | () | () | () | () |
| 1000 | () | () | () | () |
| 1500 | () | () | () | () |
| 2000 | () | () | () | () |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
A common lines approach for ab-initio modeling of cyclically-symmetric molecules
Gabi Pragier and Yoel Shkolnisky
Department of Applied Mathematics, School of Mathematical Sciences, Tel-Aviv University, Israel
[email protected], [email protected]
Abstract
One of the challenges in single particle reconstruction in cryo-electron microscopy is to find a three-dimensional model of a molecule using its two-dimensional noisy projection-images. In this paper, we propose a robust “angular reconstitution” algorithm for molecules with -fold cyclic symmetry, that estimates the orientation parameters of the projections-images. Our suggested method utilizes self common lines which induce identical lines within the Fourier transform of each projection-image. We show that the location of self common lines admits quite a few favorable geometrical constraints, thus allowing to detect them even in a noisy setting. In addition, for molecules with higher order rotational symmetry, our proposed method exploits the fact that there exist numerous common lines between any two Fourier transformed projection-images of such molecules, thus allowing to determine their relative orientation even under high levels of noise. The efficacy of our proposed method is demonstrated using numerical experiments conducted on simulated and experimental data.
*[inlinelist,1]label=(),
1 Introduction
Cryo-electron microscopy is a method for determining the three-dimensional structure of a molecule from its two-dimensional projection-images [3]. The method consists of generating projection-images of copies of the investigated molecule, where each copy assumes a random unknown orientation before being imaged. Formally, if we denote by the electrostatic potential of the molecule, and consider a rotation matrix
[TABLE]
then the projection-image is given by the line integrals of along the beaming-direction . That is,
[TABLE]
In this work, we focus on molecules that have an -fold, , rotational symmetry about some unknown axis. Such molecules are referred to as molecules with symmetry. We assume that the cyclic symmetry order of the underlying molecule is known from prior knowledge, or may be inferred from rotational invariants computed by spherical harmonics (see e.g., [7]). Intuitively, molecules with symmetry “look exactly the same” when rotated by (we subsequently denote by the set ), radians about their axis of symmetry. Mathematically, it means that the electrostatic potential function of any such molecule satisfies
[TABLE]
where represents a rotation of radians about the unknown axis of symmetry, and is the largest for which (3) holds. Since rotating the three-dimensional coordinate system has no effect on the three-dimensional structure of the molecule, it may be assumed without loss of generality that the axis of symmetry coincides with the -axis. Thus, the matrix which satisfies (3) may be written as
[TABLE]
and any three-dimensional coordinate system for the molecule must be set such that it keeps the axis of symmetry aligned with the -axis (with possibly flipping its orientation). That is, the degrees of freedom in setting the three-dimensional coordinate system of the molecule consist of in-plane rotations about the -axis, and an in-plane rotation of radians about either the -axis or the -axis.
Finding the three-dimensional structure of the molecule amounts to recovering the unknown electrostatic potential function given only the projection-images . Typically, this is done by first finding the rotation matrices , followed by a standard tomographic inversion algorithm e.g., [5, 11].
A fundamental limitation of cryo-electron microscopy is the handedness ambiguity [14], whereby the best one may expect is to recover either the set or the set where , not being able to distinguish between the two. Indeed, given any molecule whose electrostatic potential function is , consider the molecule whose electrostatic potential function is given by . For any , the projection-image of the molecule is given according to (2) by
[TABLE]
where . Noting that , and by changing the variable to we have
[TABLE]
As such, both and are consistent with the same set of projection-images , yet the reconstructed model using the latter set is biologically infeasible, and the true set may therefore only be determined by visual examination of the reconstructed model, possibly exploiting other known structural information. Furthermore, an important property of molecules with symmetry is that any projection-images , are identical. Indeed, from (3)
[TABLE]
Thus, by letting , and integrating over , it follows from (2) that for any ,
[TABLE]
Equations (5), (6), and (8) show that for each in the forward imaging model (2) there exists an equivalence class , with , , and , such that for any . The matrices are due to the ambiguity stated in (7) and so may be chosen independently for each image , while and are common to all images, with being some in-plane rotation matrix about the -axis (the axis of symmetry). In other words, there are no “true” rotation matrics that need to be recovered, as may be implied from (2), but rather we need to recover for each any member of the equivalence class . In summary, in light of the above, given projection-images of a molecule with symmetry, the goal is to recover rotation matrices , where , , and , such that (2) is satisfied for all for some .
The rest of this paper is organized as follows. In Sections 2 and 3 we review the projection slice theorem [10] and further describe the properties of common lines and self common lines which are induced by this theorem. In Section 4 we review some previous work. In Section 5 we present an outline of our proposed method. Next, in Section 6, we describe how to estimate the third row of each rotation matrix . A procedure that assures that all these estimates correspond to a single hand (see Section 1 which points out the inherent handedness ambiguity in cryo-electron microscopy) is presented in Section 7. In Section 8 we then describe a method to determine the remaining first two rows of each . Then, in Section 9 we report some numerical experiments that we conducted on simulated and experimental datasets that show the efficacy of our proposed method. Finally, in Section 10 we present some conclusions and possible extensions of this work.
As the cases of and symmetry exhibit special geometry which is advantageous from the computational point-of-view, we describe in Appendix A an algorithm for estimating the third row of each rotation matrix for molecules with such symmetries. This algorithm may replace the algorithm of Section 6 for such symmetries, as it is much faster and is empirically observed to be more robust to noise. Unfortunately, generalizations of the geometry derived in Appendix A for with are of little practical use.
2 Common lines
The projection slice theorem [10] is a key theorem underpinning many of the methods for recovering the rotation matrices of a given set of projection-images. It states that the two-dimensional Fourier transform of any projection-image is equal to the restriction of the three-dimensional Fourier transform of the electrostatic potential function to the central plane , whose normal coincides with (the third column of of (1)). Mathematically, if we denote by the three-dimensional Fourier transform of , and denote by the two-dimensional Fourier transform of the projection-image , then the projection slice theorem states that
[TABLE]
As such, since any two central planes and which are not parallel intersect along a common axis, it follows that any two Fourier-transformed images and have a pair of lines (one line in each image) on which their values agree. Such pairs of lines are called “common lines”. Hereafter we refer to and simply as “images”. The following lemma shows that any two images of a molecule with symmetry have pairs of common lines (see Figure 1 for an illustration of the case ).
Lemma 2.1**.**
For any two images and , such that , there exist angles , , such that for any ,
[TABLE]
Proof.
Denote by , the unit vector in the direction of the common axis between and , namely
[TABLE]
Also denote for every
[TABLE]
That is, the angles and express in the frame of reference of and in the frame of reference , respectively. Now, by (7), (9) and (12), for and ,
[TABLE]
∎
The relative orientation between any two central planes and is expressed algebraically by (henceforth regarded as the relative orientation between and ). In addition, denoting by , the acute angle between the central planes and , we get that is the “-convention” Euler angles parameterization [17] of . Specifically,
[TABLE]
where
[TABLE]
denote the matrices that rotate vectors by radians about the -axis and -axis, respectively. Finally, we mention in passing that any and may be recovered from the entries of using (see [15])
[TABLE]
3 Self common lines
A “self common line” is a common line between any two images and with . Thus, similarly to (10), for any there exist angles (the subscripts will be clarified shortly), such that for any
[TABLE]
Since by (8), any two such images and are identical, it follows that
[TABLE]
That is, any image has pairs of identical lines (regarded henceforth as self common lines as well). Similarly to (12), the angles and may be expressed in terms of the unit vector of (11) in the direction of the common axis between the central planes and as
[TABLE]
The following lemma shows that the direction of the -th and the -th pairs of self common lines in any image coincide.
Lemma 3.1**.**
For any , and for any ,
[TABLE]
The proof of Lemma 3.1 is given in Appendix B.1. Figure 2a and Figure 2b illustrate Lemma 3.1 for and symmetries, respectively.
Corollary 3.2**.**
By (17) and (19) it follows that for any , , and ,
[TABLE]
In addition, as is conjugate-symmetric (since is real-valued (2)),
[TABLE]
where denotes complex conjugation. As a result,
[TABLE]
We therefore see from Lemma 3.1 and Corollary 3.2 that the pairs of self common lines in any image consist in fact of different pairs. In addition, in case is even, the two self common lines that constitute the -th pair of lines are collinear. In particular, for molecules with symmetry, the sole pair of self common lines in any image consists of two collinear lines. In contrast,
- •
for molecules with symmetry (i.e., ), the two pairs of self common lines (i.e., four lines altogether) in any image coincide (i.e., there are merely two non-collinear lines), see Figure 2a,
- •
for molecules with symmetry (i.e., ), one of the pairs of self common lines in every image consists of collinear lines whereas the other two remaining pairs of self common lines coincide, see Figure 2b, and
- •
for molecules with cyclic symmetry of higher order (i.e., ), each image has in general pairs of non-collinear self common lines.
In what follows, we refer to any relative orientation , as the self relative orientation between and . By (14), any such self relative orientation is parameterized by the ordered triplet which by Lemma 3.1 is equal to . As such,
[TABLE]
where, by (16),
[TABLE]
4 Previous work
The method of angular reconstitution [4, 18] recovers the orientations corresponding to a given set of projection-images in a sequential manner. Specifically, it first determines the common lines between arbitrary three images. This establishes a coordinate system, and the orientations of all remaining projection-images are then recovered one after the other by using the common lines between any given image and those three images. However, due to the sequential nature of angular reconstitution, the method critically depends on correctly detecting the common lines between the first three images. As a result, when the input projection-images are noisy, this positioning might be wrong, which would lead to errors when inferring the relative orientations of the remaining central planes.
In [13], a non-sequential common-lines-based method was suggested, which estimates simultaneously the orientations of all projection-images. As such, contrary to the method of angular reconstitution, it does not suffer from accumulation of errors and is therefore robust even when the input images are noisy. However, this method may not be applied to projection-images of a molecule with symmetry, since there is no way to guarantee a consistent choice of common lines, and so for each pair of images and we can only estimate for an unknown , whereas to apply the method in [13] all must be the same. Unfortunately, such a consistency in all cannot be guaranteed using the method of [13].
In [2], a method which successfully handles molecules with symmetry was proposed. The method is based on determining for every two images and both of their two common lines. Unfortunately, generalizations of this method to for are of little practical use, as it requires to detect all common lines (Lemma 2.1) of every two images, which is impractical when the input images are noisy. Nevertheless, the method presented in the current paper is in some sense an extension of [2], but with two key differences. First, detecting all common lines between each pair of images is replaced by a maximum-likelihood-type procedure that directly estimates their relative orientation, from which their common lines can be easily computed. Second, our method takes advantage of self common lines, which do not exist in the case of symmetry.
5 Outline of our method
In this section, we provide an outline of our method for recovering the orientations of a given set of projection-images of a molecule with symmetry with . Figure 3 depicts a flowchart of the method.
In what follows, we denote by the third row of . Given projection-images with corresponding unknown rotations (to be estimated), the first step of our method, denoted as “relative viewing directions estimation” (Section 6), consists of estimating for each pair of projection-images and , , one of the two matrices or , where , using common lines and self common lines of and . The two matrices and are indistinguishable using common lines due to the handedness ambiguity discussed in Section 1. We thus denote the estimated matrix by . The second step, which is denoted as “handedness synchronization” (Section 7), enforces that either all estimates have a spurious , or none have at all. Once these estimates all correspond to a single hand, the next step, denoted as “viewing directions estimation”, consists of forming a symmetric matrix whose -th block of size is given by . That is,
[TABLE]
Depending on the output of the handedness synchronization step described above, the factorization
[TABLE]
yields at once either the estimates of the third rows of all rotation matrices , or the estimates of all -multiplied third rows . The next step is to form rotation matrices where the third row of each is set to be equal to the estimate for either or . The first two rows of each are set arbitrarily (so that ). In Lemma 8.1 below, we prove that for any there exist and such that
[TABLE]
where is a rotation matrix that rotates vectors by an angle of about the -axis (see (15)). As a result, since each may be replaced by either or for some , it suffices to recover the in-plane rotation angles . The step which consists of recovering all these angles is denoted by “in-plane rotation angles estimation” (Section 8). Finally, the last step is “orientations estimation” in which all rotation matrices (or ) are formed using (27).
We next present a method for estimating the set of all relative viewing directions which may applied to molecules with symmetry with (a full treatment of symmetry may be found in [2]). In addition, we describe in Appendix A an alternative method that is applicable only to molecules with either or symmetry. While the method of Section 6 may be applied to molecules with or symmetry, the alternative method of Appendix A utilizes the underlying geometry that is induced by such molecules. As a result, it produces better results in practice and is much faster.
6 Relative viewing directions estimation for symmetry with
Loosely speaking, estimating is based on inspecting for each pair of images all possible pairs of rotation matrices (discretized in some proper manner, as explained later on in Section 9), and finding the pair that induces common lines and self common lines which are most correlated. Such an approach is advantageous for the following two reasons:
By Lemma 2.1, any two images have pairs of common lines, and by Corollary 3.2, any image has pairs of (non collinear) self common lines. Therefore, for any two images, the degree by which a given pair of candidate rotation matrices induces the correct relative orientation of the two images, may be ascertained with greater confidence as increases. Specifically, if most of the common lines and self common lines that are induced by a given pair of candidate rotation matrices are highly correlated, then it is likely that these candidates correspond to the true rotation matrices. 2. 2.
Recall that any may be replaced by where is arbitrary. In light of that, for each , if we write the third column of in its spherical coordinates representation, i.e., for some and , then by a direct calculation, we get that for any ,
[TABLE]
Thus, since for any there exists such that , it follows that instead of considering the set of “all possible” candidate rotations for each image, it suffices to only consider the set given by
[TABLE]
where is the mapping of any vector in to its azimuthal angle in its spherical coordinates representation. For large this significantly restricts the search space of the candidate rotation matrices.
In light of reason (2) above, throughout the remaining of this section, we assume without loss of generality that (see (29)). The following lemma, whose proof is given in Appendix B.2, and the corollary that follows will be used in the sequel.
Lemma 6.1**.**
For any , , and such that ,
[TABLE]
Corollary 6.2**.**
By applying (30) to any , , with , we get that for any ,
[TABLE]
At first glance, it would appear that the rotation matrices may be estimated by searching for each pair of images for the pair of rotation matrices that induce pairs of lines whose values are most correlated. That is, since the values along common lines and self common lines have perfect correlation over all other pairs of lines, it follows from (10) and (22) that the pair of rotation matrices attains the maximal value of , given by
[TABLE]
with denoting the real part of , where each ray in each of the images in (32) is normalized to have norm equal to one, and where (in accordance with (16) and (24)),
[TABLE]
However, for any , the maximum of (32) over is not necessarily unique, and—depending on and —might be attained by other pairs of rotation matrices besides . For example, by (33), any pair of rotation matrices such that
[TABLE]
and
[TABLE]
would also maximize (32). Nevertheless, while not every maximizer of (32) is necessarily equal to , we observed empirically (using extensive simulations) that all such maximizers satisfy (34) and (35). As a result, it follows from (31) that
[TABLE]
where and are the third rows of and , respectively. That is, any outer product is invariant to permutations of the common lines in and . Thus, for any , we choose an arbitrary pair which maximizes (32), and obtain an estimate for which, due to the inherent handedness ambiguity satisfies . In a similar vein, it follows from (34) and (31) that, for any , any maximizer of (32) yields
[TABLE]
Thus, for any , any of the pairs where , induces an estimate for using (37). In practice, however, due to self common lines misidentification, any such estimate of may contain some error. Thus, choosing any one of the estimates to be the single estimate for is sub-optimal. Furthermore, averaging over all estimates doesn’t make any sense, since due to the handedness ambiguity, independently of other estimates. Instead, since any is a matrix of rank-, we set to be equal to the estimate that is closest to a rank- matrix. Specifically, for every we compute the estimates , and for every such estimate we find (using SVD) its three singular-values , and set where
[TABLE]
The procedure for finding all estimates and which, due to the inherent handedness ambiguity, satisfy and is summarized in Algorithm 1.
7 Handedness synchronization
At this stage, we have determined the estimates , of all relative viewing directions, where for each estimate either or independently of other estimates. In this section, we describe the “handedness synchronization” step, where the task is to manipulate these estimates so that either for all , or for all . Once this procedure is completed, we can form the matrix from (25) and infer the third rows of all rotation matrices (or all -multiplied third rows ) using (26). Handedness synchronization is done in two steps. The first step synchronizes the estimates . The second step synchronizes each of the remaining estimates , with the synchronized estimates from the first step. The reason for separating the synchronization of from that of will be clarified below.
7.1 Step : Synchronizing the estimates
We employ a procedure similar to the one described in [13]. Specifically, the task of synchronizing the set of estimates may be reduced to the task of partitioning this set into the following two disjoint sets
[TABLE]
Indeed, once all estimates are partitioned into and , we can choose either one of the sets (does not matter which one), and replace each estimate in it by . As a result, since , we get that either for all , or for all , as needed. We now describe how we obtain such a partition (39).
Let us define the “handedness graph” to be the undirected graph whose set of nodes consists of all estimates , that is
[TABLE]
and whose set of edges consists of the undirected edges between all triplets of estimates , , and (hence each triplet forms a “triangle”), that is,
[TABLE]
The weight of each edge is set to either or as explained below. For each , we consider the three estimates , and , along with the “triangle” they form in the graph. The goal is to set the weight to all edges of the triangle whose incident nodes correspond to estimates that belong to the same set in (39), and to set the weight of all other edges to . A crucial observation is that for any ,
[TABLE]
As such, we check which of the following expressions
[TABLE]
equals the zero matrix, and assign the weight of each of the three edges in the corresponding triangle as illustrated in Figure 4. The expressions listed in (43) allow to decide for each triplet , , which nodes belong to the same set in (39) and which belong to different sets. For three nodes the only two possibilities are that either all nodes belong to the same set (case 1 in (43)), or that one node belongs to one of the sets and the remaining two nodes belong to the other set (cases 2, 3, 4 in (43)). Thus, by (43) we determine the partition (39) “locally” for each triplet of nodes. For example, if all three estimates belong to the same set in (39) (i.e., either all belong to , or all belong to ), then the first expression in (43) is equal to the zero matrix, meaning that the weights of all three edges in the corresponding triangle in the graph are set to (as per Figure 4a). Indeed, in case , , and belong to , then since , we get
[TABLE]
and in case all these three estimates belong to , then since , it also follows that
[TABLE]
Once we have determined the local assignment of each triplet of nodes (via the weights on the edges between the vertices corresponding to the nodes), we obtain a global partition into the sets in (39) similarly to [13]. Specifically, we define the weighted adjacency matrix of , also denoted by , as the matrix whose entries are given by
[TABLE]
We then calculate the eigenvector that corresponds to the leading eigenvalue of the matrix . As was shown in [13], this eigenvalue has multiplicity one and its corresponding eigenvector encodes the set membership of the estimates. Specifically, if then belongs to one of the sets of (39), and if then belongs to the other set of (39). As such, by -conjugating all estimates in either one of the sets we are guaranteed that either for all , or for all .
Notice that, in practice, the estimates are computed from noisy projection-images, and thus for many triplets of estimates none of the four expressions listed in (43) might be equal exactly to the zero matrix. Thus, instead, we search for the expression that is as close as possible to the zero matrix. Specifically, we minimize
[TABLE]
over , subject to the constraint that , where each possible triplet corresponds to one of the four expressions in (43), and denotes the Frobenius norm.
7.2 Step : Synchronizing the estimates
The second step of handedness synchronization consists of synchronizing each of the estimates so that if the previous step (described in Section 7.1) resulted in every satisfying , then the goal is to enforce that for every . Otherwise, if the output of the previous step is such that every satisfies , then the goal is to enforce that for every . Recall, however, that it is unknown which of the above two possible outputs was obtained. Nevertheless, since for any , and since , it follows that for any ,
[TABLE]
As such, we can in principle synchronize every by choosing an arbitrary such that , and reset as follows:
[TABLE]
Indeed, if (which we cannot tell), then by (46), case (48a) occurs if and therefore should not be altered. Similarly, if then by (46), case (48b) occurs if (since ) so that by assigning we indeed end up having , as needed. The case of is analogous in light of (47).
In practice, however, since each of the above estimates is computed from noisy images, it might be that neither (48a) nor (48b) hold. In addition, it is desirable to synchronize each based on all estimates such that rather than only using a single such estimate. Thus, instead, each estimate is set according to the majority-vote over all . Specifically, let us denote by the sign function (which equals if and equals otherwise). Then, for every we reset to be in case that
[TABLE]
Once this second step is completed, we are guaranteed that all estimates are synchronized, i.e., either for all , or for all . As such, we then construct the matrix of (25), factorize it as in (26), and obtain all third rows (or ). The procedure for handedness synchronization is summarized in Algorithm 2.
Note that synchronizing the estimates is done separately from synchronizing the estimates . Synchronizing based on (42) involves triplets of indices , and so is the construction of the graph in (44). On the other hand, the synchronization of is based on pairs of indices (see (48a) and (48b)). So while it may be better to synchronize and simultaneously, it is currently unclear how to combine the pairwise information required for synchronizing into the triplets stricture required for constructing .
8 In-plane rotation angles estimation
At this stage, all third rows (or ) of the rotation matrices (or ) have been obtained. In this section, we describe a procedure to determine the remaining first two rows in each of these rotation matrices. The following lemma shows that any two such rows are determined by a single parameter, namely an in-plane rotation angle about the -axis (the axis of symmetry).
Lemma 8.1**.**
Let and be any two rotation matrices with identical third rows. Then, for any , there exist a unique and a unique such that
[TABLE]
where (given by (15)) is the matrix that rotates vectors by an angle about the -axis.
The proof of Lemma 8.1 is given in Appendix B.3. In light of Lemma 8.1, we next form rotation matrices by setting the third row of each to be equal to the estimate of the third row of (see Section 7), and by arbitrarily setting the first two rows of each (while ensuring that ). As a result, due to Lemma 8.1, for any , there only remains to recover the in-plane rotation angle , and form the matrix . In principle, all such angles may be determined in a sequential manner. However, we instead employ a more robust approach in which all are determined in a single step. Specifically, as we next show, the task of determining all angles may be reduced to finding for every , the relative in-plane rotation angles
[TABLE]
Indeed, by the definition of , there exist unique , , such that
[TABLE]
Thus, by letting be the matrix whose -th entry is equal to , it follows that
[TABLE]
As such, is a hermitian rank- matrix whose factorization is given by
[TABLE]
from which all angles may be retrieved using
[TABLE]
Since is a rank- matrix, its factorization (54) yields the eigenvector where with is arbitrary and unknown. As a result, (55) actually yields the angles , for some arbitrary, unknown . This poses no problem, since for any , recovering is just as good, as we have the degree of freedom of applying any single in-plane rotation about the -axis to all rotation matrices (see the paragraph following (4)). In light of (52)–(55), we describe below how to determine all of (51). We then form the matrix above, obtain its factorization (54) using SVD, and recover all in-plane rotation angles using (55).
By applying Lemma 8.1 to any two rotation matrices and , we get that there exist and unique angles , such that for any ,
[TABLE]
where the last equality follows by denoting and using the fact that . As such,
[TABLE]
where the second equality used (51) and the fact that for any . Given and , exactly one of the angles among , , lies in . Thus, by (51) and (57), is the only angle in which satisfies
[TABLE]
Equation (58) is the set of all possible relative orientations between and . Thus, if we know in (58), then the common lines between and induced by these relative orientations will perfectly correlate (see (10)). Therefore, for any , we set the angle to be the maximizer over all of
[TABLE]
where each ray in each image in (59) is normalized to have norm equal one, and where (in accordance with (16)), for any ,
[TABLE]
As a result, any of (51) may be obtained by optimizing of (59) over all . The procedure for determining all in-plane rotation angles is summarized in Algorithm 3.
Based on Algorithms 1, 2, 3 and on Section 5, the end-to-end algorithm for recovering all rotation matrices is summarized in Algorithm 4. Note that the fact that we assumed in (4) without loss of generality that the axis of symmetry coincides with the -axis, means that any reconstructed volume that is based on the output rotation matrices of Algorithm 4 will have its axis of symmetry aligned with the -axis as well.
9 Numerical experiments
We implemented Algorithm 4 in Matlab and tested it on both simulated and experimental projection-images. Section 9.1 provides some of the implementation details of Algorithm 4, and Section 9.2 analyzes its time and space complexity. Section 9.3 describes the experiments conducted using noisy simulated projection-images of the three-dimensional density map EMD- [20], which has symmetry. Section 9.4 presents results for the Trimeric HIV- envelope glycoprotein [22] which has symmetry. Section 9.5 focuses on the Human HCN1 hyperpolarization-activated cyclic nucleotide-gated ion channel [24] which possesses symmetry. Finally, Section 9.6 reports results for the GroEL protein [9]. Technically, GroEL has both a -fold cyclic symmetry about some axis, as well as a -fold cyclic symmetry (about a different axis which is perpendicular to the above-mentioned axis). Thus, strictly speaking, GroEL has (dihedral) symmetry. We nevertheless applied Algorithm 4 to it while taking into account only its -fold cyclic symmetry (i.e., we treated the molecule as ). The code of all algorithms presented in this paper is available as part of the ASPIRE software package [21].
9.1 Implementation details
All tests were executed on a dual Intel Xeon X5560 CPU (12 cores in total), with 96GB of RAM running Linux and an nVidia GTX TITAN GPU. Whenever possible, all cores were used simultaneously, either explicitly using Matlab’s parfor, or implicitly, by employing Matlab’s implementation of BLAS, which takes advantage of multi-core computing. Some loop-intensive parts of the algorithm were implemented in C as Matlab mex files. We next describe the discretization of the set of (29). As implied by Lemma 8.1, any rotation may be written as where
1 is a rotation matrix whose third row is given by for some angles and , and whose first two rows are set arbitrarily, and
2 is the matrix (15) that rotates a vector by an angle about the -axis.
As such, we sampled evenly-spaced points , and formed a matrix from each such point (we found experimentally that the resulting number of rotations is adequate).
9.2 Complexity analysis
Strictly speaking, the computational complexity of Algorithm 4 is cubic in the number of images due to Algorithm 2. However, in practice, the running time of Algorithm 4 is governed by Algorithm 1, which is quadratic in both the number of images as well as in the size of of (29). As such, the computational complexity of Algorithm 4 is where is the number of rotations in the discretization of the set . The space (storage) complexity of Algorithm 4 is .
As of timing, using a simulated set of projection-images of a molecule with symmetry, it took seconds to compute all relative viewing directions , seconds to resolve the handedness synchronization, seconds to estimate the in-plane rotation angles, and seconds in order to reconstruct the density map.
9.3 Simulated noisy data
We first tested Algorithm 4 on several datasets of simulated noisy images. Specifically, we generated four sets of projection-images, with . The projection-images in each set were generated from the three-dimensional density map EMD- [20] (which possesses symmetry) available in the Electron Microscopy Data Bank (EMDB) [23]. The orientation of each projection-image was drawn uniformly at random from the uniform distribution on , and the size of each projection-image was pixels. We generated four different copies of each of these sets of projection-images by corrupting the projection-images in each set by an additive Gaussian white noise with , where SNR (signal-to-noise ratio) is defined as the ratio between the energy (variance) of the signal and the energy of the noise.
We then applied Algorithm 4 to each of the resulting sixteen sets of projection-images. For each such set we obtained an estimated set of rotations, which we denote by , and we reconstructed the volume using the noisy projection-images and the estimated set of rotations. We next describe how we assessed the degree by which each of these sets of estimated rotations differs from the true rotations . To this end, we sampled the Fourier-transform of each projection-image along the direction vectors , where , and then lifted these direction vectors to , i.e., we defined , . We then computed for each pair the angle between and , that is
[TABLE]
which measures the error in the estimation of the three-dimensional position of the Fourier ray . In addition, the quality of each of the reconstructions was assessed by the resolution obtained using the criterion of the Fourier shell correlation (FSC) curve [19] with respect to the reference density map EMD- available in [23]. Table 1 lists the median angular error of (over all and ) for each of the sets along with the obtained resolutions. The table illustrates the robustness of Algorithm 4 to noise as the number of images increases.
9.4 Trimeric HIV-1 envelope glycoprotein (
We next tested Algorithm 4 on the Trimeric HIV- envelope glycoprotein dataset. The dataset consists of raw particle images provided in the EMPIAR- dataset [22] from the EMPIAR archive [6]. The raw particle images are of size pixels, with pixel size of . We processed the raw particle images using the ASPIRE [21] software package as follows. First, all images were phase-flipped (in order to remove the phase-reversals in the CTF), down-sampled to size of pixels (hence with pixel size of ), and normalized so that the noise in each image has zero mean and unit variance. We next used the class-averaging procedure in ASPIRE [21] to generate class averages from the raw particle images, where each image was averaged with its most similar images (after proper rotational and translational alignment). Next, we sorted the class averages according to their contrast (i.e., according to the standard deviation of the pixel values of each average). The input to Algorithm 4 was the class averages with the highest contrast. A sample of these class averages is displayed in Figure 5.
Next, we applied Algorithm 4 to estimate the rotation matrices that correspond to the class averages . Then, instead of reconstructing the three-dimensional density map using merely the pairs , we made full use of the fact that the underlying molecule is symmetric by applying the reconstruction to the pairs (see (8)). Figure 6a displays the reconstructed density map, and Figure 6b displays the reference density map EMD- available in [1] in [23]. The renderings of all volumes in this section were generated using USCF Chimera [12]. The quality of the reconstruction was assessed using the Fourier shell correlation (FSC) curve [19], implying that the resolution of the model estimated by our algorithm is equal to according to the criterion (Figure 6c).
9.5 Human HCN1 hyperpolarization-activated cyclic nucleotide-gated ion channel ()
Next, we applied Algorithm 4 to class averages of the Human HCN1 hyperpolarization activated channel which possesses symmetry. The class averages were generated from the particle images provided in the EMPIAR- dataset [24]. This dataset comprises of raw particle images of size pixels, with pixel size of . First, the raw particle images were phase-flipped, down-sampled to size of pixels, and normalized so that the noise in each image would have zero mean and unit variance. To examine the consistency of Algorithm 4, the raw projection-images were randomly split into two groups of images each, and the class-averaging procedure in ASPIRE [21] was used to generate class averages from each of the two groups independently. The class averages were generated by averaging each raw image with its most similar images (using resulted later on in inferior results). The input to subsequent steps were the top (highest contrast) class averages from each set. A sample of class averages is displayed in Figure 7.
Next, we applied Algorithm 4 to each of these two groups of class averages and estimated the rotation matrices corresponding to each group. We then reconstructed the two density maps using the class averages and the corresponding estimated rotation matrices, while considering the symmetry in the reconstruction process as was described in Section 9.4. The consistency of the reconstructions from the two groups of the data was first assessed using the criterion of the FSC curve [19], and was found to be equal to (Figure 8a). In addition, we compared (using the criterion) the reconstructions against the reference density map which was reconstructed from the same dataset as described in [8], and found the resolution to be equal to (Figure 8b). Figure 9a displays a two-dimensional rendering of a density map generated by Algorithm 4 (only the reconstruction of the first group is shown), and Figure 9b displays a two-dimensional rendering of the reference density map [8].
9.6 GroEL protein ( symmetry)
As was already mentioned above, strictly speaking, the GroEL protein [9] has (dihedral) symmetry, meaning that it has both a -fold cyclic symmetry as well as a -fold cyclic symmetry. We nevertheless applied Algorithm 4 to it while taking into account only its -fold cyclic symmetry. To this end, we used once again the class-averaging procedure in ASPIRE [21] to generate class averages from the raw particle images, where each image was averaged with its most similar images. A sample of class averages is displayed in Figure 10. We then picked the top (highest-contrast) class averages, estimated the set of corresponding rotation matrices using Algorithm 4, and reconstructed the density map. The resolution was found to be (using the criterion). Three different views of the reconstructed density map are shown in Figure 11.
10 Discussion and future work
In this paper, we proposed a method for finding the orientations that correspond to a given set of projection-images of a cyclically-symmetric molecule. In addition, we described the inherent geometry that underlies such molecules, as well as the way this geometry is expressed in their projection-images. We further demonstrated the efficacy of our proposed method by providing some numerical results using simulated and experimental datasets.
A typical pipeline for reconstructing a three-dimensional volume from a dataset of raw particle images consists of first generating a low-resolution model of the molecule from a subset of the dataset, which is then refined to a high-resolution model using the entire dataset. The reason for breaking the reconstruction process into these two steps is that all high-resolution refinement algorithms are based on non-convex optimization schemes (such as the EM-algorithm or stochastic gradient descent) which must be initialized properly in order not to converge to a molecule which is inconsistent with the data. Thus, the first step of the reconstruction pipeline (known as ab-initio reconstruction) should generate a reliable low-resolution model of the molecule. The algorithm presented in this paper addresses this step of the pipeline, and has been shown using experimental data to produce reliable low-resolution models for three different datasets. The resolutions achieved by our algorithm are significantly lower than that of the reference models since our algorithm was applied to only part of the data (in the form of a few hundreds or thousands of class averages), and not to the entire raw data set. Nevertheless, the obtained resolutions are consistent with the required resolutions at the ab-initio modeling step. In particular, high-resolution refinement algorithms typically low-pass filter the ab-initio model, and so any reconstruction whose resolution is better than 30 Å is typically sufficient.
Obviously, all existing software packages include some functionality for ab-initio modelling. However, all of them are based on some local optimization and have no mathematical guarantees. The algorithm presented in this paper is the first to be specifically designed to the geometry of the problem.
A natural future research is to extend the work to other symmetry groups. We are currently at the final stages of devising an algorithm for molecules with symmetry. As it turns out, the geometry of molecules with symmetry is completely different from that of molecules with symmetry. The reason is that molecules with symmetry have three perpendicular symmetry axes, and three corresponding generators of the symmetry group. The resulting algorithm is thus completely different than the one in the current paper. A preliminary analysis of the algorithm suggests that it may be extended to symmetry with .
Once we derive an algorithm for molecule with symmetry for , there remain three symmetry groups to be handled – (tetrahedral), (octahedral), and (icosahedral) symmetries. Those enjoy very high-order symmetry, that should be advantageous due to the large number of common lines between any two images and within each image (self common lines). These high-order symmetry groups are currently under investigation.
Appendix A Relative viewing directions estimation for or symmetry
In this section, we describe an alternative method for estimating the set of all relative viewing directions which is applicable only to molecules with either or symmetry. The following lemma, whose proof is given in Appendix B.4, is central to the proposed method.
Lemma A.1**.**
For any and for any and ,
[TABLE]
Similarly, for any and ,
[TABLE]
As a result of Lemma A.1 we have the following:
- •
For , since the only self relative orientations for any are and , and since and vice versa, it follows that in order to recover and it suffices to determine either one of these two self relative orientation for each . To recover , it is also required to determine a single, arbitrary, relative orientation .
- •
For , we shall later show that recovering may be easily avoided for any . Thus, since for any the only remaining self relative orientations besides are and , it follows that in addition to recovering a single, arbitrary, relative orientation , recovering either one of the remaining two self relative orientations or for each is sufficient in order to recover any and any .
Applying Lemma A.1 to molecules with symmetry with requires the estimation of more than just a single self relative orientation per image, and was found to be not robust in practice, and so the method of this section may be applied to either or symmetry. The advantage of this method is that it provides more accurate results in practice than the method of Section 6 and is also significantly faster.
A.1 Estimating self relative orientations
We next describe a robust procedure for determining, for both and , and for every , an estimate such that
[TABLE]
By (23), is parameterized by and is parameterized by . We next show that for any . As a result, since by (62a)–(63b) we are oblivious as to which of the self relative orientations in (64) corresponds to, the angles and in the aforementioned parameterizations may be freely interchanged. By the projection slice theorem
[TABLE]
Thus, since , it follows that . As a result, since both of these angles are acute, it follows that indeed (which we subsequently denote by ). In order to recover the angles and , for which the values along the lines they subtend in are conjugate equal (see (22)), let us define for any the mapping by
[TABLE]
where each ray in is normalized to have its norm equal to one. By (22), the two angles and subtend lines in whose Fourier transforms agree up to conjugation. As such, both and are solutions of the optimization problem
[TABLE]
The constraint in (68) is needed because any ray through the origin in is conjugate symmetric, and therefore, any two collinear lines would otherwise maximize (68). Moreover, for this constraint also guarantees that the collinear lines that constitute the second pair of self common lines do not maximize (68). For otherwise, it would lead to estimating instead of the desired for (62a) and (63a), or for (62b) and (63b).
The main incentive to use self relative orientations is the ability to estimate them in a robust manner due to the following two properties:
The domain of each mapping defined in (67) may in fact be restricted to a narrower range of angles. Specifically, we show in Lemma A.2 below that for it holds that , and for it holds that . As such, when the input projection-images are noisy, constraining the optimization problem (68) to these narrower ranges of angles increases the probability of detecting and . 2. 2.
We show in Lemma A.3 below that, for both cases and , each of the angles may be computed directly from . This is in sharp contrast to relative orientations in general, in which the common lines with a third arbitrary central plane are required [16] in order to determine any such angle . In particular, when the input images are noisy, the common lines with the third image might be misidentified, leading to a wrong estimation of .
Lemma A.2**.**
For any ,
[TABLE]
Lemma A.3**.**
For any ,
[TABLE]
The proof of Lemma A.2 is given in Appendix B.5, and the proof of Lemma A.3 is given in Appendix B.6. Based on (68), on Lemma A.2, and on Lemma A.3, the procedure for determining for molecules with either symmetry or symmetry is summarized in Algorithm 5.
A.2 Estimating relative orientations
In light of Lemma A.1, we next describe how to determine for each of the cases or , and for every , a single relative orientation where is arbitrary, and may be different for each . To this end, for every two images and , we first determine using normalized cross correlation a single common line between these images. We then find the acute angle between the underlying central planes using the voting scheme [16]. Finally, using (14) we find an estimate for a relative orientation of the central planes, which due to the handedness ambiguity corresponds to either or for some unknown .
A.3 Local handedness synchronization
At this stage, any two estimates and (obtained by Algorithm 5) satisfy
[TABLE]
for some unknown . In light of (63a) and (63b), we therefore set
[TABLE]
which guarantees that for every . Similarly to (71), any estimate (obtained in Section A.2) satisfies
[TABLE]
for some unknown . However, in order to find an estimate for using either (62a) or (62b), it is essential that
1 2 either all three estimates , , and have a spurious , or none do at all.
In other words, for every , the task is to manipulate , , and so that they correspond to one of the sets
[TABLE]
followed by setting
[TABLE]
as per (62a) and (62b). By so doing, it follows that whenever one of the first two sets in (74) is obtained, and whenever one of the last two sets in (74) is obtained. The task of obtaining for every either one of the four sets in (74) is referred as the “local handedness synchronization” of the estimates, and will be addressed next. Once this task is completed, we are guaranteed that for every .
A crucial observation is that both matrices and are rank-. In addition, if an estimate has a spurious , e.g., if , then since it follows that . Also, if , then since , we get that and vice versa. Thus, for the case of (i.e., ), we examine for every , which of the following expressions yields a rank- matrix, each obtained by -conjugating a subset of , and choosing either (expressions - in (76)) or (expressions - in (76)).
[TABLE]
For example, the first expression in (76) would yield a rank- matrix in case and in addition either all three estimates have a spurious in them, in which case , or when none of these estimates have a spurious in them, in which case . The fifth expression in (76) would yield a rank- matrix in the same cases as the first expression, only that . As another example, consider the case where only has a spurious in it. Then, the second expression in (76) would yield a rank- matrix in case (otherwise, if , then the sixth expression would prevail).
Similarly, due to (71) and (73), in order to use (62a) and (62b) for the case of (i.e., ), we need to examine which of the following eight expressions
[TABLE]
yields a rank- matrix. However, it can be shown (see Appendix C) that these eight expressions are in fact equivalent to
[TABLE]
We thus inspect which of the expressions in (78) results in a rank- matrix. Note that in practice, due to misidentification of (self) common lines, it might be that none of the above expressions yields a rank- matrix. Therefore, we choose the expression that is closest to be rank-. Specifically, we first compute the three singular values of each of the expressions of (76) for , or of (78) for . Then, we choose the expression such that
[TABLE]
and decide accordingly whether or not to transpose , whether or not to -conjugate , and whether or not -conjugate . Finally, we apply (75) to obtain .
The procedure for finding (for molecules with either symmetry or symmetry) all estimates and which, due to the inherent handedness ambiguity, satisfy and , is summarized in Algorithm 6. This algorithm may replace Algorithm 1 for molecules with either or symmetry.
Appendix B Proofs
B.1 Proof of Lemma 3.1
Let and let . By (11),
[TABLE]
Since , (79) is equal to
[TABLE]
Next, since the cross-product is an anti-symmetric operation, (80) is equal to
[TABLE]
where the first equality in (81) is because the vector -norm is rotation invariant, and therefore the denominator may be written as
[TABLE]
and the last equality in (81) is due to (11). From (79)–(81) we get that
[TABLE]
Thus, from (18),
[TABLE]
from which it follows that
[TABLE]
∎
B.2 Proof of Lemma 6.1
By (4),
[TABLE]
Next, by using the assumption that , we get that
[TABLE]
As such,
[TABLE]
and therefore, the right hand side of (82) is equal to , as needed.
∎
B.3 Proof of Lemma 8.1
We first prove that given any such and , there exists a unique angle such that
[TABLE]
To this end, let us denote the three rows of by , and the three rows of by . Since , it follows that and , as well as and . Thus, by direct calculation
[TABLE]
for some . Since is closed with respect to matrix multiplication, it follows that as well. As such, there exists a unique angle such that
[TABLE]
Finally, by right multiplying (85) by we get that (since ),
[TABLE]
which proves (83). We next prove (50). Let , and let be the unique number such that , and further define . By construction, . In addition,
[TABLE]
where (86) follows from (83) and from the fact that , (87) is because for any , and (88) uses the definition of above. Finally, the uniqueness of follows from the uniqueness of , the uniqueness of and from (86)–(88).
∎
B.4 Proof of Lemma A.1
We first prove (62a) and (62b). Let and let . Since , we get that for any . In addition, since , it follows that is a generator of . As a result,
[TABLE]
Thus, for any and for any ,
[TABLE]
Since for any it holds that , applying Lemma 6.1 using to the right hand side of (90) yields
[TABLE]
where the first equality used the fact that for any the third column of is equal to , and therefore . This proves (62a) and (62b). We next prove (63a) and (63b). To this end, for both and ,
[TABLE]
where (91) follows from (89), and (92) follows from Lemma 6.1 using which indeed satisfies for any .
∎
B.5 Proof of Lemma A.2
We need first the following two lemmas.
Lemma B.1**.**
For any , , and ,
[TABLE]
where and are defined in (11).
Proof.
Let , , and . By (11)
[TABLE]
where the second equality uses the fact that the vector -norm is rotation invariant, and therefore
[TABLE]
The nominator in (94) may be simplified by using Lagrange’s identity which holds for any , so that
[TABLE]
where (96) uses the fact that , and (97) is because . As for the denominator in (94), since the magnitude of the cross-product is given by the sine of the angle between its arguments, and since , it follows that
[TABLE]
Plugging-in (97) and (98) into (94) yields (93) which completes the proof. ∎
Lemma B.2**.**
For any , let be the representation of in spherical coordinates for some and . Then,
[TABLE]
Proof.
Let . For any ,
[TABLE]
where (100) is due to (18) and because , and (101) is the application of Lemma B.1 with . By (4), we get by a direct calculation that for any ,
[TABLE]
from which it follows that
[TABLE]
We first prove (99a) (i.e., the case ). To this end, (100)–(101) become
[TABLE]
Since it follows that
[TABLE]
As such, (103) reduces to
[TABLE]
Next, applying (102) with and gives
[TABLE]
and plugging this in (104) yields (99a).
We next prove (99b). Applying (102) with and , and with and yields
[TABLE]
Plugging this in (101) yields
[TABLE]
which proves (99b). ∎
We are now ready to prove Lemma A.2.
Proof.
We begin by proving (69a). Consider the function where
[TABLE]
Then, in light of (99a), it suffices to show that for any . As it may readily be verified, is periodic with period , and therefore is periodic with the same period. In addition, for any , see Figure 12a. As a result, it suffices to show that for any . To this end, by the chain-rule,
[TABLE]
for any . As such, is non-decreasing in . Thus, since is continuous, it follows that for any ,
[TABLE]
which proves (69a), see Figure 12a. Similarly, in order to prove (69b), consider the function where
[TABLE]
In light of (99b), it suffices to show that for any . Since is periodic with period , it follows that is periodic with the same period. In addition, for any , see Figure 12b. As a result, it suffices to show that for any . To this end, by the chain-rule
[TABLE]
Thus, since for any , we conclude that is non-decreasing in , and since is continuous, it follows that for any ,
[TABLE]
which proves (69b), see Figure 12b. ∎
B.6 Proof of Lemma A.3
Let . We begin by proving (70a) (i.e., the case where ). It holds that (in fact for any and not just for ),
[TABLE]
where the first equality in (109) is due to (18), the second equality in (109) is because , and (110) is the result of applying Lemma B.1 with . For , (109)–(110) become
[TABLE]
Next, by the projection slice theorem, it follows that
[TABLE]
where we have used the fact that (see text after (66)). Thus, (111) may be written as
[TABLE]
and solving (113) for yields (70a).
We next prove (70b) (i.e., the case where ). Let , and let be the representation of in spherical coordinates for some and . On the one hand, by the projection slice theorem,
[TABLE]
Next, since , we get by a direct calculation that
[TABLE]
and therefore by a direct calculation it follows that
[TABLE]
Thus, from (114) and (115) we get that
[TABLE]
On the other hand, by (99b) in Lemma B.2
[TABLE]
Equation (70b) now follows from (116) and (117).
∎
Appendix C Justification for the expressions in (78)
We prove that the expressions listed in (78) are equivalent to the expressions in (77). To this end, notice that since , it follows that for any ,
[TABLE]
As such, if for example, and , or for example, and , then setting , and in (77), we get using (118) that
[TABLE]
which is (twice) the first expression in (78). Alternatively, if for example and , or for example and , then by (118), setting and in (77) yields (twice) the fifth expression in (78). Other cases are similar and correspond to cases where either has a spurious , or has a spurious , or both have a spurious , in which case setting in (77) (respectively) either , or , or would yield each of the remaining expressions listed in (78).
Acknowledgments
This research was supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement 723991 - CRYOMATH) and by Award Number R01GM090200 from the NIGMS.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] Bartesaghi, A, Merk, A., Borgnia, M. J., Milne, J. L. S., and Subramaniam, S. Prefusion structure of trimeric HIV-1 envelope glycoprotein determined by cryo-electron microscopy. Nature structural and molecular biology , 20(12):1352–1357, 2013.
- 2[2] X. Cheng. Random Matrices in High-dimensional Data Analysis . Ph D thesis, Princeton University, 2013.
- 3[3] Frank, J. Three-Dimensional Electron Microscopy of Macromolecular Assemblies: Visualization of Biological Molecules in Their Native State . Oxford, 2006.
- 4[4] Goncharov, A. B. Integral geometry and three-dimensional reconstruction of randomly oriented identical particles from their electron microphotos. Acta Applicandae Mathematicae , 11:199–211, 1988.
- 5[5] Herman, G. T. Fundamentals of Computerized Tomography: Image Reconstruction from Projections . Springer, London, UK, 2nd edition, 2009.
- 6[6] Iudin, A., Korir, P. K., Salavert-Torres, J., Kleywegt, G. J., and Patwardhan, A. Empiar: a public archive for raw electron microscopy image data. Nature Methods , 13(5):387–388, 2016.
- 7[7] Z. Kam. The reconstruction of structure from electron micrographs of randomly oriented particles. J. Theor. Biol. , 82(1):15–39, 1980.
- 8[8] Lee, C. and Mac Kinnon, R. Structures of the Human HCN 1 Hyperpolarization-Activated Channel. Cell , 168(1-2):111–120.e 11, 2017.