Methods for Cryo-EM Single Particle Reconstruction of Macromolecules having Continuous Heterogeneity

Macromolecules change their shape (conformation) in the process of carrying out their functions. The imaging by cryo-electron microscopy of rapidly-frozen, individual copies of macromolecules (single particles) is a powerful and general approach to understanding the motions and energy landscapes of macromolecules. Widely-used computational methods already allow the recovery of a few distinct conformations from heterogeneous single-particle samples, but the treatment of complex forms of heterogeneity such as the continuum of possible transitory states and flexible regions remains largely an open problem. In recent years there has been a surge of new approaches for treating the more general problem of continuous heterogeneity. This paper surveys the current state of the art in this area.


Introduction
Over the past few years, the combination of cryoelectron microscopy (cryo-EM) imaging and singleparticle analysis has been established as the method of choice for studying the structure of large protein complexes at atomic or near-atomic resolution [1]. Its recent success has been enabled by advances in detector technology, sample preparation techniques and the availability of advanced image processing software packages. The two other major structuredeterminations methods are X-ray crystallography, which requires the sample to be crystallized, and NMR, which is useful only with relatively small proteins.
Cryo-EM single-particle analysis (we will denote this simply as cryo-EM) involves the imaging of individual copies (called particles) of a macromolecular structure. Through the computational processing of 10 4 to 10 6 of such particle images, 3D density maps can be obtained through single-particle reconstruction. To density maps of sufficient resolution, atomic structures can then be fitted, with ∼ 4Åbeing the worst resolution for successful ab-initio fitting. Because the macromolecules are suspended in solution before they are rapidly frozen, particle images are likely to reflect a more native conformation, but also may contain frozen instances of flexibility or conformational variation. Therefore, one of the promises of cryo-EM is that researchers will be able to construct a complete picture of all the possible conformations of the imaged structures.
Conformational changes are key to the function of many macromolecular machines. The molecular motors dynein and kinesin undergo cyclical changes as a chemical reaction (hydrolysis of ATP) drives a mechanical stepping motion that moves cargo along a microtubule filament. Glucose transporters allow cells to take up this nutrient through a conformational cycle that enforces the transport of one or two Na + ions with each glucose molecule. DNA replication is carried out by the replisome, a large combination of molecular machines that, in stepwise fashions, unwind the double-stranded DNA and synthesize new complementary DNA strands.
While existing methods excel at reconstructing clearly-defined discrete conformations from cryo-EM data, the problem of reconstructing continuouslyvarying conformations of a macromolecule is where most state of the art methods fall short. This is the continuous heterogeneity reconstruction problem. Fortunately, there has been a surge in the effort devoted to developing computational tools for continuous heterogeneity reconstruction, and this is the focus of this survey article. Our aim is to highlight the defining characteristics of each method and the conceptual overlaps and differences between them.
We note that there are many ideas and approaches for continuous heterogeneity and many papers intro-duce multiple new ideas and combine multiple approaches. For brevity and clarity, we cluster together different works and omit important implementation details. This paper covers the conceptual families of ideas and does not make specific recommendations about software to use. Some of the work surveyed is theoretical or less accessible to the user. Where available, we included links to some of the software that may be more accessible.

Complementary Surveys
There are a number of recent surveys covering other aspects of the cryo-EM pipeline.
A general review of the computational challenges and the main components in the analysis pipeline are available in [2]. A comprehensive description of the mathematical aspects of the problem, focusing on homogeneous reconstruction and validation, are available in [3].
Discrete heterogeneity is discussed in [4]. A survey of earlier work on continuous heterogeneity, with a comprehensive survey of normal mode analysis (Section 6) is available in [5]. In [6], the focus is on the interpretation of the energy landscape resulting from the heterogeneity analysis using likelihood-based methods. An up-to-date overview of the full cryo-EM pipeline, from sample preparation to the latest reconstruction methods, including time-resolved cryo-EM is available in [7].
The recent reviews [8] and [9] focus on machine learning approaches to cryo-EM. Specifically, the former gives an overview of machine learning algorithms used in each step in the cryo-EM pipe-line, from preprocessing and particle picking to 3D reconstruction and post-processing, while the latter is a thorough survey of deep generative modeling techniques for 3D reconstruction.
In this review, we summarize the state of the art methods for analyzing continuous heterogeneity in cryo-EM. Our aim is to sort the main families of ideas in the area and convey some of the main ideas of each technique.

Outline of the Paper
We first discuss a simplified image formation model that most cryo-EM reconstruction methods assume, as well as discrete heterogeneous reconstruction and multi-body refinement, in Section 2.
In Section 3 and Section 4, we describe manifold learning approaches to continuous heterogeneous re-construction, specifically applied to particle images and reconstructed volumes respectively. In Section 5 and Section 6, we discuss linear models based on covariance estimation and normal mode analysis respectively.
In Section 7, we discuss nonlinear models for continuous heterogeneous reconstruction which conceptually fit into the "hyper-molecules" framework, including traditional and deep learning algorithms. In Section 8, we discuss inference methods that bypass the traditional treatment of latent variables.
Finally, we conclude the article with a discussion in Section 9.

