Electron tomography (ET) offers a unique capacity to image biological structures in situ. However, the resolution of ET reconstructed tomograms is not comparable to that of the single-particle cryo-EM. If many copies of the object of interest are present in the tomograms, their structures can be reconstructed in the tomogram, picked, aligned and averaged to increase the signal-to-noise ratio and improve the resolution, which is known as the subtomogram averaging. To date, the resolution improvement of the subtomogram averaging is still limited because each reconstructed subtomogram is of low reconstruction quality due to the missing wedge issue.

In this article, we propose a novel computational model, the constrained reconstruction model (CRM), to better recover the information from the multiple subtomograms and compensate for the missing wedge issue in each of them. CRM is supposed to produce a refined reconstruction in the final turn of subtomogram averaging after alignment, instead of directly taking the average. We first formulate the averaging method and our CRM as linear systems, and prove that the solution space of CRM is no larger, and in practice much smaller, than that of the averaging method. We then propose a sparse Kaczmarz algorithm to solve the formulated CRM, and further extend the solution to the simultaneous algebraic reconstruction technique (SART). Experimental results demonstrate that CRM can significantly alleviate the missing wedge issue and improve the final reconstruction quality. In addition, our model is robust to the number of images in each tilt series, the tilt range and the noise level.

The codes of CRM-SIRT and CRM-SART are available at https://github.com/icthrm/CRM.

Supplementary data are available at Bioinformatics online.

Electron tomography (ET) is a unique technology to image biological structures, such as an entire cell, in their native environment (Frank, 2006). In ET, the object of interest is tilted for different degrees of rotations and a series of images is taken, which is then used to reconstruct the 3D structure of the object through reconstruction methods, such as back-projection algorithms based on inverse Radon transform. In the ideal case, if the tilt angles cover the full range of projections, a perfect reconstruction can be obtained. However, in practice, the tilt range of the projections is almost always limited within ±40° to ±70°⁠, which causes the missing wedge issue and leads to remarkable artifacts. Compared with other technologies, such as single-particle cryo-EM, this has severely limited the application of ET in high-resolution imaging.

The limited range projection issue is not only a key problem in ET but also widely exists in various fields, such as computed tomography, radial physics, ultrasonic and biomedical imaging (Harris et al., 2017; Shimizu et al., 1998; Wu et al., 2003; Zhu et al., 2005). A large number of studies have been carried out to suppress the effect of missing wedge (Carazo, 1992; Chen and Förster, 2014; Deng et al., 2016; Han et al., 2017; Kovacik et al., 2014; Kupsch et al., 2016; Leary et al., 2013). Almost all of these works attempt to introduce additional constraints representing priori knowledge about the specimen to narrow down the solution space. However, since once the information is lost, it will never be recovered, the priori knowledge is more useful on reducing artifacts than recovering lost information.

If many copies of the object of interest are present in the tomograms, their structures can be reconstructed separately, averaged together to increase the signal-to-noise ratio (SNR) and improve the resolution, which is known as the subtomogram averaging. In subtomogram averaging, the entire tomogram is firstly reconstructed from a series of projections, and then the ultra-structures (sub-volumes) containing the object of interest are picked from the reconstructed tomogram, aligned and averaged to obtain the final subtomogram. This strategy has been shown to effectively increase the SNR and improve the resolution of the subtomograms. However, the reconstruction of each subtomogram still suffers from the missing wedge issue. Additionally, since the object of interest in ET is imaged in the highly noisy native environment, the surrounding background consists of not only noise but also various artifacts and other objects, which poses further challenges to subtomogram averaging.

Subtomogram averaging is a quite new technology with enormous potential (Briggs, 2013). To date, several attempts have been made to increase the reconstruction quality of the subtomograms. The geometry calibration has been extensively studied to provide more accurate averaging operation (Bartesaghi et al., 2012; Bharat and Scheres, 2016; Xu and Alber, 2012; Xu et al., 2012; Zhao et al., 2018). Other works try to change the information interpretation in tomographic reconstruction and improve the consequent subtomogram averaging (Chen and Förster, 2014; Deng et al., 2016; Han et al., 2019; Kunz and Frangakis, 2014). The effect of the missing wedge has also been considered in the process of alignment and averaging, by applying an angle-limited mask to the subtomogram to compensate the ‘missing wedge bias’ (Förster et al., 2005; Galaz-Montoya and Ludtke, 2017; Schmid and Booth, 2008). Despite the contribution of these works, the core idea remains the same, i.e. the final structure is obtained by averaging the structures of all the subtomograms.

