1 Introduction
CryoEM has been named Method of the Year 2015 by the journal Nature Methods due to the breakthroughs that the method facilitated in mapping the structure of molecules that are difficult to crystallize. CryoEM does not require crystallization necessary for Xray crystallography, and unlike NMR it is not limited to small molecules. One of the additional great opportunities in cryoEM, which has been noted, for example, in the surveys accompanying the Nature Methods announcement [1, 2, 3], is to overcome heterogeneity in the sample: in practice many samples contain two (or more) distinct types of molecules (or different conformations of the same molecule); methods like Xray crystallography and NMR, which measure ensembles of particles, have a difficulty distinguishing between these different types. Many existing algorithms for the analysis of cryoEM experimental data address the problem of recovering a finite number of distinct structures in heterogeneous samples. However, ”[large macromolecular assemblies] tend to be flexible, and although classification methods have come a long way when applied to biochemical mixtures or welldefined conformational states, the continuous motions seen for certain samples will challenge classification schemes and set a limit to achievable resolution, at least for a number of years.”[2] The purpose of this manuscript is to propose an approach to mapping continuously heterogeneous structures, and to demonstrate its applicability to cryoEM.
The idea behind our approach is to generalize the tomography and cryoEM problems of recovering a function over , representing a 3D object, to the problem of recovering of a function over , where captures the topology of the heterogeneity. We refer to these more general functions as hypervolumes, hyperobjects or hypermolecules. For example, loosely speaking, some types of continuous heterogeneity can be represented as functions over ; tomography and cryoEM can be generalized to the problem of recovering such 4D objects. Indeed, the concept of 4D tomography (4DCT) has been used in CT scans of patients, whose bodies change periodically with time as they breath (see, for example, [4]), but it is restricted to the case where both the orientation of each projection and its phase (time) in the breathing cycle are known. We discuss general topologies of heterogeneity, the case of unknown directions and “time,” and new approaches to the representation and regularization of the problem.
To illustrate the idea, we consider the heterogeneous object whose level sets are depicted in Figure 1. This object has a one dimensional heterogeneity variable in the range , so for each value of in this range, there is a 3D object that is different from the 3D object at a different value of . However, the objects are continuous, in the sense that the object at is very similar to the object at if the difference
is small. This idea is somewhat similar to a short 2D video where the frame at each moment is generally relatively similar to adjacent frames. In movies, this similarity allows us, for example, to compress the video, because the video can be represented using much less information than a collection of independent frames.
Going back to our heterogeneous object, suppose that we have a 3D movie of this object as we change the value of the parameter in time (examples available at http://roy.lederman.name/cryoem), and suppose that we focus on a single “voxel” or coordinate in 3D space as time passes and measure the mass at this point at every moment. Suppose that at none of the mass of the object is positioned at our point in space, but as time passes, one of the “ears” passes through our point and than exits before
. Then, we would see a curve, to which we could fit some polynomial; we could fit this polynomial even if we had only sample points in time, and we could use it to interpolate the value of our voxel at times when we did not sample the voxel. The mass at our voxel can be represented using the expansion
(1) 
Looking at all the voxels in our 3D movie, we can fit such curves to each voxel, and obtain a representation of our object (discretized in space); for every sample point in space, we have the polynomials with the coefficients of that particular point:
(2) 
A very high degree polynomial would allow us to capture a very high frequency function, but the expansion is typically truncated (or the high frequency are weighted down) to reflect smoothness in the model and also due to the limited number of samples that are available, and for practical computational reasons. Now consider how a 3D object is typically represented when there is no heterogeneity; typically the representation is equivalent to some linear combination of coefficients so that
(3) 
with some spatial basis functions. In cryoEM, the objects are often represented as samples at gridpoints in the Fourier domain, so that the index encapsulates the three indices of points in the Fourier domain. The expansion is truncated to a finite number of frequency samples, and often some frequencies are penalized so that they would have lower amplitude; this reflects a preference to smoother objects. Combining our discussion of the object and its variability over “time” we obtain a generalization of a 3D object to a “3D movie,” or hyperobject:
(4) 
One of the main difficulties in cryoEM is that the 2D images of the object are given without information about the directing from which each image was taken, and the challenge is to reconstruct the 3D object from these images without knowledge of these direction; various algorithms have been constructed for this difficult task. In the case of heterogeneous objects, the challenge is to reconstruct hyperobjects from images in the absence of both the direction of each image and the value that reflects the “time in the 3D movie” or the version of the object which has recorded in that image. Various algorithms have been constructed for the case of discrete heterogeneity, where there are several distinct objects (a collection of independent “3D still scenes” as opposed to a continuous “3D movie”). Continuously heterogeneous objects, such as the object in Figure 1, are often treated as if they were distinct independent objects; this does not capture the continuous nature of the model and does not take advantage of this property to improve the reconstruction. Very loosely speaking, this would be analogous to averaging a movie in 1 second windows, then presenting these onesecond averaged windows in random order.
We submit that 1) capturing the continuous nature of hyperobjects contributes to the reconstruction and understanding of the objects, 2) the problem of reconstructing without knowing the direction and is analogous to the problem of reconstructing without knowing the directions (we discuss some of the differences that do indeed arise), and 3) our general approach can be applied in a wide range of reconstruction algorithms, based on this analogy. Indeed, the reconstruction of distinct objects (discrete heterogeneity) can be viewed as a special case of our approach.
We propose several directions for implementing the idea of using continuity or approximate continuity, and discuss in more detail one such approach based on Equation (4). We argue that different variations are useful for different models (for example, spacial bases for local variability in a larger object). Furthermore, we argue that the approach can be applied to rather arbitrary topologies of heterogeneity: the line segment in our example, the cyclic pattern in 4DCT, two dimensional surfaces etc. We demonstrate our approach in a basic prototype as proofofconcept, but emphasize that these ideas can be implemented in any of the many approaches applied to problems like cryoEM. The results of the reconstruction using our prototype are demonstrated in Figure 2.
This manuscript is organized as follows. In Section 2, we summarize some standard results used in this manuscript. In Section 3, we briefly review a simplified model of cryoEM and reformulate some of the tools used in cryoEM algorithms in a way that we find useful for generalization. In Section 4, we discuss the generalization of tomography and cryoEM to heterogeneous objects. In Section 5, we present a crude algorithm which we implemented to investigate one version of the ideas presented in this manuscript. For the sake of completeness we briefly review some of the new ideas used in this implementation, but argue that this is certainly not the only way to implement the idea of hyperobject reconstruction. Our preliminary results are presented in Section 6. A Brief summary and discussion of future work is presented in Section 7. Brief conclusions are presented in Section 8.
Additional examples and video visualization of the results are available at http://roy.lederman.name/cryoem.
Remark 1 (Terminology: “representation”).
Our use of the term “representation” in the context of this manuscript is very different from the context in which we use the term in [5]. However, we have not found a better term that would avoid this confusion. In this manuscript “representation” is a way of expressing a function or a problem, typically an expansion of a function in some basis, whereas in [5] it is a technical representation theory term. These two works are independent and largely unrelated on a technical level; the conceptual relation between the two is the motivation to treat heterogeneity as “just another variable.”
2 Preliminaries
the Fourier transform of the function in the variable 

