# Deep Generative Modeling for Volume Reconstruction in Cryo-Electron
Microscopy

Claire Donnat  Axel Levy  Stanford University, Department of Electrical
Engineering, Stanford, CA, USA  LCLS, SLAC National Accelerator Laboratory,
Menlo Park, CA, USA  Frédéric Poitevin  LCLS, SLAC National Accelerator
Laboratory, Menlo Park, CA, USA  Ellen Zhong  Massachusetts Institute of
Technology, Computer Science and Artificial Intelligence Lab, Boston, MA, USA
Nina Miolane

###### Abstract

Recent breakthroughs in high-resolution imaging of biomolecules in solution
with cryo-electron microscopy (cryo-EM) have unlocked new doors for the
reconstruction of molecular volumes, thereby promising further advances in
biology, chemistry, and pharmacological research. Recent next-generation
volume reconstruction algorithms that combine generative modeling with end-to-
end unsupervised deep learning techniques have shown promising preliminary
results, but still face considerable technical and theoretical hurdles when
applied to experimental cryo-EM images. In light of the proliferation of such
methods, we propose here a critical review of recent advances in the field of
deep generative modeling for cryo-EM volume reconstruction. The present review
aims to (i) unify and compare these new methods using a consistent statistical
framework, (ii) present them using a terminology familiar to machine learning
researchers and computational biologists with no specific background in cryo-
EM, and (iii) provide the necessary perspective on current advances to
highlight their relative strengths and weaknesses, along with outstanding
bottlenecks and avenues for improvements in the field. This review might also
raise the interest of computer vision practitioners, as it highlights
significant limits of deep generative models in low signal-to-noise regimes —
therefore emphasizing a need for new theoretical and methodological
developments.

## Introduction

Figure 1: Acquisition of 2D cryo-EM images (2D projections) from 3D
biomolecular volumes.

High-resolution reconstruction of molecular volumes from single particle
images has the potential to facilitate new breakthroughs in our ability to
understand fundamental biological mechanisms and engineer macromolecular
function [7, 52]. In this context, cryo-electron microscopy (cryo-EM) has
fostered a revolution in structural biology by allowing the imaging of
biomolecules in solution at atomic resolution [8, 9]. However, the estimation
of these molecules’ 3-dimensional (3D) volume from cryo-EM data continues to
pose a formidable challenge. In this setting, observations are limited to the
raw 2D projections of molecules (also called particles) relative to an
incoming electron beam, while their 3D orientation and position (jointly
called poses) are unknown — see Figure 1. Reconstructing molecular volumes
therefore also requires recovering a number of hidden variables such as each
particle’s 3D orientation. The difficulty of this task is further compounded
by a combination of factors, including the variability in the shape of any
given molecule (also referred to as structural “heterogeneity”), the non-
linear physics of the data acquisition process, as well as extremely low
signal-to-noise ratios — concepts formalized in the image formation model
below.

##### Image Formation Model.

The process of image formation in cryo-EM involves several physical phenomena,
including pairwise interactions between atoms, interactions between the
electron beam and the molecule’s electrostatic potential, or microscope
effects. We refer the reader to Dill et al.[10], Kohl and Reimer[11], and
Vulovic et al.[12] for in-depth descriptions of these phenomena. Nonetheless,
in most cases[13, 12], each image $X_{i}$ in a dataset of $n$ images of single
particles can be modeled as a random sample from the following generative
model:

$X_{i}=\text{PSF}_{i}*(t_{i}\circ\Pi_{2D}\circ
R_{i})(V^{(i)})+\epsilon_{i},\qquad\text{with}\quad i=1\cdots n.$ (1)

Here, $R_{i}$ is a 3D rotation representing the 3D orientation of the volume
$V^{(i)}$ with respect to the direction of the electron beam. The oriented
volume is subsequently “pierced through” by the electron beam and projected
onto the detector — an operation represented in Equation (1) by the
2D-projection operator $\Pi_{2D}$. The variable $t_{i}$ represents the 2D
translation of the projected volume with respect to the center of the image.
The effect of the microscope’s lens is modeled through the convolution $*$ of
the 2D projection by an image-dependent operator $\text{PSF}_{i}$ called the
Point Spread Function (PSF) of the microscope whose parameters can depend on
the image. Finally, additional noise $\epsilon_{i}$ is introduced in the
observed image, and typically assumed to be Gaussian with zero mean and
variance $\sigma_{i}^{2}$. Note that the underlying volume $V^{(i)}$ is
allowed to depend on $i$. This allows us to account for “conformational
heterogeneity”, a concept whereby a molecule does not necessarily exist in a
single state, but rather, that its volume corresponds in fact to one of
several stable geometries that can be achieved by this molecule.

An equivalent generative model can be formulated via the Fourier transform of
Equation (1) which we present in Appendix Appendix B: Differentiability of the
Generative Models. In Fourier Space, the convolution $*$ of the projected
volume by the PSF becomes a computationally lightweight element-wise matrix
multiplication $\odot$ between the 2D Fourier transform of the projected image
and that of the PSF (known as the Contrast Transfer Function, or CTF).
Operating in Fourier space is thus common in many cryo-EM volume
reconstruction algorithms.

##### Challenges of Cryo-EM and Deep Generative Modeling.

Whether in image or Fourier space, the image formation model described by
Equation (1) gives us the first clue about the difficulty of the
reconstruction problem. Only the recovery of the set of possible 3D molecular
volumes $V^{(i)}$ (called the conformational landscape) is relevant - yet the
image also depends on a number of additional unknown nuisance variables, such
as the pose ($R_{i}$, $t_{i}$) or the microscope-dependent $\text{PSF}_{i}$.
As noted by Singer et al. [14], these nuisance variables do not bring any
information of biological relevance but may impact the quality (and
difficulty) of the reconstruction. Reconstructing molecular volumes thus
becomes a highly non-convex optimization problem, putting algorithms at risk
of being overly sensitive to initialization [17] or converging to one of many
local minima [19]. The difficulty of the task is also compounded by the high
levels of radiation damage, structured (non-white) noise and altogether
remarkably low signal-to-noise ratios that cryo-EM images usually exhibit.
Given these difficulties, reconstruction algorithms rely on using copious
amounts of data, but often struggle to process the wealth of incoming
images[15], to the point where “the cost of image analysis can exceed 500,000
CPU hours on large, expensive computer clusters”[16] for a single experiment.
In this context, recent efforts have turned to unsupervised deep learning for
cryo-EM reconstruction, i.e. approaches allowing significant speed-ups through
the use of gradient methods and GPUs. Given their potential to advance the
field by addressing the challenges mentioned here, we propose here a critical
review of deep generative modeling for cryo-EM volume reconstruction.

##### Related Work.

Several reviews have already begun surveying challenges and advances in cryo-
EM reconstruction. The reviews by Singer et al.[14] and Bendory et al.[20]
provide a complete description of cryo-EM reconstruction, but focus on
mathematical foundations of general computational methods, rather than
specifically on deep learning approaches. Reviews by Si et al.[22], Ede et al.
[26] and Wu et al. [27] describe the pervasiveness of deep learning methods
along all steps of the cryo-EM pipeline, without specialising to high-
resolution volume reconstruction. By contrast, our review is a deep dive into
the most recent deep generative models for cryo-EM reconstruction.

##### Objective and Contributions.

Our objective is to shed light into the similarities and differences among
recent state-of-the-art, deep-generative reconstruction methods, which we
classify according to (i) their parametrization of the generative model
(Section 1) and (ii) the inference tools deployed to fit this generative model
(Section 2). This unification of recent works along a consistent statistical
framework allows us to highlight trends, outstanding challenges, and avenues
for improvements in the field (Section 3). We draw insights from our
experiences with these approaches to provide an account of the hurdles and
potential difficulties that arise in their deployment to real cryo-EM data.
Finally, through this review, we hope to catalyze deep learning advances by
providing machine learning practitioners and computer vision experts a
thorough overview of the challenges that are unique to cryo-EM.

##  1 Generative Modeling for Cryo-EM

The objective of cryo-EM imaging algorithms is to produce a 3D reconstruction
of a given molecule from a dataset of images $\\{X_{i}\\}_{i=1\cdots n}$,
where each image corresponds to a “2D projection” of the molecule at a
different (unknown) orientation and position (Figure 1). A fundamental hurdle
to this objective lies in the fact that each molecule has its own unknown
conformation (or shape). Methods that account for this heterogeneity are
called heterogeneous reconstruction methods, and can typically aspire to
higher resolution reconstructions — often at the cost of more involved and
expensive computations. Conversely, homogeneous reconstruction methods neglect
this shape variability and usually represent the conformational landscape as a
unique molecular volume. This section reviews how choices both in accounting
for conformation heterogeneity and in parametrizing the volume yield different
formulations of the cryo-EM image formation model.

###  1.1 Conformation Variable $z$.

Heterogeneous reconstruction methods introduce an additional variable $z_{i}$
for each image $i$ within the formation model of Equations (1), which we call
the conformation variable. Depending on whether conformation heterogeneity is
modeled through a discrete number of states or as a continuous variable, the
conformational landscape can be encoded as a discrete family of volumes
$\mathcal{V}=\\{V_{z},\leavevmode\nobreak\ z\in\\{1,\ldots,K\\}\\}$ (discrete
heterogeneity), or as a continuous family
$\mathcal{V}=\\{V(z),\leavevmode\nobreak\ z\in\mathbb{R}^{L}\\}$ for some
integer $L$ (continuous heterogeneity)[1]. In both cases, the family is
indexed by the variable $z$. We use the notations $V_{z}$ and $V(z)$
interchangeably. Homogeneous reconstruction can in fact be taken as the
special case where $\mathcal{V}$ only comprises a single volume, so that $K=1$
or $L=0$ (i.e. forcing $z=0$). We now write $V^{(i)}=V(z_{i})$ in Equation
(1).

##### Interpretation of the Conformation Variable

From a statistical mechanics perspective, the conformation variable $z_{i}$
encodes the location of any given single particle along the conformational
landscape[10]. For example, if $z_{i}\in\mathbb{R}$, $z_{i}$ can be used to
sort conformations along a “conformation path”, that is, a sequence of small
transformations that would interpolate two main preferred, dynamically stable
states (or conformations) for the molecule. When continuous, the dimension $L$
of this cursor variable $z$ could in principle take any value between 0 (no
heterogeneity) and $O(N)$, with $N$ the number of atoms in the molecule.
However, two factors tend to drastically limit the number of dimensions of
$z$. First, most of the main global dynamics of a molecule are captured by a
few collective variables associated with its low-frequency motion, effectively
averaging out a lot of the effects of the high number of degrees of freedom
associated with faster motions [2]. Second, imaging conditions often reduce
the ability to resolve the remaining motions, thus reducing the effective
dimensionality of $z$[3]. In other words, limits in the imaging technology
itself restrict the dimension of the variable $z$. In the case of discrete
heterogeneity, $z\in\\{1,..,K\\}$ is an index of minimum energy wells
(conformations) in the conformational landscape. Imaging conditions also
reduce the ability to resolve too many metastable states, thereby restricting
practitioners to choose a low value for $K$.

##### Discrete vs Continuous Conformational Heterogeneity: Pros, Cons and
Discussion.

Discrete heterogeneity has a rich history in cryo-EM. Popularized by the
Class3D[13, 30] extension of RELION[29] , it offers the advantage of
delivering readily interpretable results: a set of $K$ volumes, representing
$K$ main stable states of the molecules. Discrete heterogeneity is thus
particularly adequate in certain (common) scenarios where the conformation
landscape has local energy minima that produce distinct states. However, one
of the main drawbacks of this method consists in the necessary selection of
the number $K$ of appropriate conformations. Theoretically, this could be done
by cross-validation. In practice, due to the significant computing costs that
cross-validation implies, $K$ is chosen in an ad hoc fashion by the
experimenter and rarely motivated by strong quantitative arguments (see
Haselbach _et al_[4] for a rare example).

Consequently, many recent methods have turned to a continuous representation
of heterogeneity which does not require specifying a number $K$ of
conformations. This representation is also often deemed to be closer to the
underlying biology, as molecules do not exist as finite/discrete sets of
shapes. Rather, a more realistic analogy is to think of molecules as random
samples from the equilibrium distribution over their conformational space[10].
However, while a continuous representation could be more scientifically
relevant, it remains to be determined how accurate the reconstruction of the
conformational space by the space indexed by $z$ truly is. This latter point
will be critical to address for heterogeneous reconstruction methods to become
more quantitative and directly comparable to other measures from biophysicists
and biochemists. We discuss in section 3 the challenges of assessing the
precision of such approaches, which probably constitutes one of the main open
questions in the field. Additionally, despite its initial appeal, continuous
conformation heterogeneity comes with significant theoretical and practical
caveats. From a physics perspective, it is still unclear whether the full
landscape (at room temperature) is sufficiently well sampled by cryo-EM to
justify modeling conformations with a continuous rather than discrete
distribution: the sample preparation process in cryo-EM — and most
specifically the grid-freezing step— affects the distribution of conformations
which might not reflect the heterogeneity of conformations at room
temperature[5]. From a statistical perspective, using a continuous
distribution necessitates the generative model to be able to sample from the
full conformation landscape, a requirement that is itself a considerable
challenge for large molecules: the strong constraints, e.g. on on bond lengths
and torsion angles, make up for a complex, non-convex landscape that is
difficult to sample from. Despite these caveats, Table 3 shows that continuous
heterogeneity is gaining traction amongst the most recent reconstruction
advances.

###  1.2 Molecular Volume $V(z)$

The cryo-EM reconstruction problem can thus be understood as recovery of the
underlying conformational landscape $\mathcal{V}$ and the corresponding
probability distribution. The next critical step thus consists in finding a
judicious parametrization for each volume $V(z)\in\mathcal{V}$. This requires
choosing first an “output space” (image space vs. Fourier space, inducing real
vs complex values), second, an “encoding style” (reference-free vs. reference-
based) and third, and “input domain” (continuous vs. discrete).

####  1.2.1 Defining the output space: image space or Fourier space

Equations (1)-(9) show that the image formation model can be described
equivalently in image space or Fourier space. Thus, each volume within the
family of conformations $\mathcal{V}$ can be described either in terms of its
pixel intensities or its Fourier coefficients. In either case, the volume
$V(z)$ associated to the conformation variable $z$ is defined on an input
domain $\Omega\subset\mathbb{R}^{3}$ (the space of coordinates) and outputs
values in an output space that is either $\mathbb{R}$ for pixel intensities
representing the electron scattering potential of the molecule, or
$\mathbb{C}$ to encode the amplitude and phase of the Fourier coefficients.

##### Image versus Fourier space: Pros, Cons and Discussion.

From a practical standpoint, the choice of the output space is guided by the
set of properties and constraints that the analyst wishes to use to guide
volume reconstruction. Historically, the Fourier approach has been preferred.
As summarized by Punjani et al. [32], working in Fourier space has the
benefits of (a) reducing the computational cost of the image formation model
(see discussion of Equation 9), and (b) allowing closed-form maximum
likelihood reconstructions when molecules’ orientations and positions are
known. However, recent methods such as 3DFlex[32] have favored image space,
where constraints (e.g. smoothness of the deformation, conservation of energy,
etc.) are more interpretable and where operations such as interpolation and
deformation of the molecule’s density map are more naturally parametrized —
whereas the same operations require a careful treatment in Fourier space. For
example, interpolation in Fourier space can introduce unwanted artifacts. As
highlighted in Table 3, image space computations constitute a promising and
increasingly popular avenue for future developments in cryo-EM reconstruction.

####  1.2.2 Defining an encoding: reference-free or reference-based volume

The next step lies in the choice of an “encoding" for the volume $V(z)$. Cryo-
EM analysts typically have two choices: (i) using a reference-based
parametrization, which encodes the conformation landscape through its
deviation $\Delta V(z)$ from a reference conformation $V_{0}$, such that
$V(z)=V_{0}+\Delta V(z)$; or (ii) using a reference-free parametrization which
directly describes each $V(z)$, for instance as a set of atomic coordinates or
a low-dimensional embedding, but with no notion of “reference” conformation.

