# Robust ab initio solution of the cryo-EM reconstruction problem at low
resolution with small data sets

Aaditya V. Rangan and Leslie Greengard Center for Computational Mathematics,
Flatiron Institute and Courant Institute, NYU

(July, 2023)

###### Abstract

Single particle cryo-electron microscopy has become a critical tool in
structural biology over the last decade, able to achieve atomic scale
resolution in three dimensional models from hundreds of thousands of (noisy)
two-dimensional projection views of particles frozen at unknown orientations.
This is accomplished by using a suite of software tools to (i) identify
particles in large micrographs, (ii) obtain low-resolution reconstructions,
(iii) refine those low-resolution structures, and (iv) finally match the
obtained electron scattering density to the constituent atoms that make up the
macromolecule or macromolecular complex of interest.

Here, we focus on the second stage of the reconstruction pipeline: obtaining a
low resolution model from picked particle images. Our goal is to create an
algorithm that is capable of ab initio reconstruction from small data sets (on
the order of a few thousand selected particles). More precisely, we propose an
algorithm that is robust, automatic, and fast enough that it can potentially
be used to assist in the assessment of particle quality as the data is being
generated during the microscopy experiment.

##  1 Introduction

In single-particle cryo-electron microscopy (cryo-EM), a purified preparation
of particles (proteins or macro-molecular complexes) is frozen in a thin sheet
of ice, with each particle held in an unknown orientation. Using a weak
scattering approximation, the image obtained can be viewed as a noisy
projection of the particle’s electron scattering density, which we denote by
${{F}}({\boldsymbol{x}})$, with optical aberrations modeled by what is called
a contrast transfer function (CTF) [1, 2]. The task at hand is to take the
individual particle images, determine their orientations, and recover a high
resolution characterization of ${{F}}({\boldsymbol{x}})$. This is typically
accomplished through the Fourier slice theorem [3, 4], which states that the
two-dimensional Fourier transform of each projection image is a slice through
${\hat{F}}({\boldsymbol{k}})$, the three-dimensional Fourier transform of
${{F}}({\boldsymbol{x}})$ (modulated by the CTF).

The specific slice obtained from a single image is determined by the
orientation of the individual particle. With many particles at many distinct
orientations, ${\hat{F}}({\boldsymbol{k}})$ is well-sampled and
${{F}}({\boldsymbol{x}})$ can be obtained through the inverse Fourier
transform. Thus, if the orientations are known, the reconstruction problem is
linear and easily solved. The orientations, however, are not known, and the
reconstruction problem is both nonlinear and non-convex. In practice, not only
are the orientations unknown, but the particles must be correctly centered in
order to invoke the Fourier slice theorem. The unknown displacement vector
must also be learned for each particle image. We postpone a proper discussion
of the mathematics to section 2 and the appendices (Supporting Information).

From an experimental perspective, the initial stages of cryo-EM involve some
hours of imaging, resulting in the collection of multiple micrographs, each
containing a number of particle images. After sufficiently many micrographs
are collected, the next steps involve particle picking (identifying individual
particles in the micrographs), and the construction of a low- to medium-
resolution molecule for further analysis and refinement. The latter step
typically involves class averaging - finding particles that appear to have the
same orientation and taking the average in order to improve the signal-to-
noise ratio, but some algorithms proceed from the raw particle images
themselves (as we do here).

Despite recent advances in sample preparation, imaging and analysis, the cost
of data collection can still be quite significant, and the availability of
microscope time is a limited resource. In order for the picked-particle images
to be useful, they must meet certain criteria: each should contain a single
well-centered particle with a minimum of clutter, and the full set must cover
a wide variety of viewing angles, sufficient to determine the full three-
dimensional structure. Assuming that these criteria are met, these images
(hopefully the least noisy and cleanest projection views) can then be used in
downstream analyses. If the sample or the images are not of sufficient
quality, it may be necessary to terminate the current experiment and repeat
with a new sample, hoping to produce higher-quality particle images the next
time.

To reduce the time and cost of such experiments, there are many quality-
control tasks that would be useful to carry out as early as possible. These
include:

  1. 1.

determining whether any picked-particle image contains a reasonably well-
centered and isolated molecule.

  2. 2.

estimating the viewing angles of the images (and checking that the range of
viewing angles is sufficient to reconstruct the molecule).

  3. 3.

estimating the displacements/shifts required to more precisely center the
images.

  4. 4.

estimating the image quality and the number of particles that will be needed
to reconstruct the desired molecule at the desired resolution.

  5. 5.

estimating the molecular structure itself at low-resolution (which feeds back
to the estimates 1-4)

In current practice, many of these issues are addressed after carrying out
class averaging, which typically requires dozens of particles per class and
tens of thousands of picked particles, although a number of reconstruction
algorithms work directly on the raw data.

In this paper, we describe an algorithm which can serve as a complement to
more standard strategies. Our algorithm operates on a relatively small set of
picked-particle images (say, a few thousand), and attempts to directly
reconstruct (without class averaging) a low-resolution version of the imaged
molecule, along with estimates of the viewing-angles, shifts and quality of
the images involved. As we demonstrate below, our method is significantly
faster, more accurate and more reliable than many existing de-novo
reconstruction algorithms. We hope that it is robust enough to permit the
assessment of data quality as the experiment is being carried out. Moreover,
our algorithm runs quickly and automatically with no ‘tuning’ required,
allowing for multiple small pools of approximately 1000 picked-particles to be
independently assessed without supervision (each corresponding to $\sim 30$
micrographs). This would help set the stage for more robust statistical
validation techniques such as bootstrap or jack-knife estimation.

###### Remark 1

We should emphasize that the starting point of our experiments here is a
collection of picked particle images (obtained from the EMPIAR database), not
the raw movies or entire micrographs that are sometimes deposited as well.
Particle picking is itself a non-trivial task. Restricting our attention to
molecules for which particles are curated allows us to focus exclusively on
the problem of reconstruction from a small data set.

At its core, our algorithm is a modification of the simplest expectation-
maximization method. We use the framework of classical iterative refinement -
alternating between a reconstruction step and an alignment step. In the
reconstruction step, we use the current estimate of the viewing angles for
each particle/image and their displacements to approximate the molecule. In
the alignment step, we use the current estimate of the molecule to better fit
the viewing angles and displacements for each image. We will refer to this
approach as alternating minimization (AM).

It is well-known that standard AM does poorly in this setting; an
initialization with randomized viewing angles typically fails to converge to a
useful low-resolution approximation of the molecule. To address these issues
and to guide AM to a better structure, we introduce two changes: (i) we
compress and denoise the individual images by computing their two-dimensional
Fourier transform on a polar grid and focusing on carefully selected rings in
the transform domain with maximal information content, and (ii) we increase
the entropy of the viewing-angle distribution used during the reconstruction
step (in a manner to be described in detail). These modifications both reduce
the computational cost and improve the robustness of AM in the low-resolution
low-image-number regime.

In the numerical experiments below, we show that our modified method
(alternating minimization with entropy maximization and principal modes or
‘EMPM’) performs surprisingly well on several published data sets. Our
reconstructions are often on par with the best results one could expect from
an ‘oracle-informed’ AM (i.e., standard AM starting with the ‘ground-truth’
viewing-angles and displacements for each image). Moreover, our results are
similar to (or better than) what one might expect from a reconstruction using
‘optimized’ Bayesian inference (i.e., using a ground-truth informed estimated
noise level). Importantly, our strategy is efficient, running about twenty
times faster than standard Bayesian inference or the de-novo reconstruction in
the widely used, open source package Relion [5, 6, 7]. We also make direct
comparisons with two other open source methods: the neural network approach in
cryoDRGN [8] and the recent low-resolution reconstruction pipeline in Aspire
[9]. We have not made a direct comparison with the stochastic gradient descent
approach in the GPU-based cryoSPARC package [10].

In the remainder of this paper, we briefly summarize the EMPM method and
demonstrate its performance on a variety of data sets. The details of our
method, as well as some additional motivating examples, are deferred to the
Supporting Information. It is currently implemented for multi-core CPUs, but
we believe it should be straightforward to implement on GPUs as well, and to
integrate into any of the existing software packages for cryo-EM.

##  2 Overview of the EMPM iteration

Our general goal is to solve the following problem: given a set of $N_{A}$
picked particle images (2D projections) from a handful of micrographs (each
with their own CTF), estimate the unknown viewing-angles $\tau_{j}\in SO(3)$
and displacements $\delta_{j}\in\mathbb{R}^{2}$ for each image, as well as the
underlying volume ${\hat{F}}({\boldsymbol{k}})$, where
${\boldsymbol{k}}\in\mathbb{R}^{3}$ denotes a point in the Fourier transform
domain. In addition to determining the volume ${\hat{F}}({\boldsymbol{k}})$,
we would also like an estimate of the viewing-angle distribution $\mu_{\tau}$
and displacement distribution $\mu_{{\boldsymbol{\delta}}}$ (taken across the
images). If the estimates of the volume ${\hat{F}}({\boldsymbol{k}})$ and the
distributions $\mu_{\tau}$, $\mu_{{\boldsymbol{\delta}}}$ are accurate
estimates of the ground truth, they will clearly be useful for quality
control, determining (for example) whether the experimental conditions are
suitable for larger scale data collection and subsequent high resolution
reconstruction.

The main tool we employ is a form of alternating minimization (AM), which
carries out some number of iterations of the following two step procedure:

Reconstruction:

Given the estimated viewing-angles $\\{\tau_{j}^{{\text{\tiny estim}}}\\}$ and
displacements $\\{{\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}\\}$ for
$N_{A}$ 2d picked-particle images $\\{{\hat{A}}_{j}\\}$, reconstruct the 3D
molecular volume ${\hat{F}}^{{\text{\tiny estim}}}$ (using the Fourier slice
theorem and least squares inversion of the Fourier transform).

Alignment:

Given the reconstruction of the 3D molecule ${\hat{F}}^{{\text{\tiny
estim}}}$, update the estimates $\\{\tau_{j}^{{\text{\tiny estim}}}\\}$ and
$\\{{\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}\\}$.

In a standard AM scheme, the reconstruction step solves for
${\hat{F}}^{{\text{\tiny estim}}}$ using a least-squares procedure, and the
alignment step assigns the $\\{\tau_{j}^{{\text{\tiny estim}}}\\}$ and
$\\{{\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}\\}$ using some kind of
maximum-likelihood procedure (see sections A.18 and A.22 in the Supporting
Information). While there are many more sophisticated strategies for molecular
reconstruction (see, for example, [10, 8]), we focus on AM for a few reasons:
first and foremost is its simplicity, transparency and ease of automation. AM
has modular steps, few parameters, and is easily interpretable.

At the same time, however, standard AM does have serious limitations. First,
it can be computationally expensive, particularly if a broad range of
frequencies and/or translations are considered in the alignment step. Second,
AM can fail to converge to the correct solution, especially during the initial
iterations, when the initial guess is far from the true molecule.

To address these two issues, we propose a modified version of AM which we
refer to as alternating minimization with entropy-maximization and principal
modes (EMPM). The two main modifications we introduce are:

Maximum entropy alignment:

We use both the traditional maximum likelihood alignment strategy, which finds
the best estimate for each image from the current molecular reconstruction,
and its opposite. That is, we also use a ‘maximum-entropy’ alignment strategy:
for each uniformly chosen viewing angle, we find the experimental image which
best matches the corresponding projection from the current structure. The
subsequent reconstruction step makes use of this uniform distribution of
viewing angles whether or not it corresponds to the true distribution. Use of
this counter-intuitive step is important in stabilizing our low resolution
reconstruction from a small set of images (see section A.22 in the Supporting
Information).

Principal-mode projection:

Each image is transformed onto a polar grid in Fourier space. The data on each
circle in the transform domain at a fixed modulus $|{\boldsymbol{k}}|$ will be
referred to as a Fourier ring. Rather than considering all Fourier rings when
comparing 2D projections and image data, we consider only a handful of
carefully-chosen combinations of such rings, which we refer to as principal
modes. This yields significant compression of the data (and computational
efficiency) while retaining information useful for alignment (see section A.13
in the Supporting Information).

In the EMPM pipeline, we first estimate the principal-modes using the images
themselves. Using these ‘empirical’ principal modes, we proceed with AM,
alternating between maximum-likelihood and maximum-entropy alignment. Once a
preliminary volume estimate has been obtained, we use that estimate to
recalculate the principal modes. We then restart AM using the new principal
images and, again, alternate between maximum-likelihood and maximum-entropy
alignment until we have a final volume estimate. The algorithm is summarized
below (see Supporting Information for a more detailed description):

  1. 1.

Use the experimental images ${\hat{A}}_{j}$ to calculate ‘empirical’ principal
modes $U^{{\text{\tiny empir}}}$.

  2. 2.

Compute the principal images, defined to be the projection of the images onto
principal modes, $P(U^{{\text{\tiny empir}}},{\hat{A}}_{j})$.

  3. 3.

Set random initial viewing angles $\tau_{j}^{{\text{\tiny empir}}}$ and set
initial displacements ${\boldsymbol{\delta}}_{j}^{{\text{\tiny empir}}}$ to
$0$.

  4. 4.

Carry out 32 iterations of EMPM on the principal images.

    1. (a)

Reconstruct ${\hat{F}}^{{\text{\tiny empir}}}$ using $\tau_{j}^{{\text{\tiny
empir}}}$ and ${\boldsymbol{\delta}}_{j}^{{\text{\tiny empir}}}$.

    2. (b)

Use maximum-likelihood alignment to update $\tau_{j}^{{\text{\tiny empir}}}$
and ${\boldsymbol{\delta}}_{j}^{{\text{\tiny empir}}}$.

    3. (c)

Reconstruct ${\hat{F}}^{{\text{\tiny empir}}}$ using $\tau_{j}^{{\text{\tiny
empir}}}$ and ${\boldsymbol{\delta}}_{j}^{{\text{\tiny empir}}}$.

    4. (d)

Use maximum-entropy alignment to update $\tau_{j}^{{\text{\tiny empir}}}$ and
${\boldsymbol{\delta}}_{j}^{{\text{\tiny empir}}}$.

  5. 5.

Use ${\hat{F}}^{{\text{\tiny empir}}}$ to calculate ‘estimated’ principal
modes $U^{{\text{\tiny estim}}}$.

  6. 6.

Set initial viewing angles $\tau_{j}^{{\text{\tiny
estim}}}=\tau_{j}^{{\text{\tiny empir}}}$ and displacements
${\boldsymbol{\delta}}_{j}^{{\text{\tiny
estim}}}={\boldsymbol{\delta}}_{j}^{{\text{\tiny empir}}}$.

  7. 7.

Carry out 32 iterations of EMPM on the new principal images $P(U^{{\text{\tiny
estim}}},{\hat{A}}_{j})$.

    1. (a)

Reconstruct ${\hat{F}}^{{\text{\tiny estim}}}$ using $\tau_{j}^{{\text{\tiny
estim}}}$ and ${\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}$.

    2. (b)

Use maximum-likelihood alignment to update $\tau_{j}^{{\text{\tiny estim}}}$
and ${\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}$.

    3. (c)

Reconstruct ${\hat{F}}^{{\text{\tiny estim}}}$ using $\tau_{j}^{{\text{\tiny
estim}}}$ and ${\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}$.

    4. (d)

Use maximum-entropy alignment to update $\tau_{j}^{{\text{\tiny estim}}}$ and
${\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}$.

Note that, each time we update the displacements
${\boldsymbol{\delta}}_{j}^{{\text{\tiny empir}}}$ and
${\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}$, we bound the maximum
permitted displacement magnitude by ${{\boldsymbol{\delta}}_{\max}}$
(typically ${{\boldsymbol{\delta}}_{\max}}\sim 0.1$ in our dimensionless
units). In our studies with small amounts of moderately noisy experimental
data, we have observed that the EMPM strategy is more effective than Bayesian
inference or maximum-likelihood alignment alone - even with probabilistic
alignment where a distribution of viewing angles is assigned to each image.
Roughly speaking, we believe this is due the fact that the incorrect viewing
angles assigned in the maximum-entropy step are washed out by destructive
interference, while the correctly assigned ones reinforce the true structure.
One simple model problem which illustrates these points is discussed in the
Supporting Information (see sections B and C). It is worth repeating, however,
that we only propose this approach for obtaining a low-resolution initial
guess. Existing refinement methods are still crucial for high resolution.

Once we have run our EMPM algorithm, we have the estimated volume
${\hat{F}}^{{\text{\tiny estim}}}({\boldsymbol{k}})$ as output. Using
${\hat{F}}^{{\text{\tiny estim}}}({\boldsymbol{k}})$ as an approximation of
${\hat{F}}({\boldsymbol{k}})$, we can readily compute other quantities that
may be of interest, such as estimates of the true viewing angle distribution
in the data set and other types of correlations. We will use these quantities
in our numerical experiments below to compare the results of EMPM to other
methods, such as Relion’s de-novo molecular reconstruction and full Bayesian
inference.

###  2.1 Templates, correlations, and ranks

For each centered molecular reconstruction ${\hat{F}}$, and set of Euler
angles $\tau=(\alpha,\beta,\gamma)$, we define the inner product ${\cal X}$
and normalized inner product $\hat{\cal X}$ of the reconstruction with the
“ground-truth” model ${\hat{F}}^{{\text{\tiny true}}}$ by

$\displaystyle{\cal X}(\tau;{\hat{F}},{\hat{F}}^{{\text{\tiny true}}})$
$\displaystyle=$ $\displaystyle\langle{\cal
R}(\tau)\circ{\hat{F}},{\hat{F}}^{{\text{\tiny true}}}\rangle$ (1)
$\displaystyle\hat{\cal X}(\tau;{\hat{F}},{\hat{F}}^{{\text{\tiny true}}})$
$\displaystyle=$ $\displaystyle\frac{{\cal
X}(\tau;{\hat{F}},{\hat{F}}^{{\text{\tiny
true}}})}{\|{\hat{F}}\|\cdot\|{\hat{F}}^{{\text{\tiny true}}}\|}\,.$ (2)