Image Formation Model and Homogeneous Reconstruction
In this section, we discuss a simplified model for image formation that is the basis of all the approaches in this survey. For simplicity, we begin this discussion with a model for the tomographic projections of a single conformation in the homogeneous case. While each specific method may be variations of this model and contain additional details, this model provides a useful baseline. Let V (r) represent the electrostatic potential of the molecule of interest in a specific conformation at location r ∈ R 3 . Throughout this survey, we will use the terms volume and density interchangeably to refer to V .
The particle images I i (i = 1, . . . , M ) are given by the following linear model of the image formation process: for i = 1, . . . , M , where R i is a 3D rotation operator corresponding to the orientation of the volume V , P is the 2D projection operator and T i is a 2D translation operator corresponding to the offset of the projected volume with respect to the center of the image. Both R i and T i are specific to each particle image I i and are unknown in advance. C i is the contrast transfer function (CTF) operator applied to the projected image; the operator can be different for each particle image, and in our simplified discussion let us assume that C i is known. For simplicity, we assume that the particle images are cropped in advance from the larger micrograph produced by the microscope.
Lastly, Gaussian noise η i is applied to each image. In reality, the noise introduced in cryo-EM images is shot noise [10], but the Gaussian assumption is commonly used in the reconstruction literature. While the details of the noise model depend on the specific reconstruction method, it is sufficient for our discussion to assume that the noise parameters are estimated in advance.
Remark 2.1 Cryo-EM algorithms make extensive use of the Fourier transform of images and density functions for efficient representation and computation. While the distinction between representations and operations in the Fourier domain and in the spatial domain is very important in implementations, it is not specific to the study of heterogeneity. For brevity, we omit the detailed discussion of Fourier vs. spatial domain implementations in different algorithms and focus on the main ideas that are specific to the study of heterogeneity.
In the case of homogeneous reconstruction, the usual approach is to solve the maximum-a-posteriori (MAP) problem given the set of particle images {I i }: where the log-posterior distribution is given by: Here, the likelihood function p(I i |V, φ i ) is Gaussian and determined by the image formation model (1), and the rotation and the translation R i , T i of each particle image are paired into one pose variable φ i whose assumed joint distribution is denoted by p(φ i ). From a Bayesian point of view, p(φ i ) and p(V ) play the role of priors (of the pose and the volume respectively), and in some methods (e.g., RELION), they are adjusted iteratively during reconstruction.

Discrete Heterogeneity
The traditional approach to the heterogeneity problem is to perform 3D classification, where each particle is assigned to one of a small number of different, optimized reference volumes (in some algorithms, a probability of assignment to each references volume is considered). This works best when the underlying macromolecule has a small number of discrete states, and these are distinguishable in the set of single-particle images.
The main assumption here is that the reference volumes obtained are a good representation of the different conformations of the macromolecule, which requires the removal of the outliers before reconstruction and that each class contains enough particle images for the reconstruction to be possible. Usually the number of reference volumes, or classes, is specified in advance, but it can change during optimization if some classes are not very populated [11].
First introduced in [12], this method consists of maximum likelihood or maximum-a-posteriori optimization for K different volumes, or K classes, {V 1 (r), . . . , V K (r)}. Then, the log posterior distribution (3) becomes: The K volumes are usually initialized as low dimensional reconstructions from K randomly drawn subsets of the particle images, and the optimization can be performed, for example, with the expectationmaximization algorithm [12,13,14] or with stochastic gradient descent or its variants [15,11]. 1

Multi-body Extension of Traditional Analysis
A step up in complexity is the multi-body refinement [16,17,18] 2 . Here, the volume is assumed to consist of a small number of rigid components that move relative to each other and are identical across the particle images. This method starts with applying standard homogeneous reconstruction techniques to obtain a consensus reconstruction of the bulk of the volume. Then, the user specifies a set of 3D masks that define regions of the volume containing individual components that move with respect to the bulk. Each of these components is treated like a rigid body. Lastly, independent homogeneous reconstruction is performed for each component using the particle images where, for each component reconstruction, the CTF-affected 2D projections of the other components have been subtracted from the particle images. The subtraction is done either before the separate components reconstruction step [16,17] or during the reconstruction process in an iterative fashion [18]. In the latter approach, the relative orientation of each component is refined for each particle image and updated at every iteration before subtraction, which leads to improved results.
One obvious disadvantage of this approach is that it relies on the user's previous knowledge of the structural domains of the molecule and, therefore, it is susceptible to human bias. In addition, multi-body refinement is limited to rigid variability and has difficulty at the interface between the moving components.

Conformation Space
Given a set of particle images, the objective of traditional homogeneous reconstruction is to recover a single volume which, according to the forward model (1) (or a variation of it), generated the particle images. When there is conformational heterogeneity that a single homogeneous reconstruction model cannot capture, homogeneous reconstruction algorithms recover some average structure (or even worse corrupted structures) Similarly, the main underlying assumption in discrete heterogeneous reconstruction is that the particles cluster around a small number of distinct dominating structures, which the reconstruction process aims to recover.
In both cases, an important step of the data preprocessing is the removal of outliers that do not fit with any of the structures that one wants to recover. Moreover, information from particles on the boundary between different "discrete" structures, which may represent rare or less likely conformations, is lost when they are assigned to one of the volumes or classified as an outlier.
The strength of cryo-EM for structure determination of biological samples is that it is able to capture the sample and its conformations without being constrained to a crystal lattice. Cryo-EM observes individual molecules rather than crystal ensemble averages. Therefore, avoiding averaging of these conformations through reconstruction of one or several discrete volumes is a key motivation for continuous heterogeneous reconstruction.
While the output of homogeneous or discrete heterogeneous reconstruction consists of a small number of discrete volumes, the output of continuous heterogeneity analysis should capture the full space of conformations of the molecule. This often takes the form of a low-dimensional manifold. The specifics vary significantly between different techniques. Throughout this survey, we will refer to this manifold as the conformation space or latent space representation. Given such a manifold and a coordinate vector on it, one should be able to recover a volume corresponding to the conformation at those coordinates. In Figure 1, we illustrate examples of volumes along a trajectory in the conformation space for the calmodulin protein.