##### Reference-based versus reference-free: Pros, Cons and Discussion

If the column “Reference Volume" of Table 3 reflects the historical popularity
of reference-free encodings, the most recent methods relying on deep-learning
seem to have favored a reference-based approach. For instance, E2GMM[33] first
learns a reference $V_{0}$ called the “neutral representation" which then
serves in a reference-based encoding of $V(z)$ to further refine the
reconstruction by accounting for conformational variability. In AtomVAE[17] ,
Rosenbaum et al. uses a $V_{0}$ called a “base conformation" described as a
set of atom coordinates, obtained from an auxiliary method (such as an
homogeneous reconstruction or a set of atom coordinates predicted by AlphaFold
[37]). The existing reference acts as a statistical prior on the molecular
volume, thereby further constraining and guiding the recovery of the
conformation landscape. By contrast, 3DFlex[32] uses a reference volume
$V_{0}$, called a “canonical density", which is learned jointly with the
conformational heterogeneity. This has the advantage of foregoing the need to
split the pipeline in sequential steps, while allowing to borrow strength from
the joint estimation of all parameters.

Constraining the conformation recovery using a reference offers significant
advantages for ensuring the success (and convergence) of these methods given
the non-convexity of the problem. This template can be either learned (ab
initio methods), or chosen from existing data (refinement methods — more on
this in Appendix Appendix A: Ab-initio vs Refinement Reconstructions). The
general agreement across all methods consists in tackling this hierarchically,
starting with parameters which have the strongest impact on the signal, such
as defocus or pose, and gradually focusing on those whose effect is more
subtle, such as local deformations. As such, biasing the solution $V$ towards
a reference $V_{0}$, such that $\Delta V(z)=0$ implies $V(z)=V_{0}$, can
provide an interesting way of ensuring a more reliable and consistent — but
potentially biased — solution. Depending on the optimization method used, this
can in fact be critical to the success of the pipeline: Rosenbaum et al.[17]
report that adopting a reference template and warm-starting their algorithm is
indispensable to ensure the recovery of good conformations. However, because
they fundamentally bias conformations $V(z)$ to “hover” around $V_{0}$, the
success of such methods necessitates a reliable $V_{0}$. This can also incur
higher computational costs, since such methods typically require running a
first reconstruction method. This explains the interest for alternative,
reference-free methods: three out of the six heterogeneous methods in Table 3
allow to recover molecular volumes without any prior template. The extent to
which these reference-free methods are likely to succeed on real-images still
remains to be characterized.

####  1.2.3 Defining the input domain: discrete or continuous

The volume $V(z)$ represents a scalar 3D field (electrostatic potential, or
its Fourier transform) and is defined as a function from the input domain
$\Omega\subset\mathbb{R}^{3}$ to an output space $\mathbb{R}$ (or
$\mathbb{C}$). We now describe the parametrization of $\Omega$ and distinguish
two cases, depending on whether the volume is defined as a discretized or as a
continuous scalar field.

##### Discretized Domain and Explicit Parametrization.

The first class of approaches models the electrostatic potential as a discrete
3D map. In this case, the function $V(z)$ is defined on a discretized subspace
(a grid) of $\mathbb{R}^{3}$, namely $\Omega=\\{1,\ldots,D\\}^{3}$, where $D$
represents the length of the 3D voxel grid or frequency grid. $V(z)$ is
explicitly parametrized by the values it takes at each location (or voxel) of
$\Omega$. This choice is also called an explicit parametrization, a term that
will become clear in the next paragraph. Using a vectorial formalism, the
vector $V(z)$ corresponds to voxels’ intensity values, with
$V(z)\in\mathbb{R}^{D^{3}}$ or $\in\mathbb{C}^{D^{3}}$. In this case, the
resolution of the reconstructed volume is fixed by the choice of the
granularity of the grid. However, the vectorial formalism would imply that
$V(z)$ becomes an infinite-dimensional vector when it is represented
continuously (see next paragraph). For this reason, we prefer to use a
functional formalism and define the volume $V(z)$ as a function (not as a
vector), whether it is modeled as a discrete or continuous field. Discrete
domains are adopted by methods like RELION[29] -Refine3D and RELION[29]
-Class3D[29] — which associate voxels with corresponding intensities in
Fourier space —, and like CryoPoseNet[34]or 3DFlex[32] in image space.

##### Continuous Domain and Implicit Parametrization.

The second class of methods model the volume $V(z)$ as a continuous field,
i.e. as a function on a continuous domain ($\Omega=\mathbb{R}^{3}$ or
$\Omega=[-0.5,0.5]^{3}$). The domain $\Omega$ is infinite, and one cannot
explicitly maintain in memory the values that $V(z)$ takes on $\Omega$. The
solution is then to adopt an explicit parametrization for $V(z)$ using
parameters $\theta\in\Theta\in\mathbb{R}^{p}$. Depending on whether or not
these parameters have a physical meaning (e.g. centroids of pseudo-atoms), the
function $V(z)$ can be encoded:

(i) Using Neural Networks.

Some methods use neural networks to represent $V(z)$ as a (real or complex)
function of a 3D position vector. The parametrization is called “implicit"
because the values of $V(z)$ are not stored in memory; instead, the
practitioner can “query" the neural network by inputing any location
$x\in\mathbb{R}^{3}$ and receiving a value for $V(z)$ at $x$. In this case,
the parameters $\theta$ — i.e. the weights of the neural network — do not have
a physical meaning. Examples of this approach include CryoDRGN[23] and
CryoAI[18], both operating in Fourier space using a reference-free volume
encoding.

(ii) Using Gaussian Mixtures.

Other approaches constrains the volume $V(z)$ by modeling the source of the
electrostatic potential: its individual atoms or pseudo-atoms. Indeed, at a
granular level, the molecular volume can be approximated by a mixture of $N$
Gaussian functions (called scattering form factors[11]) of the form:

$V_{z}(x)=\sum_{j=1}^{N}A_{j}\exp\left(-\frac{||c_{j}-x||^{2}}{2\sigma_{j}^{2}}\right),$
(2)

where $x\in\mathbb{R}^{3}$ represents a 3D position, and
$c_{j}\in\mathbb{R}^{3}$ are the 3D coordinates of the $N$ individual atoms or
pseudo-atoms. The parameters $A_{j}\in\mathbb{R}$ and
$\sigma^{2}_{j}\in\mathbb{R}$ describe how each (pseudo-)atom contributes to
the electrostatic potential. In practice, these approaches always implement
conformational heterogeneity, and do so through a continuous conformation
variable $z\in\mathbb{R}^{L}$ that passes through a neural network to output
$c_{j}$, and possibly $A_{j}$ and $\sigma^{2}_{j}$. This approach also models
$V_{z}$ as a continuous field, as defined by Equation (2), but the parameters
defining each volume ($\theta=\\{c_{j},A_{j},\sigma^{2}_{j})\\}_{j=1,...,J}$)
now have a physical meaning. Among this general class of methods, works differ
in whether $A_{j},\sigma_{j}^{2}$ are assumed to be known, and in the
interpretation given to the variable $c_{j}$. E2GMM[33] use a conformation
variable $z$ that encodes the coefficients $c_{j},A_{j},\sigma_{j}^{2}$ and
defines the $c_{j}$ as coordinates of “coarse grained atoms" (reference-free).
CryoFold[35] assumes $A_{j}=A$ and $\sigma_{j}=\sigma$ known and fixed while
using $c_{j}$ to represent “groups of atoms”. AtomVAE[17] also assumes
$A_{j}=A$ and $\sigma_{j}=\sigma$ known and fixed, models the $c_{j}$ as the
coordinates of the atoms, and uses the conformation variable $z$ to encode
heterogeneous deviations $\Delta c_{j}.$

##### Discretized and Continuous Domains: Pros, Cons and Discussion.

Contrary to the discretized domains, approaches using continuous domains
potentially allow to achieve sharper, enhanced resolutions (within the Nyquist
limit), as any coordinate of $\mathbb{R}^{3}$ can be fed to $V(z)$. Moreover,
within continuous approaches, pseudo-atomic methods effectively add
constraints to $V(z)$ by modeling it as a mixture of Gaussians, and even more
so when assuming a reference conformation $V_{0}$. The increasing availability
of folded protein shapes — traditionally from the Protein Data Bank[36] and
more recently through the advent of AlphaFold[37] — have indeed enabled access
to relatively reliable atom coordinates of reference conformations $V_{0}$,
that can enrich the recovery of the molecular volume. We also note that
reference-based representation such as that proposed in AtomVAE[17] and
CryoFold[35] are more amenable to the inclusion of molecular dynamics
information to the volume reconstruction process.

##  2 Inference

We now turn to the description of the inference methods used in deep
generative modeling for cryo-EM reconstruction. These methods recover the
volume $V$ by finding optimal parameters $\theta$, conformation variables
$z_{i}$ and nuisance variables $(\text{PSF}_{i},t_{i},R_{i})$ of the
generative model (Eq 1). In this section, $\theta$ collectively denotes the
parameters that describe the conformational landscape as a function of $z$,
and the parameters of the function $V_{z}:x\to\Omega$ that associates a
position $x$ to an output intensity. We refer to the conformation variable and
poses jointly as the “hidden variables” and denote them as
$H_{i}=(z_{i},\text{PSF}_{i},R_{i},t_{i})$. For the sake of concision, this
section focuses on general inference methods, and we reserve a description of
their variations for Appendix C.4.

##### Setting Up the Inference Problem: Observed Likelihood vs Full Likelihood

In the context of deep generative modelling for cryoEM, the cornerstone of
inference is simply the observed likelihood $p_{\theta}(x)=p(x|\theta)$
associated with each image $x$. This likelihood is computed from the
generative model in Equation (1) (or its Fourier counterpart -(9)), which we
seek to maximize as a function of $\theta$. However, the generative model
depends on hidden variables $H_{i}=(\text{PSF}_{i},t_{i},R_{i},z_{i})$. In
most cases, the optimization of the full likelihood of each observation
$p(x_{i},h_{i},\theta)$ would be quite simple, if only the $H_{i}$ were
observed. Thus, given $n$ observed images $x_{1},...,x_{n}$, one solution
could be to jointly recover the parameters $\theta$ and hidden variables $H$
(considered here as fixed quantities, as opposed to random variables) that
maximize the log-likelihood
$\ell(X,\theta)=\sum_{i=1}^{n}\log(p_{\theta}(x_{i},h_{i}))$. Mathematically,
this requires solving the following optimization problem:

$\displaystyle\theta^{*},H^{*}=\text{argmax}_{\theta,H}\sum_{i=1}^{n}\log(p_{\theta}(x_{i},h_{i}))$
(3)

It is in fact a classical exercise in statistics to show that in this case, as
the number of estimated variables grows with the number of data points, the
estimate of $\theta$ is no longer guaranteed to converge to the real
underlying value as $n$ goes to infinity
($\lim_{n\to\infty}\mathbb{E}[\theta^{*}]\neq\theta^{\text{true}}$). We thus
have to resort to strategies that treat hidden variables as random variables,
and that fit the parameters $\theta$ based on the “observed likelihood”
$L(X,\theta)$. In this case, the objective becomes:

$\displaystyle\theta^{*}=\text{argmax}_{\theta}L(X,\theta)\quad\text{ where
}\quad L(X,\theta)=\sum_{i=1}^{n}\log
p(x_{i}|\theta)=\sum_{i=1}^{n}\log\int_{h_{i}}p(x_{i},h_{i}|\theta)d\mu(h_{i})$
(4)

where $d\mu(h)=p(h)dh$ is the probability measure associated to the hidden
variables $H$. However, this marginal likelihood requires an integral over all
possible values of $H_{i}$. This quantity is difficult to compute directly, or
in statistical terminology, “intractable”. Consequently, the crux of the
optimization pipeline is to find a way to effectively approximate it.

Figure 2: Maximization of the log-likelihood $\theta\rightarrow L(X,\theta)$
in $\theta$ by maximizations of a series of lower bounds:
$\mathcal{L}(q^{(0)},X,\theta)$, $\mathcal{L}(q^{(1)},X,\theta)$, etc. The
$\theta^{(t)}$s across iterations $t=0,1,2,...$ are represented by colored
dots and correspond to successive maxima of the lower bounds. Left: The lower
bounds are tangent to $\theta\rightarrow L(X,\theta)$, which is realized when
$q$ is the posterior of the hidden variables. Right: The lower bounds are not
tangent to $\theta\rightarrow L(X,\theta)$, but show a “gap" that corresponds
to the KL divergence between $q$ and the posterior of the hidden variables,
see Equation (5).

###  2.1 Unifying Inference Methods

Since the observed log-likelihood $L(X,\theta)$ in the objectives of Equations
(4)-(12) is intractable, optimization is usually performed by targeting a
proxy for $L(X,\theta)$, called the Evidence Lower Bound (ELBO). For the sake
of clarity and concision, we highlight here the common statistical thread of
cryo-EM reconstruction methods leveraging deep generative modeling, that all
use an ELBO-based optimization and refer the reader to Appendix Appendix C:
Inference and ELBO Computations for further discussion on their variations.

##### Evidence Lower Bound (ELBO)

The trick behind the Evidence Lower Bound (ELBO)consists in proposing a series
of distributions $q^{(0)},...,q^{(t)}$ for the hidden variables $H$, and
maximizing a series of “easily” computable lower-bounds
$\mathcal{L}(q^{(0)},X,\theta),...,\mathcal{L}(q^{(t)},X,\theta)$ for
$L(X,\theta)$ in an iterative fashion — see Figure 2. By iteratively
maximizing these lower bounds with respect to $\theta$, the true likelihood
$L(X,\theta)$ also increases. The hope is that the value of $\theta$ obtained
through their maximization will be close to the value realizing the maximum of
$L(X,\theta)$, if the lower bounds are tight enough — i.e. for small “gaps" in
Figure 2.

The lower bounds $\mathcal{L}(q,X,\theta)$ are found by showing that, for any
probability distribution $q_{i}$ on the variables $h_{i}$, the observed log-
likelihood can be written as the sum of two terms (derivations provided in
Appendix Appendix C: Inference and ELBO Computations):

$\displaystyle
L(X,\theta)=\mathcal{L}(q,X,\theta)+\sum_{i=1}^{n}\text{KL}(q_{i}(h_{i})\parallel
p_{\theta}(h_{i}|x_{i}))=\sum_{i=1}^{n}\Big{[}\mathcal{L}_{i}(q_{i},x_{i},\theta)+\text{KL}(q_{i}(h_{i})\parallel
p_{\theta}(h_{i}|x_{i}))\Big{]}$ (5)

where $KL$ is the Kullback-Leibler divergence (KL)defined as $KL(q\parallel
p)=\int q(x)\log\frac{q(x)}{p(x)}dx$, and the terms $\mathcal{L}_{i}$ write:

$\mathcal{L}_{i}(q_{i},x_{i},\theta)=\int_{h_{i}}q_{i}(h_{i})\log
p_{\theta}(x_{i}|h_{i})dh_{i}-\text{KL}\left(q_{i}(h_{i})\parallel
p_{\theta}(h_{i})\right).$ (6)

The divergences $\text{KL}(q(h_{i})\parallel p_{\theta}(h_{i}|x_{i}))$ in
Equation (5) are always non-negative. Thus, for any joint distribution
$q=\\{q_{i}\\}_{i=1\cdot n}$, the function $\mathcal{L}(q,X,\theta)$ provides
a valid lower-bound to $L(X,\theta)$ (see Figure 2), called the Evidence Lower
Bound (ELBO):

$\forall q,\forall\theta,\quad\mathcal{L}(q,X,\theta)\leq L(X,\theta).$

The lower-bounds $\mathcal{L}(q^{(t)},X,\theta)$ are proxies for
$L(X,\theta)$, that - in contrast to $L(X,\theta)$ \- can be computed and
maximized in $\theta$.