Here $(\alpha,\beta)$ denotes the viewing angle of an image with $\alpha$ the
azimuthal angle and $\beta$ the polar angle. $\gamma$ is the final in-plane
rotation (see section A.6). ${\cal R}(\tau)\circ{\hat{F}}$ is the volume
obtained from ${\hat{F}}$ using the rotation (the element of $SO(3)$) defined
by $\tau$. We define the volumetric correlation ${\cal
Z}({\hat{F}},{\hat{F}}^{{\text{\tiny true}}})$ by

${\cal Z}({\hat{F}},{\hat{F}}^{{\text{\tiny true}}})=\max_{\tau}\hat{\cal
X}(\tau;{\hat{F}},{\hat{F}}^{{\text{\tiny true}}})\text{.}$ (3)

This measurement ${\cal Z}({\hat{F}},{\hat{F}}^{{\text{\tiny true}}})$
represents the optimal correlation (over all alignments) between the two
volumes ${\hat{F}}$ and ${\hat{F}}^{{\text{\tiny true}}}$.

In order to compare our results with Relion’s de-novo reconstruction, we will
also need to calculate the ‘masked’ correlation ${\cal Z}^{{\text{\tiny
mask}}}({{F}},{{F}}^{{\text{\tiny true}}})$, defined as:

$\displaystyle{\cal X}^{{\text{\tiny mask}}}(\tau;{{F}},{{F}}^{{\text{\tiny
true}}})$ $\displaystyle=$ $\displaystyle\langle{\cal
R}(\tau)\circ{{F}}\cdot{{G}},{{F}}^{{\text{\tiny true}}}\rangle$ (4)
$\displaystyle\hat{\cal X}^{{\text{\tiny
mask}}}(\tau;{{F}},{{F}}^{{\text{\tiny true}}})$ $\displaystyle=$
$\displaystyle\frac{{\cal X}^{{\text{\tiny
mask}}}(\tau;{{F}},{{F}}^{{\text{\tiny
true}}})}{\|{{F}}\cdot{{G}}\|\cdot\|{{F}}^{{\text{\tiny true}}}\|}$ (5)
$\displaystyle{\cal Z}^{{\text{\tiny mask}}}({{F}},{{F}}^{{\text{\tiny
true}}})$ $\displaystyle=$ $\displaystyle\max_{\tau}\hat{\cal X}^{{\text{\tiny
mask}}}(\tau;{{F}},{{F}}^{{\text{\tiny true}}})\text{,}$ (6)

where ${{G}}$ is a sharp spherical mask (i.e., the indicator-function of a
ball centered at the origin) with a radius determined by the support of
${{F}}^{{\text{\tiny true}}}$.

We will also need to define inner products and correlations of single images
and volumetric slices. For this, given a Fourier space volume ${\hat{F}}$, we
generate noiseless 2D projections (in the Fourier domain) that incorporate the
CTF, which we refer to as templates. We denote these by
${\hat{S}}(\tau,{\boldsymbol{\delta}};CTF;{\hat{F}})$, where $\tau$ is the
specific set of Euler angles and ${\boldsymbol{\delta}}$ is the displacement.
When assuming the displacement is zero, we will write these templates as
${\hat{S}}(\tau;CTF;{\hat{F}})$, Formally, we have

$\displaystyle{\hat{S}}(\tau;CTF;{\hat{F}})$ $\displaystyle:=$ $\displaystyle
CTF({\boldsymbol{k}})\odot\left[{\hat{Sl}}\circ{\cal
R}(\tau)\circ{\hat{F}}\right]({\boldsymbol{k}})\text{.}$ (7)

where $CTF({\boldsymbol{k}})$ is the given pointwise value of the CTF in
Fourier space and ${\hat{Sl}}$ denotes taking the equatorial slice in the
Fourier domain. These CTF-modulated templates are used to calculate the inner
products:

$\displaystyle{\cal
X}(\tau,{\boldsymbol{\delta}};{\hat{A}}_{j};CTF_{j};{\hat{F}})$
$\displaystyle=$ $\displaystyle\langle\ {\hat{S}}(\tau;CTF_{j};{\hat{F}})\ ,\
{\cal T}(-{\boldsymbol{\delta}},{\boldsymbol{k}})\odot{{\hat{A}}_{j}}\
\rangle\text{,}$ (8)

where ${\cal
T}(+{\boldsymbol{\delta}},{\boldsymbol{k}})=e^{i{\boldsymbol{\delta}}\cdot{\boldsymbol{k}}}$
is the Fourier signature of displacement by the vector
${\boldsymbol{\delta}}$. These are used, in turn, to calculate the
correlations:

$\hat{\cal
X}(\tau,{\boldsymbol{\delta}};{\hat{A}}_{j};CTF_{j};{\hat{F}})=\frac{{\cal
X}(\tau,{\boldsymbol{\delta}};{\hat{A}}_{j};CTF_{j};{\hat{F}})}{\|{\hat{S}}(\tau;CTF_{j};{\hat{F}})\|\cdot\|{\hat{A}}_{j}\|}\text{.}$

For convenience, we generally reduce this to a function of the viewing
direction $\hat{\boldsymbol{k}}=\left(\alpha,\beta\right)$, the first two
components of $\tau$, alone:

${\cal
Z}(\alpha,\beta;{\hat{A}}_{j};CTF_{j};{\hat{F}})=\max_{\gamma,{\boldsymbol{\delta}}}\hat{\cal
X}(\tau,{\boldsymbol{\delta}};{\hat{A}}_{j};CTF_{j};{\hat{F}})\text{.}$ (9)

The parameters $\gamma$ and ${\boldsymbol{\delta}}$ are associated with a
particular choice of $\alpha$ and $\beta$ implicitly through the
$\operatorname*{arg\,max}$. When the context is clear, we will often abuse
notation and refer to the correlations in (9) by ${\cal
Z}(\alpha,\beta;{\hat{A}}_{j})$.

Finally, for each image ${\hat{A}}_{j}$ and a set of $N_{S}$ templates defined
by the viewing directions $(\alpha_{l},\beta_{l}),\ l=1,\dots,N_{S}$, we can
order the values

${\cal Z}(\alpha_{l},\beta_{l};{\hat{A}}_{j};CTF_{j};{\hat{F}}),\
l=1,\dots,N_{S}.$ (10)

This yields what we will refer to as the template rank for the image.

Conversely, for each template defined by $(\alpha_{l},\beta_{l})$, one can
order the values

${\cal Z}(\alpha_{l},\beta_{l};{\hat{A}}_{j};CTF_{j};{\hat{F}}),\
j=1,\ldots,N_{A}.$ (11)

This defines the image rank for that template. Using the terminology of linear
algebra, we can define the $N_{S}\times N_{A}$ matrix ${\cal Z}$ with

${\cal Z}_{lj}={\cal Z}(\alpha_{l},\beta_{l};{\hat{A}}_{j}).$

Template-ranking then corresponds to sorting the rows of ${\cal Z}$ while
image-ranking corresponds to sorting the columns.

##  3 Results

We will consider a number of examples in this paper, drawn from the electron-
microscopy image archive (EMPIAR) - namely, trpv1, the capsaicin receptor (the
transient receptor potential cation channel subfamily V member 1)
[EMPIAR-10005, emd-5778], rib80s, the Plasmodium falciparum 80S ribosome bound
to emetine [EMPIAR-10028, emd-2660], ISWINCP, the ISWI-nucleosome complex
[EMPIAR-10482, emd-9718], ps1, the pre-catalytic spliceosome [EMPIAR-10180],
MlaFEDB, a bacterial phospholipid transporter [EMPIAR-10536, emd-22116],
LetB1, the lipophilic envelope-spanning tunnel B protein from E. coli
[EMPIAR-10350, emd-20993], TMEM16F, a Calcium-activated lipid scramblase and
ion channel in digitonin with calcium bound [EMPIAR-10278, emd-20244], and
LSUbl16dep, L17-Depleted 50S Ribosomal Intermediates, using an average of the
published major classes as a reference [EMPIAR-10076]. To limit the number of
picked particles, we carry out our numerical experiments on subsets of
$N_{A}=1024$ images, spanning a small number $N_{\text{\tiny CTF}}$ of
micrographs.

We set ${K}=48$ as the bandlimit for our low-resolution reconstruction,
corresponding to a ‘best’ possible resolution of about $20$Å. With these image
pools, we construct a variety of different ‘reference volumes’ for
${\hat{F}}({\boldsymbol{k}})$.

Ground truth:

First, we center and restrict the published volume to the band-limit ${K}$.
This defines our ‘ground truth’ volume ${\hat{F}}^{{\text{\tiny
true}}}({\boldsymbol{k}})$.

Oracle:

Next, we measure the correlations ${\cal
Z}(\tau;{\hat{A}}_{j};CTF_{j};{\hat{F}}^{{\text{\tiny true}}})$ for a
sufficiently sampled set of Euler angles $\tau$ and displacements and use
maximum-likelihood to align the images to ${\hat{F}}^{{\text{\tiny
true}}}({\boldsymbol{k}})$. This assigns ‘ground-truth’ viewing angles
$\tau_{j}^{{\text{\tiny true}}}$ and displacements
${\boldsymbol{\delta}}_{j}^{{\text{\tiny true}}}$ to each image. More
precisely, we start by calculating ${\cal Z}$ using a small displacement disc
of ${{\boldsymbol{\delta}}_{{\text{\tiny local}}}}=0.03$, and then iterate
multiple times, with ${{\boldsymbol{\delta}}_{\max}}=0.25$ until the
correlation between the reconstructed and true volumes plateaus (this usually
takes less than $8$ iterations). We then take the final viewing angles and
displacements to define $\tau_{j}^{{\text{\tiny true}}}$ and
${\boldsymbol{\delta}}_{j}^{{\text{\tiny true}}}$. Given these ground-truth-
based viewing angles and displacements, we use a least-squares reconstruction
to compute ${\hat{F}}^{{\text{\tiny Oracle}}}({\boldsymbol{k}})$.

Oracle-AM (OAM):

Starting with the ground-truth viewing angles and displacements, as well as
${\hat{F}}^{{\text{\tiny Oracle}}}({\boldsymbol{k}})$, we apply standard AM,
iterating until the resulting volume converges (this, again, usually requires
less than $8$ iterations). We define this final volume to be
${\hat{F}}^{{\text{\tiny OAM}}}({\boldsymbol{k}})$. We expect the volume
${\hat{F}}^{{\text{\tiny OAM}}}$ to be close to the best reconstruction one
could realistically achieve with standard AM.

Oracle-Bayesian (OBI):

Next, we apply Bayesian inference reconstruction (as described in [7, 11]).
This reconstruction depends on (i) the initial starting volume, (ii) the
estimated noise level $\sigma_{{\text{\tiny BI}}}$, (iii) the resolution ${K}$
and range of displacements ${{\boldsymbol{\delta}}_{\max}}$ considered, and
(iv) the error tolerance used when calculating the integrals required for
estimating the posterior described in [7]. We fix ${K}=48$ and
${{\boldsymbol{\delta}}_{\max}}=0.10$ and set the error tolerance to the same
value used in the rest of our numerical experiments (typically
$\epsilon_{{\text{\tiny tol}}}=0.01$). To create a reference volume associated
with the ground truth, we use ${\hat{F}}^{{\text{\tiny true}}}$ as the initial
starting volume. The resulting volume – denoted by ${\hat{F}}^{{\text{\tiny
OBI}}}({\boldsymbol{k}};\sigma_{{\text{\tiny BI}}})$ – will depend on the
noise-estimate $\sigma_{{\text{\tiny BI}}}$. Note that, when
$\sigma_{{\text{\tiny BI}}}\rightarrow 0$ we will have
${\hat{F}}^{{\text{\tiny OBI}}}({\boldsymbol{k}};\sigma_{{\text{\tiny
BI}}})\rightarrow{\hat{F}}^{{\text{\tiny OAM}}}({\boldsymbol{k}})$, and when
$\sigma_{{\text{\tiny BI}}}\rightarrow\infty$ we will have
${\hat{F}}^{{\text{\tiny OBI}}}({\boldsymbol{k}};\sigma_{{\text{\tiny BI}}})$
converge to a spherically-symmetric distribution. In our experience the
quality of the reconstructed volume depends quite sensitively on the choice of
$\sigma_{{\text{\tiny BI}}}$. To build a good noise-estimate, we scan over
several values of $\log(\sigma_{{\text{\tiny BI}}})\in[-5,\ldots,0]$, choosing
the $\sigma_{{\text{\tiny BI}}}$ for which ${\hat{F}}^{{\text{\tiny
OBI}}}({\boldsymbol{k}};\sigma_{{\text{\tiny BI}}})$ has the highest
correlation with ${\hat{F}}^{{\text{\tiny true}}}$ after the first iteration.
Fixing $\sigma_{{\text{\tiny BI}}}^{{\text{\tiny opt}}}$ to be this optimal
value, we define ${\hat{F}}^{{\text{\tiny
OBI}}}({\boldsymbol{k}}):={\hat{F}}^{{\text{\tiny
OBI}}}({\boldsymbol{k}},\sigma_{{\text{\tiny BI}}}^{{\text{\tiny opt}}})$
after iterating to convergence. Thus, ${\hat{F}}^{{\text{\tiny OBI}}}$ is
close to the best reconstruction possible using Bayesian inference, assuming
that the optimal estimated noise-level $\sigma_{{\text{\tiny BI}}}$ is known
ahead of time.

In addition to the above reference volumes, we also create several de-novo
molecular reconstructions. RDN and EMPM below are oracle-free, while for
Bayesian inference we use the optimal noise level estimate
$\sigma_{{\text{\tiny BI}}}^{{\text{\tiny opt}}}$ from above.

Bayesian Inference (BI):

We run the Bayesian inference reconstruction described in [7, 11] to produce
${\hat{F}}^{{\text{\tiny BI}}}({\boldsymbol{k}})$. As mentioned above, this
reconstruction depends on (i) the initial starting volume, (ii) the estimated
noise-level $\sigma_{{\text{\tiny BI}}}$, (iii) the resolution ${K}$ and the
range of displacements ${{\boldsymbol{\delta}}_{\max}}$ considered, (iv) the
error-tolerance used when computing the integrals in the posterior, and (v)
the number of iterations applied. We fix the noise level $\sigma_{{\text{\tiny
BI}}}^{{\text{\tiny opt}}}$ to be the value obtained in the Oracle-Bayesian
method above. We fix the resolution at ${K}=48$, set
${{\boldsymbol{\delta}}_{\max}}=0.10$, set the error-tolerance to the same
$\epsilon_{{\text{\tiny tol}}}$ used in the other methods, and set the number
of iterations at $16$. For the initial starting volume we assign each of the
available images a viewing angle randomly and uniformly, and then use a least-
squares reconstruction to generate a ‘random’ initial starting volume. While
random, this starting volume has roughly the same diameter and bandwidth as
the true molecule. We use the same procedure when initializing our strategy
below.

DRG:

We run the open-source neural-network-based algorithm cryoDRGN described in
[8], with the suggested default parameters for ab initio reconstruction.

ASP:

We run the open-source software package ASPIRE [9], which uses a common-lines-
based algorithm, with the suggested default parameters for ab initio
reconstruction.

RDN:

We run Relion’s de-novo molecular reconstruction algorithm (RDN) - a variant
of Bayesian inference - to produce ${\hat{F}}^{{\text{\tiny RDN}}}$. (For
speed, RDN uses only subsets of images in its internal iterations.) We use the
parameters listed in section D in the Supporting Information. We set the
parameter \--particle_diameter to be commensurate with the true molecular
diameter, rounded up to the nearest $50$Å.

EMPM:

We carry out the reconstruction using the strategy described in the Methods
section above, and discussed in more detail in the Supporting Information to
produce ${\hat{F}}^{{\text{\tiny EMPM}}}({\boldsymbol{k}})$. In addition to
setting the maximum frequency ${K}$ and displacement
${{\boldsymbol{\delta}}_{{\text{\tiny local}}}}$, many of the computed
quantities require a numerical tolerance $\epsilon_{{\text{\tiny tol}}}$ (see
section A). We fix ${K}=48$, and scan over $\epsilon_{{\text{\tiny
tol}}}\in[0.001,0.010]$ and ${{\boldsymbol{\delta}}_{{\text{\tiny
local}}}}\in[0.01,0.10]$. Our strategy begins with the same random initial
viewing-angles for each image as used for Bayesian inference. As we shall see
below (see, e.g., Fig. 3), our results are quite robust, provided that
$\epsilon_{{\text{\tiny tol}}}\leq 0.01$ and
${{\boldsymbol{\delta}}_{{\text{\tiny local}}}}\geq 0.01$ or so.

For trpv1, as well as rib80s and ISWINCP, the data sets are quite clean. The
first $N_{A}=1024$ picked-particle images are well-centered, devoid of clutter
and of high quality. Consequently, many of the reconstructions we attempted
were successful and RDN reconstruction was also able to achieve ${\cal
Z}^{{\text{\tiny mask}}}>0.5$. For ps1, LSUbl16dep, MlaFEDB, LetB1, and
TMEM16F, however, the individual images are not quite as clean. In the case of
ps1 and TMEM16F, many of the images are particularly noisy. For LSubl16dep
(and to a lesser extent ps1), a notable fraction of the images are poorly
centered. For MlaFEDB (and to a lesser extent TMEM16F), many of the images
appear to contain fragments of other particles. For LetB1, the viewing-angle
distribution for the images is highly concentrated around the equator. Our
reconstructions for these molecules are not as successful. Nevertheless, we
still recover volumes that are comparable in quality to Bayesian inference
with an optimal estimated noise-level, and usually better than RDN
reconstruction (all at a fraction of the computational time).

Figure 1:  In this figure, we plot level sets from the volumes
${\hat{F}}^{{\text{\tiny true}}}$, ${\hat{F}}^{{\text{\tiny Oracle}}}$,
${\hat{F}}^{{\text{\tiny OAM}}}$, and ${\hat{F}}^{{\text{\tiny OBI}}}$, as
well as one of the ${\hat{F}}^{{\text{\tiny EMPM}}}$ reconstructions (left to
right, respectively). From the top down, the panels correspond to trpv1,
rib80s, ISWINCP and ps1.

