Introduction to Cryo-EM Refinement.

Formally, the FT of 3D map $$ V$$ to be reconstructed in cryo-EM refinement can be defined as the maximizer of the penalized log marginal likelihood function (3) $$ {\displaystyle \sum _{i=1}^N\log P\left(X_i|V\right)-J\left(x\right)}$$, where $$ P\left(X_i|V\right)$$ is the marginal probability (see SI Appendix, Log marginal likelihood for derivation) of $$ X_i$$ given map $$ V$$ by marginalizing over all hidden variables, such as orientation and translation (pose), and $$ J\left(x\right)$$ is the penalty for $$ x$$, which serves to reduce overfitting and guarantee the feasibility of solution.

As 3D molecular maps are both sparse and smooth, in order to incorporate these priors into refinement, a mathematical formulation for them must be developed. Conventionally, the smoothness of a function is associated with the norm of its gradient, and sparseness is referred to as the number of zeros in the values of function (20). In the following subsections, we will formulate different smoothness priors and reveal their differences. The key equations illustrating the effects of the traditional smoothness restraint and our smoothness restraint are Eqs. 2 and 5, respectively.

Traditional Smoothness Prior.

The traditional method (5) enforces smoothness by applying a quadratic restraint on the magnitudes of FTs based on the assumption that they are distributed according to Gaussian. Since the traditional method is an instance of Wiener filtering (21), the restraint strength depends on the SNR. The 3D map reconstructed by the traditional method can be defined as the maximizer of

where $$ r=\sqrt{h^2+k^2+l^2}$$ is the modulus of the spatial frequency vector of the Fourier coefficient $$ V_{hkl}$$, $$ \overline{N\left(r\right)}$$ is the average of weights $$ N\left(\left[h,\hspace{0.17em}k,\hspace{0.17em}l\right]\right)$$ in Eq. 8, and $$ \text{SNR}\left(r\right)$$ is the SNR estimated by gold-standard FSC (22).

To understand the effect of the smoothness restraint of the traditional method, we consider the role of the restraint in the gradient ascent iteration, which is of the form

