Search tips
Search criteria 


Logo of hhmipaabout author manuscriptssubmit a manuscriptHHMI Howard Hughes Medical Institute; Author Manuscript; Accepted for publication in peer reviewed journal
J Struct Biol. Author manuscript; available in PMC 2007 December 19.
Published in final edited form as:
PMCID: PMC2147720

A maximum likelihood approach to two-dimensional crystals


Maximum Likelihood (ML) processing of transmission electron microscopy images of protein particles can produce reconstructions of superior resolution due to a reduced reference bias. We have investigated a ML processing approach to images centered on the unit cells of two-dimensional (2D) crystal images. The implemented software makes use of the predictive lattice node tracking in the MRC software, which is used to window particle stacks. These are then noise-whitened and subjected to ML processing. Resulting ML maps are translated into amplitudes and phases for further processing within the 2dx software package. Compared with ML processing for randomly oriented single particles, the required computational costs are greatly reduced as the 2D crystals restrict the parameter search space. The software was applied to images of negatively stained or frozen hydrated 2D crystals of different crystal order. We find that the ML algorithm is not free from reference bias, even though its sensitivity to noise correlation is lower than for pure cross-correlation alignment. Compared with crystallographic processing, the newly developed software yields better resolution for 2D crystal images of lower crystal quality, and it performs equally well for well-ordered crystal images.

Keywords: maximum likelihood, electron crystallography, protein structure, single particle, 2dx


Electron crystallography of two-dimensional (2D) crystals is a commonly used technique to obtain high-resolution three-dimensional (3D) structures of proteins (for a recent review, see Renault et al., 2006). The technique was developed mainly by Henderson and co-workers (Henderson and Unwin, 1975) and one of its first applications led to an atomic model of bacteriorhodopsin (Henderson et al., 1990). Many other structures have been solved at a resolution that allowed an interpretation with an atomic model. The data collected from 2D crystals come in two flavors: Images and electron diffraction pattern. Electron diffraction data provide intensities of diffraction spots, much like X-ray diffraction of 3D crystals. These can be measured, and their square root gives the amplitude components of the protein structure. The images of the 2D crystals can be Fourier transformed, and both, amplitudes and phases of the calculated diffraction spots can then be measured by Fourier extraction. This is in contrast to the situation in X-ray crystallography, where phases are not directly observable. Therefore, to collect a complete data set from 2D crystals by electron crystallography, both electron diffraction and imaging are usually performed.

Images also offer another advantage over X-ray crystallography. If disorder is present in a crystal and limits the resolution in a diffraction pattern, the crystal distortions in an image of the crystal can be corrected computationally. This procedure is commonly referred to as “unbending” and significantly improves the resolution attainable by Fourier extraction from 2D crystal images. For the unbending process, a reference image containing only a small number of unit cells is generated either from a filtered version of the image itself, or by projecting an already existing 3D reference structure using the program MAKETRAN (Kunji et al., 2000). A cross-correlation map is then calculated between the reference and the image to determine the location of each unit cell in the crystal. Using the autocorrelation function of the reference, the correlation map can be searched for peaks and unit cell locations are recorded. The MRC program QUADSEARCH (Crowther et al., 1996) performs such a search and exploits the a-priori knowledge of the approximate location of the peaks. It uses an iterative refinement process to predict the peak location in the cross-correlation map, and then searches for the actual local maximum in that map within a limited radius of the predicted location. A coordinate list of the identified lattice nodes is then generated, which is used by the MRC program CCUNBEND to correct the crystal distortions in the image by one of two methods: A better ordered crystal image is generated by either smoothly warping the image so that the unit cells fall onto a perfect lattice, or by creating a discrete montage of unit cells placed at the crystallographically determined grid points. The Fourier transform of the unbent image is then evaluated to obtain values for the amplitudes and phases of the crystal projection map in that image.

The unbending procedure described above has been applied very successfully in many cases and is in common use, but it also has some limitations. First, in its current implementation, it cannot correct for in-plane rotational disorder. Second, the use of correlation functions for alignment is prone to error when the signal is weak, due to the increased chance of mistaking a noise peak for the correlation peak produced by the signal. The signal in the correlation map depends on the size of the reference area used to find the unit cell locations. A larger reference area will produce stronger signal peaks but will also be less sensitive to short-range disorder in the crystal. To detect and remove short-range disorder, a smaller reference area containing one unit cell or even a single molecule would be desirable.

In the case of strong irregularities in a badly ordered 2D crystal, or variations among the unit cells due to sample heterogeneity, single-particle image processing can be applied to the recorded 2D (pseudo-) crystal images. The goal of the single particle processing is similar to the 2D crystal unbending procedure. In both cases alignment and averaging of individual molecules or their assemblies is done to enhance the signal. Sass et al., 1989, have combined phases from correlation averaging (CA) with amplitudes from electron diffraction, obtaining a 3.5Å projection map of porin. Schultz et al., 1993, imaged negatively stained 2D monolayer crystals by tomographic single-axis tilt series, calculated a 3D reconstruction of the unit cell via single particle methods, and combined two molecules from the unit cell with non-crystallographic symmetry to fill the missing cone. Sherman et al., 1998, applied multivariate statistical analysis to 2D crystal images. Stahlberg et al., 1998, used single particle methods to detect and correct for the non-crystallographic orientation of the photosynthetic reaction center within well-ordered 2D crystal images of the surrounding light-harvesting-complex I. Tahara et al., 2000, could significantly improve the resolution of a Na+/K+-ATPase projection map by applying single particle processing methods to 2D crystal images and allowing for sample heterogeneity.

Here, we apply a maximum likelihood (ML) approach to the single-particle processing of 2D crystal images. The ML processing was introduced for the processing of images of non-crystalline material (single particles), and was shown to have superior performance at low signal-to-noise ratios (SNR, variance ratio of signal over noise) compared with correlation-based alignment (Sigworth, 1998). Combining the single particle processing with an ML approach can therefore lead to an improvement over the currently used unbending process. We utilize a whitening filter to make the ML method applicable to the real data which have nonwhite noise. We discuss here the application of ML to 2D crystals and compare its performance with the traditional correlation-based unbending. Furthermore, we show how the contrast transfer function (CTF) of the electron microscope can be included in the processing.