Manifold Learning on Particle Images
One of the first ideas for characterizing continuous heterogeneity in cryo-EM was based on the observation that a continuous space of volumes could be represented mathematically as a manifold of volumes; similarly, the space of particle images, which are the projections of different molecules in different viewing directions, form a manifold of images [19]. Manifold learning techniques are usually based on organizing observations (particle images) based on a norm of the difference between every two observations. The study of cryo-EM data as a manifold of particle images presents several challenges. For example, in the simplest form of manifold learning, each variable such as the viewing direction, in-plane rotations, in-plane translations, and CTF, introduces additional dimensions to the problem, making the analysis impractical even before considering the heterogeneity. Furthermore, the technical properties of the difference between particle images lead to mathematical difficulty in identifying the observed manifold with the product SO(3) × τ of viewing directions and heterogeneity space. Finally, the high level of noise in cryo-EM make the difference between particle images noisy. A sequence of followup work introduced additional ideas to make the analysis practical [20,21,22,23,24,25] 3 by assuming known viewing directions and analyzing the manifold of conformations as viewed from each direction using the diffusion map algorithm [26]; the maps are then aligned across all different directions either by using nonlinear Laplacian spectral analysis (NLSA) [24] or a novel algorithm ESPER [25].
Experiments performed using simulated data in [25] show that the conformation manifold is not constructed well for certain projection directions where the range of the conformation variability is small. The manifold embeddings are also affected negatively by insufficient samples for each state in the range and low SNR, effects that are equally visible for both PCA and Diffusion Maps. To circumvent these issues, ESPER [25] rotates the manifolds' coordinates (eigenfunction realignment) so that each dimension of the heterogeneity variable is captured in a separate eigenbasis of the embedding. Then, a subspace partitioning procedure is performed to assign particles to each state along the trajectory representing the continuous heterogeneity and produce a 2D movie for each viewing direction. The direction of each movie is determined using optical flow and belief propagation using the method in [23], after which the individual states along the trajectory can be reconstructed using standard homogeneous reconstruction methods.
The manifold approach in this section is characterized by sidestepping the modeling of the volumes themselves during the analysis. Instead, this approach aims to build a map of conformations for each viewing direction. Once the map is available, one selects particle images from many viewing directions of roughly the same conformation and reconstructs that conformation from the particle images using traditional reconstruction algorithms. In the following sections, we will discuss approaches that rely on explicit modeling of the volumes.

Manifolds of Volumes
Many of the challenges presented in the previous section can be avoided by applying manifold learning techniques directly to volumes, which is the approach we discuss in this section. Clearly, cryo-EM experiments produce tomographic projections, not volumes. Therefore, the approach in this section of not selfcontained in most cases. One work on manifolds of volumes in [27] assumes that the volumes are given. A later method in [28] generates the volumes by performing 3D classification of reconstructed volumes from a number of bootstrapped subsets of the data.
Once the volumes are available, both methods propose a way of obtaining the low dimensional embedding corresponding to the manifold of conformations. In StructMap 4 [27], the volumes are aligned using rigid body alignment as well as flexible alignment using their normal mode representations (which will be discussed in Section 6), resulting in a dissimilarity matrix containing pair-wise distances between volumes. Then, multi-dimensional scaling is applied to a distance matrix obtained from the dissimilarity matrix. In AlphaCryo4D 5 [28], the low dimensional manifold embedding is obtained by first extracting features from the volumes using an autoencoder neural network, and then applying t-SNE [29] to the volumes and their features.
Once the low dimensional manifold of conforma-tions is obtained, one can analyze the conformational landscape, for example by performing clustering on the embedding, or high-resolution homogeneous refinement based on specific regions of the manifold. Since this approach relies on volumes that are produced by some other methods, in inherits some difficulties in generating such volumes from other methods. That being said, many of the methods in the remaining sections involve some components of dimensionality reduction or manifold learning on the latent conformation space, which are conceptually related to the methods in this section and the previous section. In particular, a different take on manifold learning of volumes, presented in [30] combines manifold learning techniques with covariance estimations that we discuss in the next section.

Principal Volumes and Linear Models for Volumes
The approaches in this section and the remaining sections incorporate generative models that describe the space of conformations directly. This is in contrast to Section 3, which maps the conformation space without creating an internal representation of the space of possible volumes, and in contrast to Section 4, where the manifold is not inherently a generative model (with some exceptions such as [30] that integrates ideas from this section and the previous section). One approach to representing the densities is based on linear combinations of volumes. In this representation, we have a reference volume, which we denote X 0 , and K principal volumes, which we denote by {X k } K k=1 . The density V m := V (m, ·) of a particular conformation is provided as a linear combination of the principal volumes, where the coefficient b m k is the weight that the k-th principal volume receives in the m-th conformation; there is one conformation (and one set of coefficients) corresponding to each particle image. The coefficients b m k are different for each conformation and can be positive or negative, and the principal volumes are often chosen to be orthogonal to each other, such that X k (x)X n (x)dx = 0 if k = n (excluding X 0 ). Once the principal volumes are obtained, the standard practice is to visualize each principal mode separately. The expression for the movie of the k-th principal mode as a function of "time" τ (negative or positive) is If the volumes in different conformations were given (and aligned), a natural optimal choice of principal volumes X k is obtained by a standard Principal Component Analysis (PCA) of the densities. This idea was first introduced in [31] and further developed in [32], where the volumes are obtained by homogeneous reconstruction using resampled subsets of the particle images.
Remarkably, the authors of [33] and [34] demonstrate that it is possible to estimate the 3D covariance matrix and principal volumes directly from the 2D particle images under mild assumptions. This surprising fact is evident from the expressions for the covariance in Fourier domain, extending the Fourier slice theorem. Specifically, in [34], appropriate esti-matorsX M and Σ M for X 0 and the covariance matrix Σ 0 are defined using least-squares optimization problems. These estimators are shown to be consistent, i.e. they converge to µ 0 and Σ 0 as the number of particle images M tends to infinity. Finding µ M and Σ M then involves solving two linear systems. The work in [35,36,37] improves the solvers, both in terms of scalability and generality. The authors of [33] propose a probabilistic PCA approach, where the principal volumes are estimated directly from the data without first computing the covariance matrix Σ 0 . This method is improved in [38] and in [39] (implemented in CryoSPARC), which allow reconstruction of volumes of higher resolution 6 .
The interpretation of heterogeneity based on the modes in (6) has been demonstrated in several examples [35,38,39]. However, as discussed in more detail in [40], the covariance method has limited applicability to high-resolution approximation of large continuous conformational variability. The problem, in a nutshell, is that high-resolution large continuous conformational variability does not behave like Equation (6). Instead, one would require a larger number of high-resolution principal volumes and to identify specific combinations of coefficients b m k in Equation (5) that represent valid conformations.
A refined and potentially more interpretable analysis is proposed in [30]. Each particle image can be best explained as a tomographic projection of a particular volume that is a specific linear combination (equation (5)) with particular coefficients b m k . This translates each particle image to an approximate volume, which can then be used for manifold learning on volumes (Section 4). The algorithm in [30] further extends the manifold-learning approach to produce refined basis volumes, called spectral volumes, that are compatible with the recovered manifold.
A key limitation of this approach to the analysis of heterogeneity is the prerequisite of known viewing directions for each particle image. Large conformational heterogeneity makes it challenging to define and compute consistent viewing directions.