##### Inference Methods Based on an ELBO.

While the ELBO holds for any $q$, some choices are more judicious than others.
In fact, the goal is to select an optimal $q$, such that the gap between
$\mathcal{L}(q,X,\theta)\leq L(X,\theta)$ is small: this will insure that the
maximization of $\mathcal{L}(q,X,\theta)$ with respect to $\theta$ yields
estimates $\theta^{*}$ that are also appropriate (and close to the true
optimum ${\theta}^{\text{true}}$) for maximizing $L(X,\theta)$ — see Figure 2
(right). Inference methods in cryo-EM subsequently differ in the choices of
the distributions $q_{i}^{(t)}$ for each $i$ and at each iteration $t$,
thereby yielding different lower bounds $\mathcal{L}(q,X,\theta)$:

(i) Using the posteriors given current parameters (EM algorithm):

Computing the posteriors $p_{\theta}(h_{i}|x_{i})$ using the current estimated
value $\theta^{(t)}$ of $\theta$ allows choosing
$q_{i}^{(t)}(h_{i})=p_{\theta^{(t)}}(h_{i}|x_{i})$ for each $i$ at iteration
$t$ — see Figure 5. The inequality:

$\mathcal{L}_{i}(p_{\theta^{(t)}}(h_{i}|x_{i}),X,\theta)\leq L_{i}(X,\theta),$

becomes an equality for $\theta^{(t)}=\theta$. This makes the lower-bound
$\mathcal{L}(q,X,\theta)$ tangent to $L(X,\theta)$ at $\theta=\theta^{(t)}$:
progressively
maximizing $\mathcal{L}(q,X,\theta)$ with respect to $\theta$ will induce
convergence to a local maximum of $L(X,\theta)$ in $\theta$, as seen in Figure
2 (left). This is the strategy adopted by Expectation-Maximization
(EM)algorithm (more details in Appendix Appendix C: Inference and ELBO
Computations). The EM is an iterative algorithm which consists of two steps.
In the first step (called the expectation step), given current parameters
values $\theta^{(t)}$, we compute the posterior
$q_{i}^{(t)}(h_{i})=p_{\theta^{(t)}}(h_{i}|x_{i})$ to plug into our ELBO. In
the second (the maximization step), $\theta^{(t+1)}$ is taken to be the value
of $\theta$ that maximizes the ELBO. This sequence of two steps is usually
repeated until convergence. Cryo-EM methods adopting this approach are given
in the first column of Table 4. As explained in Appendix Appendix C: Inference
and ELBO Computations, while the EM algorithm does not have any convergence
guarantees, it nonetheless guarantees to increase the likehood at each step.

(ii) Approximating the posteriors given current parameters (Variational EM
algorithm):

In certain cases, the choice of $q_{i}^{(t)}(h_{i})$ as the posterior
$p_{\theta^{(t)}}(h_{i}|x_{i})$ is neither computationally attractive nor
feasible. In this case, we might prefer approximating each posterior by
finding its “best approximation” $q_{i}^{*}$ within a family of functions
called variational family $\mathcal{Q}$. Cryo-EM reconstruction methods
consider two choices that include approximating the posteriors by (i) their
“mode", i.e. the value $\hat{h}_{i}$ of $h_{i}$ that maximizes
$p_{\hat{\theta}^{(t)}}(x,h)$. In this case, each $q_{i}$ effectively becomes
a Dirac distribution at $\hat{h}_{i}$; or (ii) or a general distribution
$q_{i}$ within a family $\mathcal{Q}$: $q_{i}$ is for example a Gaussian
distribution – see Figure 5. Cryo-EM methods adopting this approach are given
in the last two columns of Table 4, within the “non-amortized" subcolumns.

##### Exact or Approximate Posteriors: Pros, Cons, Discussion.

The EM algorithm, that uses exact posteriors, holds several advantages: it is
simple and stable, since all updates can only improve the observed log-
likelihood. However, it is also potentially slow: the rate of convergence is
known to be linear with rate proportional to the fraction of information about
$\theta$ in $L(\theta,X)$ [41]. Variational EM algorithms can be faster; yet
they potentially loose in accuracy as their ELBOs do not provide tight lower-
bounds to the log-likehood $L(X,\theta)$ (Figure 2, right). As a result, we do
not have any guarantee that they converge to an (even local) maximum of
$L(X,\theta)$.

###  2.2 Introducing Amortized Inference

Figure 3: Comparison of generative reconstruction methods. Volume represented
in image space (resp. Fourier space) are outlined in blue (resp. orange).
CryoDRGN[24] is similar to cryoDRGN2[25], except rotations and translations
are given by an upstream homogeneous reconstruction. CryoVAEGAN[42] does not
explicitly store a representation of the volume and therefore does not appear
in the figure. VG = Voxel Grid ; NN = Neural Network ; GM = Gaussian Mixture.

While potentially more computationally attractive than the original EM
algorithm, Variational EM requires solving $n$ optimization problems to find
an approximate posterior $q_{i}$ for each image $i$ in $1,...,n$. This is
computationally expensive, as the number of $q_{i}$ to estimate increases as
the number of images $n$ increases. Consequently, recent methods have resorted
to using an additional approximation called Amortized Inference (AI), which
collapses the $n$ optimizations problems into one. Instead of finding the best
corresponding $q^{*}_{i}$ for each $i$ , AI optimizes the parameters $\xi$ of
a function $\text{Enc}_{\xi}$ that predicts the parameters of the distribution
$q_{i}^{*}(h_{i})$ when given $x_{i}$ as input, i.e.:
$\text{Enc}_{\xi}(x_{i})\simeq(\mathbb{E}[h_{i}],\text{Var}[h_{i}]))$, where,
in this example, the variational family $\mathcal{Q}$ is chosen to be the set
of Gaussian distributions. In other words,instead of solving $n$ separate
problems, Amortized Inference predicts the parameters of the posterior of
image $i$ using the observed image as input. The function $\text{Enc}_{\xi}$
is traditionally called an encoder. More details — including a description of
updates performed in AI — can be found in Appendix Appendix C: Inference and
ELBO Computations. Cryo-EM methods using this approach are shown in the
“amortized" subcolumns of Table 4.

##### Implementation of Amortized Inference with Variational Autoencoders

In cryo-EM reconstruction, amortized inference is deployed in the context of
variational autoencoders, denoted VAEs. VAEs are deep architectures that model
the parameters of the variational family $\text{Enc}_{\xi}$ described above as
a neural network with weights $\xi$ — therefore leveraging the expressivity of
this class of functions to get an optimal (amortized) variational
approximation. The entire VAE pipeline thus consists of two steps: an encoder,
which is simply a neural network with weights $\xi$ corresponding to the
function $\text{Enc}_{\xi}$ described above, and a decoder, which allows to
create “mock samples” that will then be compared with the observed ones based
the generative model with parameter $\theta$ chosen in Section 1. Here, the
decoder is almost entirely dictated by the process described in Equation1 and
goes beyond the conventional fully connected networks or convolutional neural
networks used in image processing.

The learnable parameters $\xi$ and $\theta$ of the encoder and the decoder are
fitted through stochastic gradient descent via backpropagation of the ELBO
through the neural network. Compared to traditional cryo-EM reconstruction
methods leveraging the EM algorithm, variational Autoencoders can be
interpreted as extending the E-step (encoder) and the M-step (decoder) of the
EM algorithm. The VAE architectures of cryo-EM reconstruction methods using
amortized inference are given in Figure 3. In this review, we have also
included for comparison purposes a non-variational version of this procedure
(the 3DFlex[32] auto-encoder). Here, the authors consider the hidden variables
as non-random variables, but add a fix amount of gaussian noise to regularise
the embeddings. While the final loss is therefore adapted, this is essentially
a VAE where the variance is fixed, while only the mean is learned.

##### Amortized Inference: Pros, Cons, Discussion

Amortized inference is faster than its non-amortized counterparts, but adds an
additional error (called the amortization error). We observe that several
methods use amortized inference, but often to estimate one hidden variable:
e.g. only the rotation $R$ or only the conformation variable $z$. Table 4
classifies the reconstruction methods by the type of inference chosen for each
variable within $h_{i}=(R_{i},t_{i},\text{PSF}_{i},z_{i})$ and indeed, we note
that this choice does not have to be consistent across all hidden variables.
Many methods “mix and match” inference techniques, using for example a
variational EM for the hidden rotation variable $R_{i}$ and a VAE for the
conformation variable $z_{i}$. Moreover, it becomes apparent from Table 4 that
(variational) autoencoders are the most common type of approaches implemented
for cryo-EM reconstruction.

While deep generative methods for cryo-EM volume reconstruction can be unified
with the framework described above (as well as with traditional Expectation
Maximization approaches), we observe that each of them has its own
specificities or “implementation tricks". They differ, for example, in their
choice of variational family, or loss function that adapts the ELBO to
facilitate convergence of the optimization procedure, see Appendix Appendix C:
Inference and ELBO Computations. These testify to the difficulties encountered
in training these algorithms in the context of cryo-EM images with low signal-
noise ratios.

##  3 Discussion

Given the wide number of options to reconstruct molecular volumes from cryo-EM
images, it is natural to ask: which reconstruction method is in fact the most
promising? In this last section, we focus on the need for establishing a set
of metrics and benchmark tasks that can be used to quantitatively compare the
performance of these methods. Starting with a review of the tools currently
available, the first take-away of this section is the urgent need for new
metrics and benchmarks. The evaluation of these methods’ performance is
currently difficult and inherently limited. We nonetheless highlight, as a
second take-away, promising features in current developments, which, in our
opinion, these developments pave the way for future improvements in cryo-EM
reconstruction.

###  3.1 Assessing Reconstruction Performances: Need for New Metrics

Performance metrics can be categorized in two classes: (a) those that assess a
method’s ability to provide good spatial resolution (i.e. distinguishing
different atoms), and, in the case of heterogeneous methods, (b) those that
assess a method’s ability to provide good conformation resolution (i.e.
distinguishing different conformations).

####  3.1.1 Assessing Spatial Resolution

##### Resolution of discretized reconstructions (3D maps).

When the reconstructed volume is parametrized as an explicit 3D map, the most
widespread measure used to evaluate its spatial resolution is the Fourier
Shell Correlation (FSC)[50]. As described by Singer et al.[14], this quantity
measures the correlation over a 3D shell between two reconstructed volumes:

$\displaystyle FSC_{k}(U,V)=\frac{\sum_{s\in
S_{k}}U_{s}V_{s}^{*}}{\sqrt{\sum_{s\in S_{k}}|U_{s}|^{2}\sum_{s\in
S_{k}}|V_{s}|^{2}}}.$ (7)

Here, $S_{k}$ is the set of Fourier voxels in a spherical shell at distance
$k$ from the origin, and $U$ and $V$ are the Fourier transforms of the 3D
volumes that we compare. Typically, $U$ and $V$ correspond to two independent
reconstructions on separate halves of the dataset, in which case, the
criterion for a method to be deemed to perform well is for the two
reconstructed volumes to be similar. The method’s resolution is then defined
as the highest resolution for which $U$ and $V$ agree “enough”. This is
precisely what the FSC (Equation 7) captures: the FSC is close to 1 when the
two maps are close. This is usually the case for small $k$, as low-frequency
signal is strong, but the FSC generally decays to zero as the signal-to-noise
decreases. The result is often plotted as a curve, with axis $x=k$. The
resolution of the reconstruction corresponds to the maximum value of $k$ such
that $FSC_{k}\geq 0.143$ — a criterion chosen to match resolution criteria
used in X-ray crystallography[49]. For synthetic datasets where a ground-truth
volume is available, the FSC is measured between the reconstruction and the
ground-truth; in which case the resolution criterion correspond to the maximum
value of $k$ such that $FSC_{k}\geq 0.5$.

##### Resolution of continuous reconstructions.

Methods that represent the volume as a continuous field are relatively new,
and it might be worth reassessing appropriate metrics for evaluating spatial
accuracy in this case.

(i) Implicit Parametrizations

While interpolation between image pixels and map voxels is necessary in the
discrete case, both for projection and for backprojection, implicit
representations of the volume (e.g., through an neural network) enable
sampling without interpolation during training. It would be interesting to
investigate whether this provides a benefit in terms of reconstruction
quality. We do not expect implicit representations to suddenly _unlock_
information, since the information content is determined by the discrete
nature of the images and their pixel size, but they might provide new
actionable ways to implement prior informations about the volume, such as
smoothness and stereochemistry, that would result in reconstructions of higher
quality.

(ii) Atomic Parametrizations

Parametrizations of the volume with atomic models represents an opportunity to
revisit the notion of spatial resolution. The traditional measure of
similarity between two atomic models that only differ in the cartesian
coordinates $U$ and $V$ of their $N$ constituting atoms (using a consistent
orientation of the molecule for $U$ and $V$) is the Root Mean Square Deviation
(RMSD). This quantity is defined as:

$\displaystyle RMSD(U,V)=\sqrt{\frac{1}{N}\sum_{i=1}^{N}|U_{i}-V_{i}|^{2}},$
(8)

However, it was soon recognized that this metric had a very narrow range for
interpretability [57]: it is a global measure of similarity, which is not
suited to capture the local — but meaningful— changes in protein structure. To
complicate the matter, measuring the RMSD between two atomic models assumes
that they are both in the same reference frame, which might not always be
defined. To improve the sensitivity of the metric, atomic models are routinely
reduced to features with desired properties, e.g. vectors of internal
coordinates independent of the reference frames. For example, reducing the
atomic model to its backbone dihedral angles or to a list of atomic contacts
has been shown to yield better clustering of conformations [58]. The
development of new established metrics to evaluate these models is thus an
important avenue of development.

####  3.1.2 Assessing Conformation Resolution: An Ill-Defined Problem

While not flawless, the evaluation of spatial resolution is a relatively well-
characterized task. By contrast, evaluating conformation heterogeneity is a
more ambiguous problem. To evaluate the quality of the reconstructions
allowing continuous heterogeneity, methods such as 3DFlex[32] or CryoDRGN[23]
perform a post-hoc analysis of the recovered latent space, showing the
flexible deformation that are induced by sweeping through the space of
possible $z$s and visually inspecting the corresponding deformations. However,
proper objective and quantitative measures of conformation heterogeneity
remain to be established: there currently exists no standardized measure or
gold-standard task to evaluate how well a method is able to capture it.

We could design a new metric, inspired by the high-level idea of the FSC.
Using two halves of the dataset to infer two independent continuous
distributions of conformations, we evaluate whether the distributions agree
using a metric such as the Wasserstein distance - modulo change of coordinate
system for the conformation variable $z$. In the case where a ground-truth
conformation is available for each image (e.g. in simulations), the inferred
distribution could be compared to the true distribution. In the case of
methods able to generate one 3D volume for each image in the dataset, one
could consider a hierarchical clustering approach where depth in the hierarchy
tree corresponds to the conformational resolution. In more concrete terms, for
all resolution $k$, the FSC between each volume pairs would be measured and
the resulting distance matrix used for clustering. Data points that fall
within the same clusters would be indistinguishable at $k$ while images that
would fall in different clusters would correspond to conformations that differ
by at least $k$. The development of such metrics is key to make sustainable
advances in next-generation cryo-EM reconstructions.

###  3.2 Quantitative Comparison of Performances: Lack of Common Benchmarks

Beyond the need for new performance metrics that are better adapted to new
advancements in the field, it is most certainly the lack of established
benchmarks that, to this day, make reconstruction methods very hard to
compare. Such benchmarks are dramatically needed, as we cannot rely on
statistical theory since the convergence properties of estimations relying on
(amortized) variational inference are not completely characterized. In fact,
the quantitative assessment of the methods’ relative performance has yet to
overcome three main hurdles:

(i) Lack of benchmark datasets:

Current methods are developed and tested on a wide range of synthetic and
experimental datasets that differ in the nature of the biomolecule being
imaged, the dataset size, image size and associated resolution — with very
little overlap across methods - see Table 2. There is unfortunately no
“MNIST”[54] or “Imagenet”[55] for cryo-EM. Most methods resort to evaluating
their performance on synthetic data, yet no cryo-EM simulator acts as a
standard to generate simulated images in a unified way. Synthetic datasets
vary in the realism of the image formation model used for simulation, e.g. in
the noise model, the signal-to-noise ratio or the distribution of nuisance
variables (e.g. poses). Subsequent experiments are typically performed on real
“in house” data — but there too, the important diversity within the
characteristics of these evaluation datasets therefore makes the comparison of
these methods a strenuous task.

(ii) Lack of benchmarking procedures:

Reconstruction methods vary in the complexity of the task that they set out to
accomplish, assuming more or less nuisance variables (such as poses or PSF) to
be known – see Table 3 and Figure 3. This makes it difficult to compare
methods on a fair ground. We need to establish modular benchmarking procedure
that would enforce a fair comparison of reconstruction performances, eg,
testing the recovery of the pose, volume or conformations, with other nuisance
variables being known and fixed.

(iii) Lack of benchmark codebase and infrastructure:

Finally, reconstructions methods are not necessarily publicly accessible, are
implemented across different programming languages, and/or are tested on
different software or hardware. Creating a codebase that re-implements these
methods for a proper evaluation using a unified infrastructure would
unfortunately represent a gigantic implementation effort. Currently, this lack
of codebase poses a significant hurdle in the accessibility and comparison of
the methods: it is currently impossible to disentangle the effect of their
proposed statistical learning problem, their programming language,
implementation tricks, or software infrastructure.

###  3.3 Qualitative Comparison of Performances

Despite the hurdles associated with performing quantitative comparisons, we
propose a qualitative evaluation of the different methods based on both
published results and our personal experience. This allows us to highlight
promising directions — to the least, in the authors’ opinion — for further
developments.

####  3.3.1 Accuracy

Despite encouraging accuracy results, some methods seem to face considerable
challenges when applied to real cryo-EM images, as they have not been properly
vetted and stress-tested in experimental conditions [39, 35, 17, 34] – see
Table 2. We consider the lack of results on experimental data as a proxy for a
limited applicability in the context of real signal-noise ratios regimes. In
order to be adopted by cryo-EM practitionners, these methods will need to
overcome the signal-noise regime and demonstrate the accuracy reported in the
papers on a larger set of (benchmark) datasets.

Despite the difficulty of the task and lack of standardized benchmarks, recent
developments in deep generative modeling have shown impressive promise in
overcoming the current computational and accuracy bottlenecks in all three
following directions:

(i) Poses:

Poses are important nuisance variables that have the potential to damage the
reconstruction, if incorrectly predicted. Accuracy of the predicted rotation
is measured on synthetic datasets with a mean/median square error (MSE)
against the corresponding ground-truth. Historically, preference was given to
methods that did not use amortized inference for the rotation estimation (e.g.
CryoSPARC[16]or CryoDRGN[23] ), as they outperformed their amortized
counterparts predicting rotations with an encoder (e.g. CryoPoseNet[34]and
CryoAI[18]): AtomVAE[17] was for instance one of the only methods using
amortized inference for the recovery of the poses, and reported difficulties
in the joint training of poses and conformation — highlighting the difficulty
of accurate amortized inference in this setting. However, the accuracy gap
between methods is closing: CryoAI[18]now showcases an rotation accuracy at
the same order of magnitude compared to CryoSPARC[16]and CryoDRGN[23] on a
real dataset. This was facilitated by the theoretical insights drawn from
CryoAI[18], who show that amortized inference techniques tend to get stuck in
local minima where the predicted molecule contains unwanted planar symmetries
due to their projections on a 2D surface. The solution that they propose to
alleviate this problem is to use the symmetrized loss:

$\ell_{\text{sym}}=\sum_{i}\min\\{||{X}_{i}-\text{PSF}_{i}*(t_{i}\circ\Pi_{2D}\circ
R_{i})(V^{(i)})||^{2},||R_{\pi}({X}_{i})-\text{PSF}_{i}*(t_{i}\circ\Pi_{2D}\circ
R_{i})({R}_{\pi}(V^{(i)}))||^{2}\\}.$

where $R_{\pi}$ is a rotation with angle $\pi$. This has recently opened the
door to significant gains in accuracy in the prediction of the poses, allowing
for the first time pose estimation to be done through amortized inference. We
anticipate that it is through such developments and theoretical insights that
reconstruction algorithms will be able to fully leverage amortized inference
for rotation prediction, hereby providing significant speed-ups.

(ii) Volumes:

Methods based on (pseudo-)atomic volume parametrizations - E2GMM[33] ,
CryoFold[35] and AtomVAE[17] \- do not compare themselves to their
counterparts, probably due to the fact that they were published concurrently
in 2021 and/or do not use the same definition of “pseudo-atoms" that are
respectively: means of 3D Gaussian distributions, residues or actual atoms. As
a consequence, we do not comment on them. For methods generating 3D maps, it
has been reported by Punjani & Fleet [32] that amortized inference can
translate into resolution loss. Yet, recent methods such as CryoDRGN[23] and
CryoAI[18]publish examples of reconstructed volumes as 3D maps whose
resolution is visually comparable to the ones obtained by CryoSPARC[16], on
down-sampled imaged. Our opinion is that amortized methods can reach near-
atomic resolution reconstructions, but this has not been demonstrated yet. If
so, we expect them to replace traditional reconstruction pipelines in the long
run, since they offer the promise to be significantly faster and to tackle
much larger datasets - a desired feature to enable sufficient sampling of the
conformational landscape.

(iii) Conformations:

Methods based on pseudo-atomic volumes parametrizations do not provide
examples of conformation trajectories that allow us to compare them. For
methods generating 3D maps, CryoDRGN[23] and 3DFlex[32] seem to be some of the
most promising approaches, as they seem to allow greater resolutions in the
recovered trajectories, based on our personal visual assessment of the
examples of conformation trajectories shown in the corresponding papers. This
remains to be confirmed by a quantitative assessment over a larger number of
conformation trajectories.

####  3.3.2 Reproducibility

Adoption of these methods by practitioners will require their reproducibility,
or robustness to different initializations, implementation tricks or choice of
hyperparameters:

(i) Initialization:

The non-convex nature of the problem puts it at very high-risk of being non
robust and sensitive to initialization — this is a phenomenon sometimes
referred to as “Einstein from noise”[47] (also described in Singer et al.
[14]). Luckily, most of the current methods show encouraging signs of
robustness to perturbations. In our experience, CryoDRGN[23] seems consistent
for different random initializations of the neural model when fixing the
poses: the conformation space does not seem to be vastly affected. This
robustness can however be challenged by extremely low signal and/or
heterogeneous datasets, in which case certain conformations can go missing.

(ii) Tricks:

The inference methods presented often make use of additional implementation
tricks (e.g. warm-starting with a known conformation), and specific
regularization schemes: e.g. AtomVAE[17] suggests starting with an initial
phase of pose-only training, which, once realized, ensures that the further
joint learning of poses and volume is successful. Both AtomVAE[17] and
CryoFold[35] regularize the recovered structure by penalizing bond lengths,
but the impact of the regularization yet remains to be properly characterized,
and in particular, its potential to frustrate the optimization landscape. The
importance of tricks and regularizations, and combinations thereof, is still
ill-understood and would require an in-depth analysis, as it hints towards a
difficult optimization landscape for this method, and its sensitivity to
initial conditions.

(iii) Hyperparameters:

Choosing hyperparameters such as the dimension of the latent space in
algorithms such as CryoDRGN[23] induces more or less regularization: too small
and it regularizes the model too much; too large leads to underfitting of the
3D model. CryoDRGN[23] usually sets it to $d=8$, but, given how heterogeneity
arises, this is necessarily molecule dependent. The field will need — to the
least— rule-of-thumb guidelines on how to choose these hyperparameters if
these methods are to be adopted by practitioners.

The robustness of these new methods needs to be confirmed on a wider set of
datasets, including datasets with high levels of noise. The fact that they
rely on user-defined implementation tricks and hyperparameters might not be an
obstacle to their adoption, as conventional methods such as RELION[29] or
CryoSPARC[16]also do.

####  3.3.3 Efficiency

Our last axis of comparison is computational efficiency: both in time and
memory requirements. First, if we take the size of the datasets used in
experiments as a proxy for efficiency, then 3DFlex[32] , CryoAI[18], E2GMM[33]
, and CryoDRGN[23] seem to be able to process remarkable amounts of
information. Additionally, we offer our own practical experience by-way of
rule of thumb. With datasets typically of more than 100GB, training times can
take up to 10 hours (including the required pre-processing steps) for methods
like RELION[29] — that is, for a run that has little hyperparameter tuning
compared to alternative deep learning methods. Newer methods like CryoDRGN[23]
hold great promise in terms of reconstruction: however, such sophisticated
methods can further benefit from gains in efficiency, both from the
computational side and in terms of memory requirements. Efficient updates of a
model’s parameters can thus be seen as a current computational bottleneck and
offers an interesting avenue for future research.

##  4 Conclusion

This review provides a critical comparison of recent cryo-EM reconstruction
methods that are based on deep generative modeling, focusing on explaining
their relative advantages or drawbacks. We have unified, compared and
contrasted existing methods through their parametrization of the volume, as
well as through the optimization procedure chosen to recover this volume and
associated hidden variables. While the use of amortized inference is crucial
to make inference tractable in this high-data, high-dimensional setting, there
seems to be much room for improvement and research on methods allowing both
faster and better inference. On a practical side, we note that recent methods
suffer from a lack of benchmarks which severely impedes their comparison and
development. From our practical experience, beyond a necessity for benchmark
datasets, we also highlight a severe need for the development of a diagnostic
toolbox tailored to the analysis of cryo-EM data. Current methods rely on a
set of choices and hyperparameters that raise a number of questions for the
practitioner: have I chosen my hyperparameters adequately? Is this choice
going to impact the accuracy of the recovery? Is there any physical or
biological meaning or interpretation in the distance between the latent space
of conformation variables? How does error on pose or PSF affect the rest of
the volume recovery process? There is therefore a pressing need for more in-
depth and systematic quantitative comparison of these methods.

## Acronyms

AI

    Amortized Inference
CTF

    Contrast Transfer Function
ELBO

    Evidence Lower Bound
EM

    Expectation-Maximization
FSC

    Fourier Shell Correlation
GAN

    Generative Adversarial Network
KL

    Kullback-Leibler divergence
PSF

    Point Spread Function
RMSD

    Root Mean Square Deviation
VI

    Variational Inference

## Glossary

collective variables

    The conformational space of a molecule can be defined by features that are functions of the coordinates of the molecule’s atoms, often referred to as collective variables since most interesting features tend to capture global collective motions of the atoms
conformational landscape

    Molecules exhibit shape variability as they thermally diffuse on their high-dimensional energy potential surface, also known as their conformation landscape. Indeed, each atom of the molecule is associated with 6 degrees of freedom, three for their position in space and three for their momentum. The joint distribution over all these degrees of freedom is directly related to the pairwise interaction between atoms which add up to the potential energy of the molecule. The likelihood of a conformation is dictated by the potential energy value associated with this conformation and thus by the relative position of atoms. One could think of the potential energy function as a surface in high-dimension with basins corresponding to stable conformations separated by high energy intermediates. The overall energy scale and relative depths of the basins is dictated by thermodynamics observables such as the temperature; Low temperature will exacerbate differences in energy and differentially increase the likelihood of conformations found in basins relative to higher energy intermediate, for example
heterogeneity

    Heterogeneity refers to structural variability between the imaged particles, whether from compositional variability in the sample or conformational dynamics associated with the molecule. Heterogeneous reconstruction algorithm aim to model the associated with the target molecule
cryo-EM

    Cryogenic Electron Microscopy
generative model

    A generative model is a statistical model of the joint probability distribution for the data and the underlying quantity of interest. It is a process that describes how data is generated, using probability distributions. By sampling from this model, we create new data points
particles

    In cryo-EM, each copy of an imaged molecule is referred to as a _particle_
poses

    The unknown 3D orientation and 2D translation of the particle within a given projection image
single particle images

    Single Particle Imaging refers to imaging modalities where data, _albeit_ corrupted, is collected about individual particles, as opposed to ensemble averaging imaging modalities such as X-ray crystallography or NMR