Maximum likelihood for 2D crystals

The application of the ML approach to single particle images has been described previously and employs the iterative expectation maximization algorithm to maximize the likelihood function (Dempster et al., 1977). Since the images we will process as single particles will be centered on unit cells excised from a 2D crystal, it is reasonable to limit the alignment of each particle to an in-plane angle and translation. Individual unit cells may also suffer from out-of-plane tilts, for example due to undulations in the crystal (cryo-crinkling, Vonck, 2000). To perform a full 3D alignment, two additional angles would have to be considered. The computational load associated with the expectation maximization algorithm would be very high in this case and some approximations have to be made (Scheres et al., 2005). However, application of the ML approach to experimental data from 2D crystals (see below), shows that very good results can also be obtained when limiting the alignment to in-plane transformations. We will, therefore, limit our discussion to this case and follow Sigworth’s implementation of the maximization algorithm (Sigworth, 1998). Briefly, we assume that we have a set of N images X = {Xi;i = 1…N} and corresponding transformation parameters Φ = {[var phi]i;i = 1…N} describing how these images are related to the underlying structure A. We further assume that the noise in each image follows approximately a Gaussian distribution and is uncorrelated. This assumption will be discussed further below when considering experimentally observed data. We can write for each image i


where R = {Ri;i = 1…N} are independent Gaussian noise images, and σ is the standard deviation of the noise. The negative sign in front of the transformation [var phi]i indicates that it describes the transformation needed to bring image Xi into register with structure A. A transformation [var phi]i consists of [var phi]i = (qxi,qyi,qαi) where (qxi,qyi) describes an in-plane translation and qαi a rotation applied to each image. Furthermore, we assume that each image is approximately centered on a point (ζx,ζy), and that the deviation from this point can be described by a Gaussian distribution with standard deviation (σx,σy) (for isotropic lattice distortions, σx = σy). For images of single particles, the distribution of in-plane rotations is usually assumed to be uniform. The in-plane angular distribution of unit cells within a 2D crystal, on the other hand, will be quite narrow. It will be represented by a Gaussian distribution centered on the value ζα with a standard deviation σα. The goal of the ML approach is to maximize the probability P(X|Θ) that the set of images X is observed assuming a set of model parameters Θ = (A, σ,ζx,ζy,ζα,σx,σy,σα). Instead of maximizing P(X | Θ), its logarithm L(X | Θ) is maximized. Assuming independent variables associated with each image, P(X | Θ) is the product of probabilities for each image Xi, and we can write for the logarithm


The probabilities P(Xi | Θ) can be written as integration over the product of the probability density P(Xi | [var phi],Θ) for each image given the model parameters Θ and a transformation [var phi], and of the probability density f([var phi] | Θ) for a transformation [var phi] given the model parameters Θ:




Using the assumptions made above for the noise in the images and the distribution of transformation parameters, we can now write for the probability density f([var phi] | Θ) for the transformation parameters




(Sigworth, 1998). M is the number of pixels in an image. By demanding that the derivative of L(X | Θ) with respect to each model parameter vanishes, an iterative expectation maximization procedure can be established for obtaining improved estimates of each model parameter, given the parameters Θ(n) (Sigworth, 1998):


Noise whitening

The formulation in the previous section does not take into account the CTF of the electron microscope and is based on the assumption of independent Gaussian noise. This implies independence of neighboring pixels, an assumption that does not hold in many experimental situations (Scheres et al., 2005). Uncorrelated noise gives rise to a flat power spectrum, which is usually referred to as a white spectrum. As an example, we calculate the power spectrum for the image stack of unit cells from a 2D crystal of bacteriorhodopsin. The images in the stack were padded to four times their original size to achieve finer sampling of the power spectra, and their power spectra were added for averaging. The average spectrum is shown in Figure 1, together with a line scan. The power spectrum exhibits Thon rings (Thon, 1971) and an amplitude falloff and therefore, the assumption of independent noise does not hold. A noise whitening procedure is needed before maximum likelihood estimation can proceed. For single particles, a noise whitening filter has been constructed from circularly-averaged power spectra of a pure noise region (Sigworth, 2004). Based on the noise model described below, we obtained a similar whitening filter for 2D crystals.

Figure 1
(A) Averaged power spectra from 10,029 tiles of a cryo-EM image of a bacteriorhodopsin 2D crystal. Pixel size was 1.15Å and tile size 100 pixels. Tiles were padded into 400 pixel wide squares, Fourier transformed, and the intensities were averaged. ...

The power spectrum in Figure 1 can be described as a composition of four components: The first arises from the crystal and can be identified as the spot pattern visible in the spectrum (Figure 1A). This will be referred to as the signal S since it contains the information that we are trying to extract with the maximum likelihood procedure. The second component is the ring pattern (Thon rings) that shows oscillations characteristic for the CTF (Figure 1A and 1B). This will be referred to as the image background Bi since it is originating from the entire sample including the non-crystalline parts, such as the carbon support film and/or ice. The third component arises from the shot noise of the electrons forming the image. It lacks the CTF modulation but is still affected by a detector modulation transfer function (MTF). This will be referred to as sample-independent background Bs. The fourth component Be is the MTF-independent background, arising from detector noise (film scanner electronics, or CCD camera electronics). The detector-convoluted Bs together with Be contribute to the background, indicated by the red line in Figure 1B. We can write for the power in our spectrum


k is a point in the power spectrum in Figure 1A. The envelope function ENVi contains all the factors affecting the optical properties of the electron microscope, such as spatial and temporal coherence, as well as any reduction in contrast due to sample movement and charging (see, for example, Zhu et al., 1997). Assuming independence of S and Bi, we can write for the power spectrum of the background in the image