However, incorporating the image information of these subtomograms into a consistent structure is a typical information fusion task. The information residing the subtomograms is relatively intermediate compared with the ones underlying these subtomograms’ corresponding projections. A direct access to the information underlying the original projection images of these subtomograms may achieve better data usage and produce clearer final reconstruction of the ultrastructure.

In this article, we propose a novel constrained reconstruction model (CRM) to better recover the information from the multiple subtomograms to compensate for the missing wedge issue in each of them. CRM is designed to produce a refined reconstruction in the final turn of subtomogram averaging after alignment, instead of the traditional averaging. Here, we first formulate the averaging method and our CRM as linear systems, and demonstrate the difference between the two models. We then prove that the solution space of CRM is no larger than that of the traditional method. To solve the linear system of CRM, we propose a sparse Kaczmarz algorithm to solve the formulated CRM, and further extend the solution to the simultaneous algebraic reconstruction technique (SART). Finally, the effectiveness and robustness of CRM have been well studied and proved in both 2D and 3D cases.

Reconstruction algorithms can be categorized into direct Fourier inversion (Penczek et al., 2004; Scheres, 2012), iterative algorithms (Andersen and Kak, 1984; Gilbert, 1972; Gordon et al., 1970) and back-projection methods (De Rosier and Klug, 1968). Here, we give a brief introduction to the algebraic reconstruction technique (ART) (Gordon et al., 1970), which is a class of iterative algorithms that have been widely used in computed tomography (Andersen, 1989; Beister et al., 2012; Prakash et al., 2010).

Although the Kaczmarz’s algorithm has demonstrated success in ET reconstruction, it suffers from issues of numeric stability and computational efficiency. To achieve better estimation and computational efficiency, various methods have been proposed, such as the simultaneous iterative reconstruction technique (SIRT) (Gilbert, 1972) and the SART (Andersen and Kak, 1984). SART is one of the most commonly used reconstruction methods. In contrast to the Kaczmarz’s algorithm, SART only updates the solution once per copy, after solving all the projections of a copy.

Since this work focuses on the final estimation of the ultrastructure in subtomogram averaging, we assume that the projections have been posed in the correct geometry. Without loss of generality and for the simplicity of presentation, we discuss the 2D condition first and then straightforward to generalize techniques to the 3D case.

We first prove that the solution space of CRM is no larger, and in practice much smaller, than that of the traditional method.

In the discrete condition, Radon transform can be reduced to a linear system. Figure 1 illustrates the main difference between the traditional subtomogram averaging method and the proposed CRM. In the traditional method, each subtomogram is reconstructed separately, picked and aligned. Then all of them are averaged together to obtain the final reconstruction. However, the missing wedge issue still exists in each reconstruction and is thus inherited to the final estimation. In contrast, CRM constrains all the subtomograms together, by which it can effectively combine the information from different copies to compensate for the missing wedge issue. We give the formulation of both the traditional method and CRM in linear systems as follows:

Illustration of the difference between the averaging-based method and the proposed constrained reconstruction model (CRM). (A) The idea structures and surroundings supported to be observed. (B) The workflow of averaging-based method. (C) The workflow of constrained reconstruction

For the ith copy, we collect the measured angular projections pi(θij,s),j=1,…,m in a vector bi⁠. Define vector fi=[hiT,giT]T and matrix Ai=[Di,Bi]⁠, where vector hi represents the pixels covered by function hi, vector gi represents the pixels covered by function gi, matrix Di represents the projection coefficients operated on hi and matrix Bi represents the projection coefficients operated on gi⁠. The reconstruction problem for the ith copy is to obtain fi by solving Aifi=bi⁠.

Here, the estimation {hi} is directly replaced by h. Therefore, the reconstruction is based on all the projections simultaneously. We give the following theorem to guarantee the feasibility of efficiently solving the linear system in (6). 

Theorem 1. The solution space of the CRM in (6) is no larger than that of the averaging method in (5).

By this, we have proved that the solution space of CRM is no larger than that of the averaging method.□

Here, we should note that the main difference between dim(Aavg) and dim(ACRM) comes from nhi−rank(Di)⁠. In real-world datasets, the range of tilt angles is limited and the sampling is discrete, which makes nhi≫rank(Di) and results in the missing wedge issue. Therefore, in practice, dim(ACRM) should be much smaller than dim(Aavg) and thus leads to better estimations.

