Introduction

Transmission electron microscopy of biological specimens at cryogenic temperatures (cryoEM) has become a widely used method to image biomolecules at high resolution, both in solution and within the cell. To retrieve the high-resolution signal, the cryoEM images have to be corrected for the contrast-transfer function (CTF) of the microscope. Common parameters used to describe the CTF include an astigmatic defocus, the spherical aberration of the objective lens, and if appropriate, a phase shift introduced by a phase plate. These parameters are commonly estimated by fitting the Thon ring pattern (Thon, 1971) in the power spectrum of micrographs to a modeled power spectrum. The program CTFFIND4 (Rohou & Grigorieff, 2015) has been developed for this task and the model and conventions to describe the CTF are widely adopted in the field.

A limitation of CTFFIND4 is that it considers the whole imaged sample to be at the same objective defocus, which is a reasonable assumption for flat and thin samples, as is common in single-particle cryoEM. However, the increased thickness of cryoEM samples of cells may introduce additional modulations in the Thon ring pattern (Tichelaar et al., 2020) that can lead to errors in the CTF modeling when not accounted for. Furthermore, samples of cells are often tilted with respect to the optical axis of the microscope, either unintentionally due to thinning methods such as cryogenic focused ion beam (FIB) milling, or intentionally during electron cryo-tomography imaging. In both cases the effects are strongest at high-resolution, where the Thon rings are more tightly spaced.

Here we describe new features of CTFFIND5 that can fit the modulations of the Thon ring patterns and determine sample thickness and tilt using an extended CTF model with additional parameters. This not only increases the fidelity of the fit, as Thon rings at higher resolution can now be fitted reliably, but also gives valuable insight into the geometry of the sample that can aid the experimentalist.

Methods

Tilt estimation algorithm

Tilt estimation in CTFFIND5 follows a strategy that is similar to the implementation in CTFTILT (Mindell & Grigorieff, 2003). The tilt axis direction ϕ and tilt angle θ are determined by fitting Thon ring patterns locally, calculated from 128 x 128 pixel tiles that form a regular grid covering the micrograph (Fig. 1a). In this model, ϕ has a positive value ranging from 0° to 360° to describe the angle to the X-axis of the micrograph. It is assumed that the defocus variation across the sample can be described by a tilted plane. Fits are evaluated using correlation coefficients between modeled CTFs and Thon ring patterns.

Figure 1:

Initially, the micrograph pixel size is adjusted (binned) by Fourier cropping to match the resolution limit of the fit set by the user and the micrograph is cropped to be square in order to speed up computation. A power spectrum is calculated from this binned and cropped image, a smooth background is calculated using a box convolution (Mindell & Grigorieff, 2003) and subtracted, the power spectrum is further binned to the tile size (128 x 128 pixels), and the fit of the tilted Thon ring patterns across the micrograph is initialized by fitting this highly binned power spectrum with a non-astigmatic CTF. This fit is then refined using a two-dimensional CTF with astigmatism. Rough values for the tilt axis and angle are then determined in a systematic search in 10° and 5°, respectively, using the locally fitted Thon ring patterns to score each pair of tilt axis and angle, followed by local refinement of tilt axis, angle and average defocus.

Finally, an average tilt-corrected power spectrum is calculated for diagnostic purposes and to allow the determination of a fit resolution. The tilt correction is designed to remove most of the Thon ring blurring due to the defocus variation across the image. To minimize ring blurring, the power spectrum from each tile is adjusted according to its local average defocus, Δf_average, by magnifying it by a factor m with

Since Δf_local will assume values across the image that are both smaller and larger than Δf_average, m will assume values smaller and larger than 1. The magnification / demagnification of the power spectrum compensates for the contraction / expansion of the Thon rings due to the local defocus change and produces approximately constant Thon ring patterns that can be averaged without losing the pattern (Fig. 1b). The compensation will have a small error if the spherical aberration is not zero. However, this error is sufficiently small to not visibly affect the Thon rings in the average.

Verification of tilt estimation using tilted aquaporin crystals

To test the robustness and accuracy of the new fitting algorithm, the defocus and sample tilts of aquaporin 2D crystals (Murata et al., 2000) were estimated using a search range from 5000 Å to 50000 Å and a 100 Å step, low and high resolution limits of 30 Å to 5 Å, respectively, and a box size for the final power spectrum of 512 pixels. The estimated tilt angle θ and axis direction ϕ were compared with the values obtained by 2D crystallographic processing (Mindell & Grigorieff, 2003).