where CTF is given by


with (Scherzer, 1948)


k is the magnitude of k. The two constants w1 and w2 are given by the percentage of amplitude contrast, W, in the image:


Usually, the value for W ranges from 0.03 (Toyoshima and Unwin, 1988; Toyoshima et al., 1993) to 0.14 (Smith and Langmore, 1992) for proteins embedded in ice, and 0.19 (Zhu and Frank, 1994) to 0.35 (Erickson and Klug, 1971) for proteins embedded in stain (uranyl acetate). λ is the electron wavelength and Cs is the spherical aberration coefficient of the objective lens. The defocus Δf is given by (Henderson et al., 1986; Mindell and Grigorieff, 2003).


where DF1 and DF2 are the two defocus values describing the defocus in two perpendicular directions in an image when astigmatism is present, αast gives the angle between the first direction (described by DF1) and the x-axis, and αk is the angle between the direction of the scattering vector k and the x-axis. For a particular image, DF1, DF2 and αast can be determined using a fitting procedure that fits a calculated CTF to the observed Thon rings (Mindell and Grigorieff, 2003). Note that in Eqs. (14) and (16) a positive value for the defocus indicates an underfocus.

The power spectrum modulation is mainly due to the first term in the sum in Eq. (12), which includes the oscillations of the CTF. Therefore, a whitening filter Hw is obtained by averaging the power spectrum along the path of the thon rings


where k1(ki) = k1(k) indicates the pixels ki in reciprocal Fourier space that are on the same Thon ring as pixel k. N is the number of Fourier pixels along one Thon ring that are averaged. k1(k) is a function that is constant for all k lying on the same Thon ring. In the case of no astigmatism it simply returns the magnitude of k. For the more general case, k1(k) is given by


where k2 is the resolution of the pixel k, (i.e.: k2 = |k|). Eq. (18) is derived in the Appendix.

During calculation of equation (17), Fourier pixels are excluded from the averaging if their power deviates from the average by more than the standard deviation of the power. The latter was calculated from the Fourier pixels for which k1(k) is a constant value, and this value was refined by four times iteratively excluding pixels with an amplitude exceeding the standard deviation of the included pixels. Therefore, strong reflections from the crystal that contain predominantly signal, are mostly excluded from the averaging and the whitening filter reflects mainly the background modulation. The Fourier transform of the background with whitened power spectrum is


where F(k) is the pre-whitening Fourier transform. Eq. (19) includes a multiplication by Sign[CTF(k)] which equals +1 for CTF≥0 and −1 for CTF<0. The multiplication does not change the power spectrum. However, it performs a phase flip according to the CTF that is important when points with different values of the CTF are being averaged. This situation can arise, for example, when there is astigmatism present in an image. Any rotation of an image relative to another image from the particle stack would then result in the averaging of image terms that do not have the same CTF value.

Figure 2 shows the average power spectrum after whitening. As the line scan in Figure 2B shows, the power is now approximately constant. Figure 3 shows that the distribution of densities in both the original the whitened images closely match a Gaussian profile. This shows that the whitening does not significantly alter the density distribution of pixels in an image. The background in the whitened image can now be treated as uncorrelated Gaussian noise and, therefore, satisfies the condition of the likelihood formalism outlined above. Apart from the background term, the whitened image also contains the crystal. However, the SNR of the whitened image is only about 1/170 (see below), making the contribution of the crystal negligible.

Figure 2
(A) Average power spectrum as in Figure 1, but calculated from the noise-whitened image. (B) Average power along the red line, as in (A).
Figure 3
Histogram of the grey-values in the noise-whitened image (red line) that was used for Figure 2. For comparison, the histogram of the same image before whitening is shown (black dashed line). After appropriate scaling, the two histograms are virtually ...

The whitening process amplifies the noise level, particularly at high resolution and near CTF zeros. Therefore, the algorithm presented here requires a larger number of particles to obtain the same overall SNR in the final reconstruction, compared with other algorithms that do not require whitening.

Contrast transfer function

The ML procedure described above can be applied to the stack of unit cells excised from the 2D crystal after whitening. The whitened stack is not suitable for determination of the projection structure of the crystal, however, since it does not correct for the CTF of the electron microscope. To obtain a CTF-corrected projection structure of the crystal, a second stack of unit cells has to be generated that contains unaltered amplitudes and phases. Upon convergence of the ML procedure for the stack with the whitened crystal data, an ML estimate of the unaltered crystal can be calculated using Eq. (7) by replacing the data set X with CTF-corrected data set X derived from the unaltered crystal. In evaluating Eq. (7), it is important that the probability function γ is calculated using the whitened images.

If S is the signal component in the ML estimate of the unaltered crystal, an estimate of the crystal structure S can be obtained by “dividing” by the CTF. However, because the CTF contains regions of small values and zeros, a straight division is not possible. It is customary in these cases to use the well-known Wiener inverse filter (Freiden, 1975). S is then given by