Several methods can be applied to solve the optimization problem, such as least square estimation (LSE) and the Kaczmarz’s algorithm. Because LSE needs to solve the inversion of the entire matrix in Eq. (8), it is computationally costly, especially when a large number of copies exist. On the other hand, Kaczmarz’s algorithm is an iterative algorithm for solving linear equation systems. Thus it is also applicable to Eq. (8) without large changes. However, a direct implementation of the Kaczmarz’s algorithm is still costly in terms of runtime and memory. Here, we first propose a novel sparse Kaczmarz’s algorithm to solve the problem, then extend the solution to SART and finally propose the overall framework to solve the CRM.

Considering the linear system in Eq. (8), let ACRM=[A1,A2,…,AN]T,b¯=[b1,b2,…,bN]T, Di=[Di,1,…,Di,m]T,Bi=[Bi,1,…,Bi,m]T⁠, and bi=[bi,1,…,bi,m]T⁠, where Ai is the ith component of ACRM, bi is the ith component of b¯, Di,l is the lth row of Di, Bi,l is the lth row of Bi⁠, and bi,l is the lth row of bi⁠. Here, we iteratively solve the linear system with the Kaczmarz’s algorithm.

Algorithm 1 outlines the proposed sparse Kaczmarz’s algorithm to solve the CRM. In Algorithm 1, hk is the finally reconstructed subtomogram for the kth iteration by Eq. (17) and gik is the finally calculated background for the ith copy in the kth iteration by Eq. (17). It should be noted that the derivation above does not introduce any additional assumptions or constraints. Therefore, sparse Kaczmarz’s algorithm will give the same solution as the original Kaczmarz’s algorithm.

Sparse Kaczmarz’s algorithm for CRM

Require: measured bi,i=1,…,N, ε0

Ensure: h

 iteration number k←1⁠, copy index i←1

 hk=1←0,g1k=1←0,g2k=1←0,…,gNk=1←0

 hold←0

 while  k≤MAX_ITR  do

  for  i≤N  do

   hk←hold

   update hk, gik by solving Dih+Bigi=bi with Eq. (20)

   hold←hk, i←i+1

  end for

  k←k+1

 end while

 return  h←hold

Comparing to the Kaczmarz’s methods, SART is more efficient and is thus widely applied in ET reconstruction. The main difference between the Kaczmarz’s algorithm and SART is that the former needs to update the solution in each projection image of each copy, whereas the latter only updates the solution once per copy after solving all the projections of a copy.

That is, the object of interest (h) of one copy is used to initialize the other copy, and vice versa. The background (g) of a copy is initialized by the same background from the previous iteration. It is straightforward to generalize the above derivation to a linear system with N copies and to other iterative methods, such as SIRT. Here, we apply the sparse SART in the overall framework to solve CRM. Figure 2 illustrates the main idea of the sparse SART and the overall workflow, which is summarized below:

Illustration of the workflow of the sparse SART to solve the constrained reconstruction model. For a copy in a certain iteration, the background is initialized by the background of the same copy from the previous iteration, and the object of interest is initialized by the object from the previous copy in the same iteration

Let iteration number k = 1 and copy index i = 1;

For the ith copy, initialize the background of this copy by the background of the same copy from the previous iteration, and initialize the object of interest, hi, by the previous copy hi−1⁠;

Solve the ith reconstruction with the iterative procedure described in Eq. (19), and i=i+1⁠;

Repeat Step 2 and Step 3 until all the copies have been updated;

Go to next iteration (⁠k=k+1⁠) until convergence or the maximal number of iterations is reached.

We first prove the effectiveness of our method on two datasets: the widely-used Shepp–Logan head phantom (Shepp and Logan, 1974) and the ribosome structure downloaded from EMDB-3489 (the center slice is used) (James et al., 2016) (Fig. 3). The sizes of the two images are 320 × 320 and 200 × 200 pixels, respectively. The input datasets of projections were simulated as follows. We got a number of copies of the seed image. Each image was stained by various types of surroundings around the object of interest, including the black blocks and noise. Since for each stained image, we could only obtain a limited range of projections, we generated a series of limited tilt angles and projected the image according to these tilt angles, resulting in a series of Radon transformed images.

Seed images for the two datasets. (A) The Shepp–Logan phantom. (B) The center slice of the ribosome structure (EMDB 3489)

The results demonstrated below are based on the sparse SART solution to CRM. The results generated by the simple averaging method and the averaging in the Fourier space with a mask considering the missing wedge (masked averaging in Fourier) (Schmid and Booth, 2008) are also demonstrated as a comparison. For the simple averaging method, the copies are added and averaged in the real space; for the method of masked averaging in Fourier, the copies are transformed into the Fourier space, weighted with a mask according to the missing wedge, averaged in the Fourier space with the synthesized weight, and then transformed back into the real space. Furthermore, our framework can be easily migrated to other reconstruction techniques (e.g. SIRT version in Supplementary Materials).