Normal Modes and Predefined Spaces of Conformations
In this section, we discuss a direct representation of the positions of atoms (or pseudoatoms) in each conformation. Let r l = (x l , y l , z l ) be the position of the l-th atom (or a pseudoatom) in a molecule, and let r = (r 1 , r 2 , ...) ∈ R 3N be the concatenation of the coordinates of N atoms. The structural variability of molecules can be simulated and studied computationally by Normal Mode Analysis (NMA) [41]. NMA produces a linearized model for small perturbations around a reference location of the atoms. The best reference structure is an (already determined) atomic model of the molecule, where the way in which the protein is folded will distinguish rigid (e.g. hydrogen-bonded) from non-rigid regions. An alternative but much less powerful approach treats a low-resolution envelope of the molecule (say, as determined by cryo-EM) as an elastic structure. In this model-free NMA, larger pseudo-atoms are used for a coarse representation of the molecule.
Analogously to Section 5, the positions of atoms r n in any conformation can be approximated by a linear combination of the normal modes {∆ k } K k=1 added to the reference positions r 0 : where the coefficients b m k determine the amplitude of the perturbation in the m-th conformation and each set of coefficients {b m k } K k=1 corresponds to a particle image. The coefficients b m k are different for each conformation, and they can be positive or negative. While the normal modes {∆ k } K k=1 are determined by the physics of the molecule, the coefficients {b m k } K k=1 are chosen so that the conformations fit the cryo-EM particle images. The starting point of NMA is defining a potential function of the energy of the reference reconstruction, which is usually considered to be a low energy conformation.
There are two distinct approaches to obtain the reference structure: a model-based approach, where the reference structure is obtained from the atomic model of the molecule of interest, and a model-free NMA, where larger pseudo-atoms are used for a coarser representation of the molecule (obtained, for example, from a volume determined using cryo-EM). The quality of the normal modes varies with the quality of the model.
Then, the normal modes are the eigenvectors of the Hessian of the potential function, with the normal modes corresponding to the large eigenvalues describing large, collective motions of atoms and the modes corresponding to small eigenvalues describing more localized motions. Once the normal modes are approximated from an approximate model of the molecule, cryo-EM comes into the heterogeneity analysis to find the that best combination of normal modes to explain each particle image.
In [42], a simple potential function that models the distances between pairs of close atoms as harmonic oscillators is introduced and is shown to accurately describe collective motions of atoms in a molecule as well as more complex and computationally expensive potential that explicitly model bond lengths, bond angles and dihedral angles. A web server allowing a user to compute the normal modes of a protein based on this idea is described in [43], and the same idea enables recent approaches like HEMNMA 7 [44] to scale to full cryo-EM datasets.
In early work in the context of cryo-EM, NMA has been used for flexible fitting of high resolution structures from X-ray crystallography to low-resolution maps from cryo-EM [45] and for finding new conformations of a previously determined structure [46]. These ideas are taken further in HEMNMA [44], a method for computing the conformational space of a reference structure by fitting the normal mode coefficients b m k that best explain each observed particle image in a cryo-EM dataset. This works by iteratively refining the normal mode coefficients and pose 7 HEMNMA and DeepHEMNMA are distributed as part of the ContinuousFlex plugin for the Scipion software package https://github.com/scipion-em/ scipion-em-continuousflex. coordinates (rotation angles and translations) for each particle image.
In [47] and [48], global collective motions described by NMA are combined with local atomic displacements given by Hamiltonian Monte Carlo sampling or molecular dynamics simulations. In DeepHEM-NMA [49], a deep neural network is trained on particle images, normal mode coefficients and pose parameters obtained from HEMNMA for a subset of the data, which then outputs the corresponding normal mode coefficients and pose parameters for the remaining data. In [50], an unsupervised learning approach is used to estimate normal mode coefficients and inplane rotations, while in [51], the coefficient of one normal mode is sampled using a deep encoder neural network together with the CTF defocus and in-plane orientation of the particle.
Subsets of particle images that have similar values of coefficients are assumed to have the same conformations; such subsets can be used to reconstruct the volume directly from the images, or the conformations can be better visualized in a lower dimensional space by performing PCA on the normal mode coefficients. An important prerequisite for the algorithms is the reference X 0 and its atomic model, which is used to compute the normal modes {∆ k } K k=1 . WarpCraft [52] models continuous heterogeneity by using normal mode analysis to combine different components of the volume, which avoids the downsides of multi-body refinement of rigid components.
A comprehensive discussion of NMA is available in [41] and in the context of cryo-EM in [5].
There are other approaches based on the physics of the molecules to provide conceptually related, but very different conformation spaces. Specifically, in [53], molecular dynamics simulations are biased with cryo-EM data by sampling a posterior determined by an energy function that combines a standard molecular dynamics energy term and one that takes into account the cryo-EM data, and in [54], CryoBIFE uses a Bayesian approach to extract the free energy profile of a molecule directly from cryo-EM images together with its uncertainty. Remark 6.1 We note that in most normal mode analysis work described in the current section the space of conformations is predetermined by a fixed the normal modes expressed in Equation (7). These normal modes are not inherently based on cryo-EM data, but rather on some external predefined model of the molecule (although that model might be de-rived from some reference homogeneous model reconstructed from cryo-EM data).