## References

  * [1] Slavica Jonić  Computational methods for analyzing conformational variability of macromolecular complexes from cryo-electron microscopy images.  _Current Opinion in Structural Biology_ DOI: 10.1016/j.sbi.2016.12.011 (2017).
  * [2] Noe, F. and Clementi C.  ollective variables for the study of long-time kinetics from molecular trajectories: theory and methods  _Current Opinion in Structural Biology_ DOI: 10.1016/j.sbi.2017.02.006 (2017).
  * [3] Katsevich, E. and Katsevich, A. and Singer, A.  Covariance Matrix Estimation for the Cryo-EM Heterogeneity Problem  _SIAM Journal on Imaging Sciences_ DOI: 10.1137/130935434 (2015).
  * [4] Haselbach et al  Structure and Conformational Dynamics of the Human Spliceosomal Bact Complex  _Cell_ DOI: 10.1016/j.cell.2018.01.010 (2018).
  * [5] Bock, Lars V. and Grubmüller, Helmut  Effects of cryo-EM cooling on structural ensembles  _BioRxiv_ DOI: 10.1101/2021.10.08.463658 (2021).
  * [6] Noé, Frank and Tkatchenko, Alexandre and Müller, Klaus-Robert and Clementi, Cecilia  Machine Learning for Molecular Simulation  _Annual Review of Physical Chemistry_ DOI: 10.1146/annurev-physchem-042018-052331 (2020).
  * [7] Ourmazd, A.  Cryo-em, XFELs and the structure conundrum in structural biology.  _Nature Methods_ DOI: 10.1038/s41592-019-0587-4 (2019).
  * [8] Nakane, T. _et al._ Single-particle cryo-em at atomic resolution.  _Nature_ 587, 152–156 (2020).
  * [9] EMDB.  Statistics.  https://www.ebi.ac.uk/emdb/statistics/emdb_resolution_year (2022).
  * [10] Dill, Ken A & et al.  Molecular driving forces: statistical thermodynamics in biology, chemistry, physics, and nanoscience.  In _Garland Science_ (2010).
  * [11] Kohl, H. & Reimer, L.  _Transmission Electron Microscopy_.  Springer Series in Optical Sciences.
  * [12] Vulović, M. _et al._ Image formation modeling in cryo-electron microscopy.  _Journal of Structural Biology_ 183, 19–32, DOI: https://doi.org/10.1016/j.jsb.2013.05.008 (2013).
  * [13] Scheres, S. H.  RELION: Implementation of a Bayesian approach to cryo-EM structure determination.  _Journal of Structural Biology_ DOI: 10.1016/j.jsb.2012.09.006 (2012).
  * [14] Singer, A. & Sigworth, F. J.  Computational methods for single-particle electron cryomicroscopy.  _Annual Review of Biomedical Data Science_ 3, 163–190 (2020).
  * [15] Kimanius, D., Forsberg, B. O., Scheres, S. H. & Lindahl, E.  Accelerated cryo-EM structure determination with parallelisation using GPUS in RELION-2.  _eLife_ DOI: 10.7554/eLife.18722 (2016).
  * [16] Punjani, A., Rubinstein, J. L., Fleet, D. J. & Brubaker, M. A.  CryoSPARC: Algorithms for rapid unsupervised cryo-EM structure determination.  _Nature Methods_ DOI: 10.1038/nmeth.4169 (2017).
  * [17] Rosenbaum, D. _et al._ Inferring a continuous distribution of atom coordinates from cryo-em images using vaes.  _CoRR_ abs/2106.14108 (2021).  2106.14108.
  * [18] Veym, A. _et al._ CryoAI: Amortized Inference of Poses for Ab Initio Reconstruction of 3D Molecular Volumes from Real Cryo-EM Images.  _arXiv_ (2022).  2203.08138.
  * [19] Boyd, S., Boyd, S. P. & Vandenberghe, L.  _Convex optimization_ (Cambridge university press, 2004).
  * [20] Bendory, T., Bartesaghi, A. & Singer, A.  Single-particle cryo-electron microscopy: Mathematical theory, computational challenges, and opportunities.  _IEEE signal processing magazine_ 37, 58–76 (2020).
  * [21] Sitzmann, V. _et al._ Implicit Neural Representations with Periodic Activation Functions.  _arXiv preprint arxiv.2006.09661_ (2020).
  * [22] Si, D. _et al._ Artificial intelligence advances for de novo molecular structure modeling in cryo-em.  _arXiv preprint arXiv:2102.06125_ (2021).
  * [23] Zhong, E. D., Bepler, T., Davis, J. H. & Berger, B.  Reconstructing continuously heterogeneous structures from single particle cryo-em with deep generative models.  _arXiv preprint arXiv:1909.05215_ (2019).
  * [24] Zhong, E. D., Bepler, T., Berger, B.& Davis, J. H.  CryoDRGN: reconstruction of heterogeneous cryo-EM structures using neural networks.  _Nature methods_ (2021).
  * [25] Zhong, E. D., Lerer, Adam, Davis, J. H. & Berger, B.  CryoDRGN2: Ab initio neural reconstruction of 3D protein structures from real cryo-EM images.  _Proceedings of the IEEE/CVF International Conference on Computer Vision_ (2021).
  * [26] Ede, J. M.  Deep learning in electron microscopy.  _Machine Learning: Science and Technology_ 2, 011004 (2021).
  * [27] Wu, J.-G. _et al._ Machine learning for structure determination in single-particle cryo-electron microscopy: A systematic review.  _IEEE Transactions on Neural Networks and Learning Systems_ 1–21, DOI: 10.1109/TNNLS.2021.3131325 (2021).
  * [28] Sharon, N., Kileel, J., Khoo, Y., Landa, B. & Singer, A.  Method of moments for 3d single particle ab initio modeling with non-uniform distribution of viewing angles.  _Inverse Problems_ 36, 044003 (2020).
  * [29] Scheres, S. H.  A bayesian view on cryo-EM structure determination.  _Journal of Molecular Biology_ DOI: 10.1016/j.jmb.2011.11.010 (2012).
  * [30] Scheres, S. H.  Chapter eleven - classification of structural heterogeneity by maximum-likelihood methods.  In Jensen, G. J. (ed.) _Cryo-EM, Part B: 3-D Reconstruction_ , vol. 482 of _Methods in Enzymology_ , 295–320, DOI: https://doi.org/10.1016/S0076-6879(10)82012-9 (Academic Press, 2010).
  * [31] Punjani, A. & Fleet, D. J.  3d variability analysis: Resolving continuous flexibility and discrete heterogeneity from single particle cryo-em.  _Journal of Structural Biology_ 213, 107702, DOI: https://doi.org/10.1016/j.jsb.2021.107702 (2021).
  * [32] Punjani, A. & Fleet, D. J.  3d flexible refinement: Structure and motion of flexible proteins from cryo-em.  _bioRxiv_ DOI: 10.1101/2021.04.22.440893 (2021).  https://www.biorxiv.org/content/early/2021/04/22/2021.04.22.440893.full.pdf.
  * [33] Chen, M. & Ludtke, S. J.  Deep learning-based mixed-dimensional gaussian mixture model for characterizing variability in cryo-em.  _Nature Methods_ 18, 930–936, DOI: 10.1038/s41592-021-01220-5 (2021).
  * [34] Nashed, Y. S. G. _et al._ End-to-end simultaneous learning of single-particle orientation and 3d map reconstruction from cryo-electron microscopy data (2021).  2107.02958.
  * [35] Zhong, E. D., Lerer, A., Davis, J. H. & Berger, B.  Exploring generative atomic models in cryo-em reconstruction (2021).  2107.01331.
  * [36] Rose, Y. _et al._ Rcsb protein data bank: Architectural advances towards integrated searching and efficient access to macromolecular structure data from the pdb archive.  _Journal of Molecular Biology_ 433, 166704, DOI: https://doi.org/10.1016/j.jmb.2020.11.003 (2021).  Computation Resources for Molecular Biology.
  * [37] AlQuraishi, M.  Alphafold at casp13.  _Bioinformatics_ 35, 4862–4865 (2019).
  * [38] Gupta, H., McCann, M. T., Donati, L. & Unser, M.  Cryogan: a new reconstruction paradigm for single-particle cryo-em via deep adversarial learning.  _IEEE Transactions on Computational Imaging_ 7, 759–774 (2021).
  * [39] Ullrich, K., van den Berg, R., Brubaker, M. A., Fleet, D. J. & Welling, M.  Differentiable probabilistic models of scientific imaging with the fourier slice theorem.  _CoRR_ abs/1906.07582 (2019).  1906.07582.
  * [40] Rohou, A. & Grigorieff, N.  CTFFIND4: Fast and accurate defocus estimation from electron micrographs.  _Journal of Structural Biology_ DOI: 10.1016/j.jsb.2015.08.008 (2015).
  * [41] Dempster, A., Laird, N. & Rubin, D. B.  Maximum Likelihood from Incomplete Data via the EM Algorithm.  _Journal of the Royal Statistical Society. Series B (Methodological)_ 39, 1–38 (1977).
  * [42] Miolane, N., Poitevin, F., Li, Y.-T. & Holmes, S.  Estimation of orientation and camera parameters from cryo-electron microscopy images with variational autoencoders and generative adversarial networks.  _CVPR conference, Workshop on Computer Vision for Microscopy Image Analysis_ (2019).  arXiv:1911.08121.
  * [43] Goodfellow, I. J. _et al._ Generative Adversarial Nets.  Tech. Rep. (2014).
  * [44] Gupta, H., Phan, T. H., Yoo, J. & Unser, M.  Multi-cryogan: Reconstruction of continuous conformations in cryo-em using generative adversarial networks.  In _European Conference on Computer Vision_ , 429–444 (Springer, 2020).
  * [45] Gelfand, A. E. & Smith, A. F. M.  Sampling-based approaches to calculating marginal densities.  _Journal of the American Statistical Association_ 85, 398–409, DOI: 10.1080/01621459.1990.10476213 (1990).
  * [46] Duane, S., Kennedy, A., Pendleton, B. J. & Roweth, D.  Hybrid monte carlo.  _Physics Letters B_ 195, 216–222, DOI: 10.1016/0370-2693(87)91197-x (1987).
  * [47] Henderson, R.  Avoiding the pitfalls of single particle cryo-electron microscopy: Einstein from noise.  _Proceedings of the National Academy of Sciences_ 110, 18037–18041 (2013).
  * [48] Kingma, D. P. & Welling, M.  Auto-Encoding Variational Bayes.  In _Proceedings of the 2nd International Conference on Learning Representations (ICLR)_ (2014).
  * [49] Marin van Heel & Michael Schatz  Fourier shell correlation threshold criteria  _Journal of Structural Biology_ 151, 250–262, DOI: 10.1016/j.jsb.2005.05.009 (2005).
  * [50] Harauz G & Marin van Heel  Exact filters for general geometry three dimensional reconstruction  _Optik_ 73, 146-156, (1986).
  * [51] Lawson, C.L. & et al.  Cryo-EM model validation recommendations based on outcomes of the 2019 EMDataResource challenge.  _Nat Methods_ 18, 156–164, DOI: 10.1038/s41592-020-01051-w(2021).
  * [52] Renaud, Jean-Paul & et al.  Cryo-EM in drug discovery: achievements, limitations and prospects  _Nature reviews Drug discovery_ 17, 471–492, DOI: 10.1038/s41592-020-01051-w(2018).
  * [53] Arjovsky, M. & Welling, M.  Wasserstein Generative Adversarial Networks.  In _Proceedings of Machine Learning Research_ (2017)
  * [54] Deng, Lie  The mnist database of handwritten digit images for machine learning research.  In _IEEE Signal Processing Magazine_ (2012)
  * [55] Deng, Jia & et al.  Imagenet: A large-scale hierarchical image database.  In _2009 IEEE conference on computer vision and pattern recognition_ (2009)
  * [56] Thiagarajan, Ganesh & Voyiadjis, George Z.  $\beta$-VAE: Learning basic visual concepts with a constrained variational framework.  In _Proceedings of ICLR_ (2017)
  * [57] Kufareva, Irina & Abagyan, Ruben  Methods of protein structure comparison.  In _Homology modeling_ (2011)
  * [58] Scherer, Martin K & Husic, Brooke E & Hoffmann, Moritz & Paul, Fabian & Wu, Hao & Noé, Frank  Variational selection of features for molecular kinetics.  In _The Journal of chemical physics_ (2019)

## Appendix A: Ab-initio vs Refinement Reconstructions

Section 1 established how cryo-EM methods target either homogenenous or
heterogeneous reconstruction, and refined these categories by drawing a
distinction between methods that are reference-based, and those that are
reference-free — that is, those that do not try to parametrize heterogeneity
as "small" deviations from a reference shape. Here, we introduce an additional
nuance between these methods, and differentiate between methods that operate
ab-initio (ie. from scratch) , or by refinement (ie, "warm-starting " the
algorithm by using an existing and external volume estimate, such as for
example a low resolution estimate). While, for the sake of clarity and
conciseness, we did not elaborate on this distinction in the main text, the
latter introduces an other (independent) axis of variation in our comparison
of existing methods. In fact, methods can target reference-free heterogeneous
reconstruction whilst warm-starting their algorithm close to a known solution
(as is done in the heterogeneous version of RELION[29] ), or vice versa,
prefer to opt for reference-based heterogeneity, with no prior knowledge of
what the reference should look like (e.g. 3DFlex[32] ).

This distinction can be understood as a way of infusing further knowledge into
the recovery. From a statistical and compute science perspective, choosing to
operate by refinement depends on whether or not we have a reasonable prior for
the shape (refinement), or if we’d rather not bias the algorithm with external
inputs (ab initio). In this high-dimensional (and highly non-convex) setting,
imposing a specific starting point for the algorithm is in fact likely to
introduce an additional bias to the solution — thereby potentially allowing it
to recover more easily interesting structures. Conversely, in the absence of
good guesses, ab initio methods prefer to use a random initialization. From a
biological perspective, this boils down to how much we trust potential
external reconstruction methods, and we think that they could be instrumental
in recovering the structure.

##### Ab-initio vs refinement: Pros, Cons, Discussions.

Figure 4 summarizes the methods performing ab-initio reconstruction versus
refinement. While “bias-free”, ab-initio methods can be challenging to fit:
due to the highly non-convex, non-linear nature of the problem, these methods
are at higher risk of inconsistencies, especially when targeting an
heterogeneous reconstruction. In fact, solutions initialised at various random
points might not necessarily converge to same point, and/or get stuck in local
minima. By contrast, warm-starting the problem can yield faster, and more
consistent convergence towards the solution. Warm-starting the solution might
also increase consistency across recovered conformation, thereby preventing
aberrations. Homogeneous reconstruction approaches can also benefit from warm-
starting strategies: in RELION [13] for instance, the hypothesised template is
converted into Fourier coefficients, that are used in the initialisation part
of the inference process. Refinement for homogeneous reconstruction holds
several advantages, especially if the homogeneous reconstruction is either
slow or very sensitive to initialization.

Yet, refinement also holds several disadavantages. This approach requires a
two-step procedure and another fragmentation of the estimation pipeline, as it
relies on a prior estimation of an initial "blunt" value by another pipeline.
Moreover, as highlighted in the review of Bendory et al [20], it also exposes
this solution to model bias: the structure recovered will be biased towards
our initial prior — potentially hindering our ability to detect flaws in this
original template or discovering new conformations. Finally, while the
literature on the sensitivity of the methods to an erroneous warm-starting is
also nonexistent, the impact of a wrong initialization could also either
severely hinder the performance of the algorithm, or push the algorithm to
hallucinate nonexistent solutions (allowing us to recover “Einstein from
noise” [47]— see discussion in part 3).

Figure 4: Works differentiated by the objective of the reconstruction
(homogeneous/heterogeneous), and its initialization (ab-initio/warm-start).

## Appendix B: Differentiability of the Generative Models

As described in the main text, an alternative (and often simpler) generative
model can also be formulated via the Fourier equivalent of Equation (1).

##### Fourier-Slice Theorem.

In Fourier Space, the convolution $*$ of the projected 2D volume by the PSF
conveniently becomes an element-wise matrix multiplication $\odot$ between the
2D Fourier transform of the projected image and that of the PSF. The latter is
better known as the Contrast Transfer Function (CTF), with
$\text{CTF}_{i}=\mathcal{F}(\text{PSF}_{i})$, so that:

$\tilde{X}_{i}=\text{CTF}_{i}\odot(P_{i}^{(R_{i},t_{i})}\circ\tilde{V}^{(i)})+N_{i},\qquad\text{with}\quad
i=1\cdots n.$ (9)

Here, $\tilde{X}_{i}$ is the 2D-Fourier transform of the observed image
$X_{i}$ ($\tilde{X}_{i}=\mathcal{F}(X_{i})$), and $\tilde{V}^{(i)}$ is the 3D
Fourier transform of the molecule (i.e, the volume) $V^{(i)}$, with
$\tilde{V}_{i}=\mathcal{F}(V_{i})$. In Equation (9), the operator
$P_{i}^{(R_{i},t_{i})}$ first extracts out of $\tilde{V}^{(i)}$ a slice
oriented by the rotation $R_{i}$ and then applies a phase shift in Fourier
space accounting for the 2D translation $t_{i}$. Finally, $N_{i}$ represents
the noise. We note that using the operator $P_{i}$ circumvents the 3D-to-2D
projection step and is therefore more computationally efficient.

##### Efficiency and differentiability of the different generative models

The methods described in this review (Eq 1 and 9) leverage the ability to do
inference with gradient-based optimization techniques. This requires the
generative models to be differentiable — that is, to allow the differentiation
of the loss associated with the output (an image) with respect to the
parameters of the model through backpropagation. As explained in Section 1,
the volume $V$ (or $\tilde{V}$ in Fourier space) can always be seen as a
function from $\mathbb{R}^{3}$ to $\mathbb{R}$ (or $\mathbb{C}$). The
parameters defining this function depend on the choice of parametrization, as
summarized in Table 1. To make sure that the model is differentiable, we need
to ensure that all the operations required in the forward model presented by
Equation 1 are differentiable. Table 1 summarizes the operations required to
“rotate”, “translate” and “project” the volume $V$ and obtain the observed
image. Those operations are:

  * •

multiplication with a matrix $R\in\mathbb{R}^{3\times 3}$,

  * •

addition of a vector $\mathbf{t}\in\mathbb{R}^{3}$,

  * •

application of an interpolation kernel,

  * •

summation (if the volume is in image space),

  * •

and analytical integration.

All these operations, if they involve the “differentiable parameters” can be
differentiated through by backpropagation. As shown in Equation (1) (resp.
Equation (9)) the last step of the generative model is the application of the
PSF (resp. CTF) which is done with a convolution (resp. element-wise
multiplication) with the projected volume. This last operation is also
differentiable, which in turn ensures that the whole generative model,
transforming a volume $V$ into an image (potentially in Fourier space), is
differentiable.