To visually evaluate the performance of our method, we set N = 3 for the first experiment, which means we have three subtomograms to be averaged. We stained the Shepp–Logan phantom with blocks and stained the ribosome structure with random noise. We limited the tilt angle range to ±50° with an interval of 3° for the experiments in Figures 4–6. We compared our CRM with the simple averaging method and the masked averaging in Fourier. Our CRM was solved by the sparse SART (Section 3.2.2) with 20 iterations and the relaxation parameter to be 0.2. For the sake of fair comparison, the reconstruction of each subtomogram in the simple averaging method and the masked averaging in Fourier was also done by SART with the same parameters.

Comparison of the reconstruction performance on the Shepp–Logan phantom dataset. The first row is seed images for the three copies (A–C) and the corresponding projections (the Radon transforms) (D). The second row is the reconstruction results of the simple averaging method with SART on the three copies (E–G) and the final averaged reconstruction (H). The third row is the masked Fourier magnitude information of the three copies (I–K) and the final weighted reconstruction (L). The fourth row is the reconstruction results of the proposed CRM with sparse SART on the three copies from the last iteration (M–O) and the final reconstruction (P)

Comparison of the reconstruction performance on the ribosome structure (EMDB 3489) dataset. The first row is seed images for the three copies (A–C) and the corresponding projections (the Radon transforms) (D). The second row is the reconstruction results of the averaging method with SART on the three copies (E–G) and the final averaged reconstruction (H). The third row is the masked Fourier magnitude information of the three copies (I–K) and the final weighted reconstruction (L). The fourth row is the reconstruction results of the proposed CRM with sparse SART on the three copies from the last iteration (M–O) and the final reconstruction (P)

Comparison of the Fourier transformed (FT) maps of the final reconstructions by the simple averaging (left), masked averaging in Fourier (middle) and CRM (right). The top row is from the results from the Shepp–Logan phantom dataset and the bottom row is the results from the ribosome structure dataset. (A–C) The Fourier transformed maps of the results of simple averaging method (Fig. 4H), masked averaging in Fourier (Fig. 4L) and CRM (Fig. 4P), respectively. (D–F) The Fourier transformed maps of the results of simple averaging method (Fig. 5H), masked averaging in Fourier (Fig. 5L) and CRM (Fig. 5P), respectively

Figure 4 shows the reconstruction results on the Shepp–Logan phantom dataset. For the groundtruth, the seed images contain clear copies of the object of interest but with various surroundings (Fig. 4A–C). Since the ranges of tilt angles are limited, the projections (the Radon transforms) can only contain limited information for the object of interest (Fig. 4D). The reconstructed structures for the three copies are shown in Figure 4E–G, from which it is clear that the object of interest was blurred due to the missing wedge effect. Figure 4H shows the refined structure gained by the simple averaging method, in which the averaged reconstruction still contains unnatural artifacts and looks very blurry. Then the masked averaging in Fourier is applied to the copies. The masked Fourier magnitude for the three reconstructed structures is shown in Figure 4I–K. Figure 4L shows the refined structure gained by a weighted averaging and inverse transform from these Fourier spectra. Finally, the CRM method is applied to the copies. The sparse SART solution to our CRM started with the first copy and used the reconstructed object of interest from the first copy to initialize the reconstruction of the second copy, which was then used to initialize the reconstruction of the third copy. It then used the reconstruction of the third copy to initialize the first copy for the next iteration and kept doing so, until convergence or the maximal number of iterations was reached (Fig. 4M–O). Figure 4P shows the final refined structure obtained by our CRM. By comparing the reconstructed structures of the three methods (Fig. 4H, L and P), we can find that both the masked averaging in Fourier and our CRM method produce much clearer details than the simple averaging method. However, the result of masked averaging in Fourier still keeps blurry lattice artifacts inside the structure, whereas the result of CRM presents a much smoother structure and looks much closer to the groundtruth phantom.

Figure 5 shows the reconstruction results on the ribosome structure from EMDB 3489. Rather than the block-stained way of adding artifacts, here we added random noise to the seed images (Fig. 5A–C). Similar conclusions can be drawn as in Figure 4. Both the masked averaging in Fourier and CRM outperforms the simple averaging method (Fig. 5H, L and P). However, the result of CRM is the most similar among the three to the groundtruth ribosome structure.