Nonlinear Models: Hyper-Molecules
In Sections 5 and 6 we discussed linear models, and in Section 3 we discussed a nonlinear approach based on manifolds of two-dimensional tomographic projections. In this section we discuss nonlinear models for the volume V . The general formulation of nonlinear models in cryo-EM, laid out in [55,56], models all the conformations of the molecule in one hyper-molecule function V (τ , r), where τ is the conformation variable (position in the conformation space) and r is a point in space (or a frequency). In the simplest case, τ := τ is a scalar and the function V (τ, r) is analogous to a movie: if we fix a point τ in "time" we obtain a single conformation V (τ, ·) analogous to a single still frame in a video. In other words, V (τ, r) is a complete description of a continuum of conformations. More generally, τ can be a high dimensional vector capturing complex continuous heterogeneity. For example, when τ is two dimensional, one can imagine a planar map where each point represents a conformation.
The hyper-molecule function V (τ , r) has been implemented in many different ways in different works.
It can be argued that a continuous model can be approximated with a sufficiently large number of discrete samples, or a sufficiently large number of classes in 3D classification in cryo-EM. This would be analogous to the sequence of individual frames that represents a movie. However, the discrete 3D classes are analogous to an unordered collection of frames, which is not as interpretable as a movie. More importantly, the relation between adjacent frames in a video reduces the amount of information required to represent the movie (indeed, compressed videos are more efficient than a collection of images representing the same frames); this is loosely translates to a fewer particle images one would need in order to recover V (τ, ·) compared to number of particle images required to recover a very large number of discretized classes. We revisit the motivation for continuous functions more formally in Section 9.2.

Orthogonal Basis Functions
The hyper-molecule V (τ , ·) can be represented in different ways. The first classic harmonic analysis approach proposed in [55,56] uses linear combinations of orthogonal basis functions P k : One practical implementation of (8) is a product of three dimensional basis functions (such as prolate spheroidal functions [57,58,56]) and one dimensional basis functions. The principle was demonstrated to work using a stochastic gradient decent like algorithm on synthetic data in [55], and using a Markov Chain Monte Carlo (MCMC) algorithm and synthetic and real data in [56]. In both cases, the viewing direction and conformational state are not assumed to be known apriori. While [55,56] argue the hyper-molecules are the natural generalization of traditional 3D volumes to the case of continuous heterogeneity, it argued that it would become increasingly difficult to rigorously generalize the traditional expectation-maximization and branch-and-bound algorithms to high-dimensional latent spaces that involve viewing direction, translation and complex conformation variables. Therefore, they propose formal MCMC algorithms and speculate about possible use of variational approximations. (8) and (5), we note that the coefficients a k in (8) are a shared part of the model, whereas the coefficients b m k in (5) are different for each individual particle image.

Remark 7.1 To contrast between equations
Remark 7.2 In fact, the linear density model in covariance approach in Equation (5) is a special case of nonlinear models. For example, one can define τ to simply be the vector of coefficients τ :

CryoDRGN
CryoDRGN 8 [59,60] introduce new ideas inspired by the success of Variational AutoEncoders (VAEs) [61] in other applications. VAEs have two components, an encoder which is optimized to provide approximations of the distribution of latent variables, and a decoder, which is optimized to reproduce particle images given the latent variables as inputs. Figure 2 is a schematic illustration of the cryoDRGN architecture. Many of the neural network based approaches discussed in the remainder of this section have an analogous architecture.

Remark 7.3
The use of a deep neural network in CryoDRGN and other methods below might create a perception that CryoDRGN is somehow pretrained on some dataset with many known structures; in fact, rather than being pre-trained, CryoDRGN fits the weights of the neural networks separately for each dataset to which it is applied.
CryoDRGN's decoder is a hyper-molecule function V (τ , r), implemented as a standard multilayer perceptron (MLP), with an interesting twist. The conformation latent variable τ (typically 8 dimensional) feeds into the MLP in the traditional way. The coordinates r = (r 1 , r 2 , r 3 ) (in the Fourier-Hartley domain) first go through positional encoding where each of the three coordinates is augmented; for example, the coordinate r 1 is augmented by the vectorr 1 with: where i = 1, . . . , N/2, and N is the size of the box containing the volume of interest. The positional encoding modification to the MLP by augmenting the input vectors has been used in transformers and vision applications; it plays an important role in fitting high-frequency components of functions [62].
In order to produce a tomographic projection from a particular viewing direction and conformation τ , CryoDRGN computes a grid of points r 1,1 , r 1,2 , ... in the Fourier-Hartley domain and evaluates V (τ , r 1,1 ), V (τ , r 1,2 ), ... at these points. The grid of values is the Hartley transform of the particle image.
Previous algorithms store or compute for each particle image I i some form of explicit latent variables such as the viewing direction R i and conformation τ i . For example, MCMC software stores these latent variables and software like RELION compares each particle image to projections of different volumes in different viewing directions. Instead of storing the explicit τ i for each particle image I i , CryoDRGN's encoder is optimized to take the image I i itself as input and compute a sample from p(τ i |I i ). Anecdotal evidence points to similarities between the encoder implementation and explicit latent variables [63].
The optimization objective for the encoder and decoder networks is the standard VAE variational lower bound of the model evidence, which includes the reconstruction error as the mean squared error between the reconstructed image and the input image. We note that the encoder and decoder networks are optimized together for each dataset; neither is pretrained on other datasets.
The original CryoDRGN implementation was demonstrated with known precomputed viewing directions. Followup work [64] demonstrated a hybrid use of explicit latent variables for viewing directions and translations, and an encoder for the conformation variables. Related work on CryoFIRE [65] uses an encoder that outputs both conformation and pose variables using amortized inference.

Adding Model Constraints
The functions used to express hyper-molecules in [55,56] and the deep learning version in [59] are very generic in their ability to express a wide variety of changing volumes; including ones that are not physically plausible. In [56], it is pointed out that complex heterogeneity presents an information problem. This can be illustrated by looking at a simplified discretized case. Suppose that a good reconstruction of a volume requires 10 5 particle images. Now suppose that the molecule has a flexible region and that the flexibility can be adequately represented using 10 states, each requiring 10 5 particle images for a good reconstruction (we will see in the discussion that continuous models may require fewer particle images, but the principle is similar). For convenience, let us assume that the particle images are conveniently equally distributed among the different states and that we are able to associate each image to the correct state. Now suppose that there are 2 such flexible regions; this means that there are 100 combinations of states of the two regions, each requiring 10 5 particle images. The number of particle images grows exponentially fast with the complexity and number of regions (if not mitigated somehow). [56] proposes to enforce structure in the models that capture any of the physical properties that can be assumed about the molecule; the example presented there is explicit decomposition of the volume into regions that are allowed move independently, and it is argued that the number of coefficients describing the problem (loosely reflecting on the number of particle images) grows linearly (or even remains constant) as the number of regions grows. We note that in contrast to the multi-body approach to decomposition of the volume (Section 2.3), the components in the example in [56] are not rigid.
The following sections discuss several different approaches that enforce models on hyper-molecules.