Results for four molecules are shown in Fig. 1. Focusing on trpv1 for the
moment, we see that ${\hat{F}}^{{\text{\tiny EMPM}}}$ captures the coarse
features of ${\hat{F}}^{{\text{\tiny true}}}$, and is not too different from
${\hat{F}}^{{\text{\tiny OAM}}}$ in terms of overall quality. Even though
${\hat{F}}^{{\text{\tiny EMPM}}}$ is a coarse estimate of
${\hat{F}}^{{\text{\tiny true}}}$, the reconstruction process has recovered a
significant amount of information regarding the viewing-angles of each image.
To illustrate this point, we compare the image-template correlations (see Eq.
9) estimated from ${\hat{F}}^{{\text{\tiny EMPM}}}$ with the ‘ground-truth’
image-template correlations measured from ${\hat{F}}^{{\text{\tiny true}}}$.
For each molecule of interest, three panels are presented in Fig. 2. The left-
most subplot is a scatterplot of the image-template correlations. That is, for
each image ${\hat{A}}_{j}$ and each template $(\alpha_{l},\beta_{l})$, we
compute a point in two-dimensions with coordinates

$({\cal Z}(\alpha_{l},\beta_{l};{\hat{A}}_{j};CTF_{j};{\hat{F}}^{{\text{\tiny
EMPM}}}),{\cal
Z}(\alpha_{l},\beta_{l};{\hat{A}}_{j};CTF_{j};{\hat{F}}^{{\text{\tiny
true}}})).$

With $N_{A}=1024$ images and $N_{T}=993$ templates, this yields a data set
with about one million points. The color in the heatmap corresponds to the
log-density of the joint-distribution at that location (see adjacent
colorbar). Note that there is a significant correlation between the estimated-
and ground-truth-image-template-correlations. In the middle subplot for each
molecule, we show a scatterplot of the template ranks for each image (see
(10)), rescaled to $[0,1]$ and aggregated over all images. In the right-most
subplot, we show a scatterplot of the image ranks for each template (see
(11)), rescaled to $[0,1]$ and aggregated over all templates. Once again, the
scatterplots are presented as heatmaps, with color corresponding to log-
density.

Figure 2:  In this figure we illustrate the correspondence between estimated
and ground-truth viewing angles for the reconstructions shown in Fig 1. Three
panels are shown for each molecule of interest. The left-most subplot is a
scatterplot of the image-template correlations. The middle subplot is a
scatterplot of the template ranks aggregated over all images, and the right-
most subplot is a scatterplot of the image ranks aggregated over all
templates. All are presented as heatmaps, with color corresponding to the log-
density of the joint-distribution. (See text for further explanation.)

This view of the data sheds light on the differences between maximum-
likelihood and maximum-entropy alignment. Recall that, in standard maximum-
likelihood alignment, one assigns $\tau_{j}^{{\text{\tiny est}}}$ (for a
particular image ${\hat{A}}_{j}$) to the template with the highest template
rank. By contrast, in maximum-entropy alignment, one does the reverse:
assigning the image with highest image rank to each sampled template
direction. It is interesting to note that the estimated and ground-truth image
ranks are typically more tightly correlated than the estimated and ground-
truth template ranks. We believe that this phenomenon contributes to the
success of our strategy. Indeed, when the estimated and ground-truth template
ranks are highly correlated (as for ISWINCP), the maximum-likelihood
alternating minimization method (i.e., Bayesian inference with noise-level
$0$) succeeds, performing roughly as well as our strategy. However, when the
estimated and ground-truth template ranks are poorly correlated (as for ps1),
then maximum-likelihoood alternating minimization does very poorly, often
converging to a volume that is worse than random.

The illustrations in Figs. 1 and 2 are obtained using just one trial of
${\hat{F}}^{{\text{\tiny EMPM}}}$. We have carried out extensive simulations
and these results are fairly typical, with qualitatively similar results with
repeated runs. For illustration, we show the values of the volumetric
correlations collected over multiple trials in Fig. 3. The values of ${\cal
Z}^{{\text{\tiny mask}}}({\hat{F}}^{{\text{\tiny
EMPM}}},{\hat{F}}^{{\text{\tiny true}}})$ shown in magenta are taken over a
range of values of $\epsilon_{{\text{\tiny tol}}}\in[0.001,0.010]$, and
${{\boldsymbol{\delta}}_{{\text{\tiny local}}}}\in[0.01,0.10]$. We did not see
a large correlation between either $\epsilon_{{\text{\tiny tol}}}$ or
${{\boldsymbol{\delta}}_{{\text{\tiny local}}}}$ and ${\cal Z}^{{\text{\tiny
mask}}}$, provided that $\epsilon_{{\text{\tiny tol}}}$ is sufficiently small
and ${{\boldsymbol{\delta}}_{{\text{\tiny local}}}}$ is sufficiently large
(i.e., $\epsilon_{{\text{\tiny tol}}}\leq 0.01$ and
${{\boldsymbol{\delta}}_{{\text{\tiny local}}}}\geq 0.01$). Note that, for
these data sets, our strategy reliably produces reconstructions that are
similar in quality to ${\hat{F}}^{{\text{\tiny OAM}}}$, and superior to the
results obtained with existing open-source packages (DRG, ASP, RDN). Finally,
we remark that a higher-quality low-resolution reconstruction can help improve
the results of subsequent refinement. As an example, we use Relion’s
relion_refine to produce a medium-resolution reconstruction from the first
$8192$ images of this dataset using two different initial-volumes. First, we
set the initial volume to Relion’s de-novo reconstruction
${\hat{F}}^{{\text{\tiny RDN}}}$ (referenced in blue in Fig. 3). Second, we
set the initial volume to the EMPM reconstruction ${\hat{F}}^{{\text{\tiny
EMPM}}}$ shown in Fig. 1. Fig. 4 illustrates the Fourier shell correlations
(FSCs) between the resulting medium-resolution reconstructions and the
published ${\hat{F}}^{{\text{\tiny true}}}$. Note that the higher-quality
initial volume ${\hat{F}}^{{\text{\tiny EMPM}}}$ results in a better medium-
resolution reconstruction. Moreover, as we shall discuss below, better low-
resolution reconstructions can be useful for other quality-control tasks in
the cryo-EM pipeline.

We have focused here on a subset of the proteins listed above. Addition
results are deferred to the Supporting material.

###### Remark 2

Of note is that the EMPM method works reasonably well even when the particle
orientations are highly non-uniform (see Fig. 15). One might expect that the
entropy-maximization step (which forces images to be assigned to each viewing
angle) would cause difficulties. For de novo low- to moderate-resolution
reconstructions, however, it appears that EMPM benefits from the prevention of
spurious clustering that can develop in pure AM iteration. At higher
resolution, we do not recommend the use of EMPM, as the entropy-maximization
step would result in a loss of resolution and is no longer necessary once a
good initial model has been created.

Figure 3:  In this figure we illustrate the aligned correlation ${\cal
Z}^{{\text{\tiny mask}}}$ between our various reconstructions and the ground-
truth for trpv1, rib80s, ISWINCP, and ps1 (see text). The different horizontal
bars indicate the values of ${\cal Z}^{{\text{\tiny mask}}}$ for different
reconstruction strategies (coded by color) with different trials from the same
strategy adjacent to one another (sorted by ${\cal Z}^{{\text{\tiny mask}}}$).
At the top of each panel (centered) is the “Oracle” low resolution
reconstruction for the indicated molecule. To the right in each panel, we show
an EMPM reconstruction and to the left, we show a reconstruction obtained
using one of the other pipelines which achieves poorer correlation scores.
Note that the ${\cal Z}^{{\text{\tiny mask}}}$ values do tend to correlate
with the image reconstruction quality. trpv1 has a highly non-uniform
distribution of ground-truth Euler angles. The distributions in the other
cases are more uniform. (See Fig. 15.)
Figure 4:  In this figure we show the Fourier shell correlation (FSC)
associated with a medium-resolution reconstructions produced using Relion’s
relion_refine applied to the first $N_{A}=8192$ images from the trpv1 dataset
(i.e., the first $N_{A}$ picked-particles from the MRC_1901 file). We run this
reconstruction twice: once using as initial volume the RDN reconstruction
${\hat{F}}^{{\text{\tiny RDN}}}$ obtained using the first $1024$ images and
once using as initial volume the EMPM reconstruction ${\hat{F}}^{{\text{\tiny
EMPM}}}$. We plot the FSC curve obtained by comparing the reconstructed volume
with the published ground-truth ${\hat{F}}^{{\text{\tiny true}}}$. The blue
(lower) curve is from the RDN initialization and the magenta (higher) curve is
from the EMPM initialization.

###### Remark 3

We have implemented our strategy in Matlab, using standard (double precision)
libraries. Even for our relatively naive implementation, the total computation
time for each of our trials is $30-50$ minutes on a desktop workstation. By
comparison, the de-novo reconstruction method within Relion typically requires
$20-25$ hours on the same machine, while a trial of Bayesian inference (also
in Matlab) requires $36-48$ hours. It is also worth noting that, when
constructing ${\hat{F}}^{{\text{\tiny BI}}}$, we chose the optimal noise-
estimate $\sigma_{{\text{\tiny BI}}}^{{\text{\tiny opt}}}$. Scanning over
multiple values of $\sigma_{{\text{\tiny BI}}}$ would require additional
computation time. In summary, we believe that our strategy has the potential
to be extremely efficient, especially when properly optimized.

###  3.1 On-the-fly particle rejection

Focusing on the ps1 data set, we now consider whether our algorithm can be
used to curate particles from small samples. As shown in Figs. 2 and 11, the
estimated image-template correlations ${\cal
Z}(\alpha,\beta;{\hat{A}}_{j};CTF_{j};{\hat{F}}^{{\text{\tiny EMPM}}})$ are
often quite strongly correlated with the true image-template correlations
${\cal Z}(\alpha,\beta;{\hat{A}}_{j};CTF_{j};{\hat{F}}^{{\text{\tiny
true}}})$. Motivated by this observation, we define:

${\cal Q}_{j}^{{\text{\tiny EMPM}}}=\max_{t}{\cal
Z}(\alpha_{t},\beta_{t};{\hat{A}}_{j};CTF_{j};{\hat{F}}^{\text{\tiny
EMPM}}),\quad{\cal Q}_{j}^{{\text{\tiny true}}}=\max_{t}{\cal
Z}(\alpha_{t},\beta_{t};{\hat{A}}_{j};CTF_{j};{\hat{F}}^{\text{\tiny true}})$

to be the maximum image-template correlation for each image. We can use ${\cal
Q}_{j}^{{\text{\tiny EMPM}}}$ as a proxy for the ‘true image-quality’ ${\cal
Q}_{j}^{{\text{\tiny true}}}$. We see from Figs. 2 and 11 that ${\cal Q}_{j}$
typically lies in the interval $[0.1,0.4]$. It is reasonable to conjecture
that those images with a low estimated quality do not align well to the
reconstructed molecule and could be discarded when attempting to reconstruct
the molecule. To investigate this hypothesis, we perform a numerical
experiment with the first $1024\cdot 8=8192$ images in the EMPIAR data set,
dividing these images into $8$ batches of $N_{A}=1024$. For each batch we
produce a reconstruction, using our strategy above. The ${\cal
Z}^{{\text{\tiny mask}}}$ for these single-batch runs are shown as magenta-
dots in the ‘IND’ column of Fig. 5. The average ${\cal Z}^{{\text{\tiny
mask}}}$ for these single-batch runs is $\sim 0.67$ (magenta line), with
relatively little variation.

Next, we use the estimated image-quality ${\cal Q}_{j}$ to sort the images
into octiles (each with $1024$ images). Those with the lowest quality scores
define the first octile and those with the highest quality scores define the
eighth octile. We now perform a reconstruction for each octile separately,
with the ${\cal Z}^{{\text{\tiny mask}}}$ shown in light-pink-dots. Note that
the first few octiles produce reconstructions that are quite a bit worse than
the typical reconstructions shown in the 1K column. It is interesting to note
that, for this image pool, using exclusively the highest quality octile
results in a poor reconstruction, since the highest quality images are
concentrated around a small subset of viewing-angles. The result obtained
using all $N_{A}=8192$ images together is shown as the dark-purple dot in the
‘1st’ column. We then show successive reconstructions using different quality
cut-offs. Using the top seven octiles ($7168$ images), we obtain the improved
${\cal Z}^{{\text{\tiny mask}}}$ shown as a dark-purple dot in the ‘2nd’
column. Using the top six octiles, we get the dark purple dot in the ‘3rd’
column, and so forth. In short, excluding the worst images improves the
reconstruction quality, but excluding too many images causes the quality to
drop, presumably as a result of losing the coverage from noisier projection
directions.

Figure 5:  In this figure we illustrate the potential for improving
reconstruction quality by curating the image-pool. For this example we use the
ps1 molecule, and consider the first $1024\cdot 8=8192$ images. First we
divide these images into batches of $N_{A}=1024$, and produce multiple
reconstructions (magenta, leftmost column). The average ${\cal
Z}^{{\text{\tiny mask}}}$ for these reconstructions is $\sim 0.67$ (magenta
line). Next, we use the individual estimates of image-quality (see text) to
divide the images into 8 ‘octiles’. We perform a reconstruction using the
images within each octile (light-pink). Note that those images of poorest
quality (i.e., in the first few octiles) produce an inferior reconstruction.
For this image-pool the images with the very best quality (i.e., in the final
octile) tend to be attributed to a small subset of viewing-angles; a
reconstruction using only those images is also of rather poor quality. We then
perform a reconstruction using all the images except those below a particular
octile (dark purple). Thus, the reconstruction using all $8192$ images is
shown in dark-purple within the ‘1st’ column. The reconstruction using the
$8192-1024=7168$ images which exclude those in the first octile is shown in
dark-purple within the ‘2nd’ column. Note that excluding the worst images
improves the reconstruction-quality. As the number of excluded images grows
the number of images used for the reconstruction shrinks; the dark-purple
circles are not much higher than the light-pink circles when the octile-index
is large.

###  3.2 Molecules with poor performance on small data sets

While our results can be considered moderately successful for many of the
molecules, our reconstructions for LetB1, LSUbl17dep and MlaFEDB using $1024$
particles are typically quite mediocre. For LetB1 and LSUbl17dep, different
runs of EMPM yielded reconstructions that were quite variable. We suspect that
this inconsistency stems from some underlying heterogeneity in the data set
itself (see [12]). For MlaFEDB, on the other hand, different trials of EMPM
yielded reconstructions ${\hat{F}}^{{\text{\tiny EMPM}}}$ that were not close
to ${\hat{F}}^{{\text{\tiny true}}}$, but were often quite close to one
another. This consistency across trials suggests that the structures observed
in the typical reconstruction might actually exist in the image pool we are
using. A recurring motif in the reconstructions is that of a cylindrical shape
with a disjoint structure off to one side. In this case, we suspect that the
image pool contains artifacts that are at least partially responsible for the
mediocre performance of EMPM (as well as BI, RDN, DRG, ASP) for this molecule.

##  4 Discussion

As demonstrated above, the EMPM iteration (a modification of alternating-
minimization) can produce a reasonable low-resolution reconstruction from a
small set of raw images. As seen in Figs. 3 and 12, it is more robust than
many of the currently used de-novo molecular reconstruction techniques, and at
least as accurate as Bayesian inference with an optimally chosen estimated
noise level. We did not make an exhaustive comparison with all available
packages, since our main purpose here is to provide a new approach which is
only relevant in the low resolution regime. We believe that the maximum
entropy step of EMPM accounts for its success, effectively preventing spurious
clustering of assigned viewing-angles. It should also be noted that the
principal mode reduction described here can be used to accelerate image
alignment and volume registration in more sophisticated and higher resolution
reconstruction algorithms.

EMPM also produces useful estimates of image quality and viewing angle
distribution – even when the image pool only contains a few thousand images.
These estimates could play a role in online quality-control: identifying poor
images to be discarded or possibly terminating imaging sessions that seem to
be yielding poor quality data.

Because EMPM is efficient (our MATLAB implementation requires less than 5% of
the computation time of RDN and Bayesian inference), it should make cross-
validation more compelling. Our reconstruction requires about 40 min., RDN
requires about 20 hrs, full Bayesian inference requires 40 hrs, cryoDRGN
requires 5 hrs and Aspire (the fastest) requires 15 min. (all on the same
laptop). Strong correlation between volumes reconstructed from different
subsets of images can reasonably be interpreted as an indication that the data
set is of high quality. On a more speculative note, the image curation
presented here uses only image-template correlations. It is possible that
image cross-correlations and the cross-correlations of the residuals of the
least-squares reconstruction step could also be valuable in settings involving
discrete or continuous heterogeneity.

Because EMPM does not require any fine-tuning of parameters, it should be
straightforward to incorporate into existing, widely-used, high-quality
software pipelines and we hope to do so in the near future. The code is
available from [13].