We further analyzed the results of the averaging method and CRM in the Fourier space. Figure 6 shows the corresponding Fourier transform (FT) maps of the two methods. It can be clearly seen that the simple averaging method has a significant number of clearly missing wedges (Fig. 6A and D), whereas the missing wedge effect has been drastically suppressed by the masked averaging in Fourier (Fig. 6B and E) and the CRM method (Fig. 6D and F). Because the masked averaging in Fourier directly operates on the Fourier weights, it seems that the results of the masked averaging have a more uniform distribution of the signals than the CRM results (Fig. 6E vs. Fig. 6F). However, the ones of the CRM have a better signal density (Fig. 6C vs. Fig. 6B).

The results of the FT analysis are consistent with the conclusions above, which explain the superior performance of the CRM framework over the simple averaging method. In the following, we would make a further comparison between the CRM framework and the masked averaging in Fourier.

We evaluated the proposed method by quantitative measures, including the Pearson’s inner-product correlation coefficient (PCC), the structural similarity (SSIM) and the peak signal-to-noise ratio (PSNR) between the reconstruction and the ground-truth in the seed image.

We measured the reconstructions of the three methods on the two datasets as shown in Figures 4H, L and P and  5H, L and P. Table 1 summarizes the comparison results. The results of the masked averaging in Fourier have very close SSIM scores to the ones of CRM. However, in general, CRM outperforms the other two methods on both datasets in terms of all the evaluation measures.

Performance comparison between the simple averaging method (AVG), the masked averaging in Fourier (M-AVG) and CRM in terms of PSNR, SSIM and PCC on the two datasets

Note: The best performance is shown in bold.

Performance comparison between the simple averaging method (AVG), the masked averaging in Fourier (M-AVG) and CRM in terms of PSNR, SSIM and PCC on the two datasets

Note: The best performance is shown in bold.

We further conducted parameter sensitivity analysis to test the robustness of our method with respect to the number of copies, the range of the tilt angles for each copy, and the noise level.

We generated datasets with different numbers of copies (i.e. the subtomograms), ranging from 3 to 49 with a step size of 2. In this experiment, we fixed the tilt angle range for each copy as −50°∼+50° with an interval of 3°⁠. Figure 7 shows the PSNR, SSIM and PCC curves with the changing number of copies.

Parameter sensitivity analysis of the number of copies on the two datasets. Here, the tilt angle range is set to from −50° to +50° with an interval of 3°⁠. The number of copies ranges from 3 to 49 with a step size of 2. (A–C) The PSNR, SSIM and PCC curves on the Shepp–Logan phantom dataset with respect to the number of copies, respectively. (D–F) The PSNR, SSIM and PCC curves on the ribosome structure dataset with respect to the number of copies, respectively

It can be seen that the reconstructions produced by CRM always have better performance than that of the simple averaging and the masked averaging in Fourier method regardless of the number of copies. When the number of copies increases from 3 to 15, the PSNR of CRM quickly improves (Fig. 7A and D). The reason is that at the beginning, a small number of copies is not sufficient to compensate for the lost information due to the missing wedge issue. When the number of copies is more than 15, there is still clear improvement in PSNR of CRM although the speed of the improvement is slower. In contrast, the PSNR of the simple averaging and the masked averaging in Fourier method cannot benefit from the increase in the number of copies, though the masked averaging in Fourier method generally behaves much better than the simple averaging. The stagnation of the quality improvement for the two averaging methods is consistent with previous findings (Bharat et al., 2015) that the quality of the tomograms themselves is a limiting factor for the resolution of the subtomogram averaging. This demonstrates that CRM can effectively integrate information from different copies to compensate for the missing wedge issue whereas the averaging method cannot do so. Similar conclusions can be drawn from the results in terms of SSIM (Fig. 7B and E). Additionally, on the ribosome structure dataset, CRM can generate reconstructions with very high structural similarities (close to 1) when the number of copies is higher than 20. For PCC, CRM can almost always generate reconstructions that are perfectly correlated (PCC≈1) with the ground-truth regardless of the number of copies. On the ribosome structure dataset, the averaging methods stop the quality improvement very early, whereas CRM can still benefit significantly with the increased number of copies (Fig. 7C and F).

Figure 8 shows selected reconstruction results of the three methods with four different numbers of copies, 3, 11, 21 and 49. Consistent with the conclusions from Figure 7, the reconstructions produced by CRM are much clearer and more similar to the groundtruth compared to the ones produced by the two averaging methods. An interesting observation is that the results produced by the masked averaging in Fourier still carry the bell-ring artifacts (Fig. 8F–H), though the specious information within the missing wedge for each copy has been masked. One possible reason is that these artifacts come from the stained backgrounds during the iterative optimization of the similarity in reconstruction for each copy, in which the mutual information within the entire system is discarded. The averaging methods achieve little visual improvement when the number of copies increases from 21 to 49. In contrast, CRM can continuously improve the quality of the reconstructions when more copies are available.