GMM Model
A Gaussian mixture model (GMM) was introduced in [66] as a different model for expressing V (τ , r) 9 . In GMM models, the volume is represented as a sum of N Gaussians centered at spatial coordinates r j (τ ) N j=1 ⊂ R 3 , which vary continuously with the conformation variable τ : where {a j (τ )} N j=1 , {σ j (τ )} N j=1 ⊂ R + are the amplitudes and widths of the Gaussians respectively and are also dependent on the conformation variable τ .
The GMM model in [66] is based an autoencoder architecture like CryoDRGN. The decoder in this model takes the latent conformation variable τ (with a default dimension four) and produces a list of centers r j (τ ) N j=1 that can then be used to evaluate the tomographic projection for any viewing direction. The encoder computes the latent conformation variable τ for any given images (with notable differences that we omit here from the CryoDRGN encoder).
The training of the neural network is performed in two steps: the encoder is first trained to match a precomputed neutral structure corresponding to τ = 0, and then both the encoder and the decoder are trained to match the variability in the data. The second step of the training is similar to that of CryoDRGN. variable φ i .

Remark 7.4
The NMA representation using pseudoatoms with locations given in (7) is a special case of the GMM representation in (12) (but not the specific algorithm discussed in this section). We can see this by letting a j (τ ) = a 0 and σ j (τ ) = σ 0 for all j = 1, . . . , N , where a 0 and σ 0 are the fixed values of the amplitudes and widths of the pseudoatoms in (7) and by letting where r 0 j ∈ R 3 is the vector of coordinates corresponding to the j-th pseudoatom in the reference structure, ∆ k j ∈ R 3 contains the entries of the normal mode ∆ k corresponding to the j-th pseudoatom, and the normal mode coefficients b m k := b k (τ m ) are functions of the conformation variable τ m . Figure 2: CryoDRGN architecture. In order to train the encoder and the decoder networks, a particle image I i is given as input to the encoder, which outputs a latent space representation (or conformation variable) τ i of the particle image. Then, the decoder network takes τ i and the pose variable φ i as input and outputs a projected imageĨ i of the estimated volume. The loss L(I i ,Ĩ i ) between the input particle image and the output image is computed and the weights of the networks are adjusted accordingly. The GMM and GMM with folding constraints models described in Section 7.4 and Section 7.5 follow a similar workflow, with the most significant difference being that the decoder outputs the atom positions in the conformation corresponding to τ i , which are then used to generate the output particle imageĨ i .

GMM with Folding Constraints
The methods described thus far make limited use of what we know about the molecules imaged in the experiment. Biomolecules are a folded chain of amino acids, and in many cases one has information about the chain and even some approximation of the folded structure. This information can be incorporated into a more nuanced model.
AtomVAE [67] 10 models the molecule at the atomics level, expressing f (τ , r) as a GMM model where each atom (other than hydrogen) is represented by a Gaussian. The architecture of AtomVAE is also based on a VAE approach, where the encoder consists of three networks: one that takes a particle image and encodes it into a general latent vector, one that takes the general latent vector and encodes the conformational landscape, and one that take the general latent vector and encodes the pose. The decoder takes the sampled latent conformation vector and the pose as inputs and outputs the translations with respect to a neutral structure, previously obtained using homogeneous reconstruction (or purely computational methods).
CryoFold [68] models the molecule at the residue level, expressing f (τ , r) as a GMM model where each amino acid is represented by two Gaussians, one for 10 The survey in [9] refers to this method as AtomVAE.
the backbone and one for the side chain. It assumes known, fixed, viewing directions and uses a VAE architecture similar to that of CryoDRGN and outputs the translations of the Gaussians with respect to a neutral structure, as done by AtomVAE.
In both methods, structural constraints related to the positions of the atoms in the backbone and the in side chains are specified in the loss function, in addition to the error in the reconstructed images. They have been demonstrated on synthetic data and they require a good initialization with a reasonable reference structure.

Deformation Based Models
It is appropriate to model many cases of continuous heterogeneity as deformations of a reference volume; this may exclude the binding and unbinding of subunits or ligands, such as elongation factors on ribosomes. Given a reference volume V 0 (r), a hyper-molecule description of the deformed volume can loosely be formulated as where g[τ ] expresses the deformation of the reference to the conformation τ (a more nuanced technical description may also account for the Jacobian of the deformation). Deformation based models attempt to fit both the reference volume V 0 (r) and the deformations g(τ ).
Deformation functions based on Zernike polynomials have been used in [69,70], which combine ideas from manifolds of volumes described in Section 4. 3DFlex 11 [71] uses a deep neural network that produces a deformation field generator g as a function of the latent conformation variable τ . 3DFlex optimizes the parameters of the deformation generator network, the canonical volume, and the latent conformation variable of each particle image in a maximumlikelihood framework. The viewing directions and the CTF parameters are assumed to be known, but can be further refined by 3DFlex.
Follow up work on the GMM method in [66] is presented in [72], which combines physical constraints and deformation models in a more scalable version of the original GMM algorithm.

Bypassing the Estimation of Latent Variables
Let us assume that the particle images are generated according to a classic generative model, where latent variables such as conformation, viewing direction, translation (and possibly CTF parameters) are selected at random and used to generate a particle image along the lines of Equation (1). Let us further assume that the distributions of these latent variables and the noise are known. In theory, this would give us a distribution of particle images that is conditional on the hyper-molecule. It was shown in [73,74,75,76] that a given distribution of particle images can be theoretically traced back to a unique volume (up to trivial symmetries). For simplicity, we limit the description in this section to the problem of homogeneous reconstruction when there is no conformational heterogeneity in the data and we will point to the way it generalizes to the heterogeneous case without delving into details.
One of the methods to compute parameters from a distribution in many statistical problems is the method of moments, where moments of a distribution are computed and used to estimate the parameters of the distribution. Indeed, back in the 1970s, Zvi Kam found that one can compute moments of the distribution of observed images and use these to reconstruct the volume [77]. The method has been 11 3DFlex is part of the CryoSPARC software package (link in footnote 1).
further developed in recent years, for example in [75,76,78]. Remarkably, this can be applied directly to micrographs, without particle picking [79]. This approach, as well as the broader area of multi-reference alignment, has been shown to extend to the case of heterogeneity (e.g., [80,81]).
Recently, machine learning approaches have emerged for approximating distributions.
Cryo-GAN [74] is a machine learning approach that directly finds a volume that generates the distribution of observed particle images. We note that, in practice, the data consists of a large but finite number of samples from the distribution, not the distribution itself. Considering the many latent variables (viewing direction, shifts, CTF), it might be surprising that this is enough for CryoGAN to recover a volume, but the idea has been demonstrated in [74]. The method has been extended to the case of heterogeneity in [82].