S(k)={S(k)SNR(k)CTF(k)[1+SNR(k)]for CTF(k)00for CTF(k)=0

where SNR(k) is the SNR of S at point k. The SNR is usually not known at every point in the spectrum. However, an estimate of the average SNR can be obtained from the crystal by comparing the integrated intensity of all diffraction peaks in a power spectrum of the crystal with the integrated intensity of all other points in the spectrum. In a typical low-dose image, such as the image of the bacteriorhodopsin 2D crystal considered below, the average SNR before averaging is about 1/50. After averaging over N unit cells in a stack (for example N=10029 unit cells from the bacteriorhodopsin crystal, see below) the SNR in the average is N times the SNR of the raw image (in our example the average SNR is about 200). If we assume that the SNR is approximately proportional to CTF2,


the average SNR within the considered resolution range, SNR¯, will be α·CTF2¯. Here, α is a proportionality constant. It follows for the SNR in the ML estimate from N images that


Using Eq. (22) we can write for the structure S




Eqs. (25) and (26) say that we can correct the observed density S by dividing by the CTF at most points. S is multiplied by the CTF only where CTF2 becomes small. The value for CTF2¯ will be close to 0.5 if the CTF is not significantly attenuated by the envelope ENVi. For a significantly attenuated CTF, the value for the average will be smaller. For example, if the envelope function at the high-resolution end of the considered resolution range attenuates the power spectrum by 90%, the value for CTF2¯ is only 0.15. For large values of SNR¯, for example 200 (see above), the action of the filter (Eq. (26)) will not depend strongly on the precise value of CTF2¯. Therefore, in the following calculations, we will approximate CTF2¯ by its upper bound of 0.5.

If the defocus does not change across the image (untilted specimen), the CTF correction factor Ω can be applied to the original image of the 2D crystal before excising the unit cells. Applying the CTF correction to the entire image or at least to larger patches of the original image before excising the unit cells has the advantage that the CTF-dependent point-spread-function (PSF) in the image is not truncated and can be completely included for the CTF correction (Rosenthal and Henderson, 2003; Philippsen et al., 2006).

A further factor affecting the contrast is the resolution-dependent amplitude falloff. To restore the signal at higher spatial frequencies, the amplitudes can be scaled by an exponential function


where a negative temperature factor B determines the amount of amplification. Note that Eq. (27) uses the definition of the temperature factor commonly found in X-ray crystallography. To relate this to the Debye–Waller temperature factor, it has to be divided by 4.

Resolution determination

In single particle image processing, the resolution of the final map can be estimated with the help of the Fourier Ring Correlation (FRC) function between the Fourier transformations of two reconstructions that are made from two subsets of the dataset. As described in Stewart and Grigorieff, 2004, this bears the risk of noise correlation, when image data of low SNR are aligned to the same reference with an algorithm suffering from strong reference bias. Therefore, we calculate the FRC between two structures obtained by independent alignment cycles with two different references. As a resolution cut-off, we used FRC = 0.33 to match the crystallographic criterion. A FRC of 0.33 corresponds roughly to a phase error of 45° (Rosenthal and Henderson, 2003). The crystallographic resolution was determined by the IQ-weighted phase residual determined during symmetry averaging, also with a cut-off of 45°.

When applied to averages calculated from a 2D crystal image, the resolution determination by direct FRC generally fails: The Fourier transformation underlying the FRC calculation assumes periodic boundary conditions in the image. As the unit cell stack generated from a 2D crystal image generally does not have window dimensions and orientations that are exactly aligned with the crystal unit cells, the periodic boundary conditions are not met. The discontinuities at the boundaries correspond to strong and sharp contrast features that will lead to spurious correlation in the FRC curve.

We have therefore opted for resolution determination of the final image after application of a circular mask with a smooth edge. In addition, our software extracts the precise crystal unit cell of the final map and re-interpolates this onto a 200 × 200 pixel wide square image, which is then Fourier transformed, and the amplitudes and phases of the Fourier pixels are extracted, following the procedure by Sherman et al., 1998. These values for amplitudes and phases are then directly compatible with the 2D crystallography package of the MRC image-processing suite, and can be used for further processing with the MRC software. Specifically, Fourier merging and lattice line interpolation performed by the MRC software for 3D reconstructions from differently tilted 2D crystal images is made possible.


Crystallographic image processing with the MRC software programs includes the program QUADSEARCH (Crowther et al., 1996), which determines a coordinate list that describes the positions of the unit cells of the 2D crystal. This list is subsequently interpolated into a much finer raster by the program CCUNBEND, which corrects crystal distortions in the real-space image. The corrected image is Fourier transformed and the amplitude and phase values of the diffraction peaks are evaluated. In QUADSERCH, the unit cell positions are found by searching for maxima in a cross-correlation map, calculated by correlating a reference image (extracted from a filtered version of the original image) with the raw image. However, instead of simply determining the locations of global maxima in this cross-correlation map, QUADSERCH can make use of the a-priori knowledge of the approximate crystal unit cell location, and predicts the location of the next unit cell by extrapolating from the last n (usually 6 or 7) identified unit cell locations. Therefore, QUADSERCH only performs a local peak search of the cross-correlation map in the neighborhood of the predicted location. This reduces the risk of mistaking a spurious noise peak for the correct correlation peak indicating the unit cell location. QUADSERCH allows also refinement of the predicted unit cell location table, by extrapolating a specific new unit cell location from the previously determined raster of neighboring unit cell locations, approaching the location in question simultaneously from different sides. This algorithm enables QUADSERCH to find the correct crystal unit cell location with a much higher precision than a simple cross-correlation peak search over the entire unit cell would allow. This allows the usage of a smaller reference particle size, since the prediction algorithm in QUADSERCH can find a correct local peak, even if this peak is only a weak local maximum within the noisy landscape of the cross-correlation map.

A single-particle approach to the processing of images of 2D crystals would have a disadvantage over the crystallographic approach if the unit cells were selected based only on the maxima in the cross-correlation map between reference and 2D crystal image, instead of profiting from the knowledge of the approximately predicted location of the next unit cell. We decided, therefore, to use the list of unit cell locations determined by QUADSERCH. The more accurate centering of the unit cells extracted at the determined locations will lead to a narrow distribution in Eq. (5) (small σxy), which acts as a significant constraint in the likelihood maximization. The different steps of processing including QUADSERCH and ML analysis are depicted in Figure 4.

Figure 4
Flow-chart of the implemented ML algorithm.

Using the environment of the 2dx software (Gipson et al., 2007a), we utilize the peak position profile determined by QUADSERCH to define the unit cell locations for further processing. The 2dx_image standard-script “Maximum Likelihood” was created, which allows the application of the ML algorithm during the standard image processing sequence. This script calls the program 2dx_ML.exe, which performs the ML calculations. All required parameters and input images for the ML processing are transmitted from the 2dx_image environment, which provides usable default values and offers an interactive help function in form of program-internal manuals, online-documentation and right-mouse-click activated pop-up information windows for the user. The results of the ML processing are channeled back into the 2dx_image graphical user interface (GUI) for further analysis or parameter refinement.

The program 2dx_ML.exe performs noise-whitening according to Eq. (19) on the entire input image, and then windows the particle stack. For the particle selection, the user has the choice of either accepting particles with a QUADSERCH profile correlation value above a given threshold, or explicitly defining the percentage of QUADSERCH particle locations (e.g. 90%) that should be chosen for usage in the ML processing. The program then estimates the ML parameters according to Eq. (46, 810), and iteratively calculates the next generation of estimates structure A(n+1), according to Eq. (7), with correction of the CTF according to Eqs. (25, 26).

For the first estimate for the iterative processing the user can choose either the average of the noise-whitened particle stack, a computer-generated random noise image, or a randomly chosen unit cell. Upper limits for the translation in form of a circular mask and for the rotational angle search range for the ML profiles can be defined. The iterative calculation of the next estimate can be done by a ML calculation with a weighting scheme according to Eq. (7). The SNR of the data can be increased by applying a 2, 3, 4, or 6-fold rotational symmetry to the calculated estimate, if such symmetry is expected to be present. Application of rotational symmetry automatically centers the result onto the symmetry center. An experimental option to apply a low-pass filter using either sharp cut-off or Gaussian profile to the estimated structure after each ML iteration is available. Application of a low-pass filter during the initial iterations can help the overall convergence of the ML procedure in some cases. The iterations are terminated when either a pre-defined number of iterations is reached or the change of the ML parameters is below a given threshold. For the final reconstruction a negative B-factor can be applied (Eq. 27) to amplify high-resolution terms. The final result is output as an image in MRC format, and can be inspected in the 2dx_image GUI. The final result is also translated into a list of amplitudes and phases following Sherman et al., 1998, for further processing in the 2dx_merge program (Gipson et al., 2007b).

For comparison with the ML results, a cross-correlation (CC) based single particle algorithm is also implemented. This CC algorithm prepares a set of rotated reference, and calculates the cross-correlation map between the set of rotated references and the particles. A peak search over all cross-correlation maps for each particle is then performed. The particles are then rotated and shifted according to the CC peak location, and subsequently averaged to produce a new reference. As for the ML algorithm, the result is then translated back into amplitudes and phases. The user can choose in the 2dx_image GUI between either the ML or the CC algorithm.


The ML algorithm was first applied to an image of a negatively stained OmpF 2D crystal recorded with a CCD camera. Figure 5 shows the results using the ML algorithm. After noise whitening and five ML iterations the OmpF trimer became visible (Figure 5A). However, further ML iterations tracked the hexagonal background pattern from the CCD camera’s fiber taper, superimposed on the crystal structure (Figure 5B, C). This is due to the enhancement of the hexagonal pattern in the power spectrum. Even though in the raw image and in the calculated power spectrum no signature of the CCD camera’s fiber optics background was observable, the hexagonal pattern became visible in the power spectrum of the whitened image.

Figure 5
ML processing of the image of a negatively stained 2D crystal of QmpF. (A–C) Intermediate maps obtained after noise whitening and 5, 15 and 20 iterations, respectively. (D–F) ML processing of the same data when using the modified whitening ...

The ML algorithm was therefore extended by a function that excludes Fourier pixels from whitening if their resolution was higher than half the Nyquist frequency and if their amplitude was higher than eight times the local average amplitude. The identified pixels were instead set to an amplitude of one, while the whitening algorithm sets the average of the amplitudes of the remaining pixels to one (see Eq. 17 and 19). This function prevents the signal in the Fourier transformation that originates from the fiber optics pattern at high frequencies from dominating the ML estimate. The modified whitening procedure removed the artifacts seen in the OmpF structure (Figure 5D–F).

We compared results using the new ML procedure with results using another algorithm based on linear cross-correlation alignment, and crystallographic image processing as implemented in the MRC software. These three procedures were applied to 2D crystal images of different SNR and crystalline order: We processed 1) a well-ordered bacteriorhodopsin image (frozen hydrated, image recorded on photographic film and digitized, with an SNR of 1/170, Grigorieff et al., 1995), 2) a partially disordered 2D crystal of the ammonium transporter AmtB (negatively stained, CCD image with an SNR of 1/80, Khademi et al., 2004), and 3) a semi-crystalline membrane patch of the tetrameric potassium channel MloK1 (negatively stained, CCD image with an SNR of 1/30, Chiu et al., 2007). To avoid amplification of noise, whitening was not applied in the CC alignment.