Selected reconstruction results with different numbers of copies. For each dataset, the first row is produced by the simple averaging method; the second row is produced by the masked averaging in Fourier; and the third row is produced by CRM. (A, E, I, M, Q, U) Reconstruction results with 3 copies. (B, F, J, N, R, V) Reconstruction results with 11 copies. (C, G, K, O, S, W) Reconstruction results with 21 copies. (D, H, L, P, T, X) Reconstruction results with 49 copies

Next, we fixed the number of copies to be 11 and changed the tilt angle range from ±30° to ±70° with a step size of 5° (for each copy, the interval of the tilt angle is fixed to 3°⁠). Figure 9 shows the PSNR, SSIM and PCC curves with the changing range of the tilt angle.

Parameter sensitivity analysis of the range of the tilt angle on the two datasets. Here, the number of copies is fixed as 11. The range of the tilt angle changes from ±30° to ±70° with a step size of 5° (for each copy, the interval of tilt angle is fixed to 3°⁠). (A–C) The PSNR, SSIM and PCC curves on the Shepp–Logan phantom dataset with respect to the tilt angle range, respectively. (D–F) The PSNR, SSIM and PCC curves on the ribosome structure dataset with respect to the tilt angle range, respectively

When the range of the tilt angle is larger, the number of projections per copy is higher and thus the missing wedge effect is less severe. Both the averaging methods and CRM use SART within each copy for reconstruction. As expected, all the three methods have better performance with the increasing range of the tilt angle. Consistent with the conclusions from Figure 7, CRM outperforms the averaging methods in terms of all the measures regardless of the range of the tilt angle. Again, on the Shepp–Logan phantom dataset, CRM can almost always generate perfectly correlated reconstructions within any tilt angle range tested here (Fig. 9C).

Understanding the effect of noise is critical to investigate the practical usefulness of a tomography reconstruction algorithm. It is related to the information theory and the condition number in linear system analysis (Fried, 1972). Here, we evaluated the effect of noise on our CRM by measuring the performance of CRM under different levels of noise and different numbers of copies.

Due to its abundance in ultra-structures, the ribosome structure dataset was chosen in this experiment. For each copy, the projection images were then generated from the stained seed image within the ±50° tilt angle range with an interval of 3°⁠. And different levels of Gaussian random noise were added to the projection images, with SNR 1, 5 and 10. For a fair comparison, both the averaging method and our CRM used SART with 40 iterations and the relaxation parameter of 0.2. Because the masked averaging in Fourier has been widely used in practice and achieves better results than the simple averaging method, here we mainly present the comparison between the CRM and the masked averaging in Fourier.

Figure 10 shows the reconstruction results of the two methods under different levels of noise and with different numbers of copies. In terms of PSNR and PCC, CRM almost always outperforms the averaging method, except for the high noise but with a small number of copies (SNR = 1 and N≤61⁠). This implies that when the noise level is very high, a very small number of copies is not sufficient for CRM to compensate for the loss in the reconstruction resolution caused by the noise. However, if a sufficient number of copies is given, CRM can effectively utilize the information in different copies regardless of the level of noise. Given a certain noise level, the performance of CRM increases quickly with the increasing number of copies, whereas the masked averaging method cannot benefit from such additional information. This phenomenon has also been observed in our previous experiments (Figs 8 and 9), which may need further efforts to explain. In terms of SSIM, both methods have comparable values. These results demonstrate that while both the averaging method and CRM can well recover the main structure of the subtomogram, CRM is robust to the noise and can take great advantage of the multiple copies to reach a much higher quality of reconstruction than the averaging method.

The performance of reconstructions under various levels of noise and various numbers of copies. AVG stands for the averaging method and CRM stands for the constrained reconstruction model. Both models were solved by SART with 40 iterations and the relaxation parameter of 0.2. PSNR stands for the peak signal-to-noise ratio between the reconstruction and the ground-truth in the seed image; SSIM stands for the structural similarity; and PCC stands for the Pearson’s inner-product correlation coefficient.