where $$ \lambda \left(r\right)=\frac{\overline{N\left(r\right)}}{\text{SNR}\left(r\right)}$$ is the damping weight for $$ V_{hkl}$$, and $$ \eta >0$$ is the learning rate. The gradient of the traditional prior $$ \lambda \left(r\right)V_{hkl}$$ can be viewed as a smoothed map in real space. By convolution theorem, let the inverse FT of $$ \lambda \left(r\right)$$ be $$ K\left({\Vert u\Vert }_2\right)$$, which is a radial function. The gradient of traditional prior $$ x\prime \left(u\right)$$ is a convolution in real space, namely, $$ x^{\text{'}}\left(u\right)={\displaystyle \underset{{\mathrm{ℝ}}^3}{\int }K\left({\Vert u-v\Vert }_2\right)x\left(v\right)dv}$$, which represents a map smoothed by the translation-invariant isotropic kernel $$ K\left({\Vert u-v\Vert }_2\right)$$. At each step, the old solution is modified by a linear combination of the gradient of log marginal likelihood function and the radial-kernel–smoothed old solution. We thus postulate that the traditional method in RELION biases the solution toward the 3D map with homogeneous smoothness across space.

Sparseness and Smoothness Priors in OPUS-SSRI.

Sparseness resembles the idea of masking in the calculation of masked FSC, where the voxels which are below a certain threshold are setting to 0. The similar effects can be achieved by restraining the sum of the absolute values of densities, namely, the $$ l_1$$ norm of the density map (7) during reconstruction. Hence, we can encode the sparseness using the $$ l_1$$ norm (7) and the smoothness of 3D map using the TV norm (8). Though these two priors can effectively guarantee both sparseness and smoothness of the map, they also introduce biases to the final solution (23). In particular, $$ l_1$$ regularization tends to underestimate true nonzero elements (23) since the corresponding soft-threshold operator shrinks the volume by a global threshold. Fan et al. discovered that nonconcave penalty can prevent true nonzero elements from being overly shrunk while preserving sparseness (23). Therefore, we employ a nonconcave penalty log norm (24) to reduce biases in the solution. The penalized log likelihood function in OPUS-SSRI has the form

where $$ \alpha $$ and $$ \beta $$ are positive, and $$ \mathit{ϵ}$$ and $$ \mathit{ϵ}\prime $$ are positive parameters to guard against the singularity of logarithm function at zero.

Optimization Methods in OPUS-SSRI.

This subsection presents the algorithm to optimize the penalized log likelihood in Eq. 3. First, the log marginal likelihood function can be optimized by the expectation–maximization method (25) (see SI Appendix, Expectation maximization for derivation). The reconstruction process alternates between the expectation step in which the distribution of pose parameters for each particle is determined and the maximization step in which the 3D map is reconstructed. Secondly, to address the nonconcavity of log norm, we approximate the logarithm function by concave function and iteratively improve the approximation (24) at each maximization step (see SI Appendix, Weighted approximation for derivation). Lastly, to average 3D maps reconstructed in consecutive maximization steps, we consider leveraging implicit gradient ascent (26), which is a widely used technique to improve the stability of optimization method. The implicit gradient ascent restrains the Euclidean distance between the new solution and the 3D map of previous maximization step $$ x^{k-1}$$. These choices yield a target function in Eq. 4 by optimizing which can improve Eq. 3 (24). Formally, at the $$ i+1$$th iteration of the $$ k$$th maximization step, let the solution obtained in the previous iteration be $$ x^i$$. We approximate Eq. 3 with the expected log likelihood (the first term in Eq. 4) and the weighted $$ l_1$$ and TV norms and define the solution at iteration $$ i+1$$ as the maximizer of the following equation:

where $$ P\left(\varphi |X_i,\hspace{0.17em}{V}_{k-1}\right)$$ is the conditional probability of the pose parameters given the observation $$ X_i$$ and the map $$ V_{k-1}$$ from the previous maximization step, and $$ \gamma >0$$ is the weight of the implicit gradient ascent penalty. The effect of log norm becomes evident after converting it to the weighted norm in Eq. 4. Reweighting each voxel yields a spatially varying threshold (Eq. 6), which can reduce the bias of prior. Next, we will demonstrate that the gradient of TV norm enforces heterogenous smoothness to the 3D map, and the $$ l_1$$ norm achieves sparseness by the soft-thresholding rule while presenting the optimization algorithm.

Though TV norm is nondifferentiable at zero, we can approximate its gradient by Nesterov smoothing (27). The approximate gradient of TV norm at a voxel $$ \left[i,\hspace{0.17em}j,\hspace{0.17em}k\right]$$ (see SI Appendix, Nesterov smoothed TV norm for derivation) is of the form

where $$ {\mathrm{\Delta }}_a$$ is a 1 × 3 vector with 1 on the $$ a$$th entry and zeros elsewhere, $$ c_{\left|i,\hspace{0.17em}j,\hspace{0.17em}k\right|}=3{\Vert \nabla x\left[i,\hspace{0.17em}j,\hspace{0.17em}k\right]\Vert }_{\mu }^{-1}+{\displaystyle \sum _{a=1}^3{\Vert \nabla x\left[\left(i,\hspace{0.17em}j,\hspace{0.17em}k\right)+{\mathrm{\Delta }}_a\right]\Vert }_{\mu }^{-1}}$$, $$ c_{\left(i,j,k\right)-{\mathrm{\Delta }}_a}^{\text{'}}={\Vert \nabla x\left[i,\hspace{0.17em}j,\hspace{0.17em}k\right]\Vert }_{\mu }^{-1}$$ and $$ c_{\left(i,j,k\right)+{\mathrm{\Delta }}_a}^{\text{'}}={\Vert \nabla x\left[\left(i,\hspace{0.17em}j,\hspace{0.17em}k\right)+{\mathrm{\Delta }}_a\right]\Vert }_{\mu }^{-1}$$, where $$ {\Vert \nabla x\left[i,\hspace{0.17em}j,\hspace{0.17em}k\right]\Vert }_{\mu }=\left(\Vert \nabla x_{(i,j,k)}^i\Vert +\mathit{ϵ}\prime \right)\max \left(\Vert \nabla x\left[i,\hspace{0.17em}j,\hspace{0.17em}k\right]\Vert ,\hspace{0.17em}\mu \right)$$, and $$ \mu $$ is the Nesterov smoothing parameter. The gradient of the smoothed TV norm at the voxel $$ \left[i,\hspace{0.17em}j,\hspace{0.17em}k\right]$$ depends on the gradients at this voxel and its neighboring voxels $$ \left[\left(i,\hspace{0.17em}j,\hspace{0.17em}k\right)\pm {\mathrm{\Delta }}_a\right]$$. Specifically, the voxel on the direction with smaller gradient has a larger contribution to the gradient of TV norm. We can write the continuous form of Eq. 5 as $$ x\prime \left(u\right)={\displaystyle \underset{{\mathrm{ℝ}}^3}{\int }K\left(u,\hspace{0.17em}v\right)x\left(v\right)dv}$$, where the kernel $$ K\left(u,\hspace{0.17em}v\right)$$ depends on the voxel $$ u$$ and its neighboring voxel $$ v$$ simultaneously. Therefore, the weighted TV norm smooths the map by a spatially varying anisotropic kernel $$ K\left(u,\hspace{0.17em}v\right)$$, which adapts to heterogenous smoothness in the 3D map.

The differentiable function with $$ l_1$$ penalty can be optimized by a joint application of gradient ascent and soft-thresholding operator (28). Denote the implicit gradient-restrained expected log likelihood $$ -{\displaystyle \sum _{i=1}^{N}P\left(\varphi |X_i,\hspace{0.17em}{V}_{k-1}\right){\Vert X_i-{\text{CTF}}_i{P}^{\varphi }V\Vert }^2-\frac{\gamma }{2}{\Vert x-x^{k-1}\Vert }_2^2}$$ as $$ h\left(x\right)$$and the weight $$ \frac{1}{\left|x_j^i\right|+\mathit{ϵ}}$$ as $$ w_j$$, and let the learning rate be $$ \eta $$ for the $$ j$$th component of $$ x$$ at $$ i+1$$ iteration. Eq. 4 can be optimized by the following equation:

where $$ x^i\text{'}=x^i+\eta \left(\nabla h\left(x^i\right)-\beta \nabla f_{\mu }\left(x^i\right)\right)$$ is the 3D map updated by gradient ascent, and $$ \text{sign}\left(\right)$$ extracts the sign of a voxel value. Eq. 6 is referred to as the soft-thresholding operator. The sparseness of the volume is preserved by the soft-thresholding operator since it sets the voxels with relatively small values to zeros. Moreover, the threshold $$ \eta \alpha w_j$$ is inversely proportional to the strength of the voxel. Therefore, the voxels with small values are more likely to be shrunk to 0 s. Since the values of electron densities for molecules are often higher than the background noises, this soft-thresholding operator can suppress noises while leaving the electron densities of molecules intact. We hence demonstrate how our $$ l_1$$ restraint yields an unbiased cleaner 3D map.

In summary, Eq. 6 is applied iteratively in the maximization step. The gradient of the TV norm enforces spatially varying smoothness in gradient ascent, while the soft-thresholding operator induced by the weighted $$ l_1$$ norm guarantees the sparseness of the 3D map. In all of our experiments, convergences were achieved in a maximization step after 100 iterations.

Back-Projection as Local Kernel Regression.

To reconstruct the 3D map, cryo-EM researchers introduced a back-projection operator, which puts the 2D FT of the image into the 3D map. As the inverse of slice operator $$ P^{\varphi }$$ (SI Appendix, Eq. S2), back-projection puts the data with 2D index $$ \left[i,\hspace{0.17em}j\right]$$ to the 3D index $$ n\left(\varphi \right)$$ according to $$ n\left(\varphi \right)=R_{\varphi }^{-1}{\left[i,\hspace{0.17em}j,\hspace{0.17em}0\right]}^T$$, where $$ R_{\varphi }^{-1}$$ is the inverse of the rotation matrix parameterized by orientation parameter $$ \varphi $$. Back-projection prompts us to create a 3D volume to store the data on nonintegral indices $$ \left[h,\hspace{0.17em}k,\hspace{0.17em}l\right]\in {\mathrm{ℝ}}^3$$ if the orientation $$ \varphi $$ is randomly sampled, while the indices of the 3D map operated in computer memory are integral, namely, $$ \left[h,\hspace{0.17em}k,\hspace{0.17em}l\right]\in {\mathrm{\mathbb{Z}}}^3$$. Therefore, the limitation of computer memory prevents us from storing the voxel of the back-projected 2D slice $$ P^{-\varphi }{X}_i$$ in the 3D map $$ V$$. Nonparametric regression can be leveraged to eliminate this discrepancy by storing weighted back-projected data on integral grids (29). Using local kernel regression (see SI Appendix, Local kernel regression for derivation), the expected log likelihood can be written as

where $$ K$$ is the kernel, $$ n_j\left(\varphi \right)$$ is a voxel back-projected with parameter $$ \varphi $$, and $$ n=\left[h,\hspace{0.17em}k,\hspace{0.17em}l\right]\in {\mathrm{\mathbb{Z}}}^3$$ is the voxel storing the kernel-weighted back-projected data. Hence, in contrast to the first term in Eq. 4, this new formulation decouples the slice operator from $$ V$$ and has variables $$ V_{hkl}$$ on integral grids.

OPUS-SSRI used Gaussian kernel, which is of the form $$ K\left(n_j\left(\varphi \right),\hspace{0.17em}n\right)=\exp \left(-\frac{{\left(h_j-h\right)}^2+{\left(k_j-k\right)}^2+{\left(l_j-l\right)}^2}{{\sigma }^2}\right)$$. We set the bandwidth $$ \sigma $$ of the Gaussian kernel to be $$ \sqrt{2}$$ in our implementation of OPUS-SSRI. For practical consideration, a Fourier coefficient $$ n_j\left(\varphi \right)$$ cannot be scattered to every voxel. Since the Gaussian kernel is quickly diminishing, we let $$ 〚\cdot 〛$$ be the round operator, and we consider storing $$ X_{ij}$$ with the index $$ \left[h_j,\hspace{0.17em}{k}_j,\hspace{0.17em}{l}_j\right]$$ on the neighboring voxels $$ \left[h_j+{\mathrm{\Delta }}_h,\hspace{0.17em}〚k_j〛+{\mathrm{\Delta }}_k,\hspace{0.17em}〚l_j〛+{\mathrm{\Delta }}_l\right]$$, where $$ {\mathrm{\Delta }}_a=\pm 1\hspace{0.17em}\text{or}\hspace{0.17em}0$$. In the original implementation of RELION (5), they chose a kernel with trilinear interpolation–like weight. Using the formalism introduced before, their kernel functions are of the form $$ K\left(n_j\left(\varphi \right),\hspace{0.17em}n\right)=w\left(h_j,\hspace{0.17em}h\right)w\left(k_j,\hspace{0.17em}k\right)w\left(l_j,\hspace{0.17em}l\right)$$ with $$ w\left(x_j,\hspace{0.17em}x\right)=\{\begin{array}{c}\Delta x+1-x_j,\hspace{0.17em}x\le x_j\\ \Delta x_j+1-x,\hspace{0.17em}x>x_j\end{array},$$ and $$ x\in \left[\lfloor x_j\rfloor ,\lceil x_j\rceil \right]$$, where $$ \lfloor \cdot \rfloor $$ and $$ \lceil \cdot \rceil $$ are the floor and ceiling operator, respectively. Hence, a Fourier coefficient is scattered to eight neighboring voxels in the original implementation of RELION (5).

Eq. 7 has a closed form solution,

with $$ N\left(n\right)=\sum _{i=1}^N\sum _{j=1}^J\sum _{\varphi }P\left(\varphi \text{|}{X}_i,\hspace{0.17em}{V}_{k-1}\right)K\left(n_j\left(\varphi \right),\hspace{0.17em}n\right){\text{CTF}}_{ij}^2$$. Ignoring regularization terms, local kernel regression estimates the value of a voxel $$ \left[h,\hspace{0.17em}k,\hspace{0.17em}l\right]$$ by an integral between the kernel $$ K$$ and neighboring experimental data $$ V_{h\prime k\prime l\prime }=P\left(\varphi \text{|}{X}_i,\hspace{0.17em}{V}_{k-1}\right){\text{CTF}}_{ij}{X}_{ij}$$, $$ V_{hkl}=\text{Π}{\displaystyle \underset{h\prime k\prime l\prime }{\int }K\left(n,\hspace{0.17em}n\prime \right)V_{h\prime k\prime l\prime }dh\prime dk\prime dl\prime },$$ where $$ \text{Π}$$ is the normalizing factor. This convolution projects $$ V_{hkl}^{\text{'}}$$ into the RKHS associated with the given kernel according to the representer theorem (9), thus enforcing $$ V_{hkl}$$ to exhibit the smoothness of the RKHS induced by that kernel.

OPUS-SSRI Implementation.

The implementation of OPUS-SSRI is based on RELION. The 3D refinement program in RELION consists of two modules, expectation and maximization. We implemented our method as a new routine in its maximization module. Therefore, when performing refinement, the expectation steps of RELION (5) and OPUS-SSRI use exactly the same settings. The gradient calculation and soft-thresholding operators are implemented with CUDA, thus allowing fast maximization.

Gold-Standard FSC.

The gold-standard FSC is the FSC between two independently refined half maps F and G (22). The gold-standard FSC of Fourier coefficients at shell $$ k$$ relates to the SNR through $$ \text{FSC}\left(k\right)=\frac{\text{SNR}\left(k\right)}{\text{SNR}\left(k\right)+1}$$, as noises inside two maps are independent (30). The frequency where the gold-standard FSC passes through 0.143 is denoted as the estimated resolution of the reconstruction. A high-resolution noise substitution method is often used to correct effects of masking in the masked FSC (30).

Model versus Map FSC.

If there exists a high-resolution atomic structural model, we can validate the cryo-EM map by comparing it to this atomic model. The first step in calculating the model versus map FSC is fitting the atomic model into the cryo-EM density map. The model map is constructed from the fitted atomic model by sampling on the same grid as the experimental map. The model versus map FSC (31) is the correlation between the FT of the model map and the FT of cryo-EM map. The point where the model versus map FSC approaches 0.143 can be regarded as the resolution of the experimental map.

Refinement Protocol.

The single-particle datasets used in this paper were obtained from either the deposited particle stack or the coordinate files. In all experiments, we built the initial maps ab initio in RELION 3.0 and refined those initial maps using the three methods to be compared. The initial map building began with one round of 2D classification in RELION 3.0. The particles belonging to the major classes were then selected to build the initial map ab initio using the 3D classification procedure in RELION 3.0. The same low-pass–filtered initial maps were subsequently supplied into the three methods, RELION 3.0, THUNDER, and OPUS-SSRI, for refinement. For the datasets with specific symmetry, the symmetry was enforced throughout the refinement process. For RELION 3.0 and OPUS-SSRI, we also used the same convergence criteria [i.e., no resolution improvement and pose changes for the last two iterations (5)]. In THUNDER, the particle grading and CTF search options were set as “True” for better results. Finally, the gold-standard FSC calculations and density map postprocessing of the refinement results of all methods were carried out in RELION 3.0. In the postprocessing step, the mask was created from the final reconstruction using all particles in the 3D refinement procedure. Using relion_postprocess (30), we obtained gold-standard FSCs and the postprocessed map from independent maps by correcting the modulation transfer function of the detector and sharpening with automatically estimated B-factors. We then compared the postprocessed map with respect to the corresponding published atomic model(s) by calculating model versus map FSC using Phenix.Mtriage (32). Before comparison, the atomic model was fitted into the postprocessed density maps reconstructed by different methods using the rigid-body fit in Chimera (33).