## References

  * [1] Jacqueline L S Milne, Mario J. Borgnia, Alberto Bartesaghi, Erin E H Tran, Lesley A. Earl, David M. Schauder, Jeffrey Lengyel, Jason Pierson, Ardan Patwardhan, and Sriram Subramaniam.  Cryo-electron microscopy - A primer for the non-microscopist.  FEBS Journal, 280:28–45, 2013.
  * [2] Yifan Cheng, Nikolaus Grigorieff, Pawel A. Penczek, and Thomas Walz.  A primer to single-particle cryo-electron microscopy.  Cell, 161:439–449, 2015.
  * [3] C. L. Epstein.  Introduction to the Mathematics of Medical Imaging.  SIAM, 2008.
  * [4] F. Natterer.  The Mathematics of Computerized Tomography.  SIAM, 2001.
  * [5] Dari Kimanius, Bjorn O Forsberg, Sjors HW Scheres, and Erik Lindehl.  Accelerated cryo-EM structure determination with parallelisation using GPUs in RELION-2.  eLife, 5:e18722, 2016.
  * [6] D. Kimanius, L. Dong, G. Sharov, T. Nakane, and Scheres S.H.W.  New tools for automated cryo-em single-particle analysis in relion-4.0.  Biochem, 478:4169–4185, 2021.
  * [7] Sjors H W Scheres.  A Bayesian view on cryo-EM structure determination.  J. Mol. Biol., 415:406–418, 2012.
  * [8] E.D. Zhong, T. Bepler, B. Berger, and J.H. Davis.  Cryodrgn: reconstruction of heterogeneous cryo-em structures using neural networks.  Nature Methods, 18:176–185, 2021.
  * [9] Aspire software, 2022.  http://spr.math.princeton.edu.
  * [10] A. Punjani, J.L. Rubinstein, D.J. Fleet, and M.A. Brubaker.  cryoSPARC: algorithms for rapid unsupervised cryo-EM structure determination.  Nat. Methods, 14:290–296, 2017.
  * [11] Sjors H W Scheres.  RELION: Implementation of a Bayesian approach to cryo-EM structure determination.  J. Struct. Biol., 180(3):519–530, 2012.
  * [12] J.H. Davis, Y.Z. Tan, B. Carragher, C.S. Potter, D. Lyumkis, and J.R. Williamson.  Modular assembly of the bacterial large ribosomal subunit.  Cell, 167:1610–1622, 2016.
  * [13] EMPM software, 2023.  https://github.com/adirangan/dir_cryoem.
  * [14] R Bracewell.  The Fourier Transform and Its Applications.  McGraw-Hill, 3rd edition, 1999.
  * [15] Zh. Zhao and A. Singer.  Rotationally invariant image representation for viewing direction classification in cryo-EM.  J. Struct. Biol., 186:153–166, 2014.
  * [16] A. Barnett, L. Greengard, A. Pataki, and M. Spivak.  Rapid solution of the cryo-EM reconstruction problem by frequency marching.  SIAM J. Imaging Sci., 10(3):1170–1195, 2017.
  * [17] Aaditya Rangan, Marina Spivak, Joakim Andén, and Alex Barnett.  Factorization of the translation kernel for fast rigid image alignment.  Inverse Problems, 36(2), 2020.
  * [18] Zh. Zhao, Y. Shkolnisky, and A. Singer.  Fast steerable principal component analysis.  IEEE Trans. Comput. Imaging, 2(1):1–12, 2016.
  * [19] A H Barnett, J F Magland, and L af Klinteberg.  A parallel non-uniform fast Fourier transform library based on an “exponential of semicircle” kernel, 2019.
  * [20] F J Sigworth.  A maximum-likelihood approach to single-particle image refinement.  J. Struct. Biol., 122(3):328–39, 1998.
  * [21] S H W Scheres, M Valle, P Grob, E Nogales, and J.-M. Carazo.  Maximum likelihood refinement of electron microscopy data with normalization errors.  J. Struct. Biol., 166(2):234–240, 2009.
  * [22] F. J. Sigworth, Doerschuk P.C., J.-M. Carazo, and S.H.W. Scheres.  An introduction to maximum-likelihood methods in Cryo-EM.  In Methods in Enzymology. Cryo-EM, Part B: 3D reconstruction, pages 263–294. Academic Press., 2010.
  * [23] Dmitry Lyumkis, Axel F. Brilot, Douglas L. Theobald, and Nikolaus Grigorieff.  Likelihood-based classification of cryo-EM images using FREALIGN.  J. Struct. Biol., 183(3):377–388, 2013.
  * [24] W. Baxter, R. Grassucci, Haixiao Gao, and J. Frank.  Determination of signal-to-noise ratios and spectral snrs in cryo-em low-dose imaging of molecules.  Journal of structural biology, 166 2:126–32, 2009.
  * [25] Aaditya V Rangan.  Radial-recombination for rigid rotational alignment of images and volumes.  Inverse Problems (submitted), 2022.
  * [26] Tamir Bendory, Ariel Jaffe, William Leeb, Nir Sharon, and Amit Singer.  Super-resolution multi-reference alignment.  Information and Inference: A Journal of the IMA, 11(2):533–555, 02 2021.
  * [27] Noam Janco and Tamir Bendory.  An accelerated expectation-maximization algorithm for multi-reference alignment.  IEEE Transactions on Signal Processing, 70:3237–3248, 2022.
  * [28] Amelia Perry, Jonathan Weed, Afonso S. Bandeira, Philippe Rigollet, and Amit Singer.  The sample complexity of multireference alignment.  SIAM Journal on Mathematics of Data Science, 1(3):497–517, 2019\.

##  Appendix A Supporting Information

Various aspects of the EMPM iteration are described here in greater detail. We
begin with the discretization of images and volumes in the Fourier domain, and
the techniques used for image alignment and volumetric reconstruction. Many of
the techniques used require some specification of a ‘global tolerance’
$\epsilon_{{\text{\tiny tol}}}$, which could be considered a user-specified
parameter. In practice, it appears to be sufficient to set it to $10^{-2}$,
corresponding to two digits of accuracy. In our experiments, we actually scan
over $\epsilon_{{\text{\tiny tol}}}\in[0.001,0.01]$ but find that the quality
of the reconstruction does not improve markedly for smaller values of
$\epsilon_{{\text{\tiny tol}}}$. In short, EMPM is not limited by the global
tolerance, but presumably by the experimental error and the structural
features of the nonconvex landscape over which we are trying to optimize.

###  A.1 Manipulating images in the Fourier domain

We use ${\boldsymbol{x}},{\boldsymbol{k}}\in\mathbb{R}^{2}$ to represent
spatial position and frequency, respectively. In polar-coordinates these
vectors can be represented as:

$\displaystyle{\boldsymbol{x}}$ $\displaystyle=$
$\displaystyle(x\cos\theta,x\sin\theta)$ (12) $\displaystyle{\boldsymbol{k}}$
$\displaystyle=$ $\displaystyle(k\cos\psi,k\sin\psi)\text{,}$ (13)

for appropriately chosen angles $\theta$ and $\psi$.

###### Definition 1

The Fourier transform of a two-dimensional function $A\in
L^{2}(\mathbb{R}^{2})$ is defined as

${\hat{A}}({\boldsymbol{k}}):=\iint_{\mathbb{R}^{2}}A({\boldsymbol{x}})\mathrm{e}^{-\mathrm{i}{\boldsymbol{k}}\cdot{\boldsymbol{x}}}\mathop{d{\boldsymbol{x}}}~{}.$
(14)

We recover $A$ from ${\hat{A}}$ using the inverse Fourier transform:

$A({\boldsymbol{x}})=\frac{1}{(2\pi)^{2}}\iint_{\mathbb{R}^{2}}{\hat{A}}({\boldsymbol{k}})\mathrm{e}^{+\mathrm{i}{\boldsymbol{k}}\cdot{\boldsymbol{x}}}\mathop{d{\boldsymbol{k}}}\text{.}$
(15)

The inner-product between two functions $A,B\in L^{2}(\mathbb{R}^{2})$ is
written as

$\langle
A,B\rangle=\iint_{\mathbb{R}^{2}}A({\boldsymbol{x}})^{\dagger}B({\boldsymbol{x}})\mathop{d{\boldsymbol{x}}}~{},$
(16)

where $z^{\dagger}$ is the complex conjugate of $z\in\mathbb{C}$. We will also
use Plancherel’s theorem [14],

$\langle
A,B\rangle=\frac{1}{(2\pi)^{2}}\langle{\hat{A}},{\hat{B}}\rangle~{},\qquad\forall
A,B\in L^{2}(\mathbb{R}^{2})~{}\text{.}$ (17)

We will ignore the factors of $2\pi$ associated with the Fourier transform
when they are not relevant.

We represent any given picked-particle image as a function $A\in
L^{2}(\mathbb{R}^{2})$, with values corresponding to the image intensity at
each location. As a consequence of Plancherel’s theorem, any inner-product
between $A$ and another image (or 2D function) $B$ can be calculated in either
real or frequency space, the latter of which is often more convenient when
aligning images to projections of a volume [15, 16, 17].

Abusing notation slightly, we will refer to $A({\boldsymbol{k}})$ and
${\hat{A}}({\boldsymbol{k}})$ in polar coordinates as:

$\displaystyle A(x,\theta)$ $\displaystyle:=$ $\displaystyle
A({\boldsymbol{x}})=A(x\cos\theta,x\sin\theta)$ (18)
$\displaystyle{\hat{A}}(k,\psi)$ $\displaystyle:=$
$\displaystyle{\hat{A}}({\boldsymbol{k}})={\hat{A}}(k\cos\psi,k\sin\psi)\text{.}$
(19)

With this notation each ${\hat{A}}(k,\psi)$ for fixed $k$ and
$\psi\in[0,2\pi)$ corresponds to a ‘ring’ in frequency space with radius $k$.

###  A.2 Rotation and translation

Using the notation above, a rotation ${\cal R}(\gamma)$ of an image $A$ by
angle $\gamma\in[0,2\pi)$ can be represented as:

${\cal R}(\gamma)\circ A(x,\theta):=A(x,\theta-\gamma)\text{.}$ (20)

Since rotation commutes with the Fourier transform, we have:

$\widehat{\left[{\cal R}(\gamma)\circ A\right]}(k,\psi)={\cal
R}(\gamma)\circ{\hat{A}}(k,\psi)={\hat{A}}(k,\psi-\gamma)\text{.}$ (21)

In this manner, a rotation of any image by $+\gamma$ can be represented as an
angular-shift of each image-ring by $\psi\rightarrow\psi-\gamma$.

Likewise, a translation ${\cal T}({\boldsymbol{\delta}})$ of an image $A$ by
the shift vector ${\boldsymbol{\delta}}\in\mathbb{R}^{2}$ can be represented
as:

${\cal T}({\boldsymbol{\delta}})\circ
A({\boldsymbol{x}}):=A({\boldsymbol{x}}-{\boldsymbol{\delta}})\text{.}$ (22)

In the Fourier domain, this action can be expressed as

$\widehat{\cal
T}({\boldsymbol{\delta}},{\boldsymbol{k}})\odot{\hat{A}}({\boldsymbol{k}})\text{,}$
(23)

where the translation kernel $\widehat{\cal
T}({\boldsymbol{\delta}},{\boldsymbol{k}})$ is

$\widehat{\cal
T}({\boldsymbol{\delta}},{\boldsymbol{k}}):=\mathrm{e}^{-\mathrm{i}{\boldsymbol{\delta}}\cdot{\boldsymbol{k}}}\text{.}$
(24)

###  A.3 Fourier-Bessel coefficients

Over the course of our molecular reconstructions, we calculate many image-
image inner products (see section 2.1). To ease the computational burden
associated with these calculations, it is convenient to use the Fourier-Bessel
basis [15, 18, 16, 17].

To define the Fourier-Bessel coefficients of an image, recall that each image-
ring ${\hat{A}}(k,\psi)$ (for fixed $k$) is a $2\pi$-periodic function of
$\psi$, which can itself be represented as a Fourier series in $\psi$:

${\hat{A}}(k,\psi)=\sum_{q=-\infty}^{+\infty}{a}(k;q)\mathrm{e}^{\mathrm{i}q\psi}\text{,}$
(25)

for $q\in\mathbb{Z}$. The Fourier-Bessel coefficients ${a}(k;q)$ of the image
ring ${\hat{A}}(k,\psi)$ are given by

${a}(k;q)=\frac{1}{2\pi}\int_{0}^{2\pi}{\hat{A}}(k,\psi)\mathrm{e}^{-\mathrm{i}q\psi}\mathop{d\psi}~{}.$
(26)

These coefficients can be represented in a more traditional fashion by
recalling that the Bessel function $J_{q}(kx)$ can be written as:

$J_{q}(kx)=\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{e}^{-\mathrm{i}kx\cos(\psi+\pi/2)-\mathrm{i}q\psi}\mathop{d\psi}~{},$
(27)

which, when combined with the definition of the Fourier transform, implies
that:

$\displaystyle{a}(k;q)$ $\displaystyle=$ $\displaystyle\iint
A(x,\theta)\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{e}^{-\mathrm{i}kx\cos(\psi-\theta)-\mathrm{i}q\psi}\mathop{d\psi}xdx\mathop{d\theta}~{}$
(28) $\displaystyle=$ $\displaystyle\iint
A(x,\theta)\mathrm{e}^{-\mathrm{i}q(\theta+\pi/2)}\frac{1}{2\pi}\int_{0}^{2\pi}\mathrm{e}^{-\mathrm{i}kx\cos(\psi+\pi/2)-\mathrm{i}q\psi}\mathop{d\psi}xdx\mathop{d\theta}~{}$
(29) $\displaystyle=$ $\displaystyle\iint
A(x,\theta)\mathrm{e}^{-\mathrm{i}q(\theta+\pi/2)}J_{q}(kx)xdx\mathop{d\theta}~{},$
(30)

which is the inner product between the original image (in real space) and a
Bessel function.

Note that, using the Fourier-Bessel representation of an image, the rotation
${\cal R}_{\gamma}$ can now be represented as:

$\displaystyle{\cal R}_{\gamma}{\hat{A}}(k,\psi)$ $\displaystyle=$
$\displaystyle{\hat{A}}(k,\psi-\gamma)$ (31) $\displaystyle=$
$\displaystyle\sum_{q=-\infty}^{+\infty}{a}(k;q)\mathrm{e}^{\mathrm{i}q(\psi-\gamma)}$
(32) $\displaystyle=$
$\displaystyle\sum_{q=-\infty}^{+\infty}{a}(k;q)\cdot\mathrm{e}^{-\mathrm{i}q\gamma}\mathrm{e}^{\mathrm{i}q\psi}\text{,}$
(33)

such that the Fourier-Bessel coefficients of the rotated image ring ${\cal
R}_{\gamma}\circ{\hat{A}}(k,\cdot)$ are given by the original Fourier-Bessel
coefficients ${a}(k;q)$, each multiplied by the phase-factor
$\mathrm{e}^{-\mathrm{i}q\gamma}$. Note that (33) easily allows for arbitrary
rotations by any $\gamma$, requiring only a pointwise multiplication by the
appropriate phase factor. This observation allows for the efficient alignment
of one image to another (see [16, 17]).

###  A.4 Image discretization

We denote by $\Omega(1)$ and $\Omega({K})$ the ball of radius $1$ and ${K}$,
respectively, in either real or frequency space:

$\Omega(1):=\\{{\boldsymbol{x}}\in\mathbb{R}^{2}\mbox{~{}such~{}that~{}}\|{\boldsymbol{x}}\|\leq
1\\}\text{,}\quad\Omega({K}):=\\{{\boldsymbol{k}}\in\mathbb{R}^{2}\mbox{~{}such~{}that~{}}\|{\boldsymbol{k}}\|\leq{K}\\}\text{,}$
(34)

and we will assume that all the images considered are supported in
$\Omega(1)$. The support constraint $A({\boldsymbol{x}})\subseteq\Omega(1)$
implies that the representation ${\hat{A}}({\boldsymbol{k}})$ will have a
bandlimit of $1$, so that ${\hat{A}}$ can be accurately reconstructed from its
values sampled on a frequency grid with spacing $\mathcal{O}(1)$ [19].

We also assume that any relevant signal within the images has a maximum
effective spatial-frequency magnitude of ${K}$; i.e., that the salient
features of ${\hat{A}}$ are concentrated in ${\boldsymbol{k}}\in\Omega_{{K}}$.
Consequently, we expect that the inversion

$A({\boldsymbol{x}})\approx\frac{1}{(2\pi)^{2}}\iint_{\Omega_{K}}{\hat{A}}({\boldsymbol{k}})\mathrm{e}^{\mathrm{i}{\boldsymbol{k}}\cdot{\boldsymbol{x}}}\mathop{d{\boldsymbol{k}}}$
(35)

accurately reconstructs $A$ to the resolution associated with ${K}$ (ignoring
high frequency ringing artifacts). When these assumptions hold we won’t lose
much accuracy when applying these transformations in a discrete setting (as
discussed below).

Note that the largest ${K}$ that we can reasonably expect to attain is related
to the number of pixels ‘$N$’ in the original image. More specifically, once
we assume that $A({\boldsymbol{x}})$ is supported in
${\boldsymbol{x}}\in\Omega(1)$, the pixel-spacing will be ${\Delta x}=2/N$,
implying the Nyquist spatial-frequency is given by $\pi/{\Delta x}=N\pi/2$.

Since we are using $N_{A}$ picked-particle images $A_{j}$ to produce a low-
resolution estimate of the imaged molecule (denoted $F$ below), we will
typically choose ${K}$ to be significantly smaller than $N\pi/2$, often in the
range of ${K}\sim 50$. Given the pixel spacing, this corresponds to a low-
resolution reconstruction with approximately $20$Åresolution. Additionally,
the fact that $x\leq 1$ implies that the Bessel coefficients ${a}(k,q)$ will
be concentrated in the range $|q|\lesssim k$, meaning that the Bessel
coefficients ${a}(k,q)$ across all $k\in[0,{K}]$ will be concentrated in
$q\in[-{Q}/2,+{Q}/2-1]$ for ${Q}=\mathcal{O}({K})=\mathcal{O}(N)$.

With the notation above, we can consider an $N\times N$ image as a discrete
set of pixel-averaged samples within $[-1,+1]^{2}$:

$A_{n_{1},n_{2}}=\frac{1}{{\Delta x}^{2}}\int_{n_{1}{\Delta
x}}^{(n_{1}+1){\Delta x}}\int_{n_{2}{\Delta x}}^{(n_{2}+1){\Delta
x}}A(x_{1},x_{2})\,dx_{2}dx_{1}$ (36)

for indices $n_{1},n_{2}\in\\{0,\ldots,N-1\\}$. We approximate the Fourier
transform ${\hat{A}}$ at any ${\boldsymbol{k}}\in\mathbb{R}^{2}$ via the
simple summation:

${\hat{A}}({\boldsymbol{k}})=({\Delta
x})^{2}\sum_{n_{2}=0}^{N-1}\sum_{n_{1}=0}^{N-1}A_{n_{1},n_{2}}e^{-\mathrm{i}{\boldsymbol{k}}\cdot{\boldsymbol{x}}_{n_{1},n_{2}}}~{}\text{,}$
(37)

where ${\boldsymbol{x}}_{n_{1},n_{2}}$ is the appropriately-chosen pixel-
center ${\Delta x}\left(n_{1}+\frac{1}{2},n_{2}+\frac{1}{2}\right)$. Because
we have assumed that the image is sufficiently well sampled (i.e., that
${\hat{A}}$ contains little relevant frequency-content above the Nyquist-
frequency ${K}$), we expect the simple sum above to be accurate.

We will typically evaluate $\hat{A}({\boldsymbol{k}})$ for ${\boldsymbol{k}}$
on a polar-grid, with $k$\- and $\psi$-values corresponding to suitable
quadrature nodes, using the non-uniform FFT (NUFFT) to compute
${\hat{A}}(k,\psi)$ at those nodes (see [19, 17]). The quadrature in $k$
corresponds to a set of ${R}$ nodes: $k_{1},\ldots,k_{{R}}$ and radial weights
$w_{1},\ldots,w_{{R}}$ so that

$\int_{0}^{{K}}g(k)kdk\approx\sum_{r=1}^{{R}}g(k_{r})w_{r}$ (38)

to high accuracy for a smooth function $g(k)$. The exact radial quadrature
nodes we select are usually those inherited from a volumetric grid (see
section A.6 below). For convenience, we will sometimes refer to the weight
function $w(k)$ as some continuous interpolant of the quadrature-weights, with
$w(k_{r})=w_{r}$. The ${Q}$ angular-nodes $\psi_{0},\ldots,\psi_{{Q}-1}$ are
equispaced in the periodic interval $[0,2\pi)$, with a spacing of
${\Delta\psi}=2\pi/{Q}$, and $\psi_{q^{\prime}}=q^{\prime}{\Delta\psi}$. These
equispaced $\psi$-nodes allow for spectrally accurate trapezoidal quadrature
in the $\psi$-direction, and we approximate the Fourier-Bessel coefficients of
each image-ring ${\hat{A}}(k,\psi)$ as follows:

${a}(k,q)\approx\sum_{q^{\prime}=0}^{{Q}-1}{\hat{A}}(k,\psi_{q^{\prime}})\exp\left(-\mathrm{i}q\psi_{q^{\prime}}\right){\Delta\psi}\text{,}$
(39)

with the index $q$ considered periodically in the interval
$[-{Q}/2,\ldots,+{Q}/2-1]$ (so that, e.g., the $q$-value of ${Q}-1$
corresponds to the $q$-value of $-1$).

###  A.5 Inner products

The inner-product between an image $A$ and a rotated and translated version of
image $B$ is denoted by:

$\displaystyle{\cal X}(\gamma,{\boldsymbol{\delta}};A,B)$ $\displaystyle:=$
$\displaystyle\langle{\cal R}(\gamma)\circ A,{\cal
T}({\boldsymbol{\delta}})\circ B\rangle\text{,}$ (40)

and (ignoring factors of $2\pi$) is approximated by:

$\displaystyle{\cal X}(\gamma,{\boldsymbol{\delta}};{\hat{A}},{\hat{B}})$
$\displaystyle\approx$ $\displaystyle\iint_{\Omega({K})}\left[{\cal
R}(\gamma)\circ{\hat{A}}({\boldsymbol{k}})\right]^{\dagger}\cdot\left[\widehat{\cal
T}({\boldsymbol{\delta}},{\boldsymbol{k}})\odot{\hat{B}}({\boldsymbol{k}})\right]d{\boldsymbol{k}}\text{.}$
(41)

The expression ${\cal X}$ can be thought of as a bandlimited inner product.

###  A.6 Volume notation

We use ${\boldsymbol{x}},{\boldsymbol{k}}\in\mathbb{R}^{3}$ to represent
spatial position and frequency, respectively. In spherical coordinates, the
vector ${\boldsymbol{k}}$ is represented as:

$\displaystyle{\boldsymbol{k}}$ $\displaystyle=$ $\displaystyle
k\cdot\hat{\boldsymbol{k}},\text{\ \ with\ \
}\hat{\boldsymbol{k}}=(\cos\phi\sin\theta,\sin\phi\sin\theta,\cos\theta)\text{,}$
(42)

with polar angle $\theta$ and azimuthal angle $\phi$ representing the unit
vector $\hat{\boldsymbol{k}}$ on the surface of the sphere $S^{2}$.

Using a right-handed basis, a rotation about the third axis by angle $\alpha$
is represented as:

${\cal R}_{z}(\alpha)=\left(\begin{array}[]{ccc}+\cos\alpha&-\sin\alpha&0\\\
+\sin\alpha&+\cos\alpha&0\\\ 0&0&1\end{array}\right)\text{,}$ (43)

and a rotation about the second axis by angle $\beta$ is represented as:

${\cal R}_{y}(\beta)=\left(\begin{array}[]{ccc}+\cos\beta&0&+\sin\beta\\\
0&1&0\\\ -\sin\beta&0&+\cos\beta\end{array}\right)\text{.}$ (44)

A rotation ${\cal R}(\tau)$ of a vector ${\boldsymbol{k}}\in\mathbb{R}^{3}$
can be represented by the vector of Euler angles $\tau=(\alpha,\beta,\gamma)$:

${\cal R}(\tau)\cdot{\boldsymbol{k}}={\cal R}_{z}(\alpha)\circ{\cal
R}_{y}(\beta)\circ{\cal R}_{z}(\gamma)\cdot{\boldsymbol{k}}\text{.}$ (45)

We represent a given volume as a function $F\in L^{2}(\mathbb{R}^{3})$, with
values corresponding to the intensity at each location
${\boldsymbol{x}}\in\Omega(1)$. We will refer to ${\hat{F}}({\boldsymbol{k}})$
in spherical coordinates as ${\hat{F}}(k,\hat{\boldsymbol{k}})$. With this
notation, each ${\hat{F}}(k,\cdot)$ corresponds to a ‘shell’ in frequency-
space with radius $k$. The rotation of any volume ${\cal
R}(\tau)\circ{\hat{F}}(k,\hat{\boldsymbol{k}})$ corresponds to the function
${\hat{F}}(k,{\cal R}(\tau^{-1})\cdot\hat{\boldsymbol{k}})$, where $\tau^{-1}$
corresponds to the vector of Euler angles
$\tau^{-1}=\left(-\gamma,-\beta,-\alpha\right)$.

###  A.7 Spherical harmonics

Using the notation above, we can represent a volume
${\hat{F}}(k,\hat{\boldsymbol{k}})$ as:

${\hat{F}}(k,\hat{\boldsymbol{k}})=\sum_{l=0}^{+\infty}\sum_{m=-l}^{m=+l}{f}(k;l,m)Y_{l}^{m}(\hat{\boldsymbol{k}})\text{,}$
(46)

where $Y_{l}^{m}(\hat{\boldsymbol{k}})$ represents the spherical harmonic of
degree $l$ and order $m$:

$Y_{l}^{m}(\theta,\phi)=Z_{l}^{m}\mathrm{e}^{+\mathrm{i}m\phi}P_{l}^{m}(\cos\theta)\text{,}$
(47)

with $P_{l}^{m}$ representing the usual (unnormalized) associated Legendre
function, and $Z_{l}^{m}$ the normalization constant:

$Z_{l}^{m}=\sqrt{\frac{2l+1}{4\pi}\times\frac{\left(l-|m|\right)!}{\left(l+|m|\right)!}}\text{.}$
(48)

The coefficients ${f}(k;l,m)$ define the spherical harmonic expansion of the
$k$-shell ${\hat{F}}(k,\cdot)$.

Given a rotation $\tau=(\alpha,\beta,\gamma)$, we can represent the rotated
volume ${\hat{G}}:={\cal R}(\tau){\hat{F}}$ as:

${g}_{l}^{m_{1}}=\sum_{m_{2}=-l}^{m_{2}=+l}\mathrm{e}^{-\mathrm{i}m_{1}\gamma}d_{m_{1},m_{2}}^{l}(\beta)\mathrm{e}^{-\mathrm{i}m_{2}\alpha}{f}_{l}^{m_{2}}\text{,}$
(49)

where $d_{m_{1},m_{2}}^{l}(\beta)$ represents the Wigner d-matrix of degree
$l$ associated with the second Euler angle $\beta$.

###  A.8 Volume discretization

As with our discretization of images, we assume that the volume
$F({\boldsymbol{x}})$ is supported in $\Omega(1)$, and that the relevant
frequency content is contained in $\Omega({K})$.

We then discretize the radial component of $\Omega({K})\in\mathbb{R}^{3}$
using a Gauss-Jacobi quadrature for $k$ built with a weight-function
corresponding to a radial weighting of $k^{2}dk$; the number of radial
quadrature nodes ${R}$ will be of the order $\mathcal{O}({K})$. On each of the
$k$-shells, we build a quadrature for $\hat{\boldsymbol{k}}$ that discretizes
the polar-angle $\beta$ using a legendre-quadrature on $[0,\pi]$, and the
azimuthal-angle $\alpha$ using a uniform periodic grid in $[0,2\pi)$, with
associated weights designed for integration on the surface of the sphere; the
total number of spherical-shell-quadrature-nodes ${T}(k_{r})$ for any
particular $k_{r}$ need only be $\mathcal{O}(k_{r}^{2})$ [16].

Each of the $k$-shells ${\hat{F}}(k_{r},\cdot)$ can be accurately described
using spherical harmonics of order $l=\mathcal{O}(k_{r})$. Thus, the number of
spherical-harmonic coefficients required for each shell
${\hat{F}}(k_{r},\cdot)$ is $\mathcal{O}(k_{r}^{2})$. The total number of
spherical-harmonic coefficients required to approximate
${\hat{F}}(k,\hat{\boldsymbol{k}})$ over $\Omega({K})$ is
$\mathcal{O}({K}^{3})$, with a maximum degree of ${L}=\mathcal{O}({K})$. The
associated maximum order will be ${M}=1+2{L}$, which is also
$\mathcal{O}({K})$.

###  A.9 Template generation

Given any volume ${\hat{F}}({\boldsymbol{k}})$, we can generate templates
efficiently using the techniques described in [16], which make use of the
spherical harmonic representation ${f}(k;l,m)$. For our purposes, we generate
the templates described in Eq. 7 for a collection of viewing directions
$(\alpha_{t},\beta_{t})$, for $t\in\\{1,\ldots,{T}\\}$. The particular choice
of the third Euler angle ($\gamma_{t}$) for each of the templates is dictated
by the strategy for template generation in [16], and does not play an
important role, since it corresponds to an in-plane rotation of the relevant
slice. We discretize the set of viewing directions to correspond with the
$\hat{\boldsymbol{k}}$ discretization on the largest $k$-shell in our
volumetric (spherical) quadrature grid: that is, ${T}:={T}(k_{{R}})$. These
viewing directions will be used to enumerate the templates (i.e.,
$N_{S}={T}$). Each of the templates are stored on a polar or Fourier-Bessel
grid identical to those used to store the images.

###  A.10 Image noise model

We assume that each of the picked-particle images ${\hat{A}}_{j}$ is a noisy
version of a ‘signal’ corresponding to a 2-dimensional projection of some
(unknown) molecular electron-density-function $F^{{\text{\tiny true}}}$:

${\hat{A}}^{{\text{\tiny signal}}}({\boldsymbol{k}})={\cal
T}(+{\boldsymbol{\delta}}_{j}^{{\text{\tiny
true}}},{\boldsymbol{k}})\odot{\hat{S}}(\tau_{j}^{{\text{\tiny
true}}};CTF;{\hat{F}})({\boldsymbol{k}})\text{,}$ (50)

where $\tau_{j}^{{\text{\tiny true}}}$ and
${\boldsymbol{\delta}}_{j}^{{\text{\tiny true}}}$ correspond to the true (but
unknown) viewing-angle and displacement associated with image ${\hat{A}}_{j}$.
In this context we say the image’s ‘true viewing direction’ is defined by the
first two Euler angles in $\tau_{j}^{{\text{\tiny true}}}$, namely the
azimuthal and polar angles $\alpha_{j}^{{\text{\tiny true}}}$ and
$\beta_{j}^{{\text{\tiny true}}}$.

We will model the noise in real space as independent and identically
distributed (iid), with a variance of $\sigma^{2}$ on a unit scale in two-
dimensional real space. This is a simple model of detector noise [20, 21, 22,
23], and does not take into account more complicated features, such as
structural noise associated with image preparation (see [24]). We assume that
the image values are the sum of the signal plus noise:

$A_{n_{1},n_{2}}=A({\boldsymbol{x}}_{n_{1},n_{2}})=A^{{\text{\tiny
signal}}}_{n_{1},n_{2}}+A^{{\text{\tiny noise}}}_{n_{1},n_{2}}\text{,}$ (51)

where the signal and noise are represented by the arrays $A^{{\text{\tiny
signal}}}_{n_{1},n_{2}}$ and $A^{{\text{\tiny noise}}}_{n_{1},n_{2}}$,
respectively. Because each entry in $A^{{\text{\tiny noise}}}_{n_{1},n_{2}}$
corresponds to an average over the area-element ${\Delta x}^{2}$ associated
with a single pixel, we expect the variance of each $A^{{\text{\tiny
noise}}}_{n_{1},n_{2}}$ to scale inversely with ${\Delta x}^{2}$. That is to
say, we’ll assume that each element of the noise array $A^{{\text{\tiny
noise}}}_{n_{1},n_{2}}$ is drawn from the standard normal distribution ${\cal
N}\left(0,\frac{\sigma^{2}}{{\Delta x}^{2}}\right)$, with zero mean and a
variance of $\sigma^{2}/{\Delta x}^{2}$. With these assumptions, downsampling
the image by averaging neighboring pixels will correspond to an increase in
the effective pixel-size and a simultaneous reduction in the variance of the
noise associated with each (now larger) pixel.

Given the assumptions above, ${\hat{A}}$ will be modeled by:

${\hat{A}}({\boldsymbol{k}})={\hat{A}}^{{\text{\tiny
signal}}}({\boldsymbol{k}})+{\hat{A}}^{{\text{\tiny
noise}}}({\boldsymbol{k}})\text{,}$ (52)

where ${\hat{A}}^{{\text{\tiny noise}}}$ is now complex and iid. Because the
noise-term $A^{{\text{\tiny noise}}}$ is real, the noise-term
${\hat{A}}^{{\text{\tiny noise}}}$ will be complex, with the conjugacy
constraint that ${\hat{A}}^{{\text{\tiny
noise}}}(+{\boldsymbol{k}})={\hat{A}}^{{\text{\tiny
noise}}}(-{\boldsymbol{k}})^{\dagger}$.

###  A.11 Variance Scaling

We define ${\hat{\sigma}}^{2}=\frac{\pi^{2}\sigma^{2}}{{\Delta x}^{2}}$ as the
variance on a unit scale in two-dimensional frequency space. We expect that
the noise term ${\hat{A}}^{{\text{\tiny noise}}}({\boldsymbol{k}})$ integrated
over any area element ${\Delta k}^{2}$ will have a variance of
$\frac{{\hat{\sigma}}^{2}}{{\Delta k}^{2}}$. Recalling our polar quadrature
described above, we typically sample ${\boldsymbol{k}}$ along ${R}$ radial
quadrature nodes $k_{1},\ldots,k_{{R}}$ and ${Q}$ angular quadrature nodes
$\psi_{0},\ldots,\psi_{{Q}-1}$, with radial and angular weights $w_{r}$ and
${\Delta\psi}$, respectively. In this quadrature scheme, each quadrature node
${\boldsymbol{k}}_{rq}=(k_{r},\psi_{q})$ is associated with the area element
$w_{r}{\Delta\psi}$, which approximates the ‘$kdk$’ integration weight.
Consequently, we expect that the noise term ${\hat{A}}^{{\text{\tiny
noise}}}({\boldsymbol{k}}_{rq})$ evaluated at $(k_{r},\psi_{q})$ on our polar
quadrature grid will have a variance of
$\frac{{\hat{\sigma}}^{2}}{w_{r}{\Delta\psi}}$.

Given the radially-dependent variance associated with each degree of freedom
in the images, we will typically consider rescaled images
${\hat{A}}(k_{r},\psi_{q})\sqrt{w_{r}{\Delta\psi}}$, rather than the raw image
values ${\hat{A}}(k_{r},\psi_{q})$. Fro simplicity, we will often write this
as ${\hat{A}}({\boldsymbol{k}})\sqrt{w(k)}$, where we assume that
${\Delta\psi}$ is the same for every $k$-value, and $w(k)$ is some functional
interpolant of the quadrature weights with $w(k_{r})=w_{r}$. This rescaling
allows us to equate the standard vector 2-norm with the $L^{2}$ integral (up
to quadrature error):

$\sum_{r=1}^{r={R}}\sum_{q=1}^{q={Q}}\|{\hat{A}}(k_{r},\psi_{q})\cdot\sqrt{w_{r}{\Delta\psi}}\|^{2}\approx\int_{k=0}^{k={K}}\int_{\psi=0}^{\psi=2\pi}\|{\hat{A}}(k,\psi)\|^{2}kdkd\psi\text{.}$

This rescaling will allow our definition of the radial objective function
below to correspond to a standard log-likelihood (see (55) in section A.14).
In a similar fashion, the least-squares formulation we will use for volume
reconstruction corresponds to a maximum-likelihood estimate, with similar
properties holding for the inner products used for alignment (see sections
A.18 and A.22 below).

###  A.12 Radial Principal Modes

