Cryogenic electron tomography reconstructions from phaseless data

We perform three-dimensional (3D) reconstructions from simulated and real cryogenic electron tomography (cryo-ET) data. Our reconstructions are based on a nonlinear and phaseless forward model very reminiscent of a commonly used model for phase contrast x-ray tomography.


Introduction
Today cryogenic electron microscopy (cryo-EM) plays an important role in structural biology [MBB + 16,MNM05,NKLB06]. Indeed, it has the potential of providing structural information of biological particles (e.g. macromolecules and supramolecular structures), including particles that have resisted sustained efforts from the conventional techniques of x-ray crystallography and nuclear magnetic resonance (NMR) spectroscopy [VBLJ17,LFB05,Pep17,Cal15]. This unique ability of cryo-EM was recognized by a Nobel Prize for Chemistry in 2017 [Brz17].
In cryo-EM, the specimen is frozen so rapidly that the aqueous buffer forms an amorphous solid that does little or no damage to the structural integrity of the object one seeks to reconstruct (a process known as vitrification). Once the specimen is vitrified, it can be imaged with an transmission electron microscope (TEM) following the data acquisition protocols associated with either electron tomography (ET) [Fra06] or single particle analysis (SPA) [CGPW15]. These techniques yield a series of TEM-images that relate to tomographic projections of the specimen under different incident angles of the probing electron beam. The 3D structure of the imaged object is then reconstructed from the measured data using techniques from computed tomography.
Traditionally, methods for tomographic reconstruction within the context of cryo-ET of biological specimen have been based on an affine forward model, with a ray transform as an essential ingredient [Haw92,VVvVR14]. The affine model can be derived as an approx imation of a more accurate nonlinear and phaseless forward operator for ET Ökt15, that has a mathematical structure closely connected to a standard model in phase contrast x-ray tomography [JL04].
Based on the afore-mentioned phaseless operator for ET, and utilizing the recently introduced generalized-SART principle [Mar18], we present a novel application of regularized Newton-Kaczmarz methods [MBK + 16, BK06] to reconstructions from simulated and experimental cryo-EM data. We also compare its performance to a few other methods, including the well-established simultaneous iterative reconstruction technique (SIRT) and weighted back projection (WBP) [Pen10,EARR15,WLL14].
This article is organized as follows. In section 2, we define our forward operator for electron tomography. Section 3 constitutes an introduction to regularized Newton-Kaczmarz methods and the generalized SART-principle. In section 4, we present the result of our numerical experiments. Finally, in section 5, we give a conclusion and outline possible future extensions of this work.

Model for TEM-imaging
Following the in-depth treatment of [Ökt15] and [FÖ08] (see also [KR08,HK96,Fra06]), we provide a brief outline of a nonlinear model for electron tomography.

Illumination
We start out with the assumption that the electron source has provided us with a monochromatic plane wave, with space-dependent part given as where ω ∈ S 2 defines the TEM optical axis and k denotes the wave number of the incoming electron.

Specimen interaction
As sketched in figure 1, we assume that the specimen resides in a slab Ω ⊂ R 3 bounded by the two parallel planes Γ in and Γ out := Γ in + tΩ (where t > 0 is the slab thickness), and is characterized by its electrostatic potential U : Ω → R 0 . According to the projection assumption, the electron experiences a phase-shift as it travels through the specimen. More precisely, if the scattering operator T sc is the operator that maps U to the electron wave ψ out at Γ out , then it is given by 4 : Here, C in := exp(ikt), P is the is the Ray transform, i.e.
and σ := me k 2 , where m is the electron mass, e is the absolute value of the elementary charge and is the reduced Planck's constant. Inelastic scattering is accounted for in the aforeoutlined framework by replacing the real-valued electrostatic potential with a complex-valued scattering potential U : Ω → C 0 := {a + bi | a, b 0}.
Remark 2.1. In the standard model for electron tomography, which we do not apply, one utilizes in addition to the projection assumption also the weak phase object assumption, which amounts to linearizing the exponential function in (2).