Verification of tilt estimation using tilt series

Lamellae prepared by FIB milling usually exhibit a pre-tilt with respect to the grid surface due to the stage tilt in the FIB instrument. In the microscope, the direction of this pre-tilt will generally not line up with the goniometer tilt axis. For the alignment of a tomogram recorded from such a lamella, the relative orientation of these two axes will have to be determined, together with the precise amount of pre-tilt. We wrote a new cisTEM (Grant et al., 2018) program, called fit_tilt_model, to read the tilt angles and axes determined for each image in a tomographic tilt series and fit them to a model incorporating a pre-tilt and a single tomographic tilt axis. Using a rotation matrix R₀ to represent the pre-tilt and rotation matrices to represent the tomographic tilt angles and axis read from the microscope, the overall sample orientations are given by

R₀ and are calculated from the tilt angles θ and axes ϕ as

In CTFFIND5, both tilt axis and angle are defined in the clockwise direction, with the angle of the axis measured from the x-axis. This may be different from the definition used by the microscope. To ensure consistency with the widely accepted angular convention in the cryoEM field, all the θ and ϕ used in this manuscript refer to an anti-clockwise direction, with ϕ measured from the x-axis.

Using the tilt information obtained with CTFFIND5, we now have a set of rotation matrices Rⁱ, and together with the rotation matrices read from the microscope, , we can calculate a set of pre-tilt estimates from equation (2). To determine the best overall pre-tilt R₀, we determine the plane-normal vectors of the sample by applying to the vector [0,0,1](z-coordinate along the beam direction), followed by calculating their mean as the normal vector of the best overall pre-tilt estimate.

By calculating the root mean squared deviation of the normal vectors , outliers can be identified and excluded to further refine. The pre-tilt can then be determined as:

To generate more reliable defocus and tilt estimates, the defocus search range and resolution fitting range can be adjusted according to the experimental tilt range and image quality. For our cryoEM samples, the low and high resolution limits were set to 50 Å to 10 Å, respectively, and the defocus search interval was set to be between ± 10000 and ±20000 Å from the nominal defocus set during data collection.

Sample thickness estimation

In CTFFIND5 we implemented a new CTF_t model function, based on the CTF function implemented in CTFFIND4 (Rohou & Grigorieff, 2015) and extended by the formula described by (McMullan et al., 2015):

where χ denotes the phase-shift as a function of the electron wavelength λ, the spatial frequency vector |g|, the objective defocus Δf, the spherical aberration C_s, the additional phase shift Δφ, and the fraction of amplitude contrast ω₂. The modulation of the CTF due to sample thickness t is described by the function ξ:

If a user requests sample thickness estimation, the program will first fit the CTF model function as implemented in CTFFIND4 and the “goodness of fit” resolution will be used as an estimate of the frequency g of the first node of the CTF_t function, with t given by:

If the option “Brute-force 1D fit” is selected, CTFFIND5 will further refine t and Δf by calculating the normalized cross-correlation between the radial average of the power spectrum (corrected for astigmatism, as described in) and CTF_t, searching systematically for the best combination of t in the range of 50-400 nm in 10 nm steps, and Δf in the range of ±200 nm from the previously fitted value, also in 10 nm steps.

Finally, if the option “2D-refinement” is selected, CTFFIND5 will optimize t, Δf₁, Δf₂, and ω using the same conjugate gradient algorithm used in CTFFIND4 and the normalized cross correlation between CTF_t and the 2D power spectrum as a scoring function.

After the optimal values for t and Δf have been obtained the “goodness of fit” crosscorrelation is recalculated using CTF_t, with a frequency window that is 1.5 time larger than in CTFFIND4 to avoid the drop-off in the node regions of CTF_t.

Verification of sample thickness estimation using Lambert-Beer’s law

We used 655 micrographs collected from one lamella of ER-HoxB8 cells (dataset Lamella_EUC 1 from (Elferich et al., 2022)). For each micrograph we calculated , where I was the sum of all pixels in the illuminated area of the movie and I₀ was the average of this sum for 45 micrographs collected over vacuum with the same energy filter settings. This value is expected to have a linear relationship with the thickness of the sample consistent with Lambert-Beer’s law (Yan et al. 2015; Rice et al. 2018):