In these sections, we briefly review the notion of radial principal modes for
image compression, described in [25]. The premise is simple: given an image
${\hat{A}}(k,\psi)$, not all the frequency rings ${\hat{A}}(k,\cdot)$ are
equally useful. By selecting the appropriate linear combinations of frequency
rings, we can compress the images while retaining information useful for
alignment. We choose the coefficients of this linear combination by solving an
eigenvalue problem, which is why we refer to these linear combinations as
principal modes.

This strategy has two benefits. The first is to make the overall computation
more efficient; rather than retaining all ${R}\approx 50$ radial quadrature
nodes, we can often compress the majority of relevant information into a
smaller number of principal modes (typically 10-16). The second benefit is
that, by retaining only those principal modes which are useful for alignment,
we explicitly avoid those degrees of freedom within the original images that
are dominated by noise and not useful for alignment. The same set of modes can
be used for the experimental images and the estimated volumes. We make use of
this compression for both the alignment and reconstruction steps of our EMPM
iteration.

###  A.13 Principal image rings and principal volume shells

Let us assume that we are given a unit vector
${\boldsymbol{u}}={[u_{1},u_{2},\ldots,u_{{R}}]}^{\intercal}\in\mathbb{R}^{{R}}$.
We define the linear combination
$[{{\boldsymbol{u}}}^{\intercal}{\hat{A}}](\psi_{q})$ as:

$\displaystyle[{{\boldsymbol{u}}}^{\intercal}{\hat{A}}](\psi_{q})=\sum_{r=1}^{{R}}u_{r}{\hat{A}}(k_{r},\psi_{q})\sqrt{w_{r}{\Delta\psi}}\text{,}$
(53)

where the rescaling factor $\sqrt{w_{r}{\Delta\psi}}$ is designed to equalize
the variances associated with the different $k$-values (see section A.11).

In a moment, we will choose our ${\boldsymbol{u}}$ to be one of the
eigenvectors of an ${R}\times{R}$ matrix, and we will refer to the linear
combination $[{{\boldsymbol{u}}}^{\intercal}{\hat{A}}](\psi_{q})$ as the
‘principal image ring’ associated with ${\boldsymbol{u}}$. Note that, based on
our noise model, we assume that the image noise is iid across image-rings.
Consequently, the variance of the noise
$[{{\boldsymbol{u}}}^{\intercal}{\hat{A}}^{{\text{\tiny noise}}}](\psi_{q})$
is equal to $\sum_{r}{\boldsymbol{u}}_{r}^{2}{\hat{\sigma}}^{2}$, which is
equal to $\|{\boldsymbol{u}}\|^{2}{\hat{\sigma}}^{2}$. Because
${\boldsymbol{u}}$ is a unit vector, this last expression is simply
${\hat{\sigma}}^{2}$, which is the same as the variance of any individual term
${\hat{A}}^{{\text{\tiny noise}}}(k_{r},\psi_{q})$ for any particular $k_{r}$.
Moreover, if we consider two orthonormal vectors ${\boldsymbol{u}}_{1}$ and
${\boldsymbol{u}}_{2}$, then
$[{{\boldsymbol{u}}_{1}}^{\intercal}{\hat{A}}^{{\text{\tiny
noise}}}](\psi_{q})$ and
$[{{\boldsymbol{u}}_{2}}^{\intercal}{\hat{A}}^{{\text{\tiny
noise}}}](\psi_{q})$ will be independent random variables, each drawn from
${\cal N}(0,{\hat{\sigma}}^{2})$.

The same strategy can be used to define a principal volume shell
$[{{\boldsymbol{u}}}^{\intercal}{\hat{F}}](\hat{\boldsymbol{k}})$ as:

$\displaystyle[{{\boldsymbol{u}}}^{\intercal}{\hat{F}}](\hat{\boldsymbol{k}})=\sum_{r=1}^{{R}}u_{r}{\hat{F}}(k_{r},\hat{\boldsymbol{k}})\sqrt{w_{r}{\Delta\psi}}\text{,}$
(54)

using the same weights as those used for the principal image rings. In most
situations, we will assume that the angular discretization ${\Delta\psi}$ is
uniform across all the image rings, allowing us to ignore the
$\sqrt{{\Delta\psi}}$ factor in the rescaling-factor above.

###  A.14 Choice of principal modes

To select the actual principal modes used above, we construct the following
objective function for the vector ${\boldsymbol{u}}$:

$\displaystyle
C({\boldsymbol{u}};F,\mu_{\tau},\mu_{{\boldsymbol{\delta}}},\mu_{CTF})$
$\displaystyle=$
$\displaystyle\int\cdots\int\left\|{{\boldsymbol{u}}}^{\intercal}{\hat{S}}(\tau,{\boldsymbol{\delta}};CTF,F)-{{\boldsymbol{u}}}^{\intercal}{\hat{S}}(\tau^{\prime},{\boldsymbol{\delta}}^{\prime};CTF^{\prime},F)\right\|^{2}\cdot$
(55) $\displaystyle d\psi
d\mu_{\tau}\left(\tau\right)d\mu_{\tau}\left(\tau^{\prime}\right)d\mu_{{\boldsymbol{\delta}}}\left({\boldsymbol{\delta}}\right)d\mu_{{\boldsymbol{\delta}}}\left({\boldsymbol{\delta}}^{\prime}\right)d\mu_{CTF}(CTF)d\mu_{CTF}(CTF^{\prime})\text{,}$

where ${\hat{S}}(\tau,{\boldsymbol{\delta}};CTF;F)$ is the usual noiseless
image template associated with viewing-angle $\tau$, displacement
${\boldsymbol{\delta}}$ and the contrast-transfer-function
$CTF({\boldsymbol{k}})$. The integral in (55) is taken over the distributions
of viewing angle, displacement and CTF function (i.e., $\mu_{\tau}$,
$\mu_{{\boldsymbol{\delta}}}$ and $\mu_{CTF}$, respectively). Given a volume
$F$, as well as distributions $\mu_{\tau}$, $\mu_{{\boldsymbol{\delta}}}$ and
$\mu_{CTF}$, the objective function $C({\boldsymbol{u}})$ is (up to an affine
transformation) equal to the negative log-likelihood of mistaking one
noiseless principal image ring (constructed using principal-vector
${\boldsymbol{u}}$) for a randomly rotated version of another (also
constructed using ${\boldsymbol{u}}$). Those vectors ${\boldsymbol{u}}$ for
which $C({\boldsymbol{u}})$ is high correspond to linear combinations of
frequencies which are useful for alignment (i.e. for which the typical
principal image rings look very different from one another). Conversely, those
vectors ${\boldsymbol{u}}$ for which $C({\boldsymbol{u}})$ is low correspond
to linear combinations of frequencies which are not useful for alignment
(i.e., for which the typical principal image rings look quite similar to one
another).

The objective function $C({\boldsymbol{u}})$ is quadratic in
${\boldsymbol{u}}$, which means that it can be optimized by finding the
dominant eigenvector of the Hessian
$\boldsymbol{C}:=\partial_{{\boldsymbol{u}}{\boldsymbol{u}}}C$, defined as the
${R}\times{R}$ matrix:

$C({\boldsymbol{u}})={{\boldsymbol{u}}}^{\intercal}\cdot\boldsymbol{C}\cdot{\boldsymbol{u}}\text{.}$

Moreover, the corresponding eigenvalue will be equal to the objective function
$C$ for that eigenvector. To actually compute the Hessian $\boldsymbol{C}$ in
practice, we use one of two strategies.

###  A.15 Empirical principal modes

At first, when we have no good estimate of either $F$ or the distributions
$\mu_{\tau}$, $\mu_{{\boldsymbol{\delta}}}$ or $\mu_{CTF}$, we simply use the
images $A_{j}$ themselves to form a Monte Carlo estimate of the integral in
(55). This boils down to the following simple approximation:

$\displaystyle C({\boldsymbol{u}})^{{\text{\tiny empir}}}$ $\displaystyle=$
$\displaystyle\frac{1}{N_{A}^{2}}\sum_{j,j^{\prime}}\iiint\left\|{{\boldsymbol{u}}}^{\intercal}{\cal
R}(\gamma){\hat{A}}_{j}-{{\boldsymbol{u}}}^{\intercal}{\cal
R}(\gamma^{\prime}){\hat{A}}_{j^{\prime}}\right\|^{2}d\psi d\gamma
d\gamma^{\prime}\text{,}$ (56)

where we have replaced the noiseless images
${\hat{S}}(\tau,{\boldsymbol{\delta}};CTF;F)$ and
${\hat{S}}(\tau,{\boldsymbol{\delta}};CTF;F)$ in the integrand of Eq. 55 with
randomly selected images ${\hat{A}}_{j}$ and ${\hat{A}}_{j^{\prime}}$ from the
original image pool, and used the empirical distribution of the images
themselves to stand in for the unkown $\mu_{\tau}$,
$\mu_{{\boldsymbol{\delta}}}$ and $\mu_{CTF}$ (after averaging over the in-
plane rotations associated with $\gamma$ and $\gamma^{\prime}$). The
associated kernel $\boldsymbol{C}^{{\text{\tiny empir}}}$ then reduces to:

$\boldsymbol{C}(r,r^{\prime})^{{\text{\tiny
empir}}}\propto\sqrt{w_{r}w_{r^{\prime}}}\cdot\left(\frac{1}{N_{A}}\sum_{j}\sum_{q}{a}_{j}(k_{r},q)^{\dagger}{a}_{j}(k_{r^{\prime}},q)-\frac{1}{N_{A}^{2}}\left[\sum_{j}{a}(k_{r},0)\right]^{\dagger}\left[\sum_{j}{a}(k_{r^{\prime}},0)\right]\right)\text{.}$
(57)

###  A.16 Estimated principal modes

Once we have a reasonable estimate of the imaged volume $F^{{\text{\tiny
estim}}}$, we can use it (along with the current estimated distributions
$\mu_{\tau}^{{\text{\tiny estim}}}$ and
$\mu_{{\boldsymbol{\delta}}}^{{\text{\tiny estim}}}$) to calculate
$C^{{\text{\tiny estim}}}$ via (55). In practice, we typically assume that
$\mu_{\tau}^{{\text{\tiny estim}}}$ is uniform and
$\mu_{{\boldsymbol{\delta}}}^{{\text{\tiny estim}}}$ is an isotropic Gaussian
centered at the origin with a standard deviation estimated using the current
${\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}$. With these assumptions the
associated kernel $\boldsymbol{C}^{{\text{\tiny estim}}}$ can be computed
analytically from the spherical harmonic coefficients ${f}(k;l,m)$:

$\displaystyle\boldsymbol{C}_{r,r^{\prime}}^{{\text{\tiny estim}}}$
$\displaystyle\propto$
$\displaystyle+\sqrt{W_{r}W_{r^{\prime}}}\cdot\left[\sum_{l=0}^{+\infty}\sum_{m=-l}^{m=+l}{f}(k_{r};l,m)^{\dagger}{f}(k_{r^{\prime}};l,m)\right]\cdot\boldsymbol{E}^{+}(k_{r},k_{r^{\prime}})$
(59)
$\displaystyle-\sqrt{W_{r}W_{r^{\prime}}}\cdot\left[{f}(k_{r};0,0)^{\dagger}{f}(k_{r^{\prime}};0,0)\right]\cdot\boldsymbol{E}^{-}(k_{r})\cdot\boldsymbol{E}^{-}(k_{r^{\prime}})\text{,}$

as ${Q}\rightarrow\infty$, where the CTF-modulated weight $W_{r}$ takes the
form:

$\displaystyle\sqrt{W_{r}}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{w_{r}}}{N_{A}}\sum_{j}CTF_{j}(k_{r})\text{,}$ (60)

and the terms $\boldsymbol{E}^{+}$ and $\boldsymbol{E}^{-}$ denote

$\displaystyle\boldsymbol{E}^{+}(k_{r},k_{r^{\prime}})=\exp\left(-\frac{\sigma_{{\boldsymbol{\delta}}}^{2}}{2}\left[k_{r}^{2}+k_{r^{\prime}}^{2}\right]\right)\cdot\exp\left(k_{r}k_{r^{\prime}}\sigma_{{\boldsymbol{\delta}}}^{2}\right)\text{,}$
(61)

and

$\displaystyle\boldsymbol{E}^{-}(k_{r})=\tilde{k}_{r}\sqrt{\frac{\pi}{2}}\exp\left(-\tilde{k}_{r}^{2}\right)\left({\cal
I}_{-1/2}(\tilde{k}_{r}^{2})-{\cal
I}_{+1/2}(\tilde{k}_{r}^{2})\right)\text{,}$ (62)

where $\tilde{k}_{r}=k_{r}\sigma_{{\boldsymbol{\delta}}}/2$, and ${\cal
I}_{q}$ refers to the modified Bessel function of first kind of order $q$.

###  A.17 Multiple principal modes

We select the principal modes
${\boldsymbol{u}}_{1},\ldots,{\boldsymbol{u}}_{H}$ to be the dominant $H$
eigenvectors of the kernel $\boldsymbol{C}$. We typically choose $H$ such that
the $(H+1)$-th eigenvalue of $\boldsymbol{C}$ is less than the global
tolerance times the first (dominant) eigenvalue of $\boldsymbol{C}$. When we
set the global tolerance to be $\epsilon_{{\text{\tiny tol}}}\sim 1e-2$ in our
numerical experiments below, this selection criterion typically corresponds to
selecting $12-16$ principal modes from $\boldsymbol{C}^{{\text{\tiny
empir}}}$, and slightly fewer ($10-14$) from $\boldsymbol{C}^{{\text{\tiny
estim}}}$.

We will refer to the matrix $U$ as the collection of $H$ principal modes:

$U=\left[\begin{array}[]{cccc}u_{1}&u_{2}&\cdots&u_{H}\end{array}\right]\text{,}$
(63)

and refer to the collection of respective principal images as
${U}^{\intercal}{\hat{A}}$:

${U}^{\intercal}{\hat{A}}=\left[\begin{array}[]{cccc}{u_{1}}^{\intercal}{\hat{A}}&{u_{2}}^{\intercal}{\hat{A}}&\cdots&{u_{H}}^{\intercal}{\hat{A}}\end{array}\right]\text{,}$
(64) $\text{such that \
}\left[{U}^{\intercal}{\hat{A}}\right](h,\psi)=\left[{u_{h}}^{\intercal}{\hat{A}}\right](\psi)\text{.}$

We use the same convention for principal volumes:

${U}^{\intercal}{\hat{F}}=\left[\begin{array}[]{cccc}{u_{1}}^{\intercal}{\hat{F}}&{u_{2}}^{\intercal}{\hat{F}}&\cdots&{u_{H}}^{\intercal}{\hat{F}}\end{array}\right]\text{,}$
(65) $\text{such that \
}\left[{U}^{\intercal}{\hat{F}}\right](h,\hat{\boldsymbol{k}})=\left[{u_{h}}^{\intercal}{\hat{F}}\right](\hat{\boldsymbol{k}})\text{.}$

We will also refer to the Bessel and spherical harmonic representations of the
principal images and principal volumes as ${U}^{\intercal}{a}$ and
${U}^{\intercal}{f}$, respectively.

###  A.18 Volume reconstruction

Most approaches to estimating molecular volumes in cryo-em involve some kind
of reconstruction process which is performed multiple times (over the course
of iterative refinement). For standard AM, the inputs to the reconstruction
process are the collection of $N_{A}$ picked particle images
${\hat{A}}_{j}({\boldsymbol{k}})$ and associated CTF-functions
$CTF_{j}({\boldsymbol{k}})$, along with estimates of the viewing angles
$\tau_{j}^{{\text{\tiny estim}}}$ and displacements
${\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}$. The output is the estimate
${\hat{F}}^{{\text{\tiny estim}}}({\boldsymbol{k}})$.

In our case, we will make use of both the standard least-squares
reconstruction, as well as a more computationally efficient approximate
Fourier inversion which we refer to as ‘quadrature-back-propagation’ (section
A.20 below).

###  A.19 Least-squares reconstruction

The standard strategy for reconstruction typically involves solving a linear
system in a least-squares sense. In real-space this least-squares system looks
like:

$\left[\begin{array}[]{c}PSD_{1}({\boldsymbol{x}})\star{\cal
T}(+{\boldsymbol{\delta}}_{1})\circ\left\\{S\circ{\cal
R}(\tau_{1})\circ\right\\}\\\ \vdots\\\
PSD_{N_{A}}({\boldsymbol{x}})\star{\cal
T}(+{\boldsymbol{\delta}}_{N_{A}})\circ\left\\{S\circ{\cal
R}(\tau_{N_{A}})\circ\right\\}\\\ \end{array}\right]\cdot
F({\boldsymbol{x}})\approx\left[\begin{array}[]{c}A_{1}({\boldsymbol{x}})\\\
\vdots\\\ A_{N_{A}}({\boldsymbol{x}})\\\ \end{array}\right]$

whereas in the Fourier domain it is written as:

$\left[\begin{array}[]{c}CTF_{1}({\boldsymbol{k}})\odot{\cal
T}(+{\boldsymbol{\delta}}_{1},{\boldsymbol{k}})\odot\left\\{{\hat{Sl}}\circ{\cal
R}(\tau_{1})\circ\right\\}\\\ \vdots\\\
CTF_{N_{A}}({\boldsymbol{k}})\odot{\cal
T}(+{\boldsymbol{\delta}}_{N_{A}},{\boldsymbol{k}})\odot\left\\{{\hat{Sl}}\circ{\cal
R}(\tau_{N_{A}})\circ\right\\}\\\
\end{array}\right]\cdot{\hat{F}}({\boldsymbol{k}})\approx\left[\begin{array}[]{c}{\hat{A}}_{1}({\boldsymbol{k}})\\\
\vdots\\\ {\hat{A}}_{N_{A}}({\boldsymbol{k}})\\\ \end{array}\right]$

or, more conveniently,