Optics
The part of the microscope between the specimen and the detector is referred to as the optics. It contains a number of lenses, but can nevertheless be modeled accurately by a single thin lens followed by an aperture (see figure 1). The effect of the optics on the electron wave is captured by the optics operator, which maps the electron wave ψ out on Γ out to ψ image on the image plane Γ image . Defining convolution on ω ⊥ as Figure 1. Schematic representation of a TEM-setup for electron tomography with a single thin lens of focal length f and an aperture in the focal plane. The incoming electron wave ψ in scatters against the scattering potential U : Ω → C (light-grey region). The operator T sc (U) in (2) maps ψ in on Γ in to ψ out on Γ out . The operator T op in (5) maps an electron wave on Γ out onto an electron wave on the image plane Γ image . Adapted from [Ökt15]. © 2017 Springer Nature. With permission of Springer.
the optics operator can be expressed as Here φ is a function satisfying |φ(x)| = 1 for all x, M is the magnification of the optical system and In the above and we have used the following version of the Fourier transform on ω ⊥ : Moreover, f is the focal length of the lens, ∆z is the defocus, r aper is the radius of the circular aperture in the focal plane and C s is the third-order spherical aberration. Defining the contrast transfer function (CTF) as we can write the optics operator as where for c ∈ R the dilation operator D c acts on a function f (ω, x) by

Detection
According to the Born rule of quantum mechanics, one has not access to (direct) measurements of the complex-valued wave-function, but rather merely to it's point-wise squared modulus. To account for this, we define the intensity operator by Typically, one linearizes the intensity operator at zero, following the weak object assumption, but in this work we retain the fully nonlinear model. The final ingredient is the detector operator that maps an intensity distribution I at Γ image to the expected detector response. It is given by a convolution with the detector response function PSF det as Here N 0 (Ω) is the image dose, i.e. the average number of incident electrons per unit area at Γ img and C gain is the overall gain that measures the average number of digital counts per electron.

Inelastic scattering and source incoherence
The model presented above may be extended so as to account for spatial and temporal incoherence. This is done by replacing the CTF in equation (10) by where the so-called envelope functions are given by Here ν m is the mean energy spread, C s is the chromatic aberration, α c is the aperture angle of the beam furnished by the condenser and where := e 2mc 2 and U acc is the acceleration voltage of the source.

Total forward operator
The total forward operator is simply the composition of the operators described above: We remark that this operator is very closely related to a common model for phase contrast x-ray tomography as presented in e.g. [JL04]. Indeed if one assumes a perfect detector response-so that T det disappears from the above expression-the models for the two imaging modes agree except for the kernels of the convolution operators involved.

Background correction and constant amplitude contrast ratio
The specimen is made up of two parts, i.e. U = U bg + U , where U bg ≡ C bg ∈ C represents the surrounding aqueous buffer and U is the part we seek to reconstruct (so it could represent e.g. a virion). Thus we consider a forward operator acting only on the part we are interested in: By propagating U bg through all sub-operators we find: where h(ω) is the path-length of the electron through the buffer. We have not explicitly used B(ω) in our reconstructions, instead data was pre-processed by a normalization procedure, such that pre-processed data may be assumed to originate from an operator with B(ω) ≡ 1. Finally, in our experiments we have made use of a common simplifying assumption in ET, namely that the real and imaginary parts of the scattering potential are proportional, i.e. there exists a constant Q ∈ R, usually referred to as amplitude contrast ratio, such that With this assumption in place we define our final forward operator as

Model for ET
So far we have treated the model for a single cryo-EM image. In ET recorded data however consist of a tilt-series of micrographs acquired at varying specimen tilt. A common tomographic acquisition geometry is parallel-beam single-axis geometry, i.e. specimen movement in between micrograph generation is restricted to rotation around a single tilt-axis. Since the electron path-length though the slab-like sample increases with increasing specimen tilt, and also because of physical limitations of the interior of the microscope, one is limited to tiltangles in some subset of [−θ max , θ max ], where a common value for θ max is 60 • . We remark that the inverse problem of reconstructing a 3D biological structure from data generated by a cryo-ET experiment is a very challenging one, the most important reason for this being the so-called dose problem. The latter refers to the fact that the electron dose-i.e. the number of incident electrons per specimen surface area-must be kept low in order to avoid severe degradation of the specimen structure. This in turn leads to highly noisy micrographs, see figure 5(a). Another important difficulty is that in ET one has a limited data-problem, due to the restrictions on the tilt-angle described in the preceding paragraph.