|  | Differentiable
---
Parameters
Domain |  | Rotation
---
Translation
| Projection
---
(in image space)
| Voxel
---
Grid
| $V(x)$
---
for $x\in\\{1,\ldots,D\\}^{3}$
| $V:\\{1,\ldots,D\\}^{3}$
---
With interpolation kernel,
$\hat{V}:\mathbb{R}^{3}$, see [39]
|  | $\sum_{x}\hat{V}\hat{(x)}$
---
$\hat{V}:\mathbb{R}^{3}\rightarrow\mathbb{R}$
| Neural
---
Network
Weights $W$ | $V:\mathbb{R}^{3}$ | $\hat{x}=R^{-1}x-\mathbf{t}$ | $\sum_{x}V\hat{(x)}$
| Mixture of
---
Gaussians
$A_{j}$, $c_{j}$, $\sigma_{j}$ |  | $V:\mathbb{R}^{3}$
---
With analytical integration,
$\int V(.,z)dz:\mathbb{R}^{2}$
| $c_{j}\rightarrow\hat{c}_{j}=Rc_{j}+\mathbf{t}$
---
$V\rightarrow\hat{V}$
| $\int\hat{V}(y,z)dz$
---
for $y\in\mathbb{R}^{2}$
Table 1: Differentiable parameters and operations realized in the generative
model, depending of the volume parametrization. $V$ can be equivalently
replace with $\tilde{V}$ in Fourier space. Using an interpolation kernel
transforms the function $V$ defined on $\\{1,\ldots,D\\}^{3}$ into $\hat{V}$
defined on $\mathbb{R}^{3}$.

## Appendix C: Inference and ELBO Computations

This appendix provides more details on the Bayesian formulation of the cryoEM
recovery problem, and on the corresponding methods. As a reminder, in the main
text, we argue that the lower bounds $\mathcal{L}(q,X,\theta)$ are found by
showing that, for any probability distribution $q_{i}$ on the variables
$h_{i}$, the observed log-likelihood can be written as the sum of two terms:

$\displaystyle
L(X,\theta)=\mathcal{L}(q,X,\theta)+\sum_{i=1}^{n}\text{KL}(q_{i}(h_{i})\parallel
p_{\theta}(h_{i}|x_{i}))=\sum_{i=1}^{n}\Big{[}\mathcal{L}_{i}(q_{i},x_{i},\theta)+\text{KL}(q_{i}(h_{i})\parallel
p_{\theta}(h_{i}|x_{i}))\Big{]}$ (10)

where $KL$ is the Kullback-Leibler divergence (KL)defined as $KL(q\parallel
p)=\int q(x)\log\frac{q(x)}{p(x)}dx$, and the terms $\mathcal{L}_{i}$ write:

$\mathcal{L}_{i}(q_{i},x_{i},\theta)=\int_{h_{i}}q_{i}(h_{i})\log
p_{\theta}(x_{i}|h_{i})dh_{i}-\text{KL}\left(q_{i}(h_{i})\parallel
p_{\theta}(h_{i})\right).$ (11)

The objective of this appendix is to derive these bounds and present the
corresponding variational techniques in greater details.

###  C.1 A Bayesian take on cryo-EM

Before we show these bounds, let us begin by a preliminary remark regarding
the model that we are trying to fit. In our formulation in the main text, the
parameters $\theta$ are considered as unknown, fixed values. Alternatively,
one could also try and incorporate the uncertainty of these parameters by
adopting a hierarchical approach, and modelling the parameters $\theta$ as
random variables with prior $p(\theta)$. In this case, inference seeks to
estimate $\theta$ through the maximum of its posterior distribution
$p(\theta|x_{1},...,x_{n})$, written in its logarithm form as:

$\displaystyle\log p(\theta|x_{1},...,x_{n})\propto\log
p(\theta,x_{1},...,x_{n})=\log p_{\theta}(x_{1},...,x_{n})+\log
p(\theta)=L(X,\theta)+\log p(\theta).$ (12)

This can increase robustness by embedding in the design of the algorithm the
inherent variability and uncertainty associated with $\theta$. In fact,
following Scheres [13], adding priors on $\theta$ can be viewed as a kind of
regularization. Such approaches — referred to as Maximum a Posteriori (MAP)
approaches — are chosen in Relion-Refine3D and Relion-Class3D[29],
FSTdiff[39], and CryoSPARC[16]. In the original RELION[13] method for
instance, the underling volume is represented through its Fourier components
$V_{l}$, whose prior is chosen to be zero-mean Gaussian distributedwith
unknown variance $\tau_{l}^{2}$.

We observe that the ML and MAP objectives in Equations (4)-(12) only differ
through the prior term $\log p(\theta)$, while the challenge of the
optimization stems froms the intractable term $L(X,\theta)$. Both approaches
traditionally employ the same optimization techniques.

###  C.2 ELBO Computations

We provide the derivations leading to Equations (10)-(11) and the definition
of the ELBO as a lower bound to the log-likelihood of the generative model. We
start by expressing the observed log-likelihood of the model as:

$\displaystyle L(X,\theta)$ $\displaystyle=\sum_{i=1}^{n}\log p(x_{i}|\theta)$
$\displaystyle=\sum_{i=1}^{n}\left(\log p(x_{i},h_{i}|\theta)-\log
p(x_{i},h_{i}|\theta)+\log p(x_{i}|\theta)\right)$
$\displaystyle=\sum_{i=1}^{n}\left(\log p(x_{i},h_{i}|\theta)-\log
p(h_{i}|x_{i},\theta)\right)\qquad\text{( Bayes rule:}\quad
p(h_{i}|x_{i},\theta)=\frac{p(x_{i},h_{i}|\theta)}{p(x_{i}|\theta)})$

In what follows, we drop the index $i$ and associated sum for convenience.
Consider any distribution $q(z)$. If we multiply both sides of the above
equation by $q$, and integrate the $z$ out, we get:

$\displaystyle\int_{z}L(X,\theta)q(h)dh$ $\displaystyle=\int_{h}\left(\log
p(x,h|\theta)-\log p(h|x,\theta)\right)q(h)dh$ $\displaystyle\implies
L(X,\theta)$ $\displaystyle=\int_{h}\left(\log p(x,h|\theta)-\log
p(h|x,\theta)\right)q(h)dh,$

since $L(X,\theta)$ does not depend on $h$. Now expanding the right hand side:

$\displaystyle L(X,\theta)$ $\displaystyle=\int_{h}q(h)\log
p(x,h|\theta)dh-\int_{h}q(h)\log p(h|x,\theta)dh$
$\displaystyle=\int_{h}q(h)\log
p(x,h|\theta)dh-\int_{h}q(h)\log\frac{p(h|x,\theta)q(h)}{q(h)}dh$
$\displaystyle=\int_{h}q(h)\log p(x,h|\theta)dh-\int_{h}q(h)\log
q(h)dh-\int_{h}q(h)\log\frac{p(h|x,\theta)}{q(h)}dh$
$\displaystyle=\int_{h}q(h)\log p(x,h|\theta)dh-\int_{h}q(h)\log
q(h)dh+KL(q(h)\parallel p(h|x,\theta))$
$\displaystyle=\mathbb{E}_{q}[\log(p_{\theta}(x,h)]-\int_{h}q(h)\log
q(h)dh+KL(q(h)\parallel p(h|x,\theta))$
$\displaystyle=\mathbb{E}_{q}[\log(p_{\theta}(x,h)]+H(q)+KL(q(h)\parallel
p(h|x,\theta)),$

by introducing the entropy of $q$ as $H(q)=-\int_{z}q(z)\log q(z)dz$, and
writing the KL-divergence between two distributions $p_{1}$ and $p_{2}$ as
$KL(p_{1}||p_{2})=\int p_{1}(x)\log(\frac{p_{1}(x)}{p_{2}(x)})dx$.

In other words, for any distribution $q$ on the latent variable $H$, the
observed marginal likelihood $L(X,\theta)$ is the sum of three terms: (a) the
expected log-likelihood, assuming that $h$ has distribution $q$, (b) the
Entropy of $q$, which can be understood of the amount of uncertainty in the
estimation of $h$ and (c) the KL-divergence between $q$ and the true posterior
$p(h|x,\theta)$, which measures how different these distributions are.

Note that, in the previous equations, terms (a) and (b) can be computed, but
term (c) — the KL divergence between the proposed distribution $q$ for $h$ and
the actual posterior distribution $p(h|x,\theta)$ is not unknown, because we
do not know $p(h|x,\theta)$. However, because the KL-divergence is always non-
negative, we know that:

$L(X,\theta)\geq\mathbb{E}_{q}[\log(p_{\theta}(x,h)]+H(q).$

Thus, introducing the ELBO term
$\mathcal{L}(q,X,\theta)=\mathbb{E}_{q}[\log(p_{\theta}(x,h)]+H(q)$, we get:

$L(X,\theta)=\mathcal{L}(q,X,\theta)+KL(q(h)\parallel p(h|x,\theta)).$ (13)

which provides Equation (10) and the right hand side of Equation (11).
Developping the ELBO term gives:

$\displaystyle\mathcal{L}(q,X,\theta)$ $\displaystyle=\int_{h}q(h)\log
p(x,h|\theta)dh-\int_{h}q(h)\log q(h)dh$ $\displaystyle=\int_{h}q(h)\log
p(x|h,\theta)dh+\int_{h}q(h)\log p(h|\theta)dh-\int_{h}q(h)\log q(h)dh$
$\displaystyle=\int_{h}q(h)\log p(x|h,\theta)dh-KL(q(h)\parallel p(h|\theta),$

which gives the left side of Equation (11).

###  C.3 Expectation-maximization (EM) Algorithm

Historically, inference in latent variable models has been achieved through
the Expectation-maximization (EM) algorithm [41]. The EM algorithm is often
introduced as a data imputation technique (see following subsection), but it
can also be understood as a dual ascent algorithm — a perspective that allows
unifying most cryoEM inference methods and that we present here.

##### EM as a Dual Ascent Algorithm

The EM algorithm can be seen as a dual ascent procedure, i.e. as a
maximization-maximization procedure, that leverages Equation (10) rewritten
as:
$\mathcal{L}(q,X,\theta)=L(X,\theta)-\sum_{i=1}^{n}\text{KL}(q_{i}(h_{i})\parallel
p_{\theta}(h_{i}|x_{i}))$. Starting with an initial guess $\theta^{(0)}$ of
the parameter, the EM algorithm performs two steps at each iteration $t$:

(a) Inference on hidden variables (E-step)

Given the current $\theta^{(t-1)}$, maximize each term
$\mathcal{L}_{i}(q_{i},\theta^{(t-1)})=L_{i}(X_{i},\theta^{(t-1)})-KL(q_{i}(h_{i})||p_{\theta^{(t-1)}}(h_{i},x_{i}))$
(by Equation 13) with respect to the distribution $q_{i}$ defined on the
hidden variables $h_{i}$. Since $\theta^{(t-1)}$ is fixed, the best choice of
$q$ is the one that minimizes the KL term
$KL(q_{i}(h_{i})||p_{\theta^{(t-1)}}(h_{i},x_{i}))$:

$\text{For each $i=1,...,n$:}\quad\quad
q_{i}^{(t)}(h_{i})=\operatorname*{arg\,max}_{q}\mathcal{L}_{i}(q,X,\theta^{(t-1)})=\operatorname*{arg\,min}_{q}KL(q_{i}(h_{i})||p_{\theta^{(t-1)}}(h_{i},x_{i}))=p_{\theta^{(t-1)}}(h_{i}|x_{i}).$
(14)

This choice $q_{i}^{(t)}(h_{i})=p_{\theta^{(t-1)}}(h_{i}|x_{i})$ corresponds
to the green distribution in Figure 5 and makes the lower bound $\mathcal{L}$
tangent to $L$ at $\theta^{(t-1)}$ in Figure 2 (left).

(b) Maximization on the model’s parameters (M-step)

Given the current $q^{(t)}=\\{q_{1}^{(t)},...,q_{n}^{(t)}\\}$, maximize
$\mathcal{L}(q^{(t)},X,\theta)$ with respect to the model parameters $\theta$:

$\theta^{(t)}=\operatorname*{arg\,max}_{\theta}\mathcal{L}(q^{(t)},X,\theta)=\operatorname*{arg\,max}_{\theta}\sum_{i=1}^{n}\mathcal{L}_{i}(q_{i}^{(t)},X,\theta).$
(15)

The argument maximum of $\mathcal{L}$ is used to update $\theta^{(t-1)}$ to
$\theta^{(t)}$ in Figure 2 (left).

These steps are iterated until convergence in $\theta$. In other words, one
first finds the "best" lower bound to $L(X,\theta)$ given $\theta^{(t)}$ by
choosing the one that is tangent to $L(X,\theta)$ at $\theta^{(t)}$. Once this
bound has been established, we subsequently maximize this bound with respect
to the parameter $\theta$. We note that to alleviate potential computational
bottleneck in the M step (respectively $\mathcal{L}(q^{(t)},X,\theta)$), the
closed form solution can be replaced by simply taking a gradient step. This
gives rise to “gradient EM" methods, which we will refer to later, as all
methods introduced in the next subsections can be formulated as “gradient"
methods too.

##### Computational Bottleneck: E-step.

One of the issues with the EM algorithm relies in the computation of the
posterior of the latent variable $h_{i}$ given $x_{i}$ and current estimate
$\theta^{(t-1)}$:

$p_{\theta^{(t-1)}}(h_{i}|x_{i})=\frac{p_{\theta^{(t-1)}}(x_{i}|h_{i})p(h_{i})}{p_{\theta^{(t-1)}}(x_{i})}=\frac{p_{\theta^{(t-1)}}(x_{i}|h_{i})p(h_{i})}{\int_{h_{i}}p_{\theta^{(t-1)}}(x_{i},h_{i})dh_{i}}.$
(16)

Unless the posterior can be computed efficiently, this computation can be
lengthy. In fact, unless the priors of the different variables are specified
using conjugacy — which would constrain the choice of prior distributions that
we could consider —, it is difficult to get closed-form updates for the
posterior. Alternatively, one can evaluate $p_{\theta^{(t-1)}}(h_{i}|x_{i})$
for each $h_{i}$ discretized on a grid by computing the integral in Equation
(16) via a Riemann sum as in RELION[29] or via Importance Sampling as in
CryoSPARC[16]. RELION[29] and CryoSPARC[16]use the EM algorithm with MAP
estimation of $\theta$, as opposed to ML estimation.

This discretization remains however a computationally intensive approach.
Consequently, despite many advances leveraging GPU computing, the integral in
Equation (16) represents the main computational bottleneck in cryo-EM
reconstruction methods, and the main reason we might want to look for other,
more efficient alternatives.

Figure 5: Possible choices for the distribution $q$ on the hidden variables
$h$, where we have dropped the indices $i$ for convenience of notations. The
choice of $q$ determines the lower bounds to $L(X,\theta)$ in Figure 2. The
inference methods used in cryoEM reconstruction can be contrasted by the
choice of $q$: as the true posterior $p_{\theta}(h|x)$ represented by the
green line, as the Dirac distribution $\delta(h)$ represented by the vertical
dashed gray line or as a distribution within a known parameterized family
$\mathcal{Q}$, such as the family of Gaussian distributions, represented by
the brown dashed line.

###  C.4 Variations on the EM algorithm

We present details on the variations on the EM algorithm discussed in the main
text.

####  C.4.1 Modal EM algorithm

The computational bottleneck observed in the EM algorithm has also motivated
the use of approximations in the E-step. The modal approximation of EM amounts
to finding a proxy for the posterior $p_{\theta^{(t-1)}}(h_{i}|x_{i})$ in the
E-step. This approximation replaces the posterior by its "mode"
$\hat{h}_{i}=\operatorname*{arg\,max}_{h}p_{\theta^{(t-1)}}(h|x_{i})$, i.e.
the most probable $\hat{h}_{i}$ for each image $x_{i}$.

##### Modal EM: The algorithm.

The E-step of the EM algorithm is replaced by its modal approximation, but the
M-step is largely unchanged:

(a) Inference on hidden variables $h_{i}$

Given current $\theta^{(t)}$, compute the modes:

$\text{For each
$i=1,...,n$:}\quad\quad\hat{h}_{i}^{(t)}=\operatorname*{arg\,max}_{h}p_{\theta^{(t-1)}}(h|x_{i})=\operatorname*{arg\,max}_{h}p_{\theta^{(t-1)}}(x_{i}|h)p(h).$
(17)

This amounts to approximating the solution of the E step (which, we’ve shown,
is the posterior of $h$ given $X$ and $\theta^{(t-1)}$) by
$q_{i}^{(t)}(h_{i})=\delta(\hat{h}_{i}^{(t)})$, the Dirac distribution at
$\hat{h}_{i}^{(t)}$, represented in gray in Figure 5. This choice for
$q_{i}^{(t)}(h_{i})$ creates a lower bound $\mathcal{L}$ that may not be
tangent to $L$ in Figure 2 (right). Since this approximation does not require
us to compute complicated posteriors, it is a simpler, but looser ELBO to the
observed likelihood.

(b) Maximization on model’s parameter $\theta$

The parameter $\theta$ is updated via:

$\theta^{(t)}=\operatorname*{arg\,max}_{\theta}\sum_{i=1}^{n}\mathcal{L}(\delta(\hat{h}_{i}^{(t)}),\theta)=\operatorname*{arg\,max}_{\theta}\sum_{i=1}^{n}p_{\theta}(x_{i},\hat{h}_{i}^{(t)}).$
(18)

The argument maximum of $\mathcal{L}$ is used to updated $\theta^{(t-1)}$ to
$\theta^{(t)}$ in Figure 2 (right), just as in the EM algorithm in Figure 2
(left).

##### Modal EM in RELION[29] , CryoSPARC[16], 3DFlex[32] , CryoDRGN[23] ,
CryoFold[35]

The EM algorithm with modal approximation is used in cryo-EM in CryoSPARC[16],
3DFlex[32] , CryoDRGN[23] and CryoFold[35] to infer the latent variable
associated with the rotation $R_{i}$. In other words, they estimate the best
rotation $\hat{R}_{i}$ for each input image $i$, relying on a branch-and-bound
optimization algorithm to address this maximization. This approach is also
used in RELION[29] to estimate the best noise variance $\hat{\sigma}_{i}$ that
is associated with each image $i$.

####  C.4.2 Variational EM algorithm

The modal approximation in the EM algorithm speeds up the E-step; yet it has
the drawback of summarizing the whole posterior distribution
$p_{\theta^{(t)}}(h_{i}|x_{i})$ by a single estimate $\hat{h}_{i}$, reducing
accuracy and leading to non-tangent lower bounds $\mathcal{L}$. Variational
Inference (VI), also known as Variational Bayes, has appeared in cryo-EM as a
compromise between the guarantees of the EM and the efficiency of modal EM
during the computation of the E-step.

Variational Inference (VI)replaces the evaluation of the posterior of the
latent variables $p_{\theta^{(t-1)}}(h_{i}|x_{i})$, by an optimizing over a
different family of candidate distributions $\mathcal{Q}$, called the
“variational family”. $\mathcal{Q}$ is typically a family with a restricted
form, so that the updates are easier to perform: for instance, if there is a
complex dependency structure between entries in our parameter $\theta$, we
might select $\mathcal{Q}$ to be a family of posteriors where the entries are
independent, so that the updates are simpler to perform: this is usually the
‘mean field’ approximation to the posterior of $\theta$. The $q$ that is
chosen in Equation 11 is the "closest" to the posterior, where “closest" is
defined in terms of the KL divergence (see Equation (19) below).

The variational family $\mathcal{Q}=\\{q_{\eta}|\eta\\}$ with parameters
$\eta$ is typically chosen to be a family of Gaussian distributions, in which
case $\eta$ represents the mean and variance, i.e.
$\mathcal{Q}=\\{q_{\eta}=N(\mu,\sigma)|\eta=(\mu,\sigma)\\}$. Other approaches
consider the variational family of Gaussian distributions with fixed isotropic
variance $\sigma_{0}$, i.e.
$\mathcal{Q}=\\{q_{\eta}=\mathcal{N}(\eta,\sigma_{0})|\eta\\}$. In modal EM,
the parametric family can be seen as a family of Dirac distributions, as in
Equation (18).

##### Variational EM as a Dual Ascent Algorithm

The E-step of the EM algorithm is replaced by VI and the M-step is unchanged.
The VI framework thus gives rise to a variational EM which performs the
following two steps at each iteration $t$:

(a) Inference on hidden variables $h_{i}$

The $n$ parameters $\eta_{i}$ corresponding to each $q_{i}$ are computed via:

$\text{For each
$i=1,...,n$:}\quad\quad\eta_{i}^{(t)}=\operatorname*{arg\,min}_{\eta}KL(q_{\eta}||p_{\theta^{(t-1)}}(h_{i}|x_{i})$
(19)

This choice for $q$ is represented in brown in Figure 5. Similarly to the
Modal EM, this creates a lower bound $\mathcal{L}$ that may not be tangent to
$L$ in Figure 2 (right).

(b) Maximization on model’s parameter $\theta$

The parameter $\theta$ is updated via:

$\theta^{(t)}=\operatorname*{arg\,max}_{\theta}\sum_{i=1}^{n}\mathcal{L}_{i}(q_{\eta_{i}}^{(t)},\theta).$
(20)

The argument maximum of $\mathcal{L}$ is used to updated $\theta^{(t-1)}$ to
$\theta^{(t)}$ in Figure 2 (right), just as in the EM and the Modal EM
algorithms.

Thus, the Variational Inference approximation in the E-step replaces the
samplings in evaluation of the true posterior, required to compute the
integral in Equation (16), by an optimization on the fewer parameters
$\eta_{i}$ parametrizing each approximate posterior $q_{\eta_{i}}$. This
optimization can be efficiently conducted with (stochastic) gradient descent.
In this sense, it also represents a solution to the computational bottleneck
of the E-step.

##### Variational EM in 3DFlex[32] and FSTdiff[39]

Variational gradient EM is used in 3DFlex[32] to infer the rotation variable
$R_{i}$, using the variational family of Gaussian distributions with fixed
isotropic variance $\sigma_{0}$.Variational EM is also used in FSTdiff[39]in
the context of homogeneous reconstruction. Here, VI happens both in Step (a),
to estimate the approximate posterior of the rotation $R_{i}$, and in Step (b)
to also estimate the approximate posterior of the volume $V$, which is a
parameter included in $\theta$. The approximate posterior of $V$ provides an
approximate measure of uncertainty on the homogeneous reconstruction.

##### Amortized Variational EM as Dual Optimization Algorithm

Amortized Inference in the E-step of variational EM leads to an amortized
variational EM. In this case, the algorithm can be written as:

(a) Inference on hidden variables $h_{i}$

The parameters $\xi$ parametrizing the encoder is computed via:

$\xi^{(t)}=\operatorname*{arg\,max}_{\xi}\sum_{i=1}^{n}KL(q_{\text{Enc}_{\xi}(x_{i})}||p_{\theta^{(t-1)}}(h_{i}|x_{i})$
(21)

which generates $n$ distributions $q_{\eta_{i}}^{(t)}$ parametrized by
$\eta_{i}=\text{Enc}_{\xi^{(t)}}(x_{i})$.

(b) Maximization on model’s parameter $\theta$

The parameter $\theta$ is updated via:

$\theta^{(t)}=\operatorname*{arg\,max}_{\theta}\mathcal{L}(q_{\eta_{i}}^{(t)},\theta).$
(22)

This method has the advantage of further improving the flexibility and
efficiency of the Variational EM algorithm, by allowing the updates to be even
more efficient by parametrizing the $q_{i}$ as a function of the input
$x_{i}$.

###  C.5 Generative Adversarial Networks

Cryo-EM reconstruction methods from Subsections C.4.1-C.4.2 have naturally
explored adaptations of the computationally expensive E-step from the EM
algorithm of Subsection C.3. As a result, they differ in their inference on
the hidden variables of the generative model. Yet, the only hidden variable of
interest is the conformation variable and associated molecular volumes. As a
result, recent works have explored methods that avoid the inference on the
“nuisance" latent variables altogether. This is the approach taken by
Generative Adversarial Network (GAN).

##### Generative Adversarial Networks as Minimax Games

A generative adversarial model (GAN)[43] is a method that estimates the
parameters of a generative model, such as the one in Equation (1), through an
adversarial process. The GAN trains a “generator" to produce images $x_{i}$
that best capture the training data distribution, while a discriminator
estimates the probability that a given image $x_{i}$ came from the training
data rather than the generator.

In other words, the generator and discriminator play the following two-player
minimax game with value function $V$:

$\displaystyle V=\min_{\theta}\max_{\phi}\mathbb{E}_{x\sim
p_{\text{data}}(x)}\left[\log\text{Dis}_{\phi}(x)\right]+\mathbb{E}_{h\sim
p(h)}\left[\log\left(1-\text{Dis}_{\phi}(f_{V_{\theta}}(h))\right)\right],$

where the generator is the cryo-EM generative model defined in Equation (1)
and the discriminator has weights $\phi$ and is denoted by
$\text{Dis}_{\phi}$. In this equation, $p_{\text{data}}(x)$ represents the
probability distribution of the images, and $p(h)$ is a prior distribution on
the hidden variables. The GAN training iterates two steps:

(a) Discriminator - Update of parameters $\phi$

, according to the gradient step:

$\phi^{(t+1)}=\phi^{(t)}+\nabla_{\phi}\frac{1}{n}\sum_{i=1}^{n}\Big{[}\log\text{Dis}_{\phi}(x_{i})+\log\left(1-\text{Dis}_{\phi}(f_{V_{\theta}^{(t)}}(h_{i}))\right)\Big{]},$
(23)

(b) Generator - Estimation of model’s parameter $\theta$

, according to the gradient step:

$\theta^{(t+1)}=\theta^{(t)}-\nabla_{\theta}\frac{1}{n}\sum_{i=1}^{n}\Big{[}\log\left(1-\text{Dis}_{\phi^{(t)}}(f_{V_{\theta}}(h_{i}))\right)\Big{]},$
(24)

where in each Step (a) or (b), the "hidden variables" $h_{i}$ are sampled
according to a prior distribution $p_{h}$. Even though a GAN iterates two
steps, including one step related to the estimation of the model’s parameters,
its framework differ from the variations of the EM algorithm in the sense that
the hidden variables are not inferred: it is enough to be able to randomly
sample from them from some prior distribution.

##### GANs in CryoGAN[38]and Multi-CryoGAN[44]

CryoGAN[38]uses a GAN to perform homogeneous reconstruction of molecular
volumes, while using uniform sampling on the rotation, translation and CTF
hidden variables. Multi-CryoGAN[44]introduces heterogeneous reconstruction
with this approach. Both CryoGAN[38]and Multi-CryoGAN[44]use the Wasserstein
variant [53] of the traditional GANs.

## Appendix D: Details on (Variational) Autoencoders

Autoencoders (AEs) and variational autoencoders (VAEs) are the main
realizations of the amortized variational inference approaches in the cryo-EM
reconstruction literature. This appendix provides additional details to link
the traditional presentations of the AEs and VAEs to the framework described
in Section 2.

###  D.1 Autoencoders

An autoencoder traditionally aims to minimize the following loss function:

$\ell(\theta,\xi)=\sum_{i=1}^{n}||x_{i}-\hat{x_{i}}||^{2}=\sum_{i=1}^{n}||x_{i}-f_{V_{\theta}}(h_{i})||^{2}=\sum_{i=1}^{n}||x_{i}-f_{V_{\theta}}(\text{Enc}_{\xi}(x_{i})||^{2},$
(25)

where $f_{V_{\theta}}$ denotes here the output of the generative model (such
as the one in Equation 1) and is usually called a reconstruction:
$\hat{x}_{i}=f_{V_{\theta}}(h_{i})$, while $\ell(\theta,\xi)$ is called the
reconstruction loss. Note that the mean square error is used to quantify the
reconstruction loss, but it can be replaced by other metrics. For example, the
binary-cross entropy is a metric traditionally used to compare images $x_{i}$
and $\hat{x}_{i}$. In this context, the goal is to perform the double
minization in $\theta$ and $\xi$ such that:

$\hat{\theta},\hat{\xi}=\operatorname*{arg\,min}_{\theta,\xi}\ell(\theta,\xi).$
(26)

By considering the decoder as latent variable model:
$\hat{x}_{i}=f_{V_{\theta}}(h_{i})+\epsilon_{i}$ with $\epsilon_{i}$ a
standard Gaussian noise, this loss corresponds to the negative log-likelihood,
informed by the encoder.

To minimize this objective, the AE takes a gradient step at each iteration
$t$, or backward pass through the network, such that:

(a) Encoder - Inference on latent variables $h_{i}$

The encoder updates its weight $\xi$ through a gradient step with learning
rate $\alpha$:

$\xi^{(t)}=\xi^{(t-1)}-\alpha\nabla_{\xi}\ell(\theta,\xi).$ (27)

(b) Decoder - Estimation of model’s parameter $\theta$

The decoder updates its weight $\theta$ through a gradient step with learning
rate $\alpha$:

$\theta^{(t)}=\theta^{(t-1)}-\alpha\nabla_{\theta}\ell(\theta,\xi).$ (28)

We note that this gradient descent is usually performed via stochastic
gradient descent, such that only a mini-batch of the data is considered at
each iteration to compute $\ell$, as opposed to the full dataset of $n$
images.

###  D.2 Variational Autoencoders

A variational autoencoder traditionally aims to minimize the following loss
function, which is the negative ELBO [48]:

$\mathcal{L}(\xi,\theta)=\sum_{i=1}^{n}\mathbb{E}_{q_{\xi}(h_{i}|X_{i}))}(\log(p_{\theta}(X_{i}|z))+KL(q_{\xi}(h_{i}|X_{i})||p(h_{i})).$
(29)

In this loss, the first term is called the reconstruction term, estimated with
one Monte Carlo sample through the so-called “reparametrization trick", and is
akin to the reconstruction loss of the AE. The second term is a KL divergence
term that is called the regularization term, as it regularizes the posterior
of the latent variable $h_{i}$ by forcing it to be close to the prior
$p(h_{i})$ of $h_{i}$ which is modeled by a standard Gaussian distribution.

We can rewrite this loss by using the functions $\text{Enc}_{\xi}$ and
$\text{Dec}_{\theta}$, and assuming that the Monte Carlo sampling performed to
compute the expectation happens by sampling a unique $\tilde{h}_{i}$ through
$q_{\xi}(h_{i}|X_{i})$ which is the convention adopted in these architectures:

$\mathcal{L}(\xi,\theta)=\sum_{i=1}^{n}||x_{i}-\hat{x_{i}}||^{2}-KL(q_{\xi}(h_{i}|x_{i})\parallel
p(h_{i}))\\\
=\sum_{i=1}^{n}||x_{i}-\text{Dec}_{\theta}(\text{Enc}_{\xi}(x_{i})||^{2}-KL(q_{\text{Enc}_{\xi}}(h_{i}|x_{i})\parallel
p(h_{i})),\\\ $ (30)

where the KL term has a closed form in terms of the output of the encoder, due
to the fact that the $q$ distribution belongs to a Gaussian family of diagonal
covariance.

By realizing that the decoder only participates in the first term of the loss
function, the VAE takes a gradient step at each iteration $t$, or backward
pass through the network, such that:

(a) Encoder - Inference on latent variables $h_{i}$

The encoder updates its weight $\xi$ through a gradient step of learning rate
$\alpha$ such that:

$\xi^{(t)}=\xi^{(t-1)}-\alpha\nabla_{\xi}\mathcal{L}(\xi,\theta).$ (31)

(b) Decoder - Estimation of model’s parameter $\theta$

The decoder updates its weight $\theta$ through a gradient step with learning
rate $\alpha$:

$\theta^{(t)}=\theta^{(t-1)}-\alpha\nabla_{\theta}\mathcal{L}(\xi,\theta)=\theta^{(t-1)}-\alpha\nabla_{\theta}\ell(\theta,\xi).$
(32)

The main difference with the autoencoder is that the latent variable $h_{i}$
is considered as a random variable, rather than a fixed deterministic value.
That is, $h_{i}$ is endowed with a parametric probability distribution
represented by $q$ — which is a Gaussian distribution with parameters output
by the encoder, so that $h_{i}\sim
N(\mu_{\xi}(X_{i}),\sigma^{2}_{\xi}(X_{i}))$ and
$\text{Enc}_{\xi}(X_{i})=(\mu_{x}i(X_{i}),\sigma^{2}_{\xi}(X_{i})$.
Considering $h_{i}$s as random variables has been shown to lead superior
reconstruction results over the autoencoder.

The VAE loss can be adapted as:
$\mathcal{L}(\xi,\theta)=\sum_{i=1}^{n}\mathbb{E}_{q_{\xi}(z|X_{i}))}(\log(p(X_{i}|z))+\beta
KL(q_{\xi}(X_{i}|Z||p(x))$, where $\beta$ is an additional hyperparameter
introduced in $\beta$-VAE[56]. Traditional VAEs have $\beta=1$, to ensure that
the actual negative ELBO is minimized. Yet, to prevent pathological issues in
the fitting of VAEs (including posterior collapse), recent work has shown that
counterbalancing the reconstruction error with the KL divergence through a
$\beta$ could yield superior results.

In the cryo-EM implementations of the VAE architectures, the ELBO loss can be
supplemented with additional terms that we name “structure losses" for now. We
will explain it in the next appendix as it depends on the structure of the
decoder (i.e. exact choice of generative model) chosen by the method.

In practive, CryoPoseNet[34], CryoAI[18]and E2GMM[33] use an autoencoder, and
E2GMM[33] additionally implements an variational autoencoder with a
variational family of Gaussian distributions with fixed isotropic variance
$\sigma_{0}$. The AE architecture is used in CryoPoseNet[34]with a traditional
L2 reconstruction loss, in CryoAI[18]with a “symmetrized" L2 reconstruction
loss and in E2GMM[33] with a tailored reconstruction loss that relies on the
Fourier ring correlation (FRC) reconstruction metric.

CryoVAEGAN[42] , CryoDRGN[23] , CryoFold[35] , and AtomVAE[17] use VAEs with a
variational family of Gaussian distributions with diagonal covariance matrix,
to respectively infer $(R_{i},\text{CTF}_{i})$ (CryoVAEGAN), $z_{i}$ (CryoDRGN
and CryoFold) and $(z_{i},R_{i})$ (atomVAE) – see Table 4.

The VAE architecture is used in CryoDRGN[23] with the negative ELBO loss, and
in CryoVAEGAN[42] with the negative ELBO loss extended with the $\beta$ hyper-
parameter described in the supplementary materials, an additional geometric
regularization term, while the reconstruction loss relies on the binary cross-
entropy as opposed to the L2 reconstruction loss. AtomVAE[17] also implements
a VAE with a modified ELBO loss that relies on important sampling (not
detailed in this review), that leverages the $\beta$ hyperparameter and
additionally includes a “structure loss" — see supplementary materials.

## Appendix E: Constraints

As described in the introduction and following the exposition by Scheres [13]
(2012), the molecule reconstruction problem is a difficult, highly non-linear
inverse problem, which makes the parametrization of the shape and associated
constraints particularly important. This parametrization amounts to impose
structure on the desired reconstructed shapes — either by leveraging domain
knowledge on the properties of molecular volumes, or by using external
information to guide the reconstruction.

From a physics/biology perspective, this can be seen as a necessary enrichment
of the cryo-EM data with either external assumptions and/or domain knowledge
on the properties of the solution to ensure a more accurate recovery. From a
statistics perspective, formulated this way, the problem rapidly takes on a
Bayesian flavour, and the objective of this step is to find the right “prior"
on the distribution of our latent variables. From a computational perspective,
this corresponds to adding additional terms to the loss/objective defining the
optimization problem, and effectively explains what we term “structure loss"
in Section 2. We detail here the additional constraints that can equip the
inference methods across the reconstruction algorithms.

###  E.1 Smoothness

Smoothness refers to the property by which a signal, or quantity of interest,
varies with "no abrupt change" over continuous regions. As explained by
Scheres (2012) [13], "because macromolecules consist of atoms that are
connected through chemical bonds, the scattering potential will vary smoothly
in space, especially at less than atomic resolution." Smoothness of the
recovered scattering potential $V$ thus appears to be a reasonable assumption,
which is implemented in different ways depending on the parametrization
chosen:

Smoothness of the 3D image.

The smoothness assumption holds in image space as the image is a projection of
the electron field — which is itself continuous. This smoothness is encoded by
the normal distribution. This also calls to mind the "ridge penalty", a
similar type of regularisation in statistics.

Smoothness of the Fourier coefficients.

The smoothness of the 3D density map translates into smoothness over
neighbouring Fourier coefficients. In a Bayesian pipeline, this is typically
parametrized by assuming that these coefficients are sampled from a normal
distribution. The RELION algorithm (Scheres [13]) is based on such a
smoothness assumption. This prior is encoded by assuming independent Gaussian
priors on the Fourier components of the signal: $V_{l}\sim N(0,\tau_{l}^{2})$.
Note here that the algorithm is not encouraging any other type of structure
(e.g sparsity through spike-and-slab prior, etc).

Smoothness of the deformation field.

3DFlex [32] directly exploits the knowledge that conformational variability of
a protein is the result of physical processes that transport density over
space. This means that mass and local geometry are preserved. As a result,
this method implements a convection operator that outputs the deformation
field. 3DFlex exploits prior knowledge of smoothness and local rigidity in the
deformation field.

Smoothness of the function over 3D coordinates in Fourier domain.

Cryo-DRGN[23] and Cryo-Fold[35] represent the volume as the function
$f:\Omega^{3}\rightarrow\mathbb{R}$ over a 3D domain. In practice, this
function is implemented by a neural network, which constrains it to be
continuous, and possibly smooth if the activation functions used by the
network are themselves smooth (sigmoids, for example).

###  E.2 Rigidity

Other constraints exploit physical properties and knowledge of the system to
constrain the reconstruction. Specifically, due to the fact that the molecules
studied in cryo-EM are frozen, we can assume that their conformational
heterogeneity will not present large variations. As a result, we can assume
that the molecular structure has some type of “rigidity". This is implemented
in practice by refraining the variables describing the shape heterogeneity
from varying excessively. Depending on the parametrization chosen to represent
the volume heterogeneity, the rigidity constraint takes different forms:

Rigidity of the deformation field.

The deformation field that parametrizes volume heterogeneity in 3D-Flex[32]
can be constrained to only generate “small" deformation, through a
regularization term. This is integrated as a “structure loss" of the deep
learning training procedure of 3D-Flex[32].

Rigidity in (pseudo)-atoms coordinates.

The (pseudo-) atoms coordinates that parametrize the volume heterogeneity in
Cryo-Fold[35] are constrained to be close to a base conformation’s
coordinates, via an $L2$ penalization on the deformations. This forms the
“structure loss" integrated in the objective function for the training
procedure.

Rigidity in deviation from (pseudo)-atoms coordinates.

Here, the volume is parametrized as a set of deviations of pseudo-atoms
coordinates, compared to a base conformation, as in Cryo-DeepMind[17]. The
rigidity of the molecular structure is thus enforced by adding a constraint
that takes the form of a L2 regularization constraining the $\Delta c_{j}$ to
be small, and effectively forming the “structure loss" added to the
optimization objective of this approach.

We also note that the methods that do not assume any heterogeneity in volumes
$V_{i}$ but rather model the volume as a unique possible conformation $V$, are
essentially relying on the rigidity assumption, assuming that variations
around $V$ are not large.

## Appendix F: Experiments

This appendix summarizes the experiments conducted in the papers cited in this
review, in Table 2. As mentioned in the main text, we observe a great
diversity of datasets, which explains the difficulty encountered in comparing
methods’ performances and accuracies.

| Biomolecule | # of images | Image size | Data Type | Noise
---|---|---|---|---|---
|  | 50,000 | 128 x 128 | Synthetic | No/Yes
| Ribosome 80S (EMPIAR-10028) | 105,247 | 90 x 90 | Experimental | Yes
| Ribosome 50S (EMPIAR-10076) | 131,899 | 90 x 90 | Experimental | Yes
CryoDRGN[23] | Protein complex | 50,000 | 64 x 64 | Synthetic | N/A
CryoFold[35] | Haemoglobin (PDB 5NI1) | 50,000 | 128 x 128 | Synthetic | Yes
| Tri-snRNP spliceosome (EMPIAR-10073) | 102,500 | 180 x 180 | Experimental | Yes
3DFlex[32] | TRPV1 ion-channel (EMPIAR-10059) | 200,000 | 128 x 128 | Experimental | Yes
CryoPoseNet[34] | E.coli adenylate kinase (PDB 4AKE1) | 9,000 | 128 x 128 | Synthetic | No/Yes
| Ribosome 50S (EMPIAR-10076) | 124,900 | N/A | Experimental | Yes
| Precatalytic spliceosome (EMPIAR-10180) | 327,490 | N/A | Experimental | Yes
E2GMM[33] | SARS-CoV-2 spike protein (EMPIAR-10492) | 55,159 | N/A | Experimental | Yes
FSTdiff[39] | GroEL-GroES protein | 40,000 | 128 x 128 | Synthetic | Yes
|  | 2, 544 | 128 x 128 | Synthetic | No
CryoVAEGAN[42] | Ribosome 80S (EMPIAR-10028) | 5,119 - 8,278 - 4,917 | 128 x 128 | Experimental | Yes
AtomVAE[17] | Aurora A Kinase (Simulated) | 63,000 | 64 x 64 | Synthetic | Yes
| 80S ribosome (PDB 3J79 and 3J7A) | 1,000,000 | 128 x 128 | Synthetic | No/Yes
| SARS-CoV-2 spike protein (PDB 6VYB) | 1,000,000 | 128 x 128 | Synthetic | No/Yes
| Spliceosome (PDB 5NRL) | 1,000,000 | 128 x 128 | Synthetic | No/Yes
CryoAI[18] | 80S ribosome (EMPIAR-10028) | 105,247 | 256 x 256 | Experimental | Yes
CryoGAN[38] |  $\beta$-galactosidase | 41,000 | 180 x 180 | Synthetic | Yes
|  $\beta$-galactosidase (EMPIAR-10061) | 41,123 | 192 x 192 | Experimental | Yes
Multi-CryoGAN[44] | Heat-shock protein Hsp90 | 100,000 | 32 x 32 | Synthetic | Yes

Table 2: Summary of synthetic and experimental cryo-EM datasets on which the
reconstruction methods have been tested. The light red background indicates
methods that perform homogeneous reconstruction, and the light yellow
background heterogeneous reconstructions.

## Appendix G: Classification of Reconstruction Methods

This last appendix provides details in the classification of the
reconstruction methods, in the form of two tables. Table 3 compares the
generative models and Table 4 compares the inference methods.

Volume Param. |  Space | Conformational Model |  Approach |  Hidden Variable |  Known Variable |  Reference Volume
---|---|---|---|---|---|---
|  |  |  | CryoPoseNet[34] |  | -
---
Rotation
-
| CTF
---
-
Translation
Free
|  | Image |  Homogeneous $z\longrightarrow V$ | CryoGAN[38] |  | -
---
Rotation
Translation
| CTF
---
-
-
Free
Discrete Representation |  Voxel Grid |  |  | Heterogeneous
---
$z\longrightarrow V$
Multi-CryoGAN[44] |  | -
---
Rotation
Translation
| CTF
---
-
-
Free
|  |  |  | Heterogeneous
---
$z\longrightarrow f(U(z_{j}),V_{0})$
3DFlex[32] |  | -
---
-
-
| CTF
---
Rotation
Translation
Yes
|  |  |  Homogeneous $z=(\mu_{V},\sigma_{V})\longrightarrow V$ | FSTdiff[39] |  - Rotation Translation |  CTF - - | Free
|  | Fourier |  Heterogeneous $z\in\\{1,..,K\\}\longrightarrow\\{V_{1},...,V_{K}\\}$ | RELION[29] |  - Rotation Translation |  CTF - - | Free
|  |  |  Heterogeneous $z\in\\{1,..,K\\}\longrightarrow\\{V_{1},...,V_{K}\\}$ | CryoSPARC[16] |  - Rotation Translation |  CTF - - | Free
|  Neural Network | Image |  | Homogeneous
---
$z=(\mu_{V},\sigma_{V},CTF_{V})\longrightarrow V_{i}(2D)$
CryoVAEGAN[42] |  | CTF
---
2D Rotation
-
| -
---
-
Translation
Free
Continuous Field |  |  |  Heterogeneous $z\longrightarrow V_{j}$ | CryoDRGN[23] |  - Rotation Translation |  CTF - - | Free
|  | Fourier |  Homogeneous $z\longrightarrow V_{j}$ | CryoAI[18] |  - Rotation Translation |  CTF - - | Free
|  |  |  Heterogeneous $z\longrightarrow\\{c_{j},A_{j},\sigma_{j}\\}_{j\in[1,N]}$ or $\\{\Delta c_{j}\\}_{j\in[1,N]}$ ; ($N$ “blobs") | E2GMM[33] |  - - - |  CTF Rotation Translation | Free or Yes
|  Gaussian Mixture | Image |  Heterogeneous $z\longrightarrow\\{c_{j}\\}_{j\in[1,N]}$; ($N$ residues) | CryoFold[35] |  - - - |  CTF Rotation Translation | Free
|  |  |  Heterogeneous $z\longrightarrow\\{\Delta c_{j}\\}_{j\in[1,N]}$; ($N$ residues) | AtomVAE[17] |  - Rotation Translation |  CTF - - | Yes

Table 3: Classification of reconstruction methods in terms of the
parametrization of the volume $V$ (or $\tilde{V}$), i.e. the model of
conformational heterogeneity with the conformation variable $z$, use of image
versus Fourier space, representation as a discrete or continuous field, and
with a reference-free or reference-based encoding. We also indicate whether
nuisance variables are hidden or assumed to be known in the reconstruction
methods cited. When not specified, the latent variable $z$ belongs to a vector
space $\mathbb{R}^{L}$.

| Posterior $p(h_{i}|x_{i})$ | Mode $\hat{h}_{i}$ | Variational $q(h_{i})$
---|---|---|---
|  | Non-Amortized |  | Amortized
---
(encoder)
Non-Amortized |  | Amortized
---
(encoder)
Distribution on $h_{i}$ | $p(h_{i}|x_{i})$ | $\operatorname*{arg\,max}_{z}p(h_{i}|x_{i})$ | $\text{Enc}_{\xi}(x_{i})$ | $q_{\eta_{i}}\in\mathcal{Q}$ | $q_{\text{Enc}_{\xi}(x_{i})}\in\mathcal{Q}$
Name | EM | Approx. EM | Am. EM / AE | Var. EM / AD | Var. Am. EM / VAE
|  | RELION[13]
---
(pose $R_{i},t_{i}$ and
conformation $z_{i}$)
| RELION[13]
---
(noise $\sigma_{i}$) (*)
| CryoPoseNet [34]
---
(rotation $R_{i}$)
| 3DFlex[32]
---
(conformation $z_{i}$)
| CryoDRGN[23]
---
(conformation $z_{i}$)
|  |  | CryoSPARC[16]
---
(rotation $R_{i}$) (*)
| E2GMM[33]
---
(conformation $z_{i}$)
| FSTdiff[39]
---
(pose $R_{i},t_{i}$ and
conformation $z$)
| CryoVAEGAN[42]
---
(2D rotation $R_{i}$
and $\text{CTF}_{i}$)
|  |  | CryoDRGN[23]
---
(pose $R_{i},t_{i}$) (*)
| CryoAI[18]
---
(pose $R_{i},t_{i}$)
|  | atomVAE[17]
---
(pose $R_{i},t_{i}$ and
conformation $z_{i}$)
Inference on $h_{i}$ (E-step) |  |  | 3DFlex[32]
---
(conformation $z_{i}$)
|  |  | CryoFold[35]
---
(conformation $z_{i}$)

Table 4: Classification of cryo-EM reconstruction algorithms with respect to
the inference on the conformation $z_{i}$ and on the unknown hidden variables
$h_{i}$. The inference step updates the posterior $p(h_{i}|x_{i}$, or its mode
$\hat{h}_{i}$, or a variational approximation $q(h_{i})$ of it — which