Except for the reconstruction accuracy, the convergence and runtime is also an important measure for a novel iterative method. Because the solution of the constrained reconstruction framework inherits the features from the algebraic reconstruction technique, its convergence depends on the internal update strategy. Considering the usage of SART in the internal loop of CRM and the exact proof of SART’s convergence (Jiang and Wang, 2003), the sparse SART solution to CRM (CRM-SART) always has a converged station.

The comparison of reprojection error under different iteration number. (A) Results from the dataset of Shepp–Logan phantom and (B) results from the dataset of ribosome structure

We have proved the robustness of our CRM under different conditions. Here, we further extend the CRM to three-dimensional cases and try to analyze its performance. The performance of CRM is compared with that of masked averaging in Fourier. To avoid the effect of the non-structural sharp cutting edge of an EMDB density map in the estimation of resolution, we use the PDB-converted structures instead of EMDB maps in our following experiment, which are produced by EMAN2’s e2pdb2mrc.py (Tang et al., 2007).

Two different structures are used, including the MscS channel structure of Escherichia coli converted from PDB 5Y4O (Yu et al., 2018) and the GTP-binding nuclear protein Ran (protein) converted from PDB 5ZPU. A modified InSilicoTEM (Vulović et al., 2013) is used to simulate the tilt series, in which the volumes of particles (subtomograms) are stained by noise and artificial membranes. To focus on the effects of the reconstruction model, the effect of the contrast transfer function is prefixed. For each structure, volumes with 200 particles under different positions are simulated. The SNR of the projections is controlled around 1:1 and the tilt range is set to −50°+50° with 5° intervals.

Since the geometry position for each subvolume and the projection parameters are known in prior during simulation, it is not difficult to obtain an averaged subtomogram by masked averaging in Fourier for the particles as well as the accurate projection parameters for each particle in our CRM. To analyze the results, both Fourier Shell Correlation (FSC) and ResMap (Kucukelbir et al., 2014) were used to compare the performance between the masked averaging in Fourier and our CRM. For each reconstruction result, a sphere mask was applied to suppress the effects of unnecessary surroundings, and the ground-truth that converted from the PDB was used as the reference. All the results were calculated with the Å/voxel setting to 1.

Figure 12 shows the reconstruction results for both the methods. For a fair comparison, both the averaging method and our CRM used SART with 40 iterations and the relaxation parameter of 0.2. Judging from the 3D structures shown in Figure 12, we can find that the masked averaging in Fourier and our CRM achieved similar structures. However, the results of CRM has a clearer boundary and detailed structures. The ResMap results and FSC curves show that the average method didn’t achieve as high resolution as our CRM for both the MscS channel structure and the GTP-binding nuclear protein Ran. According to the frequency value truncated at 0.5 (FSC0.5), for the MscS channel structure, the result of the masked averaging in Fourier has a resolution of 3.08 Å⁠, while the result reconstructed by our CRM has a resolution of 3.0 Å⁠; for the GTP-binding nuclear protein, the result of the masked averaging in Fourier has a resolution of 3.05 Å⁠, while the result reconstructed by our CRM has a resolution of 2.98 Å⁠. Furthermore, it can be found that the FSC curves generated by CRM generally have higher correlation values compared with the ones of the averaging method, which indicates that the results of CRM are much closer to the groundtruth. In details, the local resolution distributions generated by ResMap show that the results of CRM have the local resolutions mainly distributed at 2.2 Å and 3.2 Å⁠, while the results of the averaging method have a lot of local resolutions distributed at 4.2 Å⁠. Moreover, the ResMap reported a mean resolution of 3.10 Å and 2.96 Å for the masked averaging in Fourier’s results of the MscS channel structure and the GTP-binding nuclear protein, respectively. As a comparison, the reported mean resolutions of our CRM are 2.96 Å and 2.77 Å for the results of the MscS channel structure and the GTP-binding nuclear protein, respectively. Here, both the results of FSC and ResMap indicate that our CRM can gain an improvement in the reconstruction quality compared to the masked averaging in Fourier.

The 3D reconstruction results and their detailed analysis of FSC and ResMap. (A and B) The reconstruction results of the MscS channel structure converted from PDB 5Y4O. (C and D) The reconstruction results of the GTP-binding nuclear protein Ran converted from PDB 5ZPU. Here, (A and C) are results carried out by the masked average in Fourier, whereas (B and D) are results carried out by our constrained reconstruction model