Regularized Newton-Kaczmarz
Regularized Newton-Kaczmarz methods were first proposed in [BK06] as a reconstruction technique for nonlinear inverse problems with a block-structure. The approach has been applied to x-ray phase contrast tomography in [MBK + 16]. Given a problem of the form (G 1 ( f ), . . . , G N ( f )) = (g 1 , . . . , g N ) with block-forward operators G j : X → Y j between Hilbert-spaces X, Y 1 , . . . , Y N and data-chunks g j ∈ Y j , the idea is to sequentially perform Levenberg-Marquardt iterations on the different subproblems G j ( f ) = g j , i.e.
for k = 0, 1, . . . , k stop and some processing-order { j k } ⊂ {1, . . . , N}. Here, G j k [ f k ] denotes the Fréchet-derivative of the operator G j k evaluated at the current iterate f k and β > 0 is a regularization parameter. The benefit of iterations of the form (25) is that their computation is typically much cheaper than iterations on the full inverse problem.
Iterations of the form (25) generalize classical Kaczmarz-type tomographic reconstruction methods such as ART and SART [GBH70,AK84], which correspond to the case where G j are taken as blocks of a tomographic projection operator. It is typically observed that-given a suitable processing-order {j k }-Kaczmarz-methods yield a reasonable reconstruction already within a single cycle over the whole data set, i.e. after a single iteration with respect to each data-chunk g j [NW01,EHN14]. Potentially, the approach thereby enables massive performance gains compared to bulk-iterative methods that fit the full data set in each iteration.

Generalized SART-principle
Recently, [Mar18] derived a general methodology, termed Generalized SART-principle, which enables an efficient implementation of regularized Kaczmarz-iterates within the context of tomographic inverse problems. Among other settings, the approach was demonstrated for phase contrast x-ray tomography with Sobolev-regularization.
The generalized SART approach assumes that the forward operator-blocks G j are given by the composition of a ray-transform P j (U)(x) := P(U)(x, ω j ), x ∈ ω ⊥ j for a single incident direction ω j ∈ S 2 and an additional image-formation operator that acts on the projection data. From (2) and (19), we see that the electron tomography problem exhibits such a structure: the full forward operator assumes the form We consider Newton-Kaczmarz iterations (25) with operator-blocks G j := T proj • P j and a Sobolev-W 1,2 -penalty term, corresponding to the norm f 2 Using the chain-rule, the resulting update for the scattering potential is obtained as where I j denotes the measured intensity distribution under the incident direction ω j . According to the generalized SART-principle, the minimizer in (28) (for a slightly smaller search-set, see [Mar18]) can be computed via the formula: Here, u j := P j (1 Ω ) is the projection of a function that is constant one within the reconstruction domain Ω and zero outside and P * j is the back-projection operator corresponding to P j , i.e. its L 2 -adjoint.
Most importantly, the update-formula (29) requires only a single evaluation of P j k and P * j k in (29a) and (29c), whereas the optimization problem in (29b) is completely formulated in terms two-dimensional images r k , p, p k . Accordingly, it does not involve costly computations on 3D volumetric data.  For the numerical experiments of the subsequent section, a discretized version of the generalized SART-formula (29) is iterated, where the gradients are approximated by finite differences. We solve the discrete quadratic optimization problem arising from (29b) via a conjugate gradient method. Non-negativity of the reconstructed scattering potential is imposed by setting negative values to zero after each iteration. For the processing-order of the Kaczmarz-blocks, we choose the multi-level-scheme from [GG94].

Numerical experiments
Our reconstructions were carried out using Operator Discretization Library (ODL), which is a python package for numerical functional analysis [AKÖ17].

Simulated data
We used a phantom (depicted in figure 2(a)) consisting of a number balls of varying radii and contrast. It was represented by a 250 × 210 × 40 voxelized volume with voxel-size 0.5 nm, and coincides up to an additive constant with the one considered in [Ökt15 5 ] (we have subtracted the background from the latter). Corresponding micrograph data was generated using 61 tilt-angles, with uniform angular coverage in the range ±60 • . Electron microscope-and specimen parameters were as in table 1. Ideal detector response was assumed, i.e. T det was set to the identity operator on the detector plane and our simulation assumes a detector quantum efficiency (DQE) [RÖM + 11] equal to unity. In all methods that we have implemented, we have enforced positivity on the reconstructed function.
4.1.1. Noise-free data. Firstly, we consider reconstructions from noise-free data ( figure 2(b)). Generalized-SART with β = 0.001 and γ = 0 (see section 3) yielded a reasonably good reconstruction already after a single cycle, see figure 3(a). We also performed reconstructions with Landweber iterations, a technique which generalizes the well-known SIRT algorithm [NW01]. Landweber required about 1000 iterations to produce a reasonably good reconstruction and about 5000 iterations to achieve a completely artifact-free result. For comparison, we also performed reconstructions based on a linearization of our forward operator. Moreover, in order to enhance the effect of nonlinearities, we ran additional reconstructions on data generated from a phantom with twice the contrast. Certain ring-shaped artifacts persisted in reconstructions based on the linearized model, even after 10 000 iterations (see figures 3(e) and (d)). On the contrary, no such artifacts occur in the generalized-SART reconstructions, being based on the fully nonlinear model. As quantitative measures of the reconstruction qualities, we computed structural similarity index (SSIM) [WBSS04] and peak signal to noise ratio (PSNR) using functionality provided by Python package sci-kit image. Results are presented in table 2.

