Alignment of cryo-EM movies of individual particles by optimization of image translations
Direct detector device (DDD) cameras have revolutionized single particle electron cryomicroscopy (cryo-EM). In addition to an improved camera detective quantum efficiency, acquisition of DDD movies allows for correction of movement of the specimen, due both to instabilities in the microscope specimen stage and electron beam-induced movement. Unlike specimen stage drift, beam-induced movement is not always homogeneous within an image. Local correlation in the trajectories of nearby particles suggests that beam-induced motion is due to deformation of the ice layer. Algorithms have already been described that can correct movement for large regions of frames and for MDa protein particles. Another algorithm allows individual MDa protein particle trajectories to be estimated, but requires rolling averages to be calculated from frames and fits linear trajectories for particles. Here we describe an algorithm that allows for individual MDa particle images to be aligned without frame averaging or linear trajectories. The algorithm maximizes the overall correlation of the shifted frames with the sum of the shifted frames. The optimum in this single objective function is found efficiently by making use of analytically calculated derivatives of the function. To smooth estimates of particle trajectories, rapid changes in particle positions between frames are penalized in the objective function and weighted averaging of nearby trajectories ensures local correlation in trajectories. This individual particle motion correction, in combination with weighting of Fourier components to account for increasing radiation damage in later frames, can be used to improve 3-D maps from single particle cryo-EM.
The use of CMOS technology in direct detector device (DDD) cameras for electron cryomicroscopy (cryo-EM) has enabled the acquisition of exposure series ‘movies’. Movies of radiation sensitive specimens revealed that beam-induced motion blurs images [1, 2, 3]. DDD movies are typically acquired with exposures of 1 to 3 e/Å/frame on the specimen, which corresponds to 2 to 5 e/pixel/frame on the detector, depending on microscope magnification. These low exposures result in low signal-to-noise ratios (SNRs) in individual movie frames. Optimal extraction of high-resolution information from images of single particles requires alignment of movie frames, a process that is complicated by the low SNR. Fig. 1A shows the average of a 30 frame movie acquired with 2.5 e/pixel/frame, corresponding to 1.2 e/Å/frame on the specimen at 200 kV with a K2 Summit DDD (Gatan Inc). A few representative particles are circled in red. Fig. 1B shows a single frame from the movie, illustrating the low SNR of the frames. Frame alignment is complicated further by the presence of fixed pattern noise in images from errors in sensor gain normalization. Significant progress in image analysis has already been enabled by programs to perform rigid body translational alignment of entire field-of-view movie frames (currently 4000 4000 pixels for most cameras). A method introduced by Li and colleagues  decreases the weight of high spatial frequencies in images to suppress artifacts from fixed pattern noise before calculating pairwise cross-correlation functions between movie frames. The optimal frame displacement values from the cross-correlation functions are used to create a system of over-determined linear equations. Matrix algebra is then used to determined the frame-to-frame translations that best fit the data in a least squares sense. This least squares whole frame alignment method has allowed high-resolution structures to be determined for important biological macromolecules [4, 5, 6].
Cryo-EM of large particles with DDDs has shown that beam-induced motion cannot be described completely by rigid body translation of entire movie frames . Instead, these experiments suggested that the beam-induced movement of ice-embedded protein is better described by a translation of each particle in each frame. Examination of tilt pairs of images demonstrated that rotation of specimens, probably due to movement of the ice layer, does occur . However, the magnitudes of these rotations were small and will have the most significant effect on particles with large radii, like viruses. The consequences of specimen rotation can be neglected at present without limiting map resolution for particles smaller than 1 MDa. The translation of particles in frames can be written as for each particle in frame , where and are the difference in particle position between frame and frame in the and directions, respectively. If these translations are known, their inverse () can be applied to the particle images before averaging of frames to optimize the extraction of high-resolution information from the image. It is likely that accurate individual particle motion correction could extract information from images that is neglected by whole frame alignment. Despite the success of the least squares method for whole frame alignment, it was pointed out by the authors of the method that it is not reliably able to align image regions smaller than 2000 2000 pixels for movies acquired using typical conditions. As such, the least squares method is not capable of aligning regions of frames that contain individual particles in order to correct for deformation of the ice layer during imaging. A method to align individual particle images was introduced that is tightly integrated into the single particle orientation estimation framework of the program Relion and has resulted in several high-resolution structures [8, 9, 10]. For small particles this approach requires rolling averages of frames, which increases the SNR over individual frames but loses information about true trajectories. Also, the individual particle trajectories for small particles from this method include errors, and it is necessary to fit linear trajectories with constant velocities for particles, which are not necessarily a good approximation for their true trajectories. Furthermore, the approach cannot readily be used outside of the Relion software package.
Here we aim to identify the translations for movies of individual ice-embedded particles that best bring the frames into alignment for each particle, without the use of rolling frame averages or fitted linear trajectories. In order to produce a robust and computationally efficient method for correcting the effects of beam-induced movement in small regions in images, or on individual small ( MDa) particles, we pose the problem in terms of optimization. We propose an objective function based on the correlation of the Fourier transforms of individual frames with the sum of all frames. A well-established iterative optimization algorithm that makes use of partial derivatives of the objective function is then used to find the desired translation values. Once optimized, this objective function gives frame-to-frame trajectories for images of individual particles that show strong local correlation. We show that smoothing of trajectories for individual particles can be used to identify and correct beam-induced particle movement. This approach, in combination with compensation for the fading of Fourier components due to radiation damage, was implemented in a new program that we call alignparts_lmbfgs.
2 Methods and Results
2.1 Choice of objective function
Based on the observation that averages of unaligned particle frames appear blurred, a reasonable alignment for each region of the frame that contains a particle is the alignment that makes the sum of all of the frames best agree with each of the frames. Accordingly, we propose an objective function that maximizes the sum of the correlations of the Fourier transform of each shifted frame with the sum of the Fourier transforms of the shifted frames. Prior to analysis, we apply a temperature factor in Fourier space with the form to prevent fixed pattern noise from dominating the analysis . The effect of translation on the Fourier transform of a movie frame is a phase change, , in each Fourier component of the frame, written for the Fourier component of frame . The phase shifted Fourier component is given by or where . The amount of phase change is given by
where is the extent in pixels in both the and direction of the image, and and are the distance of the Fourier component from the origin in the and directions, respectively. As described above, and are the difference in particle position between frame and frame in the and directions, respectively. The Fourier transform of a sum is equal to the sum of Fourier transforms. Consequently, the Fourier component from the sum of the shifted frames of a movie with frames is given by . The unnormalized correlation between two Fourier transforms, and , is given by where denotes the complex conjugate. For the correlation between the sum image and the individual frame, these values must be summed for the Fourier components in a resolution band . It is only necessary to consider two times the real part of the expression for the correlation, because the Fourier transforms of real functions, such as images, are Hermitian, so that for every term in the correlation there is a corresponding term and adding these two terms removes their imaginary parts. In an objective function, the factor of 2 may be neglected without changing the position of the optimum, and the negative of the function can be used in order to interface with pre-existing optimization algorithms, which typically seek to minimize functions. Consequently, we can propose an objective function, , that meets the criterion described above:
where represents the set of all and values, and is the complex conjugate of , with calculated from and according to equation 1. With equation 2 as the objective function, iterative optimization methods can be used to explore the -dimensional space of frame translations to find values of and that minimize the function.
2.2 Partial derivatives of the objective function
Numerous algorithms exist for optimizing objective functions. Optimization problems can benefit greatly from the ability to analytically determine partial derivatives, or gradients, of the objective function with respect to all variables. The derivatives of and with respect to , the shift in the direction for the frame, are
Using these simplifications, the derivative of the objective function in equation 2 with respect to is
Noting that when and when , the expression simplifies further:
Similarly, the partial derivative of equation 2 with respect to is
We elected to use the limited memory Broyden-Fletcher-Goldfarb-Shanno (lm-bfgs) algorithm  to optimize the objective function in equation 2. By providing equations 2, 3, and 4 for lm-bfgs optimization, values of and were obtained for movies of V-ATPase particles in ice. Fig. 2A shows the calculated trajectories from optimization of 200 regions of 320 320 pixels in each frame. These 200 image regions were selected by template matching from the image in Fig. 1A, and contain a mixture of usable particle images and other image features. The trajectories show local correlation, even though at this stage in the analysis individual particle trajectories are not provided with any information about the trajectories of nearby particles, except for any overlap in the 320 320 pixel boxes. Close inspection of the trajectories in two regions of the micrograph (Fig. 2Bi and ii) reveals noise in the trajectories of individual particles obtained by the optimization method.
Although encouraging, the noise seen in trajectories of particles in Fig. 2Bi and ii suggests that the optimization does not show the true trajectories of individual particle images. One obvious approach to reducing noise in a trajectory is to calculate the trajectory from a larger portion of the image, thereby increasing the signal available for calculating the objective function. Unfortunately, as the size of the box used for determining particle positions increases, particles must progressively be excluded that fall too close to the edge of the image. Increasing box sizes also results in almost identical trajectories for nearby particles that may mask the local variation in movement that this technique aims to recover. Better noise removal can be achieved by using two reasonable assumptions that are neglected in the analysis presented in Fig. 2. The first assumption is that trajectories are unlikely to have sudden changes in direction, although the possibility of these changes cannot be eliminated. The second assumption is that nearby particle trajectories are correlated. Enforcing these two conditions can be used to ‘smooth’ particle trajectories to remove noise.
2.3.1 Second order smoothing
The assumption that true particle trajectories are unlikely to undergo sudden and dramatic changes in direction can be enforced by penalizing changes in and . If and are constant ( and are 0), the expected value for is . Deviation from this expected linear trajectory can be penalized by an amount . The overall penalty imposed on the objective function to encourage smoothness is then given by
where is a user selected weighting parameter. This penalty is known as second order smoothing because it penalizes finite difference approximations of the second derivatives of and with respect to , and . The penalty function described in equation 5 is added to the objective function in equation 2 to obtain the overall objective function that is optimized. The contribution to the penalty function in equation 5 from shifting of the frame when is
and consequently the first derivative of equation 5 with respect to is given by
The derivative of the smoothed objective function is therefore the sum of the values from equations 3 and 6 for the derivative with respect to , and the sum of the values from equations 4 and 7 for the derivative with respect to . Fig. 3A shows the effect of increasing values of the user set parameter for two regions on opposite sides of the micrograph (Fig. 3Ai and ii). With , the trajectories are noisy, as seen in Fig. 2. With a significant amount of noise has been removed from the trajectories. Note also that nearby trajectories appear to be correlated even though this condition has not been enforced. With , an excessively large number for these images, trajectories have been forced to become linear. Forcing trajectories to be linear is equivalent to fitting a single drift rate for each particle in the movie. The effects of different values of depend on the number of frames in a movie and the pixel values in those frames.
2.3.2 Local averaging for smoothing
Local correlation of nearby particle trajectories without the use of an increased box size can be achieved by weighted averaging after trajectories are calculated. In this approach, ‘raw trajectories’ are determined for individual particle images with or without second order smoothing. Once raw trajectories are determined, locally averaged trajectories are calculated according to
where is the smoothed displacement vector for the particle in the frame and is the original displacement vector for the particle in the frame. The weight is given by
where is the distance between the and particles and is a user set parameter that determines the extent to which the smoothing is applied. This Gaussian weighting is equivalent to the local averaging used for fitting linear trajectories in Relion . Because of the Gaussian form of equation 9, 95 % of of the weight for a particle trajectory will come from the trajectories within pixels of that particle. Fig. 3B shows the effect of increasing the parameter for two sets of nearby trajectories (Fig. 3Bi and ii), without the use of second order smoothing. With , the trajectories are noisy, as seen in Fig. 2. With a significant amount of noise has been removed from the trajectories, even though smoothness has not be enforced. With , an excessively large number, trajectories on opposite sides of the micrograph from each other have been forced to be similar. In this situation, depending on the number of particles selected in the micrograph, the method becomes a nearly rigid frame alignment. Ideal smoothing of particle trajectories comes from combining the two approaches described above. Fig. 4A shows trajectories with and . As can be seen in two enlarged regions from opposite sides of the micrograph (Fig. 4Bi and ii), individual particle trajectories appear smooth with strong local correlation but significant variation from one edge of the micrograph to the other.
2.3.3 Exposure weighting
Correcting for the movement between frames restores the phases of each Fourier component to their correct values. 2-D crystal studies have shown that diffraction spots fade according to the relationship , where is the instantaneous intensity of the diffraction spot before electron exposure, is the instantaneous intensity at exposure , and is the critical exposure at which the intensity of the diffraction spot fades to times its initial intensity [12, 13]. The use of estimates of to optimize the combination of DDD movie frames was proposed previously , with the necessary values measured from the fading of calculated diffraction spots from crystals. More recently these values were recorded for the same purpose from single particle analysis experiments . We included this exposure-weighting approach in alignparts_lmbfgs, utilizing estimates from the single particle data .
3 Characterization of the algorithm
To test the alignparts_lmbfgs algorithm, we compared the resolutions of maps calculated from a small dataset of images of 20S proteasome particles obtained with the same camera and microscope conditions described above for the V-ATPase (Latham, Benlekbir, and Rubinstein, in preparation). The resolution band used to align particle images was set at 500 to 40 Å and the temperature factor to 2000 Å. The alignment of particle images was not sensitive to most changes in the low-resolution frequency cutoff, while including higher-resolution information in the alignment by changing the high-resolution frequency cutoff and decreasing the temperature factor could introduce noise into the particle trajectories. Particle images were subjected to 2-D classification, 3-D classification, and 3-D map refinement in Relion. We compared whole frame alignment, individual particle motion correction with alignparts_lmbfgs, and individual particle motion correction including weighted averaging of Fourier components based on radiation damage. The resolutions that could be obtained (Fig. 5) were 4.17 Å for whole frame alignment, 3.82 Å for individual particle motion correction, and 3.68 Å for individual particle motion correction with compensation for radiation damage. It is important to note that when using Relion, the best resolution for a dataset is often obtained when the dataset is taken through the entire map calculation process, including 2-D classification and 3-D classification with each 3-D map here being constructed from approximately 14,500 particle images. Calculating a map from local motion corrected and exposure weighted images but selecting those images based on 2-D and 3-D classification of images from whole frame motion correction produced a map at 3.64 Å. However, calculating a map from whole frame motion corrected images that had been selected based on 2-D and 3-D classification of local motion corrected and exposure weighted images gave a resolution of 4.32 Å. The time required for the program to process images depends on the size of particle images, the number of frames in a movie, and the tolerance setting of the termination test within the lm-bfgs optimizer. Increasing second order smoothing by increasing the parameter also causes the algorithm to take more time to find the minimum in the objective function. However, a typical example consisting of particle images that are 256 256 pixels with 30 frames per movie required approximately 0.3 s per particle image when running as a single process on a single core of an Intel i7-4790K CPU with a 4.0 GHz clock rate.
For full frame alignment, the least squares algorithm proposed by Li and colleagues possesses an advantage over the approach described here, in that the frame translations are highly over-determined: a movie consisting of frames will provide equations that can be used to determine the frame translations needed to correct motion . However, while the least squares method correlates low SNR frames with other low SNR frames, alignparts_lmbfgs, or a whole frame implementation of the algorithm alignframes_lmbfgs, correlates low SNR frames with the relatively high SNR sum of frames. Consequently, the alignparts_lmbfgs method is able to work with image boxes at least as small as the 256 256 pixel boxes used here and alignframes_lmbfgs is likely able to align lower-contrast whole frames than the least squares methods. The alignparts_lmbfgs approach should behave similarly to a simplistic iterative approach where frames are averaged and individual frames are subsequently aligned to the average. Special care must be taken in this simple iterative approach to ensure that the frame is not aligned to an average where the frame has been included at a fixed position, which could bias the alignment of the frame. With alignparts_lmbfgs, the average always includes the frame with the translations for the frame that are being tested. Also, in alignparts_lmbfgs changing the translations for the frame instantaneously affects the correlation of all other frames with the sum image, while with the simple iterative approach it does not, possibly making the identification of a global optimum less robust. The simple iterative approach will also almost certainly be slower than alignparts_lmbfgs at finding the optimum alignment of frames. The alignparts_lmbfgs approach benefits from being able to incorporate the second order smoothness constraint directly into the objective function. Both algorithms could become trapped in a local alignment minimum. However, the form of equation 2 suggests that the problem may be convex and in practice local minima do not appear to cause problems. The Relion procedure  integrates estimation of particle trajectories with projection matching from a reference map of the protein complex. Both procedures attempt to regularize particle trajectories: Relion by using a running average of particle frames and fitting of a linear trajectory, alignparts_lmbfgs by introducing the second order smoothness constraint. The Relion approach has the potential advantage that projections from a refined 3D map will posses stronger signal than the sum of all frames used as a reference in alignparts_lmbfgs. The potential disadvantage of the Relion approach relative to alignparts_lmbfgs is that errors in contrast transfer function (CTF) estimation, structured noise in images from sample contamination or ice contamination, differing conformation of the protein particle in the image and map, and any other sources of inaccuracy in projection matching could affect the accuracy of trajectory estimation. Compared to the procedure introduced in Relion, alignparts_lmbfgs is much less computationally expensive.
The two different smoothing approaches, second order smoothing and local weighted averaging of trajectories, have different advantages and uses. Second order smoothing is independent of particle density in images. If images contain few particles, their trajectories should be smoothed by increasing the second order smoothing parameter . Local averaging of trajectories will have little effect for particles that are far apart from each other but can be applied effectively where there are many particles or other image features that can be aligned. For this situation, the amount of second order smoothing can be decreased somewhat. The value of the parameter that is used for aligning particle images depends on the value of the objective function in equation 2. Therefore, for the same amount of smoothing will need to be bigger for images that have larger pixel values or more frames than for images with smaller pixel values or fewer frames. Consequently, should be tuned to produce a physically reasonable trajectory for particles by testing a variety of values with a small subset of the data. In practice, we have found that for consistent imaging conditions a fixed value of produces consistent results. In comparison, the value should be set to reflect how quickly trajectories are expected to vary across the image, which will depend on the magnification at which image were obtained. At higher resolution, or for larger particles, it could be useful to interpret trajectories to include the rotations of particles in the ice layer.
5 Software and data availability
The programs alignparts_lmbfgs and alignframes_lmbfgs are available from
We thank Jianhua Zhao for providing the cryo-EM image of yeast V-ATPase used for producing figures 1 to 4 and Michael Latham and Samir Benlekbir for providing the 20S proteasome images used for the characterization of the algorithm shown in figure 5. We thank Tim Grant and Niko Grigorieff for providing data on radiation-induced fading of Fourier components ahead of publishing their manuscript and Hui Guo for characterizing an earlier version of the program. We are grateful to Alexis Rohou, Tim Grant, Niko Grigorieff, Richard Henderson, Sjors Scheres, Peter Rosenthal, and members of the Rubinstein laboratory for helpful comments on the manuscript. This work was funded by Natural Sciences and Engineering Research Council discovery grant 401724-12 (JLR) and the Canada Research Chairs program (JLR). A preprint of this manuscript was first published on arXiv.org on 24 Sept 2014 (http://arxiv.org/abs/1409.6789).
7 Author contributions
JLR conceived of the Fourier space objective function for optimizing particle movement estimates, wrote the programs, characterized their performance, and wrote the manuscript. MAB suggested the use of the gradient-based lm-bfgs algorithm, proposed the second order smoothing, and made other critical suggestions for the programs and the manuscript.
-  G. McMullan, A. Faruqi, Electron microscope imaging of single particles using the medipix2 detector, Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment 591 (1) (2008) 129–133.
-  R. Glaeser, G. McMullan, A. Faruqi, R. Henderson, Images of paraffin monolayer crystals with perfect contrast: minimization of beam-induced specimen motion, Ultramicroscopy 111 (2) (2011) 90–100.
-  A. F. Brilot, J. Z. Chen, A. Cheng, J. Pan, S. C. Harrison, C. S. Potter, B. Carragher, R. Henderson, N. Grigorieff, Beam-induced motion of vitrified specimen on holey carbon film, J Struct Biol 177 (3) (2012) 630–7. doi:10.1016/j.jsb.2012.02.003.
-  X. Li, P. Mooney, S. Zheng, C. R. Booth, M. B. Braunfeld, S. Gubbens, D. A. Agard, Y. Cheng, Electron counting and beam-induced motion correction enable near-atomic-resolution single-particle cryo-em, Nat Methods 10 (6) (2013) 584–90. doi:10.1038/nmeth.2472.
-  E. Cao, M. Liao, Y. Cheng, D. Julius, Trpv1 structures in distinct conformations reveal activation mechanisms, Nature 504 (7478) (2013) 113–8. doi:10.1038/nature12823.
-  M. Liao, E. Cao, D. Julius, Y. Cheng, Structure of the trpv1 ion channel determined by electron cryo-microscopy, Nature 504 (7478) (2013) 107–12. doi:10.1038/nature12822.
-  R. Henderson, S. Chen, J. Z. Chen, N. Grigorieff, L. A. Passmore, L. Ciccarelli, J. L. Rubinstein, R. A. Crowther, P. L. Stewart, P. B. Rosenthal, Tilt-pair analysis of images from a range of different specimens in single-particle electron cryomicroscopy, J Mol Biol 413 (5) (2011) 1028–46. doi:10.1016/j.jmb.2011.09.008.
-  S. H. W. Scheres, Relion: implementation of a bayesian approach to cryo-em structure determination, J Struct Biol 180 (3) (2012) 519–30. doi:10.1016/j.jsb.2012.09.006.
-  S. H. Scheres, Beam-induced motion correction for sub-megadalton cryo-em particles, Elife 3 (2014) e03665.
-  M. G. Campbell, D. Veesler, A. Cheng, C. S. Potter, B. Carragher, 2.8 å resolution reconstruction of the thermoplasma acidophilum 20s proteasome using cryo-electron microscopy, eLife 4 (2015) e06380.
-  R. H. Byrd, P. Lu, J. Nocedal, C. Zhu, A limited memory algorithm for bound constrained optimization, SIAM Journal on Scientific Computing 16 (5) (1995) 1190–1208.
-  P. N. T. Unwin, R. Henderson, Molecular structure determination by electron microscopy of unstained crystalline specimens, Journal of molecular biology 94 (3) (1975) 425–440.
-  S. B. Hayward, R. M. Glaeser, Radiation damage of purple membrane at low temperature, Ultramicroscopy 4 (2) (1979) 201–210.
-  L. A. Baker, E. A. Smith, S. A. Bueler, J. L. Rubinstein, The resolution dependence of optimal exposures in liquid nitrogen temperature electron cryomicroscopy of catalase crystals, J Struct Biol 169 (3) (2010) 431–7. doi:10.1016/j.jsb.2009.11.014.
-  T. Grant, N. Grigorieff, Measuring the optimal exposure for single particle cryo-em using a 2.6 å reconstruction of rotavirus vp6, eLife (2015) e06980.