where k is the apparent mean free path for inelastic scattering.

We then used CTFFIND5 to estimate the thickness t of each micrograph using the “Brute-force 1D fit” and “2D-refinement” setting, low and high resolution limits set to 30 Å and 5 Å, defocus search range set between 500 nm and 5000 nm, and low and high resolution limits for thickness estimation set to 10 Å and 3 Å. We used a “RANSAC” algorithm as implemented by the scikit-learn Python package (Pedregosa et al., 2011) to fit a linear model to the relationship of and t, while rejecting outliers. We then manually inspected every outlier of the model fit and categorized the reason for the discrepancy into “Occluded beam” (either from contamination or the edges of the lamella), “Low image signal” (in most cases exposures containing no cellular features), “Carbon/Platinum”, and “Lipid droplet” (see Fig. 4).

Verification of sample thickness estimation using tomography

Lamellae prepared from ER-HoxB8 cells were imaged using a Titan Krios 300 keV TEM controlled by SerialEM (Mastronarde, 2005). For each dataset an initial exposure was taken with a magnification of 64,000, resulting in a pixel size of 1.6 Å and an exposure of 30 e^-/Å. This was followed by the acquisition of a tilt series at a magnification of 48,000, resulting in a pixel size of 2.087 Å. A total of 35 tilt images at a tilt interval of 3° were collected from - 51° to 51°, relative to the milling angle, using a grouped dose-symmetric scheme (Hagen et al., 2017). The exposure per tilt was 3 e^-/Å, resulting in a total exposure of 105 e^-/Å.

For tomographic reconstruction, tilt movie frame motion correction was performed using SerialEM (Mastronarde, 2005), and tilt series were aligned using the IMOD software package (version 4.11, Mastronarde & Held, 2017). For coarse alignment, a high-frequency cutoff radius of 0.15 was used. A fiducial model was generated using patch tracking with patches of 450 x 450 pixels and a fractional overlap of patches of 0.33 x 0.33. High-tilt frames were omitted while generating the fiducial model. Robust fitting with a tuning factor of 1 was used for fine alignment. After computing the alignment, the fiducial model was edited by removing unreliable patches, and then alignments were re-computed. The edited models with the lowest residual mean errors and standard deviations were used for fine alignment. Tomogram positioning was used to correct the tilt angle offset. Fully aligned stacks were generated with a binning factor of 4, resulting in a tomogram pixel size of 8.3 Å. Tomograms were reconstructed using the SIRT-like filtering option in IMOD (Mastronarde & Held, 2017; Mastronarde, 1997) and manually inspected. The tomograms were back-projected along the y-axis using a homemade script, generating a small set of XZ projections. Thickness measurements on the projected central slides were performed using the display program included with the cisTEM software package (Grant et al., 2018).

CTF correction of medium magnification lamella images

The CTF of the representative medium magnification image with a pixel size of 40 Å was estimated using CTFFIND5 with the following parameters: defocus range: 1,000,000 to 4,000,000 Å; search step 50,000 Å; low and high resolution limits: 400 Å and 80 Å. We then used the program apply_ctf, included with cisTEM, to flip the phases according to the estimated CTF. We furthermore implemented the Wiener like filter described in (Tegunov & Cramer, 2019) in apply_ctf to produce the image shown in Fig. 6d.

Benchmarking CTFFIND5 runtimes

CTFFIND5 runtimes were measured using 3 representative micrographs (Table 2). As a baseline measurement, CTFFIND5 was run without estimation of tilt and sample thickness enabled. Then runtime was measured enabling either one of these option or both. Every test was repeated four times and the average and standard deviation of the last three runs are reported, to minimize the contribution of hard-drive speed. The tests were performed on a single core of an Intel Core i9-12900KF CPU.

Results

Tilt estimation

We tested the defocus correction for the Thon rings on a representative micrograph taken from a cryo-FIB milled lamella. As expected, the correction results in the observation of Thon rings at higher spatial resolution (Fig. 1c). In this example, correcting for the estimated moderate tilt of 12.3° improved the highest resolution at which a reasonable fit could be obtained from 5.9 Å to 4.6 Å. The power spectrum also appears less noisy, which can be attributed to some low-pass filtering that occurs with the interpolation of the Thon ring patterns of individual tiles to perform the defocus correction.