the Fourier transform of in spacial coordinates  
the vector rotated by 

the function rotated by , so that 
2.1 Spherical harmonics and rotations
In this subsection we summarize some of the properties of spherical harmonics.
The normalized spherical harmonic, denoted by , with integer and , is defined by the formula:
(5) 
where are the associated Legendre polynomials (see, for example [6]). The spherical harmonics are an orthonormal basis of , so that
(6) 
and the expansion of any function on the sphere in this basis is
(7) 
with the appropriate expansion coefficients .
Any arbitrary rotation of a 2D sphere can be represented by three Euler angles, we denote the rotation operator by . The expansion of a rotated function is given by the formula
(8) 
where the coefficients in the expansion are
(9) 
with the th order of the appropriate form of the WignerD matrix for the rotation (see, for example [7]).
The restriction of a function on a sphere to the “equator” is given by
(10) 
It immediately follows that the Fourier expansion of the restricted function is
(11) 
where
(12) 
2.2 Haar basis
The Haar wavelet mother function is defined by the formula
(13) 
The Haar function is defined by the formula
(14) 
for all integer . For the purpose of our discussion of a finite intervals, it is convenient to take and
and add the constant function. Higher dimensional Haar basis functions are composed as the tensor product of one dimensional Haar functions.
A truncated expansion of a function in the Haar basis is simply a piecewise constant function. In other words, the space is divided into intervals, and the function is constant in each of these intervals. However, the Haar basis can also be thought of as a multiscale/tree decomposition of an interval (or hypercube), into subintervals. The multiscale tree structure of the functions is useful in analysis.
3 Setup and reformulation
3.1 The representation of objects
An object or a “volume” is a function over (in our case, is coordinates in space). The Fourier transform of the object, which we denote by , is a function over (in our case are coordinates in the Fourier domain). For the purpose of this discussion, we assume that the objects are continuous with respect to and (this assumption can be relaxed).
In applications, the function must be discretized in some way, if only so that we can represent it in a digital computer. Furthermore, the choice of discretization reflects (implicitly or explicitly) priors or regularization of the object, which are required in order to make the tomography problem tractable.
The object is often represented using discrete samples of on a regular grid, restricted to some finite box which is sufficiently large to contain the molecule which we wish to reconstruct (e.g., a grid); points off the grid are sometimes evaluated using some continuous interpolation from the grid points (although noncontinuous voxels are also considered). Since many of operations in cryoEM are naturally represented in the Fourier domain, and as it is often assumed that (or some lowresolution version of it) is bandlimited, many cryoEM algorithms and programs use a similar regular grid in the Fourier domain to represent samples of in . In works such as [8], the Fourier transform of the object is represented as concentric shells sampled at various radii.
These different representations assume certain properties of ; the regular grid in a box in the real space assumes that the is compactly supported in real space (the choice of interpolation scheme implies additional smoothness assumptions). The regular grid in the Fourier domain implies, for example, that the function is bandlimited. The representation in concentric shells in the Fourier domain makes similar assumptions on the bandlimit of the function and is computationally convenient; each concentric shell can be represented using spherical harmonics rather than samples on the sphere, yielding a natural continuous representation that has no directional bias (unlike the regular grid).
Loosely speaking, different common representations of an object can be summarized as vectors of coefficients. For example, a vector composed of all the values sampled at the grid points of a 3D Fourier transform of an object. In standard linear representations of objects, the object (or its Fourier transform) is written as
(15) 
or
(16) 
where is the vector of coefficients, the th element of the vector, some basis functions appropriate for that representation, and their Fourier transforms.
Typically, the representation of objects is continuous, in the sense that small rotations and translations of the object, and small changes in the object, result in small changes to the coefficients in the representation. This property is not necessarily uniform: for example, some types of perturbation in the objects can have a larger effect on high frequency components than on low frequency components of some representations.
The choice of basis and the truncation of the basis to a finite number of basis functions (e.g., finite cube in space and finite resolution) restricts the functions that can be represented to certain function spaces. This serves as an implicit regularization in algorithms for the reconstruction of molecules (see below), and it is sometimes augmented by more explicit regularization or priors on the object.
Bases that are better suited to a problem allow high accuracy approximation of the relevant functions using fewer coefficients, or more efficient computation. Loosely speaking, a function that is efficiently approximated in one basis can be approximated in other reasonable bases even if the representation is less efficient and requires more basis functions. For example, a low order polynomial can be represented efficiently as a linear combination of polynomials basis functions, and it can also be approximated well by equally spaced samples with linear interpolation between the sample points. In this sense, all the bases are “equivalent.” However, an efficient basis is also an implicit regularizer which restricts the space of functions. For example, low order polynomials are smooth and have restricted oscillations, whereas the sampling scheme above has more degrees of freedom and represents functions that are not smooth and functions that oscillate more. Such restriction can be introduced in any arbitrary representation explicitly as constraints, filters, penalties or regularizers.
3.2 CryoEM
Electron Microscopy is an important tool for recovering the 3D structure of molecules. Of particular interest in the context of this manuscript is Single Particle Reconstruction (SPR), and more specifically, cryoEM, where multiple noisy 2D projections, ideally of identical particles in different orientations, are used in order to recover the 3D structure. The following formula is a simplified noiseless imaging model of SPR, for obtaining the noiseless image from a object (representing the molecule’s density):
(17) 
where and is the rotation that determines the orientation of the molecule. In other words, the model is that the molecule is rotated in a random direction, and the recorded image is the topview projection of the rotated molecule, integrating out the or axis. In the Fourier domain, the Fourier transform of an image is a slice of the Fourier transform through the origin:
(18) 
where and is the random rotation.
One of the characteristic properties of cryoEM SPR that sets it apart from other tomography techniques is that the orientation of the molecule in each image is unknown in cryoEM, whereas in other tomography techniques the rotation angles are typically recorded with the measurements.
The analysis of cryoEM images is further complicated by many additional effects, some of the most notable are:

extremely high levels of noise, far exceeding the signal in magnitude (see sample images in Figure 3),

filters associated with the imaging process (CTF),

an unknown shift of each image in the image plane,

and the discretization of the measurements.
More detailed discussions of these challenges, and various other challenges can be found, for example, in [9]. Since these aspects are studied extensively in other works, and since the goal of this work is to introduce a new generalization of previous works, we will not discuss these aspects in much detail when they do not raise too many new issues that are of particular interest in continuous heterogeneity at the level of this preliminary discussion. We note that in our experiments we use simulated data with high levels of noise.
3.3 The tomography problem in cryoEM
Suppose that we are given a set of images with the orientation of each image. Then, the reconstruction problem in cryoEM is a classic tomography problem. In the Fourier domain, the model in equation (18) suggests that the problem is to reconstruct from (very noisy) slices of the Fourier transform of the object. This problem is illposed even in the noiseless case since we have no information about the values of at points between the slices. However, additional priors and assumptions on the properties of the object make the problem more tractable.
For the purpose of our discussion of tomography, we propose a further simplified description of the tomography problem, a one frequency pixel model of tomography, where we receive scalar samples of , each at a given frequency point :
(19) 
It is easy to see that if the Fourier transform of each image in (18) is considered to be a (finite) collection of samples of at various frequency points, then the one frequency pixel model (19) describes the same model as (18), with the latter grouping together samples into images.
Taking into account the discussion of representing an object as a vector of coefficients of some arbitrary linear expansion (15), we propose a more general linear operator formulation of linear models such as (17), (18) or (19), which simply states that each (noiseless) measurement is some vector of coefficients given by a linear operation of the vector of coefficients :
(20) 
In the cryoEM imaging model, the linear operator captures the operation of rotating a molecule, projecting it and representing the resulting noiseless image using the vector of coefficients (more comprehensive practical models also take into account the CTF, shifts, discretization, noise, etc.).
3.4 Unknown orientations, and reconstruction algorithms in cryoEM SPR
One of the characteristics of the cryoEM SPR problem that sets it apart from other tomography problems, is that the orientations of each image is unknown, making the problem severely illposed.
For the purpose of this discussion, we make a simplistic distinction between two approaches to the cryoEM reconstruction problem: objectfree image alignment, and direct object estimation. The objectfree image alignment approach typically takes advantage of the projectionslice theorem to estimate the rotation
of each image with respect to others based on the intersection of slices in the Fourier domain (see, for example, [11, 12, 13, 14]). Once the rotations are recovered, the next step is the tomography problem of estimating the object given the rotations. In the current manuscript, we focus on methods of direct object estimation, often called refinement, which alternate between refining an estimate of the object, and estimating the rotation of each image with respect to the object: given an estimated object at step , estimate the most likely orientation of the image recorded in the th image, for all images:(21) 
where is the operator which produces an image from the object , rotated to orientation , and is the th recorded image. In practice, this minimization is typically performed by comparing each image to a set of computed template projections of the estimated object at some sampled values of . Given these estimated rotations, produce a new estimate of the object, for example by solving:
(22) 
possibly with some additional regularization terms. Some popular refinement software packages and algorithms follow more elaborate statistical approaches (e.g., [15, 16, 17, 18, 19, 20]
) such as the Bayesian approach, and other optimization schemes (e.g., stochastic gradient descent in
[21]). Since these approaches use similar fundamental operators and due to the limited scope of this manuscript, we will restrict our attention to conventional refinement algorithms and argue that our approach generalizes such algorithms to hypermolecules. Our approach can be used in a similar way to generalize more elaborate algorithms.3.5 Frequency marching
Several refinement tools such as RELION[18, 19] and FREALIGN [17, 22, 23] gradually increase the resolution of the estimated object as they iterate. In [8], this concept is reformulated as frequency marching: starting the iterative process with a representation of the object that uses only a small number of low frequency basis functions, and adding higher frequency basis functions to the expansion at subsequent refinement iterations.
Frequency marching highlights the fact that use of basis functions is not merely an technical implementation detail of how to store a function digitally, but rather a key tool at the heart of the algorithm.
3.6 Heterogeneity
So far, we have assumed that all the molecules being imaged in an experiment are identical copies of each other, so that all the images are projections of identical copies of the object , from different directions. However, in practice, the molecules in a given sample may differ from one another for various reasons. For example, the sample may contain several types of different molecules due to some contamination or feature of the experiment. Alternatively, the molecules which are studied may have several different conformations or states, or some local variability. The heterogeneity may be discrete (e.g., in the case of distinct different molecules), continuous (in the case of molecules with continuous variability), or a mixture of continuous elements and discrete elements.
In the heterogeneous settings, the simplified noiseless imaging model of the homogeneous case (see (17)) is generalized to
(23) 
where is the objects imaged in sample . The generalized operator formulation (20) is
(24) 
where is the vector of coefficients in the representation of the object .
Most existing algorithms and software tools treat only the case of discrete heterogeneity. In conventional refinement algorithms, discrete classes of molecules are reconstructed by modeling multiple independent objects (see, for example [24]); given multiple estimated objects , we estimate the most likely pair of class and orientation for each images, and proceed to produce a new estimate of each objects based on the the images assigned to it. For example, the generalization of (21) and (22) is
(25) 
and
(26) 
A similar generalization is used in more elaborate algorithms (e.g., [15, 18]) and in software packages such as RELION[19]. Other approaches to heterogeneity require some method of recovering the rotation of the images although the images reflect mixtures of projections of different molecules (e.g., [25, 26]). We treat heterogeneity in objectfree algorithms in [5] (with generalization to continuous heterogeneity  in preparation). Another recent independent work on objectfree reconstruction [27] proposes to iterate between estimating the orientations and estimating the class labels based on pairwise relations between images. A different perspective on objectfree algorithms for continuous heterogeneity, proposed in [28, 29], is to construct a manifold of images and study the low dimensional structures induced by rotations and heterogeneity.
Currently, the prevailing approach to continuous heterogeneity is to treat the object roughly as if it were a collection of discrete independent objects, rather than a continuum of objects. In this manuscript manuscript we investigate whether we can capture the continuous nature of the states, and even take advantage of it in reconstruction.
4 Analytical apparatus
4.1 Heterogeneous molecules  “hypermolecules”
The purpose of this section is to generalize the definition of an object to a hyperobject, which represents a heterogeneous set of objects, simply by adding a variable that identifies each object instance.
To generalize the definition of an object to an hyperobject, we define a family of objects as a function over , where is the index set or parameterization of the family of objects. The evaluation of the density of an object at the coordinates , is analogous to the evaluation of a hyperobject at the coordinates , which means choosing the object instance with the index and evaluating the density of that object at the coordinates . To illustrate the notation in the discrete case, suppose that we have only two distinct objects, then , and the object instances evaluated at are and . In the continuous case, suppose that we parameterize the hyperobjects by , then the object at evaluated at is . We note that, in general, can have various topologies (discrete, an interval, a subset of a multi dimensional space, a torus, combinations thereof, etc.). For example, in 4DCT, the topology of heterogeneity may capture the cyclic behavior of breathing (as opposed to an interval with independent ends); in cryoEM it may capture a one dimensional variability in states, or, for example, independent local variability in two different areas of the molecule.
We denote by the Fourier transform which integrate over the coordinate space , and does not interact with the parameterization , so the Fourier transform of the object indexed by , and evaluate at the frequency is .
Finally, we define a transform on the
parameterization which is analogous to .
For example, suppose that ,
and suppose that is the discrete Fourier transform.
Then, we have the “parameterwise” discrete Fourier transform of the object
evaluated at parameter frequency and coordinates ,
.
In this case, yields the average of all objects:
.
We may also apply both transforms: .
For the purpose of this discussion, we assume that the operators commute:
.
4.2 Tomography of hyperobjects
The caricature one frequency pixel model of tomography (19) describes tomography as the inverse problem of recovering the Fourier transform of a object from many noisy samples at known frequency coordinates (with some constraints or regularization to make the inverse problem tractable).
This description of the tomography problem generalizes naturally to hyperobjects: the one frequency pixel model of hyperobject tomography is the inverse problem of recovering from many noisy samples at known frequencies and values of the heterogeneity index , with some analogous constraints or regularization of the hyperobject. In the discrete case, this simply means that for each sample, we have the label of the the object which was measured, and the frequency coordinates of the sample , and we may collect all the samples of each object and proceed to process each object independently from the others.
In computed tomography (CT), 4D tomography technology[4] (4DCT) has been proposed in order to image the body of a patient in different states of the breathing cycle, so that is time or phase within a breathing cycle. In 4DCT, as in classic CT, the orientation of each image is known; in addition, the position of each image in the breathing cycle is recorded using external means. The original 4DCT reconstruction algorithms binned the images according to their position in the cycle and then reconstructed an independent objects from each discrete bin. In recent years, regularization techniques have been introduced in order to take advantage of the relation between the volumes in different phases of the cycle (see, for example, [30, 31, 32, 33]).
4.3 The representation of heterogeneous objects  the discrete case
When there are K classes of objects, they can be represented as K vectors of coefficients. Suppose that is the vector of coefficients associated with the representation of the th object, with the th element of this vector. Obviously, we can rewrite these as columns of a matrix , such that , and the th object as
(27) 
Clearly, is equivalent to used in (23), so that in the discrete case this formulation of the problem is equivalent to the classic formulation used in existing cryoEM algorithms which store multiple independent models of objects.
As a step toward the discussion below of representations of continuously heterogeneous objects, we present some nonformal motivation through the following trivial accounting exercise. Suppose that is a unitary matrix, such as the Discrete Fourier Transform matrix (this requirement can be relaxed in various ways with minor modifications to the discussion). We define the matrix , where each column is a linear combination of columns in . The th objects is now written as
(28) 
Clearly, we can choose the first column of to be constant, obtaining the the average of all molecules
(29) 
Generally speaking, we do not have a preference for the orientation of an object that we represent (with some technical exceptions). Similarly, when we represent several classes of objects, we do not have any preference for the relative orientation of one class with respect to the others. However, the trivial example above demonstrates that when we have multiple classes of objects that are somewhat similar, it may make sense to place them in similar orientations so that the average of all objects bears some resemblance to some low resolution version of the objects. The first column of can then be thought of as representing the average object, and subsequent columns can be thought of as a coefficients representing the variations from the average (e.g., SVD/PCA of the hyperobject).
Next, we consider the following accounting exercise. Suppose that we is of dimensionality ( coefficients in the representation of each of the objects). Suppose that instead of arranging the coefficients in a matrix, we concatenate them into the vector :
(30) 
with and . Then
(31) 
where
(32) 
In other words, by introducing the basis functions , we represent the multiple objects as hyperobject, with an expansion that is analogous to Equation (27).
4.4 The representation of heterogeneous objects  the continuous case
Many existing software packages represent objects as a discrete set of samples points in a 3D grid, the natural generalization of this representation is a higher dimensional grid, where some axes capture the heterogeneity. This is very similar to treating continuously heterogeneous molecules as if they were a collection of independent discrete states. The idea of this manuscript is to take advantage of relation between these states, some (but not all) of the approaches that we propose to implement this idea use implicit regularization through basis functions; the purpose of this section is to discuss the representation of hyperobjects.
Suppose that we have a family of objects, parameterized by the variable , and continuous (in the appropriate way) with respect to . In other words, for every value of , we have an objects instance of the hyperobject, and a small difference between parameters implies a small difference between the instance objects, i.e. a small (in an appropriate definition of distance).
The natural generalization of (15) to this higher dimensional case is:
(33) 
with basis functions over , or, equivalently, in the Fourier domain
(34) 
One convenient way to produce such basis function is a tensor product between some standard 3D object basis functions and basis functions appropriate for the parameter space,
(35) 
so that
(36) 
Appropriate choices of bases provide implicit regularization for reconstruction algorithm, so the basis functions should be chosen considering the types of objects and heterogeneity that we expect to model.
In fact, one such generalization of objects to hyperobjects is a natural generalization of the standard use of the sample points in the Fourier domain. A 3D object is typically represented in cryoEM algorithm using samples on grid points of the 3D Fourier transform of the object . One of the ways to formulate a simple instance of hyperobjects, is to consider the 4D hyperobject and to represent is using sample points on a 4D grid of the 4D Fourier transform of the hyperobject , where the 4D Fourier transform is the composition of the Fourier transform on the spacial axes and the Fourier transform in the heterogeneity/“time” axis (the number of sample points on each axis can be different from the number of sample points on other axes). Slightly more generally, can be replaced with some arbitrary other transform (see Section 4.1).
4.5 Unknown orientations and heterogeneity parameters, and reconstruction of hypermolecules in cryoEM
The advantage of this general form of expressing hyperobjects is that the general linear operator formulation (see (20)) applies with no conceptual change. An image is produced by the operation
(37) 
where is a vector of coefficients that represents a hyperobject (as sample points or in some other basis) and is an operator that produces an image at the correct orientation and parameter value . While this is very similar to the formulation of discrete heterogeneity (see (24)), this formulation highlights the idea of treating the heterogeneity parameter in the same conceptual way as the rotation parameter (to the extent possible).
The traditional direct reconstruction algorithm is again a generalization of (21) and (22):
(38) 
and
(39) 
What makes it different from (25) and (26) is that can be continuous (or finely discretized), and that it reconstructs the hyperobject in a way that uses relations between the object instances at different values of . For example, it is possible to have a different value of for each image (in fact, this is the realistic model since is a continuous parameter, the differences between close points would of course have to be small); traditionally, the values of would either be restricted to a discrete number of unordered values that correspond to existing classes, or this would mean that each image would be assigned to a separate class, making it impossible to reconstruct that class. At the same time, nothing would have been known about intervals where no samples had been taken. Indeed, in this generality, the reconstruction of a hyperobject appears to be severely illposed, but so is the reconstruction of a single object. Loosely speaking, the same conceptual problem exists in the homogeneous case (with somewhat different properties): each image has its own rotation so that most points in the Fourier space are never sampled more than once or twice. The reconstruction of nonheterogeneous objects is achieved using additional implicit and explicit properties, priors, regularization and penalties. We propose to generalize the existing algorithms and to apply similar principals to the generalized problem. In this section we briefly review some of the approaches for generalization that we find of particular interest for continuous heterogeneity. The same general concepts apply to more elaborate approaches and algorithms used for reconstruction of objects in the homogeneous case and in the case of discrete heterogeneity.
The following are some approaches to explicit or implicit regularization and priors. These approaches can be applied to discretized sampled states (connecting the currently independent states) or to more continuous representations of states (indeed, to some extent, these are equivalent, see Section 3.1). The thread that connects these methods is the relation between states, or “continuity.” Some of these approaches are standard constructions once the problem is stated as proposed in this paper, and some are less standard. The idea behind all these methods is to capture the relation between different “states” (heterogeneity parameter values) of the hyperobject instead of treating them as independent entities.