$\left[\begin{array}[]{c}CTF_{1}({\boldsymbol{k}})\odot\left\\{{\hat{S}}\circ{\cal
R}(\tau_{1})\circ\right\\}\\\ \vdots\\\
CTF_{N_{A}}({\boldsymbol{k}})\odot\left\\{{\hat{S}}\circ{\cal
R}(\tau_{N_{A}})\circ\right\\}\\\
\end{array}\right]\cdot{\hat{F}}({\boldsymbol{k}})\approx\left[\begin{array}[]{c}{\cal
T}(-{\boldsymbol{\delta}}_{1},{\boldsymbol{k}})\odot{\hat{A}}_{1}({\boldsymbol{k}})\\\
\vdots\\\ {\cal
T}(-{\boldsymbol{\delta}}_{N_{A}},{\boldsymbol{k}})\odot{\hat{A}}_{N_{A}}({\boldsymbol{k}})\\\
\end{array}\right]\text{.}$ (66)

In this expression, the $N_{A}$ viewing angles $\tau_{j}$ and displacements
${\boldsymbol{\delta}}_{j}$ are assumed to be given as input, as are the
image- or micrograph-specific CTF functions $CTF_{j}({\boldsymbol{k}})$ and
images ${\hat{A}}_{j}({\boldsymbol{k}})$. The unknown (to be determined) is
the volume ${\hat{F}}({\boldsymbol{k}})$.

This least-squares problem is standard in maximum-likelihood estimation, and
versions of this problem are solved in many software packages, often via
conjugate gradient iteration applied to the normal equations.

As written, the least-squares problem above has a block-structure: different
$k$ values are independent from one another. Thus, the problem can be solved
by determining the $k$-shell ${\hat{F}}(k,\hat{\boldsymbol{k}})$ from the
$k$-rings ${\hat{A}}_{j}(k,\psi)$ for one $k$ value at a time. However, the
problem stated in (66) is slightly misleading, as the residuals for each of
the $k$ values are not directly comparable due to the fact that different
$k$-shells have different variances associated with them (see section A.11).
To account for this, we rewrite (66) as:

$\left[\begin{array}[]{c}\sqrt{w(k)}\cdot
CTF_{1}({\boldsymbol{k}})\odot\left\\{{\hat{Sl}}\circ{\cal
R}(\tau_{1})\circ\right\\}\\\ \vdots\\\ \sqrt{w(k)}\cdot
CTF_{N_{A}}({\boldsymbol{k}})\odot\left\\{{\hat{Sl}}\circ{\cal
R}(\tau_{N_{A}})\circ\right\\}\\\
\end{array}\right]\cdot{\hat{F}}({\boldsymbol{k}})\approx\left[\begin{array}[]{c}\sqrt{w(k)}\cdot{\cal
T}(-{\boldsymbol{\delta}}_{1},{\boldsymbol{k}})\odot{\hat{A}}_{1}({\boldsymbol{k}})\\\
\vdots\\\ \sqrt{w(k)}\cdot{\cal
T}(-{\boldsymbol{\delta}}_{N_{A}},{\boldsymbol{k}})\odot{\hat{A}}_{N_{A}}({\boldsymbol{k}})\\\
\end{array}\right]\text{,}$ (67)

where the function $w(k_{r})$ is the radial quadrature weight. With this
formulation, the residual (calculated using the standard vector 2-norm) is (up
to a quadrature error) equal to an integral representing the log-likelihood
that the volume ${\hat{F}}({\boldsymbol{k}})$ could have produced the observed
images ${\hat{A}}_{j}({\boldsymbol{k}})$ (under the assumptions about the
noise-model described in section A.10).

###  A.20 Principal-mode least-squares reconstruction

The same general methodology above can be applied to the principal images
directly. To describe the modified least-squares system, we first define the
low-rank decomposition of the radial component of the CTF-functions. Let
$CTF(k,j)$ correspond to the radial component of the CTF-function
$CTF_{j}({\boldsymbol{k}})$. Using the discretized analog of a functional SVD
we represent $CTF(k,j)$ by

$\displaystyle CTF(k,j)$ $\displaystyle=$
$\displaystyle\sum_{h^{\prime}}u^{\prime}_{h^{\prime}}(k)\cdot\sigma^{\prime}_{h^{\prime}}v^{\prime}_{h^{\prime}}(j)\text{,}$
(68)

which we can truncate at rank $H^{\prime}$ in accordance with the global
tolerance. In the numerical experiments we present below, the image pools
typically span $30$ or so micrographs, usually resulting in $H^{\prime}$ about
$2-3$.

Using this $H^{\prime}$-rank decomposition of the CTF-array (across the image
pool), we can approximately represent the noiseless principal images as
follows:

$\displaystyle CTF_{j}({\boldsymbol{k}})\odot\left\\{{\hat{Sl}}\circ{\cal
R}(\tau_{j})\circ{\hat{F}}\right\\}({\boldsymbol{k}})$ $\displaystyle\approx$
$\displaystyle\sum_{h^{\prime}=1}^{H^{\prime}}v^{\prime}_{h^{\prime}}(j)\sigma^{\prime}_{h^{\prime}}\cdot
u^{\prime}_{h^{\prime}}(\|{\boldsymbol{k}}\|)\odot\left\\{{\hat{Sl}}\circ{\cal
R}(\tau_{j})\circ{\hat{F}}\right\\}({\boldsymbol{k}})\text{.}$ (69)

After using the radial principal modes to compress both this approximation and
the associated component of the right hand side of (67), we obtain:

$\displaystyle\sum_{h^{\prime}=1}^{H^{\prime}}v^{\prime}_{h^{\prime}}(j)\sigma^{\prime}_{h^{\prime}}\cdot{u_{h}}^{\intercal}\left[u^{\prime}_{h^{\prime}}({\boldsymbol{k}})\odot\left\\{{\hat{Sl}}\circ{\cal
R}(\tau_{j})\circ{\hat{F}}\right\\}({\boldsymbol{k}})\right]={u_{h}}^{\intercal}\left[{\cal
T}(-{\boldsymbol{\delta}}_{j},{\boldsymbol{k}})\odot{\hat{A}}({\boldsymbol{k}})\right]\text{.}$
(70)

By defining the CTF-modulated principal volume:

${\hat{G}}(k,\hat{\boldsymbol{k}};h^{\prime},{\hat{F}})={U}^{\intercal}\cdot\left[u^{\prime}_{h^{\prime}}(k)\odot{\hat{F}}(k,\hat{\boldsymbol{k}})\right]\text{,}$

we can now write the compressed least-squared problem as:

$\sum_{h^{\prime}=1}^{H^{\prime}}\sigma^{\prime}_{h^{\prime}}\left[\begin{array}[]{c}v^{\prime}_{h^{\prime}}(1)\left\\{{\hat{S}}\circ{\cal
R}(\tau_{1})\circ\right\\}\\\ \vdots\\\
v^{\prime}_{h^{\prime}}(N_{A})\left\\{{\hat{S}}\circ{\cal
R}(\tau_{N_{A}})\circ\right\\}\\\
\end{array}\right]\cdot{\hat{G}}({\boldsymbol{k}};h^{\prime},{\hat{F}})\approx\left[\begin{array}[]{c}{U}^{\intercal}\left[{\cal
T}(-{\boldsymbol{\delta}}_{1},{\boldsymbol{k}})\odot{\hat{A}}_{1}({\boldsymbol{k}})\right]\\\
\vdots\\\ {U}^{\intercal}\left[{\cal
T}(-{\boldsymbol{\delta}}_{N_{A}},{\boldsymbol{k}})\odot{\hat{A}}_{N_{A}}({\boldsymbol{k}})\right]\\\
\end{array}\right]\text{,}$ (71)

where the radial reweighting of $\sqrt{w(k)}$ in (67) is accounted for by our
definition of ${U}^{\intercal}$. Note that, when $H^{\prime}=H={R}$ the
equation above is equivalent to (67). However, when $H^{\prime}<{R}$ and
$H<{R}$ this will correspond to a projected version of (67).

Note also that, with this representation, we need not actually solve for
${\hat{F}}({\boldsymbol{k}})$. Instead, we can use the principal images (on
the right hand side) to solve the above problem (in a least-squares sense) for
the collection of CTF-modulated principal volumes
${\hat{G}}({\boldsymbol{k}};h^{\prime},{\hat{F}})$. Solving (71) can be less
computationally costly than solving (67) when the total number of degrees of
freedom within the various ${\hat{G}}$ is less than the total number of
degrees of freedom within ${\hat{F}}$ (i.e., when $H\times H^{\prime}<{R}$).
Even when this is not the case, the CTF-modulated principal volumes
${\hat{G}}$ (the solution to (71)) can still be more accurate (in terms of
comparison to the ground truth) than the solution ${\hat{F}}$ to (67). We see
this improvement most acutely when the principal mode reduction
${U}^{\intercal}$ effectively removes noise from the images ${\hat{A}}_{j}$
(e.g., if there are many frequency rings which contain little to no
information).

We remark that, when dealing with least-squares reconstruction and principal
mode reduction, ‘the diagram commutes’. That is, because the principal mode
projection ${U}^{\intercal}$ is orthonormal, the principal mode projection
${U}^{\intercal}\left\\{CTF\cdot{\hat{F}}\right\\}$ (i.e., the principal
projection of the full volume which solves the original least-squares problem
in (67)) will be the same as the solution
$\sum_{h^{\prime}}v^{\prime}\sigma^{\prime}{\hat{G}}$ to the principally-
projected least-squares problem in (71), provided of course that the
$H^{\prime}$-rank decomposition of the CTF-array is accurate (e.g., if
$H^{\prime}$ is taken to ${R}$). Thus, if $H\times H^{\prime}$ is indeed
larger than ${R}$, then instead of solving (71), we can obtain ${\hat{G}}$
more easily by just solving the original least-squares problem for ${\hat{F}}$
and then projecting down to the CTF-modulated principal volumes ${\hat{G}}$.

Note also that, due to linearity, the noiseless images associated with the
CTF-modulated principal volumes ${\hat{G}}$ will approximately equal the
noiseless principal images associated with the volume ${\hat{F}}$, with errors
controlled by the truncation $H^{\prime}$ of the CTF-values. That is to say, a
principal mode compression of (69) can be written as:

$\displaystyle{U}^{\intercal}\left[CTF_{j}({\boldsymbol{k}})\odot\left\\{{\hat{S}}\circ{\cal
R}(\tau_{j})\circ{\hat{F}}\right\\}({\boldsymbol{k}})\right]$
$\displaystyle\approx$
$\displaystyle\sum_{h^{\prime}=1}^{H^{\prime}}v^{\prime}_{h^{\prime}}(j)\sigma^{\prime}_{h^{\prime}}\cdot\left\\{{\hat{S}}\circ{\cal
R}(\tau_{j})\circ{\hat{G}}({\boldsymbol{k}};h^{\prime},{\hat{F}})\right\\}({\boldsymbol{k}})\text{.}$
(72)

This relationship can be used to align the principal images
${U}^{\intercal}{\hat{A}}_{j}$ directly to the solution ${\hat{G}}$, without
necessarily needing to reconstruct ${\hat{F}}$ (see section A.22).

#### Reconstruction via quadrature back propagation (QBP)

In our experience, an accurate approximation of the solution to the least-
squares problem above can require a significant amount of computation time
when the conjugate gradient method requires many iterations to converge. As a
more efficient (but typically less accurate) alternative, we employ a slightly
different strategy. We refer to this alternative as ‘quadrature-back-
propagation’ (QBP).

The basic observation underlying QBP is that, assuming that there was no image
noise and no CTF-attenuation, then the volume and images would solve:

$\left[\begin{array}[]{c}{\hat{Sl}}\circ{\cal R}(\tau_{1})\circ\\\ \vdots\\\
{\hat{Sl}}\circ{\cal R}(\tau_{N_{A}})\circ\\\
\end{array}\right]\cdot{\hat{F}}=\left[\begin{array}[]{c}{\cal
T}(-{\boldsymbol{\delta}}_{1},{\boldsymbol{k}})\odot{\hat{A}}_{1}({\boldsymbol{k}})\\\
\vdots\\\ {\cal
T}(-{\boldsymbol{\delta}}_{N_{A}},{\boldsymbol{k}})\odot{\hat{A}}_{N_{A}}({\boldsymbol{k}})\\\
\end{array}\right]$

with each template ${\hat{Sl}}\circ{\cal R}(\tau_{j})\circ{\hat{F}}$
corresponding to a cross-section of the volume ${\hat{F}}({\boldsymbol{k}})$.
Each of these cross-sections can be viewed as a sample of the overall volume
${\hat{F}}({\boldsymbol{k}})$, with values given by the appropriate entry on
the right-hand-side of the equation above. With sufficiently many cross-
sectional samples, a functional representation of the volume ${\hat{F}}$ can
be reconstructed via quadrature in the spherical harmonic basis via a
spherical harmonic transform. That is, we have the approximate formula

${f}(k;l,m)\approx\sum_{p}w_{p}\left[Y_{l}^{m}(\hat{\boldsymbol{k}}_{p})\right]^{\dagger}{\hat{F}}(k,\hat{\boldsymbol{k}}_{p})\
\ \forall\ \ l,\|m\|,l^{\prime},\|m^{\prime}\|\leq{L}\text{.}$

where $\\{\hat{\boldsymbol{k}}_{p},w_{p}\\}$ correspond to a quadrature rule
on the sphere. The number of such nodes is of the order $[\mathcal{O}(k)]^{2}$
on the spherical $k$-shell. (Equispaced points in the azimuthal direction and
Legendre points in the polar angle, for example, are spectrally accurate.)
Because the data points will not generally coincide with the necessary
quadrature nodes, we interpolate from the given data to the nodes in some
local neighborhood. A simple rule would be to take the local average over all
data points that lie within some distance of the node. Note that this strategy
is not iterative and is typically an order of magnitude (or so) faster than a
least-squares reconstruction.

In the absence of image-noise, the error in the reconstruction of ${f}(k;l,m)$
will come from the ‘interpolation error’ associated with interpolating from
the cross-sectional sample data to the quadrature nodes. When the image-noise
is nonzero and the interpolation regions are fixed in size, then the
approximation to ${\hat{F}}(k,\hat{\boldsymbol{k}}_{p})$ will exhibit two
sources of error. In addition to the interpolation error, there is a sampling
error due to the image noise at each cross-sectional sample-point. In
practice, we generally choosing the local neighborhood of the quadrature nodes
to be sufficiently large to contain roughly $\mathcal{O}(10^{1-2})$ sample
points. This is usually sufficient for two digits of accuracy (coinciding with
our usual global error tolerance).

The above strategy can easily be generalized to deal with image-specific CTF-
values. The QBP approach above essentially treats each sample within the local
interpolation region of each quadrature node as a noisy observation of
${\hat{F}}(k,\hat{\boldsymbol{k}}_{p})$. The average can be thought of as a
local maximum-likelihood estimate. Since different data points might
correspond to different CTF functions, we can simply compute the CTF-weighted
local averages instead. This strategy coincides with the CTF-weighting of
Bayesian inference described in [7, 11].

###  A.21 Remark regarding further template compression

If the ranks $H$ and $H^{\prime}$ are both sufficiently small (e.g., $2-4$ and
$1-2$, respectively), then the dimensionality of the space of CTF-modulated
principal templates (as measured via the numerical rank) is very small, and
usually quite a bit smaller than the actual number of principal templates
(i.e., the number of viewing-directions). This would allow for additional
compression, using a reduced set of ‘basis’-templates (which span the space of
templates) rather than the actual templates themselves. In our numerical
experiments we do not implement this compression because our $H$ and
$H^{\prime}$ are typically too large for it to be efficient.

###  A.22 Image alignment

Just as for the reconstruction step, image alignment is often performed
multiple times during molecular refinement. For standard AM, the inputs to the
alignment process are the estimated molecular volume ${\hat{F}}^{{\text{\tiny
estim}}}({\boldsymbol{k}})$, as well as the picked-particle images
${\hat{A}}_{j}({\boldsymbol{k}})$ and associated CTF-functions
$CTF_{j}({\boldsymbol{k}})$. The outputs are the estimated viewing angles
$\tau_{j}^{{\text{\tiny estim}}}$ and displacements
${\boldsymbol{k}}_{j}^{{\text{\tiny estim}}}$ for the various images.

In our case we will make use of both the standard maximum-likelihood
alignment, as well as a more stable ‘maximum-entropy’ alignment (described
below).

###  A.23 Alignment via maximum likelihood

For this, we calculate the quantities ${\cal X}$, $\hat{\cal X}$ and ${\cal
Z}$, described in section 2.1, for a collection of ${T}=\mathcal{O}({K}^{2})$
viewing directions $\alpha_{t}$ and $\beta_{t}$ corresponding to our spherical
shell quadrature grid for $\hat{\boldsymbol{k}}$ at $k=k_{{R}}\approx{K}$. For
each of these viewing-directions $(\alpha_{t},\beta_{t})$, we consider a range
of $O({K})$ values for $\gamma$ in $[0,2\pi)$. For the displacements, we
typically limit them to lie in a small disk of radius
${{\boldsymbol{\delta}}_{{\text{\tiny local}}}}\sim 0.01$ or so, taking
advantage of the low numerical rank of the set of translation-operators within
this disk [17, 25].

Given the correlations ${\cal Z}(\alpha_{t},\beta_{t};{\hat{A}}_{j})$, a
common strategy for updating the viewing angles $\tau_{j}^{{\text{\tiny
estim}}}$ and displacements ${\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}$
(referred to as ‘maximum-likelihood’ alignment), involves selecting the
$\tau_{j}^{{\text{\tiny estim}}}$ (and, as a result, the
${\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}$) for each image $j$ to
maximize ${\cal Z}(\alpha_{t},\beta_{t})$ over the ${T}$ viewing directions
for that particular image $j$. In the ideal situation where the image
formation model is accurate, this choice corresponds to the viewing-angle that
is most likely to have produced the particular image ${\hat{A}}_{j}$, given
that the molecule is indeed ${\hat{F}}^{{\text{\tiny estim}}}$. In practice we
do not define $\tau_{j}^{{\text{\tiny estim}}}$ to be the very best $\alpha$
and $\beta$ for each image, instead randomly selecting an $\alpha_{t}$ and
$\beta_{t}$ for which ${\cal Z}(\alpha_{t},\beta_{t})$ is within the top
$5^{\text{th}}$-percentile of its values for that image [16]. The results of
our numerical-experiments are quite insensitive to this value of $5\%$; any
choice from $1\%$ to $\sim 10\%$ achieves similar performance for maximum-
likelihood alternating minimization.