Discussion and Perspectives
Having surveyed the main approaches to continuous heterogeneity analysis, in this section we discuss some of the broader in these approaches pose.

Reading the Output: Explicit Models vs. Reconstruction from Images
Given the output of a continuous heterogeneous reconstruction algorithm, there are generally two approaches to obtain the volume corresponding to a conformation at a specific value of the conformation variable τ = τ * . One approach involves selecting all the particle images {I i } that have been assigned values {τ i } that are within a certain distance from τ * and feed them to a homogeneous reconstruction algorithm. This produces one volume that represents that specific region of the conformation space. Alternatively, one can directly use the explicit model generated by the algorithm. For example, the hyper-molecule models described in Section 7 have an internal representation of each state, which allows the user to generate the volume V (τ * , ·). Similarly, equation (6) gives a way to visualize principal volumes.
Some of the methods presented in this survey rely on homogeneous reconstruction to generate volumes with different conformations, as they lack an explicit representation of the volumes (the manifold of images methods in Section 3 are an example). However, most of the other methods lend themselves to both approaches to generating heterogeneous volumes. In this case, the homogeneous reconstructions from subsets of the particle images can be used as a form of validation, as the final volumes are reconstructed directly from the data: if the volumes appear to be high-resolution and biologically plausible, users tend to accept them as a successful run of the algorithm. Such a result is especially convincing when generating volumes directly from the model would result in volumes that are relatively low-resolution due to model constraints such as a limited number of principal volumes computed in the principal volume approach or a small number of Gaussians in the GMM models. Furthermore, as models get more elaborate (e.g., atomic models), there are concerns that the results present artifacts that reflect bias and error in the model, which are somewhat mitigated if one obtains homogeneous reconstructions from subsets of the data.
That being said, visualization based on models can provide higher-resolution volumes, which reflects the advantages of continuous models over discrete reconstructions. As an illustrative example, consider a molecule with one flexible region that requires 10 separate discrete models to be fully captured in the reconstruction, for which we are given 10 6 particle images. If we restrict our attention to the rigid part of the volume, applying one of the hyper-molecule algorithms would use all 10 6 particle images for the reconstruction of this region, while selecting subsets of 10 5 particle images for homogeneous reconstruction of one of the conformations leads to only this subset of the data being used to reconstruct the rigid part, which would otherwise benefit from the full dataset. Modeling considerations and the ability to refine the viewing directions taking into account flexible regions extends this argument to other less obvious cases.

Continuous vs. Discrete Models
Continuous functions can be approximated by a sufficiently large number of discrete samples, and therefore it is compelling to deduce that traditional discrete 3D classification is a good discretization of continuous heterogeneity. In this section, we point to several advantageous features of continuous models.
First, in contrast to discrete heterogeneity analysis, continuous heterogeneity analysis does not require the user to specify the number of classes to be reconstructed. This choice, which is not trivial when the underlying heterogeneity is concentrated in discrete states and influences the resolution of the inferred model, is particularly difficult when the heterogene-ity is inherently continuous. Instead, approaches for continuous heterogeneity analysis explicitly model the continuum of states, where states that are sufficiently distinct in the particle images tend to naturally appear as dense "islands" on the conformation manifold.
Second, methods for continuous heterogeneity are able to use all the available data to produce higher resolution volumes in any conformation, while discrete heterogeneity approaches only use a subset of the data for a particular conformation. The example given in Section 9.1 illustrates this aspect: while the rigid part of a volume can benefit from the full dataset, 3D classification methods only use the particle images corresponding to one specific conformation to reconstruct the volume, including the rigid part, leading to lower resolution than possible. Here, one could correctly argue that the rigid region can be reconstructed using homogeneous reconstruction separately from the heterogeneous region, as it is done in multibody refinement (Section 2.3). However, continuous heterogeneity analysis may still perform better, as it can update the viewing directions estimations (potentially yielding better estimation of the conformation variable) and it models the continuum of states more accurately than discrete models, which do not have a notion of "neighboring states". For example, the deformation model described in Section 7.6 severely restricts the space of possible conformations to a (limited) deformation of a reference volume, and therefore all the particle images from all conformations contribute high-resolution information to the reference volume.
Third, an implication of the previous argument is that continuous heterogeneity analysis has an advantage in analyzing rare conformations using the continuity between conformations. Using the deformation model as an example, when the deformation assumption is applicable, this model constructs a high resolution reference volume based on the highly populated conformations and uses the particle images from a rare conformation effectively to determine the path of deformation of the reference volume. This leads to a high resolution approximation of the rare conformation as the correct deformation of the reference volume.
One last argument regards the advantage of continuous heterogeneity analysis over multi-body analysis (Section 2.3). While it has been implemented as an extension of traditional reconstruction algorithms, multi-body analysis can also be seen as a special case of the hyper-molecule model. Specifically, is has proven to be very effective in the analysis of rigid com-ponents moving with respect to each other, but it falls short when such structural assumptions do not hold and it fails to resolve the interface regions between the rigid components, issues which are addressed in the hyper-molecule models for examples.
Finally, it is worth noting that high resolution homogeneous reconstruction and traditional 3D classification as well as multi-body analysis are available in mature software and are well understood by practitioners. On the other hand, practitioners still face the dilemma of which continuous heterogeneity software, model assumptions and parameters to use in continuous heterogeneity analysis, so it is not the aim of this article to argue that the problem has been solved or that higher resolution is always easily attainable, but to highlight conceptual advantages on which these algorithms increasingly capitalize.