Representing hyperobjects using tensor products of standard 3D object basis functions and parameterization function (space or frequency domains):
(40) for all combinations of .

A similar tensor product, with restrictions on the combinations of allowed. For example, allowing high frequency object functions (indexed by ) to be combined with more (or fewer) parameterization functions (indexed by ).

Penalties on some components, such as coefficients of high frequency heterogeneity functions.

High dimensional functions of both space/frequency and parameter, and in particular, spaces of basis functions that are chosen to better capture predicted types of structure and variability (e.g., superimposing a homogeneous global representation of an object with a localized basis that allows variability in spatially restricted certain areas of the object).

Using frames (e.g., high dimensional wavelet frames), which capture the structure and variability (possibly with sparsity constraints).

Optimal Transportation distance, Earth mover’s distance (EMD) or Wasserstein distance (see discussion below).

Sparsity in the use of basis function or in the spacial objects.

Total variation (TV) regularization in the spacial domain or a generalized form.

Priors or constraints on the distribution of images in the heterogeneity parameter space (and combinations of heterogeneity parameter and orientations).

Tree or multiscale heterogeneity structures (see below).

Hyperfrequency marching (see Section 5.2 below).

Continuity through image assignment rather than object representation: soft assignment of class and rotation, which is aware of class topology, so that an image that is assigned probability to be certain class and rotation is also assigned probability to be in nearby classes with similar rotations (up to global rotations).
We note that a source of true continuous heterogeneity in cryoEM is flexibility of the molecule. To the extent that the molecule can be modeled as a sum of smaller objects (atom or substructures), it may be useful to regularize the variability in the hyperobject using a penalty based on EMD (see, for example, [35]), so that for small , the EMD between the objects instance at parameter value and the nearby object distance at parameter value is small. Loosely speaking, EMD distinguishes between local variability in distributions of masses and more global redistributions of mass, so it may be suited to capture local variability in the location of sub units. Extensions of this treatment of distances provides, for example, methods for interpolating between instances.
The tree construction could potentially have several interesting properties. In its simplest form, the tree construction uses parameter basis functions that decompose the parameter space into “intervals,” in each interval the function is constant. Technically speaking, and in the absence of additional constraints, this construction appears to be no different than a standard decomposition to multiple discrete classes practiced in existing algorithms. Furthermore, it is essentially equivalent to using the Haar basis (see (14)) for the parameter expansion, which yields noncontinuous functions, but otherwise falls nicely into the proposed approach to expand in arbitrary bases. However, when regularization, such as highfrequency component penalties, is introduced, the different instances are tied together. Indeed, this construction does not offer a clear continuous parameterization of the heterogeneity, but, it could be less sensitive to the choice of parameterization topology and expansion (at the cost of not being adapted to using prior knowledge of the appropriate topology and expansion). This approach can be extended with additional multiscale representations; preliminary results indicate that adding nonconstant basis functions to the expansion are likely to contribute to the construction.
Due to the limited scope of this first discussion of this approach, our discussion of the first prototype algorithm is focused on implicit regularization through the choice of basis functions and frequency marching. For example, if the parameter space is , we use in the tensor product (see (36)) of standard bases used for objects, and and a family of polynomials in the heterogeneity axis (for more details, see Section 5).
4.6 A remark on ambiguity
As is wellknown, there are multiple equivalent solutions to classic cryoEM reconstruction problems due to certain symmetries. For example, classic cryoEM cannot distinguish between a molecule and its mirror image, so two mirror solutions are considered equivalent. Similarly, any solution is unique only up to rotation of the molecule. Furthermore, in the case of discrete heterogeneity, the solution is unique only up to permutations of the classes (i.e. any of the reconstructed molecule can be called “molecule number 1”), and rotation of each class independently of the others.
The reconstruction of hyperobjects is subject to some similar ambiguities. For example, a hyperobject is generally equivalent to a rotated version of itself or a reflected version. In addition, a hyperobject created by rotating the object to a different orientation at every value of the heterogeneity parameter , as define by the formula
(41) 
would be equivalent to . Typically, would be continuous with respect to the parameter space if we assume a continuous representations. Furthermore, in the absence of a unique metric on the parameter space , the hyperobject is subject to reparameterization. For example, suppose that the parameter space is one dimensional: , then the reparameterized hyperobject defined by the formula
(42) 
where is an increasing function on the interval with
, is equivalent to the original hyperobject. In some cases, it may be possible to regularize the hyperobject or introduce a metric that reduces the ambiguity in the parameterization, for example by requiring a parameterization where the samples are uniformly distributed in the parameter space, or penalizing high order heterogeneity coefficients.
We note that some ambiguities can technically yield hyperobjects that are not equivalent. For example, it is technically possible for an algorithm to use the degrees of freedom provided by the hyperobject model to capture shifts in images instead of an interesting variation (much like it is technically possible that some existing algorithm would generate two molecules that are rather similar shifted copies of one another). Such ambiguities and parameterization considerations are likely to be of interest in more elaborate implementation of the approach proposed in this manuscript.
5 Algorithms
The approach presented in this manuscript is quite general, and can be used to generalize various algorithms. The purpose of this section is to briefly present a simplified algorithm which we have implemented as proofofconcept for our approach. Since this particular algorithm is only one of many ways to implement our approach, we will present it in limited detail. We defer the more comprehensive discussion of the ideas introduced in this section to future papers.
5.1 Expansion of hyperobjects
We represent an object using concentric shells in the Fourier domain, with the function on each shell expanded in spherical harmonics, as in [8]. The th shell of an object is represented using the expansion
(43) 
where and define a position on the sphere, are spherical harmonics (see (5)), and are coefficients of the expansion. For each image we compute the Fourier expansion of concentric circles in the Fourier domain. The th circle in the Fourier transform of the image is expanded as follows:
(44) 
Unlike the implementation in [8], all our operations are computed directly with spherical harmonics using rotation operators that produce the coefficients of a rotated sphere from the coefficients of a sphere. The rotations of the sphere and the restriction of the sphere to the circles are sparse operations (see (9) and (12), respectively) in the sense that they mix very few of the coefficients in the expansion of the object.
Our preliminary implementation of hyperobjects allows one dimensional parameterization over the interval . We expand the hyperobject in a basis that is the tensor product between the basis for objects in (43) and heterogeneity basis functions which we denote by ; here we chose either normalized Legendre polynomials or Chebyshev polynomials (see, inter alia, [6]), shifted to the appropriate interval:
(45) 
5.2 “Hyperfrequency marching”
The algorithm implemented here generalizes the idea of frequency marching presented in [8], where the expansion is initially restricted a small number of shells, with later iterations allowing a growing number of shells in the expansion (larger values of allowed in later iterations). Our implementation follows a similar approach, increasing the allowed in later iterations, however, we also restrict the parameterization basis; we start with (no heterogeneity), and later increase the number of basis functions allowed in the expansion.
5.3 A simplified stochastic optimization algorithm
The algorithm implemented here is based on Stochastic Gradient Descent (SGD). An SGD algorithm for the Bayesian approach to cryoEM has been proposed in [21], our simplified implementation is an SGD version of a more traditional algorithm, with the continuous heterogeneity generalization proposed in this paper. At each iteration of the algorithm, we choose a small number of images (a “minibatch”), and estimate the orientation and heterogeneity parameter of each image. We then compute the gradient of the object that would decrease the discrepancy between the images in the minibatch and the object at the selected orientation, and make a small update to the object in that direction.
5.4 Sampling
Our algorithm requires samples the continuum of possible rotations and heterogeneity parameter values in order to generate templates at every iteration. In the current implementation, we sample both the rotations and parameter values randomly, with the exception of in plane rotations of the template, in which the computation of the cost function is accelerated using FFT as discussed in [4, 8].
6 Experimental Results
The simplified algorithm discussed above was implemented in Matlab. This simplified proof of concept does not consider shifts in the images, CTF, etc., but it does included simulated noise.
We simulated a hyperobject (a “cat”) composed of Gaussian elements in real space, where each Gaussian follows a continuous trajectory as a function of the parameter, so that we have a continuous space of objects. Examples of simulated 3D objects instances are presented in Figure 1. The hyperobject displays large scale extensive heterogeneity in the form of flexibility which is difficult to model as a combination of small number of rigid objects. We used the simplified imaging model (see (17)) to simulate pixels images in various orientations and parameter values, and added Gaussian noise to the images. The SNR in this experiment was . Examples of the simulated images are presented in Figure 4.
In our preprocessing stage, we compute the Fourier transform of the image at sample points on concentric circles in the Fourier domain, and then use FFT to obtain the image coefficients on concentric circles defined in (44).
We run the algorithm ab initio, without an initial guess. We initialize the hyperobject to a nonheterogeneous spherically symmetric object, by setting each sphere to the average of the corresponding circles across all images, i.e. we set the th order coefficient in the expansion of the th shell to the normalized average of the th order coefficients of the th circle of all the images.
We use SGD iterations to update the nonheterogeneous model object and after several steps of frequency marching we alternate between increasing the frequency limit and the number of parameterization basis functions .
For the purpose of visualization of the computed hyperobject, we choose sample values of the parameter , and then sampled the shells at quadrature points in the Fourier domain and used NUFFT (see [36]) to recover each object at regular grid points in real space. We use Matlab’s “isosurface” and “patch” to visualize level sets of the objects.
The algorithm was run on computer equipped with Intel(R) Core(TM) i74770 CPU, 32GB RAM server and a single NVIDIA GeForce GTX 980 Ti GPU with 6 GB GPU RAM for about 5 hours. We note that very similar results have been obtained even in the 12 hours of this experiment and of similar shorter experiments.
Examples of reconstructed objects are presented in Figure 2. We observe that the reconstructed objects appear to be slightly rotated one with respect to the other compared to the simulated data, due to the ambiguity discussed in Section 4.6. In fact, the relative rotations in the result may reflect a better choice than our choice in the simulation, in terms of the rate of change in the hyperobject, or in terms of the norm of the heterogeneity coefficients.
In Figure 5, we present the distribution of the pairs of trueparameter (xaxis) and estimated parameter (yaxis) assigned to images during the final steps of the refinement, before the algorithm was stopped. This current simple implementation does not regularize or reparameterize the parameter space, and has been used ab initio, without any initial estimate; while it succeeds in obtaining an appropriate parameterization in the surprising majority of the runs, in some runs, typically when using a small number of iterations or very rapid frequency marching in early stages of the run, we find some examples of inefficient parameterization, which is typically not difficult to detect. For example, in Figure 6 we find that the algorithm maps most of the parameter space to extreme points, preliminary results suggest that even this naive implementation of the algorithm gradually reparameterizes this mapping, but this process is relatively slow. This incident is relatively easy to detect by examining the marginal distribution on the yaxis (under the assumption or constraint of a more uniform distribution). Another example, shown in Figure 7 demonstrates a way in which the algorithm can use only part of the allowed parameter space. The reconstruction results are still good in the region of space that the algorithms chose to use, but it is an inefficient use of the parameter space. This case is also easy to detect and correct. Conservative optimization appear to mitigate these issues even in the current implementation; initial estimates of the objects and better regularization are also likely to mitigate these issues when they occur.
Video visualizations of the results are available at the website
http://roy.lederman.name/cryoem.
7 Summary and future work
The idea behind this work is to extend the existing tools in cryoEM to the case of continuous heterogeneity. Conceptually, we attempt to treat the heterogeneity parameter (or “class”) in the same way as the orientation. We propose various approaches for representing the heterogeneous hyperobjects and for regularizing the problem. In this manuscript, we presented one simple expansion of the hyperobject and one simple algorithm and implementation of these ideas as proofofconcept. In this section we discuss some future directions for this line of work.
We believe that various versions of this general approach have several advantages compared to the construction of independent classes: it allows to use assumptions on the topology and properties of the heterogeneity in a given problem to improve the reconstruction, it captures the topology of the heterogeneous structures and provides a natural way to reconstruct object instances at arbitrary parameter points, it combines information from different nearby object instances to produce a finegrained spectrum of states, leveraging continuity assumptions.
While many of the representation and regularization techniques proposed here can be implemented using sampled discretized values of the heterogeneity parameter, we find is particularly useful to use continuous bases in this discussion, because they offer a natural generalization of existing algorithms. In this manuscript we implemented only several continuous expansion using a tensor product of bases of functions and standard bases for the expansion of 3D objects. More elaborate bases and other continuous and discrete representations of hyperobjects are likely to be useful in capturing properties of the heterogeneity and further development of this approach.
One possible source for methods that could be adapted, is the existing body of work in 4DCT. Since 4DCT typically treats the case where both the orientation and heterogeneity parameter value (phase in the breathing cycle) are known, not all the ideas in 4DCT are directly applicable. At the same time, the investigation of unknown parameters in cryoEM may contribute to methods in CT, in cases where the breathing phase is not recorded accurately, or when there are additional heterogeneity parameters.
The examples used in this manuscript included only one dimensional heterogeneity on an interval, however the approach is general and can be used to investigate more elaborate topologies of heterogeneity. The investigation of other topologies and the sensitivity to the choice of topology, metrics, and bases is another direction in this line of work. We note that since the algorithms are relatively fast, it is possible to try multiple topologies, and we note that the tree/Haar approach is likely to be less sensitive to the precise choice of topology. However, the choice of topology is a useful implicit regularizer or prior. Another aspect of the work on topology is to develop methods of detecting local problems in the parameterizations that algorithms discover. While our discussion has been devoted mostly to continuous heterogeneity, the methods are applicable to discrete heterogeneity, for example when there are small changes between two conformations. A related matter that merits additional investigation is the case where some states in the continuum of states are far more represented than others, and the related issue of defining a metric on the space of heterogeneity.
The approach is general and we submit that it can be used in various existing cryoEM algorithms (although the actual software implementations may require significant modifications). It is possible to combine existing software with this approach in several ways, for example, by crude continuous classification of images in this approach, followed by a local reconstruction, or a reconstruction of an initial homogeneous model of the molecule, and then superimposing local basis functions in areas of variability in the model and running this type of algorithm to resolve the heterogeneity.
While the results obtained using this preliminary prototype are promising, additional work is needed to complete the implementation of this approach in an algorithm that is efficient enough to reconstruct high resolution hyperobjects in real cryoEM settings, while utilizing modest computational resources.
8 Conclusions
A framework has been presented for the tomography inverse problem in the case of continuously heterogeneous objects, where the orientation and heterogeneity parameter values are unknown. The proposed framework generalizes existing approaches for the reconstruction of molecules in cryoEM, so that, in principal, existing algorithm can be adapted to the case of continuous heterogeneity.
The approach has been demonstrated in simplified simulated cryoEM settings, using one of the proposed implementations, and a new prototype algorithm.
We are currently working on expanding our investigation of alternative representations and regularization approaches within this framework, and adapting the implementation to realworld cryoEM applications.
Acknowledgments
The authors were partially supported by Award Number R01GM090200 from the NIGMS, FA95501210317 and FA95501310076 from AFOSR, Simons Foundation Investigator Award, Simons Collaboration on Algorithms and Geometry, and the Moore Foundation DataDriven Discovery Investigator Award.
References
 [1] A. Doerr, “Singleparticle cryoelectron microscopy,” Nature Methods, vol. 13, no. 1, pp. 23–23, 2016.
 [2] E. Nogales, “The development of cryoem into a mainstream structural biology technique,” Nature Methods, vol. 13, no. 1, pp. 24–27, 2016.
 [3] R. M. Glaeser, “How good can cryoem become?,” Nature methods, vol. 13, no. 1, pp. 28–32, 2016.
 [4] D. A. Low, M. Nystrom, E. Kalinin, P. Parikh, J. F. Dempsey, J. D. Bradley, S. Mutic, S. H. Wahab, T. Islam, G. Christensen, et al., “A method for the reconstruction of fourdimensional synchronized ct scans acquired during free breathing,” Medical physics, vol. 30, no. 6, pp. 1254–1263, 2003.
 [5] R. R. Lederman and A. Singer, “A representation theory perspective on simultaneous alignment and classification,” arXiv preprint arXiv:1607.03464, 2016.
 [6] M. Abramowitz and I. Stegun, “Handbook of mathematical functions,” National Bureau of Standards, Applied Mathematics Series, vol. 55, p. 1046, 1964.
 [7] R. R. Coifman and G. Weiss, “Representations of compact groups and spherical harmonics,” Enseignement math, vol. 14, pp. 121–173, 1968.
 [8] A. Barnett, L. Greengard, A. Pataki, and M. Spivak, “Rapid solution of the cryoem reconstruction problem by frequency marching,” arXiv preprint arXiv:1610.00404, 2016.
 [9] J. Frank, Threedimensional electron microscopy of macromolecular assemblies: visualization of biological molecules in their native state. Oxford University Press, 2006.
 [10] M. Liao, E. Cao, D. Julius, and Y. Cheng, “Structure of the trpv1 ion channel determined by electron cryomicroscopy,” Nature, vol. 504, no. 7478, pp. 107–112, 2013.
 [11] M. Van Heel, “Angular reconstitution: a posteriori assignment of projection directions for 3d reconstruction,” Ultramicroscopy, vol. 21, no. 2, pp. 111–123, 1987.
 [12] A. Singer, R. R. Coifman, F. J. Sigworth, D. W. Chester, and Y. Shkolnisky, “Detecting consistent common lines in cryoem by voting,” Journal of structural biology, vol. 169, no. 3, pp. 312–322, 2010.