To test the performance of the new CTFFIND5 sample tilt estimation, we used a dataset of images of tilted aquaporin crystals that were also used to benchmark the original CTFTILT implementation (Mindell & Grigorieff, 2003; Murata et al., 2000). Table 1 compares the tilt information of the samples obtained from crystallographic analysis and the estimates obtained using CTFFIND5. Overall, the results of CTFFIND5 agree well with the aquaporin crystals information. The average discrepancy was 1.9° for the tilt axis direction and 1.5° for the tilt angle.

Table 1:

Table 2:

To test whether CTFFIND5 would be able to correctly assign tilt axis and angle for tilt series data, we analyzed two tilt series from different grids of lamellae prepared by cryo-FIB milling from mouse neutrophil-like cells (Elferich et al., 2022). We then plotted the estimated values for tilt axis and angle as a function of nominal stage tilt (Fig. 2). The estimated tilt angle shows a roughly linear relationship with the nominal stage tilt, but since CTFFIND5 reports only positive tilt angles the overall plot has a chevron-shape. The estimated tilt axis angle is approximately constant at high tilts but changes by about 180° at 0° estimated tilt, again due to the convention enforced by CTFFIND5. Notably, in both examples there is an offset of about 20° between nominal and estimated tilts, which is due to the pre-tilt of the specimen caused by FIB-milling at a shallow angle. To quantify and delineate both the tilt axis direction of the microscope and the pre-tilt of the specimen we fit all values to a model as described in Methods (Fig. 2). The fitting resulted in an estimated tilt axis angle of 178.2° and 179.8°, respectively, which is consistent with the SerialEM calibration of 178.4° and 176.3° for the stage tilt axis. The estimated pre-tilt values were 20.6 °and -21.9°, consistent with a FIB-milling angle of 20° and opposite orientation of the grids relative to the milling direction. The pre-tilt axis angles were estimated as 171.8° and 183.8°, which is consistent with the error expected from manually aligning the milling direction when inserting grids into the microscope.

Figure 2:

To estimate the accuracy of the tilt estimation in tilt series, we calculated the mean absolute difference between the tilt and axis-angle estimates and the fitted model, excluding the axis-angle estimates at tilt angles under 5°. For the first tilt series we obtained accuracy estimates of 2.08° and 2.58° for tilt and axis-angles, respectively. In the second tilt series the accuracy estimates were 3.95° and 9.47°. In both cases the accuracy was lower than for the tilted aquaporin crystals, presumably due to the relatively short exposure of each micrograph in the tilt series. However, the substantially higher mean differences in the second tilt series suggest that the accuracy is highly dependent on the quality of the underlying data.

Sample thickness estimation

Even after correcting for sample tilt we found that for FIB-milled samples we often could observe Thon-ring like modulation in the power spectrum at higher resolution than suggested by the goodness of fit estimate (Fig. 3b, top plot). These modulations are out of phase with the predicted modulations, as described by (McMullan et al., 2015) and (Tichelaar et al., 2020). We therefore implemented an extension of the CTF model as described by (McMullan et al., 2015) (Fig. 3a). For some images we found that the thickness could be well estimated by assuming that the goodness of fit resolution estimate obtained using the old model implemented in CTFFIND4 corresponds to the first node in the modulation function, according to Eq. (7). With our new model, estimated CTF parameters were very similar to those from CTFFIND4, but the fit in CTFFIND5 extended to higher resolution (Fig. 3b).

Figure 3:

In other images, mostly with defocus values under 1 μm and with a sample thickness over 200 nm, CTFFIND4 could fit the power spectrum before and after the first node using the old CTF model, with some deviations between the fit and the power spectrum (Fig. 3c).

Fitting the power spectrum with the new model in CTFFIND5 resulted in substantially different estimated CTF parameters and an improved fit, even though the goodness-of-fit estimation did not change. Based on these results we conclude that CTFFIND5 will provide more accurate CTF parameters for images of thick samples, such as those generated from FIB-milling. In addition, the fit provides a direct readout of the specimen thickness, which is important for judging specimen quality and the potential for high-resolution information that can be recovered from these images.