Outliers, Junk, and Rare Conformations
Practitioners have long used 2D and 3D classification not only to identify actual conformations, but also to identify "junk": particles that produce low quality volumes and objects that were incorrectly picked as particles are discarded. One of the caveats in this procedure is that it is difficult to determine whether high quality particle images are also discarded. Perhaps most importantly, it is plausible that particle images from rare conformations are also discarded. Preliminary informal evidence suggests that an analogous procedure is applicable to different algorithms for continuous heterogeneity analysis: junk particles appear to be "pushed" out to lower density regions of the conformation space, reducing their influence on the "valid" regions. It is therefore conceivable that some of the continuous heterogeneity analysis algorithms strike a better balance between using rareconformation particle images and discarding junk.

Interpretation of Output
The output of continuous heterogeneity analysis varies considerably between different methods. However, there are two main families of outputs. One is a model of the space of molecular structures (e.g., the principal volume expansion of V m in (5), the hyper-molecule function V (τ , r) in (8)), which does not exist in the image manifold approach in Section 3. The other family of outputs consists of the estimates of the conformation variable τ for each particle image, which does not exist in the methods in Section 8.
From both types of outputs, one can obtain a low dimensional representation of the latent conformation space of τ , where it is common to try to identify physically meaningful clusters or manifolds. For example, in [20], the density of the particles in the latent space is related to the energy landscape of the ribosome. In other work, various dimensionality reduction and clustering algorithms are applied to the latent space to identify dominating conformations. We note, however, that the latent conformation variable does not immediately correspond to a physically meaningful and interpretable quantity, it is not necessarily unique even up to trivial symmetries and can in principle be arbitrarily deformed (for some motivating examples, see [83]). Different algorithms and even the same algorithm with a different random seed can produce different latent space representations. Similarly, different post-processing dimensionality reduction and clustering algorithms (and again, even different random seeds) can produce different results. In cases of relatively simple heterogeneity, it is plausible that an expert user would be able to probe structures in the latent space together with the models or reconstructions associated with them in order to understand the mechanisms that it captures. However, complex heterogeneity and more automated tools may require more work on the interpretation of the latent space.
Mathematically, this is a problem related to the manifold metric as transformed by the measurement and reconstruction process, and while it has been studied in other fields [84,85,86,87], it still requires further investigation in a cryo-EM context. Moreover, establishing an appropriate metric for the conformation space that is physically meaningful is a problem that has not been explored yet.

Validation
Validation has always been a nuanced problem in cryo-EM. A common sanity check for homogeneous reconstruction is the Fourier shell correlation (FSC) procedure [88,89]: the data is split into two subsets and two distinct reconstructions are obtained from them. The correlations between the two volumes are computed for each frequency radius and the resolution of the final reconstruction is given by the frequency radius where the FSC curve drops below a specific value (usually 0.143). Above this resolution, the reconstruction is considered to overfit the noise. In addition, a qualitative procedure is for an expert to evaluate and confirm the plausibility of the reconstructed molecule.
There is no clear solution for validating results of continuous heterogeneity analysis. To illustrate the difficulty, consider two hyper-molecules that are generated by one method from two separate subsets of the data. Even if the solutions are equivalent, they can assign different values of τ to the same conformation, in which case it is not clear how to compare the two outputs. One idea proposed in [71] for deformation models (Section 7.6) is to compute the FSC between two reference models that the algorithm computes for two different subsets of the data. However, this idea only applies to deformation models, since other methods do not necessarily compute a reference volume. Given the wide variety of families of ideas for analysis of continuous heterogeneity and the different types of output they produce, it is even more challenging to find metrics that are applicable to all methods. Beyond the evaluation of outputs for individual problems, there is a desire to determine which tools perform best, which reinforces the need for evaluation metrics and benchmark datasets, something that is not currently available. Having said this, given the many forms of heterogeneity in applications and the fundamental differences between the algorithms for heterogeneity analysis, one should be cautious that the metrics of choice and the benchmark datasets used do not exclude some of the existing and new ideas.

Conclusions
In this survey, we covered a broad range of methods for analyzing continuous heterogeneity in cryo-EM. In addition, we highlighted the advantages of continuous heterogeneity analysis and the way they circumvent some of the tradeoffs that practitioners encounter in traditional methods involving homogeneous reconstruction, discrete heterogeneous reconstruction and multi-body analysis.
As evident in this survey, many different approaches have been proposed in recent years. Most of these ideas are still in active development and it is likely that new ideas will emerge in the coming years. Some of the algorithms, like the various members of the non-linear models in Section 7, share conceptual similarities that make it plausible that we will see convergent software packages that offer multiple and customizable models, as envisioned in [56]. At the same time, there are significant differences between the methods, making different approaches well-suited for different problems.
One of the open problems in this area is the validation of the results of the algorithms. A related difficulty is in comparing different algorithms in the absence of good metrics and benchmarks; indeed, the different approaches differ even in the form of output they produce. A silver lining in the large number of fundamentally different algorithms that are becoming available is that one could conceivably apply several different algorithms to the same dataset and examine how well their conclusions agree qualitatively, which could serve as a temporary form of validation while the area matures.
Finally, we have we witnessed in recent years the immense success of protein structure prediction algorithms such as AlphaFold [90] and RoseTTAfold [91]. While these algorithms are not currently applicable to continuous heterogeneity (and, arguably, neither to discrete heterogeneity), it is plausible that they will evolve in this direction as a next step in complexity. Conversely, as cryo-EM datasets become more massive and broader in scope, incorporating physical models as priors or constraints obtained from structure prediction software will become a necessity. Therefore, there is an opportunity for continuous heterogeneity analysis methods to bridge the gap between such software and the massive experimental datasets. Preliminary examples of this general direction can be found in [67], which integrates insights from AlphaFold, and in [54].
While we do not discuss cryo-electron tomography (cryo-ET) explicitly in this survey, many of the approaches to continuous heterogeneity in cryo-EM can be generalized to cryo-ET. One of the complications in cryo-ET is that the targets of the analysis are typically measured in situ, and do not float freely as in typical cryo-EM. The continuous heterogeneity model may capture heterogeneity in the environment instead of heterogeneity in the target (which may or may not be a desirable result, depending on the application).