For the well-ordered 2D crystal image of bacteriorhodopsin (bR, Figure 6), 3-fold symmetrization was applied during processing with each method. All processing was done with only one image of bR. For this image, as well as the images of other samples processed for this work, the effect of potential beam tilt was ignored. Using the FRC = 0.33 criterion, the ML approach yielded a resolution of 3.2Å, slightly better than that of the CC alignment (3.4 Å). The crystallographic approach produced results at a similar resolution (3.5Å, IQ-weighted phase residual during symmetrization < 45°). The FRC curves in Figure 6C show the correlation between two structures obtained by two separate ML and CC processes, using two independent reference structures. The FRC curves indicate slightly higher resolution of the ML. The ML process converged after 12 iterations with a translational standard deviation of 0.5 pixels and an angular standard deviation of 0.57°. When the alignment was done with the entire dataset and one common reference, and the FRC was calculated using two subsets of the data, the cross-correlation (CC) alignment showed a slightly higher FRC at high resolution (result not shown). This may have been due to the correlation of noise features between the two structures obtained from subsets of the data using a single reference (Grigorieff, 2000).

Figure 6
Application of the ML (A) and CC (B) algorithms to a frozen hydrated 2D crystal image of bacteriorhodopsin. 3-fold symmetry was applied during processing. (C) The Fourier ring correlation between two half-data sets of the final maps. (D) Image processing ...