Estimating the accuracy of sample thickness estimation using the Lambert-Beer law on energy filtered data

CryoEM is frequently performed using an energy filter to remove inelastically scattered electrons. The fraction of inelastically scattered electrons can be described by the Lambert-Beer law, which states that the fraction of electrons removed from the image is proportional to the thickness of the sample. The apparent mean free path for electron scattering has been experimentally determined for common cryoEM conditions (Rice et al., 2018). To test whether thickness estimation in CTFFIND5 is consistent with this method we used a dataset of 655 exposures of a lamella of ER-HoxB8 cells collected using the DeCo-LACE approach (Elferich et al., 2022). We used CTFFIND5 to estimate the thickness t of every exposure and plotted − against t (Fig. 4). Fitting the data to a linear model described in Methods (Eq. (8)), we found that 568 out of 655 exposures followed closely a linear relationship with a mean free path k of 317 nm. Manual inspection of images that did not follow this linear relationship revealed that they either contained visible ice contamination, platinum deposits, or they were collected over ice without cellular features and displayed weak Thon rings. The value of κ is consistent the value found by (Rice et al., 2018), even though our dataset was collected without an objective aperture. The x-axis intercept of the linear model was -14.1 nm, meaning that the node position systematically predicts a smaller thickness than predicted by the Lambert-Beer law. This discrepancy is further discussed in the next section. To estimate the accuracy of the sample thickness determined by CTFFIND5 we calculated the mean absolute difference to the linear model, which was 4.8 nm. These data suggest that sample thickness determination using node-fitting is an alternative to using Lambert-Beers law that has the advantage of not relying on the constant κ and the intensity I₀, both of which might not be readily available. Also, the two approaches are complementary as they rely on orthogonal mechanisms.

Figure 4:

Estimating the accuracy of sample thickness estimation using tomography

We used a dataset of seven micrographs collected from lamellae of ER-HoxB8 cells together with tilt series collected afterwards from the same locations to verify the accuracy of the thickness estimates obtained using CTFFIND5. We used CTFFIND5 to estimate the thickness (t_CTFFIND) for every location and compared it with the thickness estimated from the tomogram reconstructed from the tilt series (t_TOMO). We measured t_TOMO by manually estimating the distance between the surfaces of the lamella in three different positions.

When we plotted t_CTFFIND against t_TOMO we found that the values were highly correlated, but t_TOMO was consistently smaller than t_CTFFIND (Fig. 5). A linear fit revealed a slope of 0.95 and a y-axis intercept of 0.12 nm. This means that the CTFFIND5 thickness estimate is on average 1.05x higher than the thickness estimated by tomography. (Tichelaar et al., 2020) also report that estimating the thickness from the CTF nodes resulted in values roughly 1.1x higher than estimated by tomography. The reasons for the systematic discrepancies between thicknesses estimated by CTFFIND5 and estimates based on Lambert-Beer’s law and tomography are unclear, but since they are small and CTFFIND5 estimates lie in between the other two estimates, they will provide comparable information.

Figure 5:

CTF estimation and correction assists biological interpretation of intermediate-magnification lamella images

During data collection of cryoEM data in cells, the operator frequently relies on images taken at low magnification to select areas of interest and establish their biological context. The pixel size of these images is usually about 40 Å, with a defocus of about 200 μm. This produces strong contrast from biological membranes, but can sometimes also lead to substantial fringes near these membranes (Fig. 6a). We found that a simple CTF correction based on CTFFIND defocus estimates obtained from the overview images can reduce these fringes (Fig. 6b). A simple CTF correction can be done using the program apply_ctf, included with cisTEM, by phase flipping according to the fitted CTF (Fig. 6c). However, we found that including a Wiener filter-based amplitude correction describe by (Tegunov & Cramer, 2019) produces a more naturally looking image that might be best suited to recognize cellular features (Fig. 6d).

Figure 6:

CTFFIND5 runtimes

To gauge the ability of CTFFIND5 to provide real time feedback during cryoEM data collection we measured its runtime on three representative micrographs (Table 2). Without estimation of tilt or sample thickness CTFIND5 performed CTF estimation roughly within a second. Estimation of the sample thickness adds roughly half a second to the runtime, therefore allowing CTF estimation within a timeframe comparable to typical exposure times. Estimation of the tilt on the other hand increased runtimes substantially to the order of several minutes, due to the exhaustive search of potential tilts over hundreds of powerspectra. While these runtimes are substantially slower than cryoEM data acquisition, near real time estimation can be achieved by using multiple CPU cores.