Noisy data.
In order to add (Poisson) shot-noise to the data, a realistic value of C bg was needed. It was determined by fitting our data to data generated by TEM-simulation software described in [RÖM + 11]. As in the case of noise-free data, we compared generalized-SART (this time with β = 400 and γ = 0.95) to SIRT (see figures 4(a) and (b)). We also computed reconstructions using gradient descent on an objective function given as the sum of the L 2 -data discrepancy and the Huber-norm (with γ Huber = 0.01), and TV-regularization using a nonlinear Chambolle-Pock algorithm [CP11,Val14] (see figures 4(c)-(f)). The same figures of merit as in the noise-free case were calculated also in the noisy setting; results are presented in table 3. While generalized-SART turns out to be competitive in terms of PSNR, it did not perform very well in terms of SSIM compared to the Huber/TV regularized methods. That the latter methods managed to achieve higher SSIM-values might be due to the fact that Huber/ TV regularization produces reconstructions with sparse gradients, which clearly suits the balls phantom. Note however that both the gradient descent-and Chambolle-Pock based methods needed a very large number of iterations, whereas generalized-SART produced a good result after only three cycles.

Experimental data
The single-axis tilt-series was acquired by FEI using a conventional bright-field TEM (FEI Titan Krios with a Falcon direct electron detector), see figure 5. It consists of 81 micrographs and contains Cowpea Mosaic Virus (CPMV), bacteriophage T4, Keyhole Limpet Hemocyanin and Tobacco Mosaic Virus (TMV). We computed reconstructions based on generalized-SART from three sub-regions, and compared these to previously obtained reconstructions. The methods used in the latter include an approximate inverse method (AIM), SIRT, WBP and TV-regularization (TVreg). For details, see [Ökt15]. Electron microscope-and specimen parameters used in the reconstructions are display in table 4, while reconstructions are presented in figures 6-8. generalized-SART reconstructions were obtained using three cycles over the tilt-angles, and regularization parameters were chosen as β = 3000 and γ = 0.95. Although it is difficult to assess the reconstruction-quality for the experimental data set, the visual impression indicates that generalized-SART performs at least better than AIM, SIRT and WBP in terms of noise-suppression and contrast and is widely comparable to TVreg. In particular, while the WBP-results appear to be somewhat smeared-out horizontally due to the limited angular range (see especially 7(d)), generalized-SART (similarly as the other algebraic methods) performs quite well at suppressing such limited-data-artifacts.

Conclusions
We have implemented Sobolev-regularized Newton-Kaczmarz iterations as a reconstruction method for electron tomography, and tested it on simulated and experimental cryo-EM data. Importantly, the algorithm is derived upon the basis of a full nonlinear and phaseless TEMmodel. Therefore, it is able to outperform methods that use a linearization of the contrast, as we demonstrated on simulated data. Our comparisons to standard methods indicate that the Newton-Kaczmarz-approach may compete with state-of-the-art iterative reconstruction algorithms like SIRT and total-variation-regularization in terms of reconstruction quality. At the same time, the proposed method allows for short computation times, owing to fast semiconvergence of Kaczmarz-methods and the efficient implementation via generalized-SART schemes.
Several paths lay open for future research extending what his been presented here, including quantitative comparisons of the generalized-SART approach to other reconstructions methods based on metrics such as the Fourier Shell Correlation (FSC) [DKK12]. We would also like to assess the potential benefits in terms of reconstruction quality from further refinements of the forward model.
Finally, the missing phase information in the nonlinear TEM-model considered in this work raises the issue, whether the scattering potential can always be uniquely reconstructed. We aim to address this topic in the future by adapting uniqueness theory for x-ray phasecontrast imaging [Mar15].