[13]
A. Singer and Y. Shkolnisky, “Threedimensional structure determination from common lines in cryoem by eigenvectors and semidefinite programming,”
SIAM journal on imaging sciences, vol. 4, no. 2, pp. 543–572, 2011.  [14] Y. Shkolnisky and A. Singer, “Viewing direction estimation in cryoem using synchronization,” SIAM journal on imaging sciences, vol. 5, no. 3, pp. 1088–1110, 2012.
 [15] S. H. Scheres, M. Valle, R. Nuñez, C. O. Sorzano, R. Marabini, G. T. Herman, and J.M. Carazo, “Maximumlikelihood multireference refinement for electron microscopy images,” Journal of molecular biology, vol. 348, no. 1, pp. 139–149, 2005.
 [16] F. J. Sigworth, P. C. Doerschuk, J.M. Carazo, and S. H. Scheres, “Chapter tenan introduction to maximumlikelihood methods in cryoem,” Methods in enzymology, vol. 482, pp. 263–294, 2010.
 [17] N. Grigorieff, “Frealign: highresolution refinement of single particle structures,” Journal of structural biology, vol. 157, no. 1, pp. 117–125, 2007.
 [18] S. H. Scheres, “A bayesian view on cryoem structure determination,” Journal of molecular biology, vol. 415, no. 2, pp. 406–418, 2012.
 [19] S. H. Scheres, “Relion: implementation of a bayesian approach to cryoem structure determination,” Journal of structural biology, vol. 180, no. 3, pp. 519–530, 2012.
 [20] A. Punjani, M. Brubaker, and D. Fleet, “Building proteins in a day: Efficient 3d molecular structure estimation with electron cryomicroscopy,” IEEE Transactions on Pattern Analysis and Machine Intelligence, 2016.
 [21] A. Punjani, J. L. Rubinstein, D. J. Fleet, and M. A. Brubaker, “cryosparc: algorithms for rapid unsupervised cryoem structure determination,” Nature Methods, vol. 14, no. 3, pp. 290–296, 2017.
 [22] N. Grigorieff, “Chapter eightfrealign: An exploratory tool for singleparticle cryoem,” Methods in enzymology, vol. 579, pp. 191–226, 2016.
 [23] D. Lyumkis, A. F. Brilot, D. L. Theobald, and N. Grigorieff, “Likelihoodbased classification of cryoem images using frealign,” Journal of structural biology, vol. 183, no. 3, pp. 377–388, 2013.
 [24] M. van Heel and M. StöfflerMeilicke, “Characteristic views of e. coli and b. stearothermophilus 30s ribosomal subunits in the electron microscope.,” The EMBO journal, vol. 4, no. 9, p. 2389, 1985.
 [25] E. Katsevich, A. Katsevich, and A. Singer, “Covariance matrix estimation for the cryoem heterogeneity problem,” SIAM journal on imaging sciences, vol. 8, no. 1, pp. 126–185, 2015.
 [26] J. Andén, E. Katsevich, and A. Singer, “Covariance estimation using conjugate gradient for 3d classification in cryoem,” in Biomedical Imaging (ISBI), 2015 IEEE 12th International Symposium on, pp. 200–204, IEEE, 2015.
 [27] Y. Aizenbud and Y. Shkolnisky, “A maxcut approach to heterogeneity in cryoelectron microscopy,” arXiv preprint arXiv:1609.01100, 2016.
 [28] P. Schwander, R. Fung, and A. Ourmazd, “Conformations of macromolecules and their complexes from heterogeneous datasets,” Phil. Trans. R. Soc. B, vol. 369, no. 1647, p. 20130567, 2014.
 [29] J. Frank and A. Ourmazd, “Continuous changes in structure mapped by manifold embedding of singleparticle data in cryoem,” Methods, vol. 100, pp. 61–67, 2016.
 [30] X. Jia, Y. Lou, B. Dong, Z. Tian, and S. Jiang, “4d computed tomography reconstruction from fewprojection data via temporal nonlocal regularization,” in International Conference on Medical Image Computing and ComputerAssisted Intervention, pp. 143–150, Springer, 2010.
 [31] H. Gao, R. Li, Y. Lin, and L. Xing, “4d cone beam ct via spatiotemporal tensor framelet,” Medical physics, vol. 39, no. 11, pp. 6943–6946, 2012.
 [32] D. Kazantsev, G. Van Eyndhoven, W. Lionheart, P. Withers, K. Dobson, S. McDonald, R. Atwood, and P. Lee, “Employing temporal selfsimilarity across the entire time domain in computed tomography reconstruction,” Phil. Trans. R. Soc. A, vol. 373, no. 2043, p. 20140389, 2015.
 [33] V. P. Gopi, P. Palanisamy, K. A. Wahid, P. Babyn, and D. Cooper, “Multiple regularization based mri reconstruction,” Signal Processing, vol. 103, pp. 103–113, 2014.
 [34] D. A. Low, P. J. Parikh, W. Lu, J. F. Dempsey, S. H. Wahab, J. P. Hubenschmidt, M. M. Nystrom, M. Handoko, and J. D. Bradley, “Novel breathing motion model for radiotherapy,” International Journal of Radiation Oncology* Biology* Physics, vol. 63, no. 3, pp. 921–929, 2005.

[35]
Y. Rubner, C. Tomasi, and L. J. Guibas, “The earth mover’s distance as a metric for image retrieval,”
International journal of computer vision
, vol. 40, no. 2, pp. 99–121, 2000.  [36] L. Greengard and J.Y. Lee, “Accelerating the nonuniform fast fourier transform,” SIAM review, vol. 46, no. 3, pp. 443–454, 2004.