Furthermore, optimization of the number of tiles used, better search algorithms, or implementations employing GPUs could increase the speed to the point where real time estimation is more feasible.

Conclusion

The new features implemented in CTFFIND5 improve CTF estimation from the power spectra of cryoEM micrographs where assumptions made in its predecessor, CTFFIND4, namely a thin and untilted sample, do not hold. The tilt of the sample is estimated by fitting the CTF to the power spectra calculated from small patches across the image, similar to other software including CTFTilt (Mindell & Grigorieff, 2003), Ctfplotter (Xiong et al., 2009; Mastronarde, 2024), goCTF (Su, 2019), and Warp (Tegunov & Cramer, 2019). After estimation of the sample tilt a tilt-corrected power spectrum is produced that exhibits stronger Thon rings at higher resolution.

To take into account the modulation of the power spectra by thick samples (Tichelaar et al., 2020; McMullan et al., 2015) we fit a modified CTF model, which increases the resolution of the fitted regions of the spectra and provides a read-out of the sample thickness. While the low exposures (3-5 e^-/A2) typically used in electron cryo-tomography often preclude fitting of sample thickness from power spectra of individual images in the tilt series, we demonstrate that this works reliably for higher exposures (∼ 30 e^-/A2) typically used for 2D template matching (Lucas et al.; Rickgauer et al., 2017) and in-situ single particle analysis (Cheng et al., 2023).

While these improvements are especially relevant for in-situ samples, e.g., prepared by cryo-FIB milling, the analysis of images of purified samples recorded at lower acceleration voltages, e.g., 100 keV (McMullan et al., 2023), may also benefit since thickness-dependent CTF modulations will appear at lower resolution with longer electron wavelengths (see Eq. (6)). Per-micrograph CTF estimation can be followed by per-particle CTF refinement, as implemented in cisTEM (Grant et al., 2018), Relion (Kimanius et al., 2021), or cryoSPARC (Punjani et al., 2017). The improvements of CTFFIND5 will provide better starting values for this refinement, yielding better overall CTF estimation and recovery of high-resolution information during 3D reconstruction.

In summary, the improvements implemented in CTFFIND5 result in more accurate CTF estimation of thick and tilted samples and provide valuable information about the samples to the microscopist.

Data availability

The images of tilted aquaporin crystals were previously published (Murata et al., 2000) and are available at https://grigoriefflab.umassmed.edu/tilted_aquaporin_crystals. The untilted exposures of ER-HOXB8 cells are available at EMPIAR (EMPIAR-11063). The tomograms and tilt series from ER-HOXB8 cells have been deposited to EMDB (EMD-43419, EMD-43420, EMD-43424, EMD-43425, EMD-43427, EMD-43428, EMD-43429) and EMPIAR (EMPIAR-11854), respectively. The source code for CTFFIND5 is available at https://github.com/GrigorieffLab/cisTEM/tree/ctffind5 and binaries for most Linux distributions can be downloaded at https://cistem.org/development.

Acknowledgements

We would like to thank members of the Grigorieff lab for testing CTFFIND5 and helpful discussion. We furthermore thank Benjamin Himes for helpful discussions and code review within cisTEM. We thank Alexis Rohou for comments on the manuscript. L.K. and N.G. gratefully acknowledge funding from the Chan Zuckerberg Initiative, grant #2021-234617 (5022).

Author contributions

Conceptualization: J.E. and N.G.; Data curation: J.E., L.K. and X.Z.; Formal analysis: J.E., L.K., X.Z. and N.G.; Funding acquisition: N.G.; Investigation: J.E., L.K. and X.Z.; Methodology: J.E., L.K. and N.G.; Project administration: J.E. and N.G.; Resources: J.E. and N.G.; Software: J.E., L.K. and N.G.; Supervision: J.E. and N.G.; Validation: J.E., L.K., X.Z. and N.G.; Visualization: J.E., L.K. and X.Z.; Writing – original draft: J.E., L.K., X.Z. and N.G.; Writing - review & editing: J.E., L.K., X.Z. and N.G.;