ML processing outperformed the other two methods when applied to images of partially disordered 2D crystals of AmtB (Figure 7). The ML reconstruction showed a resolution of 12Å (FRC = 0.33). The CC alignment resulted in a visually lower-resolution map, but similar resolution at the FRC cut-off of 0.33. The ML process converged after 23 iterations with a translational standard deviation of 1.1 pixels and an angular standard deviation of 1.35°. The crystallographic processing of the same micrograph with the MRC software resulted in a projection map at 18Å resolution (IQ-weighted phase residual during symmetrization < 45°).

Figure 7
Processing of a negatively stained, partially disordered 2D crystal image of the ammonium transporter AmtB. Final maps with 2-fold averaging using the ML (A) and the CC (B) alignment approach. (C) The Fourier ring correlation between structures calculated ...

When applied to a badly ordered 2D crystal of MloK1, the single particle approaches achieved significantly better results than the crystallographic approach (Figure 8). Both the ML and CC algorithms produced structures at a resolution of about 20Å. The similar performance might be due to the high SNR, compared to the other images, produced by the high defocus used when recording the image. At high SNR values, the ML and CC algorithms tend to have similar performance (Sigworth, 1998; Scheres et al., 2005). The ML process converged after 25 iterations with a translational deviation of 3.8 pixels. The unit cells were rotated at 90° steps to detect the deviations from four-fold symmetry. The crystallographic processing of the same micrograph resulted in a blurred reconstruction of 28 Å resolution (IQ-weighted phase residual during symmetrization < 45°).

Figure 8
Image processing of a negatively stained, badly ordered 2D crystal image of a tetrameric potassium channel MloK1. (A–B) Final maps with 4-fold averaging using the ML (A), and CC (B) alignment approach. (C) The Fourier ring correlation between ...

Noise correlation produced by the ML and CC alignment methods

A major advantage of the ML approach, compared with CC alignment, is the larger convergence radius. This means that the final ML estimate is less likely to be biased by the initial reference (Sigworth, 1998; Scheres et al., 2005). Another problem often found with CC alignment is so-called over-fitting, in which the refinement procedure introduces features in the refined structure that reflect the noise in the data, rather than the signal (Grigorieff, 2000; Stewart and Grigorieff, 2004). A hallmark of over-fitting is a FRC curve that shows artificially high values at high resolution where the SNR is low. Following Stewart and Grigorieff, 2004, we designed a test to detect over-fitting in the ML and CC approach. We constructed a test dataset on the computer, consisting of 2000 images that were generated from a test pattern (Figure 9A) by random translations with a standard deviation of 5 pixels. In-plane rotation was not applied to avoid interpolation errors at high resolution. Each pattern was then overlaid with white Gaussian noise to produce a set of noisy images with a SNR of 1/200, one of which is shown in Figure 9B. Since the signal and noise components of each image were known, a detailed evaluation of the alignment performance could be done.

Figure 9
Computer simulation to detect over-fitting of noise. (A) Test pattern used to generate 2000 noisy images with a SNR of 0.005 (B). Average calculated from the noise component of the test images after one iteration of CC alignment (C), and ML estimation ...

We used the noise-free structure as the initial reference to perform one round of ML estimation, or one round of CC alignment, depending on which procedure was tested. For the ML procedure, the initial value for the standard deviation of translation was 5 pixels, matching the model parameters that were used to generate the simulated data. The performance of each procedure was evaluated by calculating the new ML estimate (or average of the CC aligned images) from only the noise component of the test images. This allowed us to detect any bias produced by the ML or CC alignment approach. A bias is expected to reveal a replica of the reference in the noise average (Stewart and Grigorieff, 2004). The results in Figure 9C, D suggest that the ML procedure exhibits over-fitting to a lesser extent than CC alignment. Over-fitting becomes significantly more apparent in the ML average after 20 iterations (Figure 9F), but still remains less prominent compared with the result from CC alignment (Figure 9E). Especially at high resolution, the FRC between the noise average produced by the ML procedure and the original pattern is much lower compared with the FRC curve calculated for the average from the CC alignment (Figure 9G).

The performance of the ML and CC alignment procedures can also be evaluated by their alignment error


where sxi and syi are the actual offset of the i-th (1 ≤ iN) image to the true reference, and sxi and syi are determined from the coordinates of peaks in the CC map calculated between the i-th image and the final map. The ML alignment showed an error of 2.8 pixels, and the CC alignment an error of 5.2 pixels.


Several recent publications reported data of crystals that diffract to less than 8Å resolution (e.g. Oling et al., 2001; Parcej and Eckhardt-Strelau, 2003; Tahara et al., 2000; Ziegler et al., 2004). In those cases, a single particle approach can lead to a substantial improvement (Tahara et al., 2000). Here we developed and implemented an adaptation of a ML algorithm for computer image processing of 2D crystal images.

The amount of noise correlation in three different alignment methods has been studied by Stewart and Grigorieff, 2004. In their experiments, the phase residual performed worst as a similarity measure, CC alignment performed better, and a newly proposed alignment using a weighted CC performed best. The ability of the new method developed by Stewart and Grigorieff, 2004, to minimize noise correlation is partly due to the weighting function. The ML algorithm can be viewed as using a weighting scheme that is related to the cross correlation profile and the in-plane transformation profile (Eq. 5). The profiles usually result in a frequency low-pass filtration, which presumably renders the ML procedure more resistant to noise over-fitting.

The new ML algorithm uses unit cells centered according to the coordinates written out by the program QUADSERCH. The distribution function describing the in-plane translations of the unit cells converged during the iterative ML processing to a narrow Gaussian-like profile with a width depending on the order of the crystal. In the examples shown here the final profiles had standard deviations of three pixels or less, but did not converge to a single point. The widened profile reflects a difference between the correlation peak locations found by QUADSERCH and unit cell locations estimated by the ML algorithm. The effect of the Gaussian-like profile is similar to the convolution of the reference image with a Gaussian filter during a conventional single-particle cross-correlation alignment. With the narrowing of the Gaussian profile during ML refinement, the influence of high-resolution terms in the images becomes stronger. The dominance of the low-resolution terms at the beginning of the refinement is presumably one of the reasons for the superior convergence radius and reduced initial reference bias of the ML procedure, compared with CC alignment. The Gaussian filter also reduces the over-fitting of noise (see above). This self-adjusting filter is similar to the self-adjusting weighting function used with the weighted CC alignment (Stewart and Grigorieff, 2004).