Let us define the function ${\text{\cal Sort}}[\cdot]$ which maps lists of
real-numbers to lists of positive integers such that, for any vector
$\boldsymbol{a}\in\mathbb{R}^{L}$, the value of ${\text{\cal
Sort}}[\boldsymbol{a}]_{l}$ is the (integer) rank of $a_{l}$ (ranging from
$1,\ldots,L$ as it stands within the entries of $\boldsymbol{a}$). As noted in
section 2.1, we can see that the template ranks are obtained by applying
${\text{\cal Sort}}[\cdot]$ to each ‘column’ of the array ${\cal
Z}(\alpha_{t},\beta_{t};{\hat{A}}_{j})$ (fixing $j$), while the image ranks
are obtained by applying ${\text{\cal Sort}}[\cdot]$ to each ‘row’ of the same
array (fixing $t$). The maximum-likelihood alignment described above depends
on the template ranks, while the maximum-entropy alignment (described below)
depends on the image ranks.

###  A.24 Alignment via maximum-entropy

The maximum-likelihood alignment described above can be thought of as a
‘greedy’ algorithm which selects the estimated viewing angles for each image
without consideration of the choices made for the other images. As a result of
this greedy strategy, the estimated viewing angle distribution (across all
$\tau_{j}^{{\text{\tiny estim}}}$) is not controlled, and can end up becoming
quite far from the true viewing-angle distribution (until the volume used to
generate the templates is sufficiently accurate).

To avoid this phenomenon, one can choose the $\tau_{j}^{{\text{\tiny estim}}}$
while constraining the distribution of estimated viewing angles in some way. A
simple special case of this strategy is to constrain the distribution of
estimated viewing angles to be uniform. We refer to such an alignment strategy
as ‘maximum-entropy’ alignment (referring to the entropy of
$\mu_{\tau}^{{\text{\tiny estim}}}$). In practice we implement this maximum-
entropy alignment as follows:

  1. 1.

Declare all the images ‘unassigned’.

  2. 2.

Define a list of $N_{A}$ uniformly distributed viewing-directions
$\hat{\boldsymbol{k}}_{j}=\left(\alpha_{j},\beta_{j}\right)$ on the sphere.
(This is typically done by cycling through a randomly-permuted list of the
${T}$ viewing-angles $(\alpha_{t},\beta_{t})$ until $N_{A}$ viewing-directions
have been selected).

  3. 3.

Step through the list of viewing directions. For each viewing-direction
(defined by a particular $\alpha$ and $\beta$), assign the image
${\hat{A}}_{j}$ that has not previously been assigned and which maximizes
${\cal Z}(\alpha,\beta;{\hat{A}}_{j})$ (over the unused image-indices $j$).
When assigning an image ${\hat{A}}_{j}$ to a particular viewing direction
$\alpha$, $\beta$, we set the in-plane rotation $\gamma_{j}^{{\text{\tiny
estim}}}$ and displacement ${\boldsymbol{\delta}}_{j}^{{\text{\tiny estim}}}$
implicitly via the $\operatorname*{arg\,max}$ (see the definition of ${\cal
Z}$ in Eq 9).

Note that the third step of this procedure uses the image ranks for each
template (rather than the template ranks for each image). In some cases the
image ranks and template ranks contain similar information (see, e.g., the
results for ISWINCP in Fig 2). However, it is more typical for the image ranks
to be more informative than the template ranks (see the other case studies in
Figs. 2 and 11).

Intuition might suggest that maximum-entropy alignment is a reasonable
strategy if the true viewing angle distribution is known to be uniform, and
not otherwise. Nevertheless, as illustrated in the main text, our EMPM
iteration still often produces results that are significantly better than
maximum-likelihood alignment, even for cryo-EM data where the true viewing
angle distribution is rather nonuniform. We study this behavior in some detail
for simplified models in sections B and C below.

##  Appendix B Multi-Reference Alignment (MRA)

Multi-Reference Alignment (MRA) is an abstraction of the cryo-EM molecular
reconstruction problem, reduced to a two-dimensional setting.

In this context the image displacements and CTF-modulation are typically
ignored, and we assume (i) the imaged molecule itself is translationally
invariant in the $x_{3}$-direction, and (ii) the imaged molecule is only
observed from the $x_{3}$-axis. Because of (i), the function
${{F}}({\boldsymbol{x}})$ will be determined by ${{F}}(x_{1},x_{2})$, and
${\hat{F}}({\boldsymbol{k}})$ will be restricted to the equatorial-plane, as
denoted by ${\hat{F}}(k,\psi)$. Additionally, due to (ii), each ‘image’
${\hat{A}}_{j}$ will be a noisy version of ${\hat{F}}(k,\psi)$, subject to an
unknown in-plane rotation $\gamma_{j}$:

${\hat{A}}_{j}(k,\psi)=[{\cal
R}(\gamma_{j})\circ{\hat{F}}](k,\psi)+\text{noise}\
={\hat{F}}(k,\psi-\gamma_{j})+\text{noise}.$

If we further assume that the true volume ${\hat{F}}$ involves only a single
$k$-shell, then each image is restricted to a single $k$-ring, and can be
considered a periodic function on $[0,2\pi]$.

The MRA problem has been well studied [26, 27, 28], and can be used to gain
insight into the behavior of many reconstruction strategies. It can be used,
for example, to investigate the sensitivity of Bayesian inference (BI) to
mismatches between the estimated and true noise levels (see Fig 6).

##  Appendix C Multi-Slice Alignment (MSA)

An ‘orthogonal’ sub-problem to MRA, which we refer to as Multi-Slice Alignment
(MSA), can be described as follows. As in MRA, we ignore displacements and the
CTF and assume that the volume is invariant in the $x_{3}$-direction. However,
in MSA, we assume that the imaged molecule is only observed from the
equatorial plane. With this assumption the images will correspond to
projections of ${{F}}({\boldsymbol{x}})$ onto one-dimensional lines in the
$x_{1}$-$x_{2}$-plane (i.e., samples of the two-dimensional Radon transform of
${{F}}(x_{1},x_{2})$). Thus, retaining the notation ${\hat{A}}_{j}$, each
image corresponds to a one-dimensional slice of ${\hat{F}}(k,\psi)$, expressed
as a function of $k\in\mathbb{R}$ for a specific (but unknown) $\psi_{j}$:

${\hat{A}}_{j}(k)={\hat{F}}(|k|,\psi_{j})+\text{noise}\text{\ if\ }k>0\text{,\
or\ \ }{\hat{F}}(|k|,\psi_{j}+\pi)+\text{noise}\text{\ if\ }k<0\text{.}$

To simplify the scenario further, we restrict each image to the ray of
positive $k$:

${\hat{A}}_{j}(k)={\hat{F}}(k,\psi_{j})+\text{noise}\quad\text{for
$k>0$}\text{.}$

If we further assume that the true volume ${{F}}$ can take on arbitrary
complex values, and that ${\hat{F}}$ involves only a single $k$-shell, then
${\hat{F}}$ corresponds to a closed-curve in $\mathbb{C}$, written as the
function ${\hat{G}}(\psi)$, parametrized by $\psi$. Each image is now simply a
point, corresponding to a noisy observation of ${\hat{G}}(\psi)$. If we assume
that the true volume ${\hat{F}}$ involves a discrete set of ${R}$ $k$-values
$\\{k_{1},\ldots,k_{{R}}\\}$, then ${\hat{G}}$ will be a closed curve in
$\mathbb{C}^{{R}}$ (once again parametrized by $\psi$), and each image will be
a vector in $\mathbb{C}^{{R}}$ corresponding to a noisy observation of
${\hat{G}}(\psi_{j})$. An example illustrating the MSA problem for ${R}=1$ is
shown in Fig 7.

#### Numerical experiments:

We believe that MSA can be a useful testbed for studying the performance of
various strategies for cryo-EM reconstruction. As a simple example, let
${R}\equiv 1$ above and define ${\hat{G}}^{{\text{\tiny
true}}}:[0,2\pi)\rightarrow\mathbb{C}$ using a fixed number of Fourier modes:

${\hat{G}}^{{\text{\tiny true}}}(\psi)=\sum_{q=-{Q}}^{+{Q}}{g}^{{\text{\tiny
true}}}(q)\cdot\exp(iq\psi)\text{,}$

with the shape of ${\hat{G}}^{{\text{\tiny true}}}$ determined by the $2{Q}+1$
coefficients ${g}^{{\text{\tiny true}}}_{q}$. In the experiments below, we
generate ${\hat{G}}^{{\text{\tiny true}}}$ randomly by drawing each of the
${g}^{{\text{\tiny true}}}_{q}$ from a (complex) Gaussian distribution with
variance $1$.

Having drawn a particular ${\hat{G}}^{{\text{\tiny true}}}$, we sample this
molecule at $N_{A}$ points, producing a collection of simulated data points
${\hat{A}}_{j}$ defined by:

${\hat{A}}_{j}={\hat{G}}^{{\text{\tiny true}}}(\psi_{j})+\epsilon_{j}\text{,}$

where the $\psi_{j}$ are drawn independently and uniformly from $[0,2\pi)$,
and each $\epsilon_{j}$ is an independent complex random variable drawn from a
normal distribution with variance $\sigma$.

Now, given the collection of $N_{A}$ images $\\{{\hat{A}}_{j}\\}$, we can try
to reconstruct ${\hat{G}}^{{\text{\tiny true}}}$ using the strategies
discussed in the main text. Due to the simplicity of MSA, much of the
apparatus described in the main text simplifies dramatically. For example,
reconstruction involves a subsampled and/or non-uniform Fourier transform, and
alignment involves calculating the distance between each image point and the
curve ${\hat{G}}$ as a function of $\psi$.

Many features of this toy problem are expected. For example, when $\sigma=0$
the images lie exactly on the curve ${\hat{G}}^{{\text{\tiny true}}}$. If we
begin the alignment process with an initial guess ${\hat{G}}^{{\text{\tiny
init}}}={\hat{G}}^{{\text{\tiny true}}}$ for the molecule, then each image
will be aligned exactly to its correct viewing-angle $\psi_{j}^{{\text{\tiny
true}}}$ (up to a global phase-shift). Subsequent reconstruction steps will
produce ${\hat{G}}^{{\text{\tiny estim}}}={\hat{G}}^{{\text{\tiny true}}}$, so
long as the number of observed points $N_{A}$ is at least as large as the
number of modes $1+2{Q}$, and the reconstruction operator is well-conditioned
(i.e., if the image-points are not too closely clumped together). Put more
simply, if $\sigma=0$ and $N_{A}\geq 1+2{Q}$ then, generically speaking, a
maximum-likelihood alternating minimization starting with
${\hat{G}}^{{\text{\tiny init}}}={\hat{G}}^{{\text{\tiny true}}}$ as an
initial-function will be stable.

However, even in this simple scenario (i.e., with $\sigma=0$), if we were to
attempt maximum-likelihood alternating minimization with a different initial
function then we are not guaranteed to converge to ${\hat{G}}$. Indeed, if
$N_{A}$ is not too large, there will typically be many different functions
with the same bandlimit of ${Q}$ that pass close to the observed points
${\hat{A}}_{j}$. Which of these functions we converge to will depend on the
details of the reconstruction process. The situation becomes more complicated
once we add image noise; when $\sigma>0$ even an initial guess of
${\hat{G}}^{{\text{\tiny init}}}={\hat{G}}^{{\text{\tiny true}}}$ will not
necessarily lead to convergence to a function ${\hat{G}}^{{\text{\tiny
estim}}}$ that is close to ${\hat{G}}^{{\text{\tiny true}}}$.

To probe this phenomenon, we perform a set of numerical-experiments. We
initialize our reconstruction with a function of the form:

${\hat{G}}^{{\text{\tiny init}}}:=\lambda\cdot{\hat{G}}^{{\text{\tiny
true}}}+(1-\lambda)\cdot{\hat{G}}^{{\text{\tiny rand}}}\text{,}$

where ${\hat{G}}^{{\text{\tiny rand}}}$ is a random function drawn the same
way as (but independent from) ${\hat{G}}^{{\text{\tiny true}}}$, and $\lambda$
is a parameter describing how close ${\hat{G}}^{{\text{\tiny init}}}$ is to
the ground truth. Given a particular ${\hat{G}}^{{\text{\tiny init}}}$, we
perform the reconstruction, iterating until the resulting
${\hat{G}}^{{\text{\tiny estim}}}$ converges. For each trial we can then
measure the error $E$ between ${\hat{G}}^{{\text{\tiny estim}}}$ and
${\hat{G}}^{{\text{\tiny true}}}$.

In order to reward trials where ${\hat{G}}^{{\text{\tiny estim}}}$ and
${\hat{G}}^{{\text{\tiny true}}}$ are close to one another, even if they are
parametrized differently, we define the error $E$ as follows: First we define
$d({\hat{G}}^{{\text{\tiny estim}}}\leftarrow{\hat{G}}^{{\text{\tiny
true}}}(\psi))$ to be the distance between ${\hat{G}}^{{\text{\tiny
true}}}(\psi)$ (i.e., a point in $\mathbb{C}$) and the curve
${\hat{G}}^{{\text{\tiny estim}}}$. We then define the error
$E({\hat{G}}^{{\text{\tiny estim}}}\leftarrow{\hat{G}}^{{\text{\tiny true}}})$
as the integral of $d^{2}({\hat{G}}^{{\text{\tiny
estim}}}\leftarrow{\hat{G}}^{{\text{\tiny true}}}(\psi))$ over $\psi$.
Similarly, we define $d({\hat{G}}^{{\text{\tiny
true}}}\leftarrow{\hat{G}}^{{\text{\tiny estim}}}(\psi))$ to be the distance
between ${\hat{G}}^{{\text{\tiny estim}}}(\psi)$ and the curve
${\hat{G}}^{{\text{\tiny true}}}$, and define the error
$E({\hat{G}}^{{\text{\tiny true}}}\leftarrow{\hat{G}}^{{\text{\tiny estim}}})$
as the integral of $d^{2}({\hat{G}}^{{\text{\tiny
true}}}\leftarrow{\hat{G}}^{{\text{\tiny estim}}})$ over $\psi$. Finally, we
define the overall error $E=\max(E({\hat{G}}^{{\text{\tiny
true}}}\leftarrow{\hat{G}}^{{\text{\tiny estim}}}),E({\hat{G}}^{{\text{\tiny
estim}}}\leftarrow{\hat{G}}^{{\text{\tiny true}}}))$ to be the larger of the
two.

Figure 8 illustrates the average $E$ (accumulated over $512$ random trials)
for three different reconstruction strategies across a range of $\lambda$ and
$\sigma$ values. On the left we show the results for Bayesian inference, where
we set the estimated noise level equal to the true noise level (i.e.,
$\sigma$) when reconstructing the molecule. Note that Bayesian inference is
quite accurate when the initial guess is accurate (i.e., when $\lambda=1$).
However, when $\lambda<1$ and there is a mismatch between the initial guess
and the true function, Bayesian inference typically performs quite poorly.

Similarly, we can apply standard maximum-likelihood alternating minimization
(AM). This is equivalent to Bayesian inference with an estimated noise level
of $0$. The results for AM are identical to those shown on the left when
$\sigma=0$. However, as soon as a little image noise is added ($\sigma>0$), AM
typically converges to a ${\hat{G}}^{{\text{\tiny estim}}}$ which is wildly
different from ${\hat{G}}^{{\text{\tiny true}}}$, with errors that are an
order of magnitude larger than those shown for Bayesian inference (which was
run with an estimated noise-level of $\sigma$). We do not display the (large)
errors observed for standard AM.

In the middle, we show the results for a simple alternating minimization
scheme using maximum-entropy alignment at each step. On the right we show the
results for our preferred EMPM strategy: namely, alternating between maximum-
likelihood and maximum-entropy alignment. Note that, while our preferred
strategy is less accurate than Bayesian inference when $\sigma$ is small and
${\hat{G}}^{{\text{\tiny init}}}={\hat{G}}^{{\text{\tiny true}}}$, it is more
accurate in this regime when $\sigma>0$ and $\lambda<1$.

Figure 9 illustrates the results of a similar numerical experiment, where the
images are drawn from a non-uniform distribution of viewing angles. For this
numerical experiment, we consider viewing angle distributions of the form
$\exp(\beta\cos(\psi))/(2\pi{\cal I}_{0})$, controlled by a parameter $\beta$.
When $\beta=0$ the viewing angle distribution is uniform, and as $\beta$
increases, the viewing angle distribution becomes more and more sharply peaked
around $\psi=0$. Results for $\lambda=0$ are shown in Fig 9. The results with
$\beta=0$ are the same (in distribution) as the $\lambda=0$ scenario in Fig 8.
Note that EMPM is typically more accurate than Bayesian inference in this
regime.

##  Appendix D Relion De-Novo parameters:

When using Relion’s de-novo molecular reconstruction, we use the following
parameters.

\--sgd_ini_iter 50 | \--sgd_inbetween_iter 200
---|---
\--sgd_fin_iter 50 | \--sgd_write_iter 20
\--sgd_ini_resol 35 | \--sgd_fin_resol 15
\--sgd_ini_subset 100 | \--sgd_fin_subset 500
\--sgd | \--denovo_3dref
\--ctf | \--K 1
\--sym C1 | \--flatten_solvent
\--zero_mask | \--dont_combine_weights_via_disc
\--pool 3 | \--pad 1
\--skip_gridding | \--oversampling 1
\--healpix_order 1 | \--offset_range 8
\--offset_step 2 | \--j 6

Additionally, we set the \--particle_diameter  to be commensurate with the
true molecule diameter (rounded up to the nearest $50$ or $100$Å).