In this article, we proposed a novel CRM for high-resolution subtomogram averaging. Our model can alleviate the missing wedge issue caused by the limited tilt angle range by effectively integrating information from multiple copies of the object of interest. We gave the first theoretical proof to guarantee that the solution space of our model is no larger, and in practice much smaller, than that of the averaging method. We then proposed a sparse Kaczmarz algorithm-based solution for our model and extended it to a novel sparse SART solution, which is used in our final framework. It should be noted that although our solution can compensate for the missing wedge effect, it does not use any prior knowledge. Our solution only focuses on the reinforcement of natural information underlying the large number of projections. Therefore, our proposal will not introduce any false information related to the presumption of prior knowledge. Experimental results demonstrate the effectiveness, robustness and generalization power of our model, which promise its huge potential for high-resolution subtomogram averaging.

The proposed CRM is focused on the final information usage in reconstruction. Here, we did not try to change the optimization of ultrastructure’s orientation and geometry in the alignment procedure. Instead, one could use the classical averaging in the parameter refinement but use the CRM in the final information fusion. Especially, the existence of ring-bell artifacts and the stagnation of the quality improvement in the averaging method may require further comprehensive study. Furthermore, though we have proposed the concept of constrained reconstruction, we do not investigate any prior information and regulation term (e.g. L1-form for sparsity and nonuniform sampling) within the model. Future study to better understand the underlying mechanism of constrained reconstruction is needed, as well as the combination with other optimization techniques. Considering the complex situation in cryo-electron microscopy, further efforts are required to apply the constrained reconstruction to practical applications.

This work was supported by the National Key Research and Development Program of China (2017YFE0103900 and 2017YFA0504702), the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under Awards No. FCC/1/1976-18-01, FCC/1/1976-23-01, FCC/1/1976-25-01, FCC/1/1976-26-01 and FCS/1/4102-02-01, and FCS/1/4102-02-01, and the NSFC projects Grant (U1611263, U1611261, 61932018, and 61672493).

Conflict of Interest: none declared.

Andersen A., Kak A. (

Google Scholar

WorldCat

Andersen A.H. (

Google Scholar

WorldCat

Bartesaghi A.  et al.  (

Google Scholar

WorldCat

Beister M.  et al.  (

Google Scholar

WorldCat

Bharat T.A., Scheres S.H. (

Google Scholar

WorldCat

Bharat T.A.  et al.  (

Google Scholar

WorldCat

Briggs J.A. (

Google Scholar

WorldCat

Carazo J.-M. (

Google Scholar

WorldCat

Chen Y., Förster F. (

Google Scholar

WorldCat

De Rosier D.J., Klug A. (

Google Scholar

WorldCat

Deng Y.  et al.  (

Google Scholar

WorldCat

Förster F.  et al.  (

Google Scholar

WorldCat

Frank J. (

Google Scholar

Google Preview

WorldCat

Fried I. (

Google Scholar

WorldCat

Galaz-Montoya J.G., Ludtke S.J. (

Google Scholar

WorldCat

Gilbert P. (

Google Scholar

WorldCat

Gordon R.  et al.  (

Google Scholar

WorldCat

Han R.  et al.  (

Google Scholar

WorldCat

Han R.  et al.  (

Google Scholar

WorldCat

Harris W.  et al.  (

Google Scholar

WorldCat

James N.R.  et al.  (

Google Scholar

WorldCat

Jiang M., Wang G. (

Google Scholar

WorldCat

Kaczmarz S. (

Google Scholar

WorldCat

Kovacik L.  et al.  (

Google Scholar

WorldCat

Kucukelbir A.  et al.  (

Google Scholar

WorldCat

Kunz M., Frangakis A.S. (

Google Scholar

WorldCat

Kupsch A.  et al.  (

Google Scholar

WorldCat

Leary R.  et al.  (

Google Scholar

WorldCat

Penczek P.A.  et al.  (

Google Scholar

WorldCat

Prakash P.  et al.  (

Google Scholar

WorldCat

Scheres S.H. (

Google Scholar

WorldCat

Schmid M.F., Booth C.R. (

Google Scholar

WorldCat

Shepp L.A., Logan B.F. (

Google Scholar

WorldCat

Shimizu K.  et al.  (

Google Scholar

WorldCat

Tang G.  et al.  (

Google Scholar

WorldCat

Vulović M.  et al.  (

Google Scholar

WorldCat

Wu T.  et al.  (

Google Scholar

WorldCat

Xu M., Alber F. (

Google Scholar

WorldCat

Xu M.  et al.  (

Google Scholar

WorldCat

Yu J.  et al.  (

Google Scholar

WorldCat

Zhao Y.  et al.  (

Google Scholar

WorldCat

Zhu N.  et al.  (

Google Scholar

WorldCat

The authors wish it to be known that, in their opinion, Renmin Han and Lun Li authors should be regarded as Joint First Authors.