The rotational distribution for the ML processing is assumed to be Gaussian with a standard deviation that depends on the quality of the crystals, for example +/− 2.5° for relatively well ordered crystals. Some crystal images, however, may be best processed by assuming a uniform distribution of rotations, i.e. equal likelihood for all angles within a specific interval. The newly implemented software allows the specification of an angular search range and step size. A large step size can be used to test for systematic deviations from nominal symmetry. For example, deviations from four-fold symmetry of the tetrameric potassium channel MloK1 can be detected by testing alignments of the unit cells in 90° steps (Chiu et al., 2007). Therefore, a prior knowledge of the distribution of translation and rotation is very important for the ML method, as also been noticed by Scheres et al., 2005.

The whitening procedure is essential for satisfying the current ML theory, and needed for the ML method to work reliably. Further development of the theory to allow also non-white noise distributions may improve the performance of the ML algorithm because it avoids amplification of noise. The described whitening filter was applied to images of entire crystals with subsequent excision of the unit cells, to avoid truncation of the PSF. CTF correction and whitening of smaller areas containing only a small number of unit cells is required for adapting the method to tilted images with different defocus areas (this work is ongoing).


Single particle processing of 2D crystal images has been considered previously for images of 2D crystals (Sherman et al., 1998; Tahara et al., 2000). Here we apply the Maximum Likelihood approach for the first time to 2D crystal images. The correlation of real-space noise of neighboring pixels in the image is eliminated by using a whitening filter. Compared with cross-correlation alignment, the ML method is significantly less sensitive to, but not entirely free of noise-correlation and reference bias. The better performance on data of low SNR makes its application to stacks of images containing unit cells excised from the crystal image possible. The ML approach has its greatest advantage over the crystallographic approach when applied to partially disordered crystals where the disorder would normally limit the attainable resolution. Based on our tests, the ML approach can yield structures at better than 3.5Å resolution, sufficient for their interpretation by atomic models.

The newly developed ML algorithm is available within the software package 2dx (


This work was in part supported by the NSF, grant number MCB-0447860, to XZ and HS, and by the National Institutes of Health, grant U54 GM-074929, to HS. NG gratefully acknowledges financial support by the National Institutes of Health, grant 1 P01 GM-62580, and by a research fellowship from the Humboldt Foundation. We thank H.-T. Chou for providing the image of a 2D membrane crystal of AmtB, and P.-L. Chiu for the MloK1 image.


In the absence of astigmatism, Thon rings are circular and the function k1(k) in Eq. 17 corresponds to the identity function k1(k) = |k|. Astigmatism produces elliptical or hyperbolical Thon rings, in which case averaging of the pixels along a Thon ring ellipse requires identifying all pixels k in reciprocal space that lie on the same Thon ring ellipse or hyperbole. This is done by the function k1(k), which assigns to each reciprocal pixel k the resolution of the corresponding pixel on the same Thon-ring in the direction of the largest defocus, which is defined as the direction of DF1 or by αast in Eq. 16. The function k1(k) can be derived by calculating the phase shift χ(k) (Eq. 14) for the Fourier pixel k in question and setting that equal to the phase shift χ(k) for the defocus in direction αast:


This can be transformed into


giving rise to Eq. 18:


The sign of the square root in Eq. A2 can be identified by testing for the non-astigmatic case DF1 = DF2 = Δf (k), which leads to


This is equal to


which resolves to k1 = k2, when the minus sign in Eq. A5 and therefore the plus sign in Eq. A2 is chosen.


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.


  • Chiu PL, Pagel M, Evans JE, Chou HT, Zeng X, Gipson B, Stahlberg H, Nimigean CM. The structure of the prokaryotic cyclic nucleotide-modulated potassium channel MloK1 at 16Å resolution. Structure. 2007 in press. [PMC free article] [PubMed]
  • Crowther RA, Henderson R, Smith JM. MRC image processing programs. J Struct Biol. 1996;116:9–16. [PubMed]
  • Dempster A, Laird N, Rubin D. Maximum likelihood from incomplete data via the EM algorithm. J of the Royal Statis Society. 1977;Series B:1–38.
  • Erickson HP, Klug A. Measurement and compensation of defocussing and aberrations by fourier processing of electron micrographs. Phil Trans R Soc. 1971;261B:105–118.
  • Freiden BR. In: Topics in Applied Physics. Huang TS, editor. Springer; Berlin: 1975. p. 177.
  • Gipson B, Zeng X, Zhang ZY, Stahlberg H. 2dx-User-friendly image processing for 2D crystals. J Struct Biol. 2007a;157:64–72. [PubMed]
  • Gipson B, Zeng X, Wouts R, Kühlbrandt W, Stahlberg H. 2dx_merge - User-friendly merging for 2D crystals. J Struct Biol. 2007b this issue.
  • Grigorieff N. Resolution measurement in structures derived from single particles. Acta Cryst D Biol Crystallogr. 2000;56(Pt 10):1270–7. [PubMed]
  • Grigorieff N, Beckmann E, Zemlin F. Lipid location in deoxycholate-treated purple membrane at 2.6 A. J Mol Biol. 1995;254:404–15. [PubMed]
  • Henderson R, Unwin PN. Three-dimensional model of purple membrane obtained by electron microscopy. Nature. 1975;257:28–32. [PubMed]
  • Henderson R, Baldwin JM, Downing KH, Lepault J, Zemlin F. Structure of purple membrane from Halobacterium halobium: recording, measurement and evaluation of electron micrographs at 3.5 Å resolution. Ultramic. 1986;19:147–178.
  • Henderson R, Baldwin JM, Ceska TA, Zemlin F, Beckmann E, Downing KH. Model for the structure of Bacteriorhodopsin based on high-resolution electron cryo-microscopy. J Mol Biol. 1990;213:899–929. [PubMed]
  • Khademi S, O’onnell J, 3rd, Remis J, Robles-Colmenares Y, Miercke LJ, Stroud RM. Mechanism of ammonia transport by Amt/MEP/Rh: structure of AmtB at 1.35 Å Science. 2004;305:1587–94. [PubMed]
  • Kunji ER, von Gronau S, Oesterhelt D, Henderson R. The three-dimensional structure of halorhodopsin to 5 Å by electron crystallography: A new unbending procedure for two-dimensional crystals by using a global reference structure. Proc Natl Acad Sci USA. 2000;97:4637–4642. [PubMed]
  • Mindell JA, Grigorieff N. Accurate determination of local defocus and specimen tilt in electron microscopy. J Struct Biol. 2003;142:334–47. [PubMed]
  • Oling F, Bergsma-Schutter W, Brisson A. Trimers, dimers of trimers, and trimers of trimers are common building blocks of annexin a5 two-dimensional crystals. J Struct Biol. 2001;133:55–63. [PubMed]
  • Parcej DN, Eckhardt-Strelau L. Structural characterisation of neuronal voltage-sensitive K+ channels heterologously expressed in Pichia pastoris. J Mol Biol. 2003;333:103–16. [PubMed]
  • Philippsen A, Engel HA, Engel A. The contrast-imaging function for tilted specimens. Ultramicroscopy 2006 [PubMed]
  • Renault L, Chou HT, Chiu PL, Hill RM, Zeng X, Gipson B, Zhang ZY, Cheng A, Unger V, Stahlberg H. Milestones in electron crystallography. J Comput Aided Mol Des. 2006;20:519–27. [PMC free article] [PubMed]
  • Rosenthal PB, Henderson R. Optimal determination of particle orientation, absolute hand, and contrast loss in single-particle electron cryomicroscopy. J Mol Biol. 2003;333:721–45. [PubMed]
  • Sass HJ, Buldt G, Beckmann E, Zemlin F, van Heel M, Zeitler E, Rosenbusch JP, Dorset DL, Massalski A. Densely packed beta-structure at the protein-lipid interface of porin is revealed by high-resolution cryo-electron microscopy. J Mol Biol. 1989;209:171–5. [PubMed]
  • Scheres SH, Valle M, Carazo JM. Fast maximum-likelihood refinement of electron microscopy images. Bioinformatics (Oxford, England) 2005;21(Suppl 2):ii243–ii244. [PubMed]
  • Scherzer O. The theoretical resolution limit of the electron microscope. Appl Phys. 1948;20:20–29.
  • Schultz P, Celia H, Riva M, Sentenac A, Oudet P. Three-dimensional model of yeast RNA polymerase I determined by electron microscopy of two-dimensional crystals. EMBO J. 1993;12:2601–7. [PubMed]
  • Sherman MB, Soejima T, Chiu W, van Heel M. Multivariate analysis of single unit cells in electron crystallography. Ultramic. 1998;74:179–99. [PubMed]
  • Sigworth FJ. A maximum-likelihood approach to single-particle image refinement. J Struct Biol. 1998;122:328–39. [PubMed]
  • Sigworth FJ. Classical detection theory and the cryo-EM particle selection problem. J Struct Biol. 2004;145:111–22. [PubMed]
  • Smith MF, Langmore JP. Quantitation of molecular densities by cryo-electron microscopy. Determination of the radial density distribution of tobacco mosaic virus. J Mol Biol. 1992;226:763–74. [PubMed]
  • Stahlberg H, Dubochet J, Vogel H, Ghosh R. Are the light-harvesting I complexes from Rhodospirillum rubrum arranged around the reaction centre in a square geometry? J Mol Biol. 1998;282:819–831. [PubMed]
  • Stewart A, Grigorieff N. Noise bias in the refinement of structures derived from single particles. Ultramic. 2004;102:67–84. [PubMed]
  • Tahara Y, Oshima A, Hirai T, Mitsuoka K, Fujiyoshi Y, Hayashi Y. The 11 A resolution projection map of Na+/K+-ATPase calculated by application of single particle analysis to two-dimensional crystal images. J Electron Microsc (Tokyo) 2000;49:583–7. [PubMed]
  • Thon R. In: Electron microscopy in materials sciences. Valdre U, editor. Academic Press; London: 1971. pp. 571–625.
  • Toyoshima C, Unwin N. Contrast transfer for frozen-hydrated specimens: determination from pairs of defocused images. Ultramic. 1988;25:279–91. [PubMed]
  • Toyoshima C, Sasabe H, Stokes DL. Three-dimensional cryo-electron microscopy of the calcium ion pump in the sarcoplasmic reticulum membrane. Nature. 1993;362:467–71. [PubMed]
  • Vonck J. Parameters affecting specimen flatness of two-dimensional crystals for electron crystallography. Ultramic. 2000;85:123–9. [PubMed]
  • Zhu J, Frank J. Accurate retrieval of transfer function from defocus series. Proceedings of the 13th International Congress on Electron Microscopy; Paris, France. 1994.
  • Zhu J, Penczek PA, Schroder R, Frank J. Three-dimensional reconstruction with contrast transfer function correction from energy-filtered cryoelectron micrographs: procedure and application to the 70S Escherichia coli ribosome. J Struct Biol. 1997;118:197–219. [PubMed]
  • Ziegler C, Morbach S, Schiller D, Kramer R, Tziatzios C, Schubert D, Kühlbrandt W. Projection structure and oligomeric state of the osmoregulated sodium/glycine betaine symporter BetP of Corynebacterium glutamicum. J Mol Biol. 2004;337:1137–47. [PubMed]