A new solution to the curved Ewald sphere problem for 3D image reconstruction in electron microscopy
AA new solution to the curved Ewald sphere problem for 3D imagereconstruction in electron microscopy
J. P. J. Chen, K. E. Schmidt, J. C. H. Spence, R. A. Kirian
Department of PhysicsArizona State UniversityTempe, AZ 85287 USA
Abstract
We develop an algorithm capable of imaging a three-dimensional object given a collection of two-dimensional images of that object that are significantly influenced by the curvature of the Ewaldsphere. These two-dimensional images cannot be approximated as projections of the object. Suchan algorithm is useful in cryo-electron microscopy where larger samples, higher resolution, or lowerenergy electron beams are desired, all of which contribute to the significance of Ewald curvature.
1. Introduction
In this paper we propose an iterative projection algorithm to reconstruct a 3-dimensional (3D)object from 2-dimensional (2D) real-space images in the case where the depth-of-field is less thanthe thickness of the sample, so that these 2D images are not projections of the sample scatteringdensity. This is also referred to as the “Ewald sphere curvature problem” in the context of cryo-electron microscopy (cryo-EM) Wolf et al. (2006); Leong et al. (2010); Russo & Henderson (2018);Zhu et al. (2018) and it becomes important for resolutions and incident wavelengths where thethickness of the sample exceeds the depth-of-focus of the lens system Spence (2013). In that case,bright-field transmission electron microscope (TEM) images no longer represent valid 2D projectionsof the object, so that conventional cryo-EM software (based, for example, on the method of filteredback-projection) cannot be used to merge the 2D images into a 3D reconstruction without lossesin resolution. Unlike some of the previous approaches discussed below, our method does not relyon application of a corrected transfer function, but builds in the relevant scattering theory fromthe outset, and inverts this from all the data using the technique of iterative projection algorithms(IPAs). A solution to the Ewald curvature problem would allow 3D reconstructions from lower beamenergies, higher resolutions, and with thicker samples, all factors of which contribute to the severityof the Ewald curvature.For a scattering vector q = k − k , the Ewald sphere is defined as the spherical surface tracedout by k under elastic scattering conditions, so that negligible energy is transferred to the sample.Here k is the incoming wavevector of the electron beam, while k is the scattered wavevector. TheEwald curvature problem arises when images are formed from a range of q -vectors over which theEwald sphere cannot be well-approximated as a plane. In that case, the Fourier-projection theorem,relating planes in reciprocal space to projections of the object in real space, as assumed in many Preprint submitted to Ultramicroscopy February 9, 2021 a r X i v : . [ phy s i c s . i n s - d e t ] F e b urrent reconstruction and data merging algorithms, cannot be used. For a crystal (where q extendsto a reciprocal lattice point, defining a Bragg condition), a curved sphere means that beams scatteredat equal angles to the direct beam on opposite sides of the diffraction pattern are no longer Friedelpairs represented by conjugate complex structure factors. The two different complex structure factors(four real numbers) are combined by the lens to form one set of interference fringes (defined by justtwo quantities) contributing to the image. In previous work by Wolf et al. (2006); Russo & Henderson(2018), a correction algorithm was described that relies on the fact that (for a crystal), Bragg beams(whose width is limited to that of a small particle) separate laterally as they progress downstream.This is seen in far-out-of-focus shadow images of small crystals, which separate into one shadow imagefor each Bragg beam. At resolution d , the separation due to defocus, D f , alone is 2 D f λ/d (otheraberrations may also contribute) Zuo & Spence (2017). This allows the interferences contributed byone side of the diffraction pattern to be separated from those due to the other. Russo & Henderson(2018) then generalize their method to non-periodic samples with continuous scattering, and applytheir corrections using an elegant segmented transfer function divided into wedges. By comparison,our method does not require the use of images recorded at large defocus or additional images whichmay introduce radiation damage. However our method proposed here is much more computationallyintensive, and may be well suited to images formed with the use of a Zernike phase plate.The conventional formulation of bright-field, weak-phase object, high-resolution transmission elec-tron microscopy imaging assumes either a sample transmission function for a phase object (basedon an eikonal approximation) using a projected potential Spence (2013), or uses the first Born ap-proximation for electron scattering. The phase shift introduced by the sample is assumed to be lessthan ninety degrees. However, a new formulation is needed for the case of a curved sphere, sincethe image is then formed from scattering on the Ewald sphere rather than onto a plane in reciprocalspace, and a projected potential cannot be used. As detailed in Section 2, we use the standard Bornapproximation for the wave function, but we do not make a far-field approximation. Therefore ourmethod gives the correct relationship between the phase of the direct and scattered radiation Lentzen(2014). As for the case of optical microscopy where the Ewald sphere curvature causes the depth-of-field to be smaller than the sample thickness (as used for the optical sectioning technique), it is thennecessary to compute wavefields on successive planes normal to the beam within the sample, andto propagate these onward to an area detector, or to the “exit face” plane across the downstreamface of the sample. Focus adjustment then makes it possible to “look inside” a semi-transparentobject albeit with limited depth-of-focus and artifacts caused by the other planes of the sample. Formodern aberration-corrected transmission electron microscopes, this depth-of-field can be as smallas a few nanometers. The validity domain of the approximations used in this paper therefore limitour method to weakly scattering objects for which multiple scattering is negligible but Ewald spherecurvature is significant. As discussed elsewhere in the context of the inversion of multiple scatteringDonatelli & Spence (2020), multiple scattering may also occur with a flat Ewald sphere for samplesconsisting of heavy elements.Recently, Gureyev et al. (2019a,b) gave a comprehensive treatment of the Ewald curvature problemand proposed a solution method they call “pattern matching tomography.” In a separate work,Ren et al. (2020) proposed an iterative strategy based on finding a 3D potential that best fits themeasured intensities with regularization, for inorganic non-periodic samples. They also incorporatethe multislice algorithm, allowing their technique to account for multiple scattering. The algorithmwe propose here is similar in nature to the method developed by Ren et al. (2020), with the difference2eing that we use “projection operators” (operations that make the minimum change to an input) toseek out a 3D potential that fits the measured intensities, as opposed to cost function minimization bycalculating gradients. Iterative procedures similar in spirit of our work have been proposed previouslyfor example by ? where they reconstruct 2D objects (exit wave functions) from a through-focus seriesof 2D images.This paper is structured as follows: Section 2 introduces the mathematical model that we will useto describe the imaging process and our criterion for the Ewald curvature to become significant.Section 3 outlines our proposed algorithm. Section 4 shows the results of simulations with imple-mentation details, before drawing conclusions in Section 5.
2. The model
The objective of our work is to reconstruct a 3D image of an object given a collection of 2D imagesformed with the object in different orientations and located at different distances from the focal plane.Assuming no multiple scattering (the first-order Born approximation), ideal plane-wave illumination,and an ideal imaging system, the measurement from the experiment can be calculated as follows.Starting with the Schr¨odinger equation as detailed in Appendix A, it can be shown that the scatteredwavefront ψ ( r ) from an incoming plane wave incident upon an object described by a potential energydistribution V ( r ) under the Born approximation is given by ψ ( x, y, z ) = e ik z (cid:20) − i πEλ (cid:90) d r (cid:48) (cid:90) dk x π (cid:90) dk y π e ik x ( x − x (cid:48) ) e ik y ( y − y (cid:48) ) e i (cid:16) √ k − k x − k y − k (cid:17) ( z − z (cid:48) ) V ( r (cid:48) ) (cid:21) , (1)where r = ( x, y, z ), z > z (cid:48) , k x and k y are the spatial frequency components along the transversedirections x and y , k = 2 π/λ , E and λ are the energy and wavelength of the incident plane wave,respectively. For z < z (cid:48) a virtual image is produced. Introducing the scaled object potential f ( r ) = πEλ V ( r ) (2)and the wavevector transfer q = ( q x , q y , q z )= k − k = ( k x , k y , k z − k ) , (3)where k = ( k x , k y , k z ) is the outgoing wavevector, k = (0 , , k ) is the incoming wavevector, k z = (cid:113) k − k x − k y , and q z = k z − k . The scattering intensity, | ψ ( x, y, z ) | , can be put into theform (cid:12)(cid:12) ψ ( x, y, z ) (cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) − i (cid:90) dz (cid:48) (cid:90) dq x π (cid:90) dq y π e i ( q x x + q y y ) e iq z ( z − z (cid:48) ) (cid:90) dx (cid:48) (cid:90) dy (cid:48) e − i ( q x x (cid:48) + q y y (cid:48) ) f ( r (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12) . (4)Equation (4) can be written as (cid:12)(cid:12) ψ ( x, y, z ) (cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) − i (cid:90) dz (cid:48) F − ⊥ e iq z ( z − z (cid:48) ) F ⊥ f ( r (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12) , (5)3here F ⊥ and F − ⊥ are the 2D Fourier transform and the inverse 2D Fourier transform, respectively,in the transverse directions x and y to the propagation axis. Expanding equation (5) gives (cid:12)(cid:12) ψ ( x, y, z ) (cid:12)(cid:12) = 1 + 2 Im (cid:18)(cid:90) dz (cid:48) F − ⊥ e iq z ( z − z (cid:48) ) F ⊥ f ( r (cid:48) ) (cid:19) + (cid:12)(cid:12)(cid:12)(cid:12)(cid:90) dz (cid:48) F − ⊥ e iq z ( z − z (cid:48) ) F ⊥ f ( r (cid:48) ) (cid:12)(cid:12)(cid:12)(cid:12) . (6)Since the Born approximation only keeps terms linear in the potential, the quadratic term in thepotential in equation (6) will be correct only at points where the cross term between the directand scattered waves is negligible. Because we are in the bright-field regime, we must either drop thesecond order term or include the second Born approximation in this cross term. In typical bright-fieldelectron microscopy, the second order term in f ( r ) is much smaller than unity and thus we chooseto ignore it here. We also drop the constant bright-field term in equation (6) under the assumptionthat it can be removed at an early stage in the data processing. Thus the 2D images obtained in abright-field electron microscope in our model are proportional to I z ( x, y ) = Im (cid:18)(cid:90) dz (cid:48) F − ⊥ e iq z ( z − z (cid:48) ) F ⊥ f ( r (cid:48) ) (cid:19) , (7)where I z ( x, y ) is understood as the continuous 3D scattering intensity evaluated at a specific valueof z . That value of z is referred to here as the “focal distance” that defines the focal plane of theimage.Equation (7) is the governing equation that we use in this paper to describe the electron microscopyimaging process. We do not incorporate the effect of the lens transfer function in this current work,however we present a derivation of it now for context and completeness. Define the scattering angle θ as the angle between k and k (see Figure 2). For θ small, q z can be approximated as q z = k (cos θ − ≈ − πλ θ k z − k ≈ − λ π (cid:0) q x + q y (cid:1) . (9)With the inclusion of spherical aberration C s , the scattering intensity becomes I z ( x, y ) = Im (cid:18)(cid:90) dz (cid:48) F − ⊥ exp (cid:16) − iχ ( q x , q y , ∆ f ) (cid:17) F ⊥ f ( r (cid:48) ) (cid:19) , (10)where ∆ f = z − z (cid:48) is the deviation from Gaussian focus and χ ( q x , q y , ∆ f ) is the familiar lens transferfunction χ ( q x , q y , ∆ f ) = λ π (cid:0) q x + q y (cid:1) ∆ f − C s (cid:18) λ π (cid:19) (cid:0) q x + q y (cid:1) = 2 πλ (cid:18) θ ∆ f − C s θ (cid:19) . (11)As mentioned above, our proof-of-principle work presented here in this paper do not consider thelens transfer function and we work with equation (7) instead of (10).4n order to simulate equation (7) on a computer, we discretize f ( r ) into a total of M
2D sectionsalong the direction of propagation, each separated by a distance ∆ z . Then we may write I z ( x, y ) (cid:39) Im (cid:32) M − (cid:88) m =0 F − ⊥ e iq z ( z − m ∆ z ) F ⊥ f ( x, y, m ∆ z ) (cid:33) (12)up to a proportionality constant ∆ z . Since in an electron microscopy experiment, many imageswith the sample in different orientations and different distances from the focal plane are typicallymeasured, we write the n th measured data image, I n ( x, y ), referred to here as a “view,” as I n ( x, y ) = Im (cid:32) M − (cid:88) m =0 F − ⊥ e iq z ( z n − m ∆ z ) F ⊥ f n ( x, y, m ∆ z ) (cid:33) , (13)where z n is the focal distance for the n th view and f n ( r ) = R n f ( r ) (14)is the scaled object potential rotated by a 3D rotation operator R n that rotates a 3D function to the n th orientation. If we define the free-space propagator as P ( ζ ) ≡ F − ⊥ e iq z ζ F ⊥ , (15)where ζ is the propagation distance, we may simplify equation (13) to I n ( x, y ) = Im (cid:32) M − (cid:88) m =0 P ( z n − m ∆ z ) f n ( x, y, m ∆ z ) (cid:33) . (16)The above expressions show that each view can be understood as the imaginary part of a coherentsuperposition of the scattering amplitudes from a series of thin 2D sections that are all propagatedto a common plane. This process is depicted schematically in Figure 1 for two example views. Theobjective of our work can then be stated as follows: reconstruct the 3D scaled object potential f ( r )given a set of N
2D views, { I n ( x, y ) } . The effect of Ewald curvature is contained in the term q z m ∆ z in the exponential of equation (13). Interms of more useful quantities, this effect of the Ewald curvature depends on three factors: (1) thewavelength of the incoming plane wave, λ , (2) the desired resolution, d , and (3) the size of the sampleitself, L . For single-particle cryo-EM we assume a roughly spherical particle, where L becomes thepropagation distance along the beam path. For a slab-shaped sample, L is the particle thickness.The relationship between the three factors, λ , d and L is derived in this section and has been shownpreviously by ? Spence (2013); ? .Consider an object of uniform density with thickness L in the z direction. Such an object has theFourier transform in the q z direction proportional to F ( q z ) = (cid:90) L/ − L/ dz e iq z z = L sin ( q z L/ q z L/ . (17)5 igure 1: Two schematic examples of how views are calculated. Wavefronts emanating from sections of the object(dashed line) separated by ∆ z and indexed by m are propagated to a common plane (thick solid line) a distance z n along the z axis from a chosen origin and coherently summed. The imaginary part of the resultant 2D image, I n ( x, y ),is what is referred to as a “view.” Note the two examples have the object in different orientations ( R n and R n +1 ) anddifferent focal distances ( z n and z n +1 ). A microscope lens collects elastic scattering amplitudes that are proportional to this Fourier trans-form. For a given view, the amplitudes are constrained to lie on the Ewald sphere, and are limitedby the maximum scattering angle that the lens can collect. The values sampled by the Ewald sphereare indistinguishable from the values sampled on a planar slice through the origin in Fourier spacewithin regions where F ( q z ) is relatively constant. The first zero of F ( q z ), in this case a sinc function,occurs at q z L/ π , so the effect of Ewald curvature is definitely significant when q z ≥ πL . (18)Writing q in terms of spherical coordinates gives, q = k (sin θ cos φ ˆ x + sin θ sin φ ˆ y + cos θ ˆ z − ˆ z )= 2 k sin( θ/ (cid:104) cos( θ/
2) cos φ ˆ x + cos( θ/
2) sin φ ˆ y − sin( θ/ z (cid:105) (19)so that equation (18) becomes L ≥ λ ( θ/ . (20)With q = | q | = 2 k sin( θ/
2) = 4 π sin( θ/ λ , (21)6 igure 2: Reciprocal space construction for determining the condition for Ewald curvature to become significant. Thecircle shows the Ewald sphere. The dashed line indicates the position of the first zero of the sinc function describedin the main text. O is the origin in Fourier space. the resolution d defined so that qd = 2 π is d = λ θ/ . (22)Relating the resolution from equation (22) to (20) then yields L ≥ d λ (23)as the condition for the effect of the Ewald sphere to be significant. Equation (23) is equivalent tothe requirement that the Fresnel number F = d / ( Lλ ) exceeds 0.5 for the effect of the Ewald sphereto be ignored, in agreement with ? . If L is interpreted as a propagation distance along the beamdirection, and d the resolution limit imposed by a lens, then 2 d /λ is the depth-of-field, i.e., therange of planes considered to be in focus at resolution d , referred to the object space. Equation (23)therefore expresses the condition that the depth-of-field is less than the sample thickness, so that“optical sectioning” (with resolution 2 d /λ reckoned along z ) is possible for 3D imaging by recordingimages for many values of z in equation (7). For a defocus value of 2 d /λ , d also gives the width ofthe first Fresnel edge fringe, or the width of zones in a Fresnel zone plate with focal length 2 d /λ .Equation (23) written in another way becomes ε = 2 d λL ≤ , (24) The phase factor in the Fresnel propagator gives L = d / (4 πλ ) which is equivalent to requiring q z L ≥ / q z L ≥ π , as is used in the derivation here for Ewald curvature to become significant. ε , which is the ratio of the depth-of-field to theobject thickness, and is equal to twice the Fresnel number. Equation (24) says: if ε is less than orequal to one, the effect of Ewald curvature is significant. Figure 2 shows how the parameters λ , d and L can be related in Fourier space at the condition where the effect of the Ewald sphere becomessignificant under our criterion. Figure 3 shows geometrically the effect on the region where the Ewaldsphere is considered “flat” when one of those three parameters is altered while the other two remainfixed.As an example, for 100keV electrons ( λ ≈ .
037 ˚A) and an object of thickness L = 60nm, theresolution at which Ewald curvature becomes relevant is d = 3 .
33 ˚A. Figure 4 shows the resolutionat which the effect of Ewald curvature becomes significant for different object thicknesses at fourdifferent beam energies. Figure 5 shows example views of a virus particle calculated at three differentincident beam energies and two different focal distances.
Figure 3: The effect of the curvature of the Ewald sphere can be ignored by choosing (a) low enough resolution d , (b)small enough wavelength λ , or (c) small enough object size L . The thick solid curve shows the region in Fourier spacewhere the Ewald curvature is considered to be not significant. The dashed line indicates the position of the first zeroof the sinc function described in the main text. O is the origin in Fourier space. igure 4: Object size L versus resolution d at which Ewald curvature is considered to be significant according toequation (24) for electrons at four different kinetic energies.Figure 5: Example views of a virus particle at three different electron kinetic energies and two different focal distances.The virus in this simulation is 29nm across. The greyscales are set to the same maximum and minimum values in allsubfigures. . Algorithm Many problems can be posed in terms of the satisfaction of multiple constraints. The solution to theoriginal problem requires that all of the constraints be satisfied. The geometric interpretation of thisis that the solution is located at the intersection of all surfaces defined by the constraints in a highdimensional space. One way to arrive at the intersection, and hence solve the original problem, isvia algorithms called iterative projection algorithms (IPAs). Operations to “project” onto individualconstraint surfaces, referred to as “projection operators”, are constructed where they make theminimum possible change to an iterate, denoted here by f , such that a specific constraint is satisfied.The boldface indicates that f is the vectorized representation of a discrete function, but can also bethe vectorized representation of a collection of functions in general, as is the case for our proposedalgorithm. The projection operators can be combined to form deterministic rules that update f insuch a way as to progressively satisfy all constraints. See Gerchberg & Saxton (1972); Fienup (1982);Bauschke et al. (2002); Elser (2003); Luke (2005); Elser et al. (2007); Marchesini (2007); Millane &Lo (2013) and references therein for more in-depth discussions on IPAs and constraint-satisfactionproblems.Typically the update rule for the IPA is formulated such that the next iterate is generated from acombination of two projection operators, denoted here by P S and P M , acting on the current iterateat the j th iteration, f ( j ) . The simplest IPA is the Error Reduction (ER) algorithm Gerchberg &Saxton (1972); Fienup (1982) where the j th iterate is updated according to the rule: f ( j +1) = P S P M f ( j ) . (25)The ER algorithm moves the iterate steadily towards a fixed-point but is unable to escape and exploreother regions of the multi-dimensional space should that fixed-point turn out not to be a solution.Such outcomes are common when the constraints involved are non-convex in the defined space. Amore effective update rule is the Relaxed-Averaged-Alternating Reflections (RAAR) algorithm Luke(2005) in which the j th iterate is evolved according to the rule: f ( j +1) = β f ( j ) − βP S f ( j ) + (1 − β ) P M f ( j ) + 2 βP S P M f ( j ) , (26)where β is a real-valued parameter of the algorithm.The above formulation deals naturally with two constraints. When there are more than two con-straints, an effective strategy was proposed by Gravel & Elser (2008) in an approach they have coined“divide-and-concur.” In this approach the iterate f contains as many copies of the discrete function f as there are constraints: this is the so-called “divide” step of the algorithm. In one of the projectionoperators each copy of f is made to satisfy a single individual constraint (hence the reason for havingone f for each constraint), the other projection operator replaces the set of all f with the average f , calculated by N (cid:80) Nn =1 f n for N equally weighted constraints: this is the so-called “concurrence”step, which enforces the requirement that all copies of f must describe the same discrete function.We shall identify the iterate in the case of N constraints with curly brackets f = { f n } where thisnotation is understood to mean the iterate f contains the set of all N f n , and more precisely, f is theconcatenated vectorization of all N discrete functions f n .10 .1. Our problem In our problem we have N constraints where each of the constraints corresponds to a single view, I n ( x, y ), recorded by the microscope. Therefore, following the divide-and-concur recipe, the iterate f consists of a set of N volumes { f n ( x, y, z ) } . The projections P S and P M for our problem then requireoperators that map a single volume to a set of N volumes, and the corresponding inverse operatorsthat take the set of N volumes and map them back to a single common volume. These operators aredefined as follows. The forward and inverse propagation operators, as introduced in Section 2, are P ( z ) = F − ⊥ e iq z z F ⊥ (27) P − ( z ) = F − ⊥ e − iq z z F ⊥ = P ( − z ) . (28)Define L mn = P ( z n − m ∆ z ) (29) L − mn = P − ( z n − m ∆ z ) , (30)so that we may write ρ mn ( x, y ) = L mn f n ( x, y, m ∆ z ) (31) f n ( x, y, m ∆ z ) = L − mn ρ mn ( x, y ) , (32)where ρ mn ( x, y ) is the m th 2D section from f ( x, y, z ) rotated by R n and propagated to focal distance z n . We further define the operators L and L − such that L f ( x, y, z ) = (cid:110) R n f ( x, y, z ) (cid:111) = (cid:110) f n ( x, y, z ) (cid:111) = f (33) L − f = L − (cid:110) f n ( x, y, z ) (cid:111) = 1 N N (cid:88) n =1 R − n f n ( x, y, z ) = f ( x, y, z ) . (34)The operator L takes f ( x, y, z ) and makes N copies of it with each copy rotated by R n to formthe set { f n ( x, y, z ) } . The operator L − takes the set { f n ( x, y, z ) } , inverse rotates each f n ( x, y, z ) inthat set and outputs the average of those N inverse rotated volumes. The action of L and L − issummarized graphically in Figure 6. Using the defined operators, the n th view from equation (16)can be written as I n ( x, y ) = Im (cid:32) M − (cid:88) m =0 ρ mn ( x, y ) (cid:33) = Im (cid:32) M − (cid:88) m =0 L mn f n ( x, y, m ∆ z ) (cid:33) . (35)Note that the information contained in the set of all M N ρ mn ( x, y ), i.e., { ρ mn ( x, y ) } , is in principleexactly the same as that contained in the set of all N f n ( x, y, z ), i.e., { f n ( x, y, z ) } , due to theunitarity of the forward and inverse propagation operators P and the rotation operator R . Inpractice however there will be some loss of information going forward and coming back due to theinterpolations required for rotations onto a Cartesian computational grid.11 igure 6: The action of L and its inverse. The 3D arrays each represent an individual f ( x, y, z ). The square bracketsdenote the set { f n ( x, y, z ) } , which is the iterate f in our iterative projection algorithm. Here we define the two projection operators, P M and P S , for our problem as follows. First, denote ρ mn ( x, y ) at the j th iteration of the algorithm by ρ ( j ) mn ( x, y ). The iterate at the j th iteration can bewritten similarly as f ( j ) and the views given by that iterate as I ( j ) n ( x, y ). The views I ( j ) n ( x, y ) will ingeneral not be equal to the measured data views I data n ( x, y ) and the projection operator P M seeks tomake them equal with the minimum amount of change to ρ ( j ) mn ( x, y ). The appropriate operation canbe shown to be (see Appendix B for more details) P Mn ρ ( j ) mn ( x, y ) = ρ ( j ) mn ( x, y ) + i M (cid:16) I data n ( x, y ) − I ( j ) n ( x, y ) (cid:17) , (36)where P Mn operates on the ρ ( j ) mn ( x, y ) derived from f n ( x, y, z ), and P M f ( j ) = (cid:110) P Mn ρ ( j ) mn ( x, y ) (cid:111) (37)is the result of operating on all N f n ( x, y, z ). Equation (37) will make the least amount of changeto the set { ρ ( j ) mn ( x, y ) } in the Euclidean sense such that the views calculated from the output of P M will equal the measured views I data n ( x, y ) for all n given any input ρ ( j ) mn ( x, y ). P M may be understoodas a “data satisfaction” projection operator that generates 3D models that agree with the measured2D views (one model for each view). This kind of problem is analogous to the situation studied by ? where an object is reconstructed from the averaged diffracted intensities from a number of object“clusters.” The “clusters” in this case are the individual sections in each 3D volume. The datasatisfaction projection operator P M results in a set of inconsistent models that need to be addressedby an additional projection operator, which we describe next.The second projection operator, P S , consists of the following steps: (1) Inverse propagate each ρ ( j ) mn ( x, y ) from the set { ρ ( j ) mn ( x, y ) } with the operator L − mn . (2) Form N
3D volumes from the setof all
M N inverse propagated sections. (3) Inverse rotate each volume to a common orientationby R − n . (4) Average the N volumes together (the concur step). (5) Apply any other constraintssuch as support, reality or positivity that may be available, collectively denoted here by the operator C . (6) Rotate each volume such that their orientation corresponds to their respective views by R n ,and then finally, (7) forward propagate the result with the forward propagation operator L mn and12edistribute the propagated sections to form an updated set of { ρ ( j ) mn ( x, y ) } . The projection operator P S as described by these steps can be written concisely as P S f ( j ) = L C L − f ( j ) (38) The number of views, N , is typically large in practice, thus the computer memory required to storethe iterate f can become an issue. We have therefore structured our algorithm in a way such thatthe IPA is done in groups. The set of all data views { I (data) n ( x, y ) } is split into G groups, each workedon by a separate IPA. The results from each group are averaged together after some number ofIPA iterations and that averaged volume becomes the new input to each IPA in the next iteration.This procedure is repeated some number of times. The result of this restructuring is that onlyapproximately N/G + 1 volumes need to be stored in memory at any given time (the +1 is for theaccumulator volume). Thus the memory demand of our algorithm can be controlled by selecting thenumber of groups G . We refer to the iterations where the output of the groups are averaged as the“outer loop”, indexed by k , and the iterations of the IPA itself as the “inner loop”, indexed as beforeby j . We denote the total number of iterations in each of these loops by K and J , respectively. Aflow chart of our overall algorithm is shown in Figure 7.A downside of this restructuring is that the operations no longer consist of making the least amountof change to the iterate to satisfy the set of all data views. The output of P M now satisfies thedata views only in their respective groups. Only in the limit of one group ( G = 1), does thealgorithm revert back to the original algorithm composed solely of proper projection operations.However this restructuring and grouping is necessary when there could potentially be hundreds ofthousands of views. Such an algorithm structure is also extremely amenable to parallelization, bothfor the individual groups of IPAs and also within the P M projection operator itself where each P Mn operation can be carried out in parallel across all n .13 …P Mn I data P Mn I A P Mn I data data … j = j +1 max inner iterations j No Yes P M P S … j = 1 ……P Mn P Mn data I ( G- A +2 I ( G- A +1 data I N data … j = j +1 No P Mn Yes P M P S j = 1No Yes … max inner iterations j max outer iterations k + k = k +1 k = 1 …… Figure 7: Flow diagram of the proposed algorithm with the iterative projection algorithm chosen to be in the errorreduction (ER) configuration for ease of illustration. Each green box contains a separate iterative projection algorithm,each running on a separate “group” of data. The variable A is the number of data views assigned to the first N − G is the total number of groups. The blue and red boxes contain the P S and P M projection operators,respectively. The iterate consists of the sets of 3D arrays shown enclosed in the square brackets. . Simulations We tested our algorithm on simulated 3D volumes. The particle that was used for the simulationswas the virus capsid of a tobacco necrosis virus (protein data bank (PDB) identification: 1TNV)determined by x-ray crystallography Bando et al. (1994). The virus capsid at the center of the unitcell was taken as the particle to be reconstructed with a size of approximately L = 29nm. Sincethe PDB file only provided the atomic coordinates of the capsid, the interior of the virus was filledwith the average electron density of the capsid for our simulations. The Python software package“reborn” Kirian et al. (2020) was used to convert the PDB file into a 3D array of real and positivevalues which we take as the scaled object potential f ( r ). The views were calculated via equation (16)with the freespace propagators, (27) and (28), implemented using the fast Fourier transform (FFT)algorithm. Rotations of the 3D volume required to compute the views, i.e., equation (14), wasimplemented via a series of 2D rotations of the planar slices of the 3D volume. Each 2D rotation wasin turn performed by three shears as described by Unser et al. (1995). The rotations were specifiedby three Euler angles, and the random rotations that were needed to generate randomly orientedviews were obtained by sampling those Euler angles.The IPA used for the reconstructions in this paper was the RAAR algorithm. Reality and positivityconstraints were applied, i.e., the operator C , introduced in equation (38), consists of setting theimaginary parts and the negative real parts of f ( r ) to zero. No other additional constraints wereapplied, in particular, the support constraint was not used. The starting 3D volume as input to thealgorithm was filled with uniformly distributed random values between 0 and 1. The progress of thereconstruction was monitored by calculating the errors e ( j,k ) = (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) G G (cid:88) g =1 (cid:80) r (cid:16) f ( j,k ) g ( r ) − f true ( r ) (cid:17) (cid:80) r ( f true ( r )) (39)and E ( j,k ) = (cid:118)(cid:117)(cid:117)(cid:117)(cid:116) G G (cid:88) g =1 (cid:80) n (cid:80) x,y (cid:16) I ( j,k ) gn ( x, y ) − I data gn ( x, y ) (cid:17) (cid:80) n (cid:80) x,y (cid:0) I data gn ( x, y ) (cid:1) , (40)where e ( j,k ) is the root-mean-squared error between the reconstructed virus at the j th and k thiterations of the inner and outer loop, and the ground truth, f true ( r ), where the mean is taken overall groups G ; similarly, E ( j,k ) is the root-mean-squared error between the set of all N views generatedby the iterate at the j th and k th iterations of the inner and outer loop, calculated according toequation (35), and the set of all data views, { I data n ( x, y ) } , where again, the mean is taken over allgroups. The errors are appended at the beginning of the next outer loop to the errors from theprevious outer loop to form the full error vectors, e and E , at the end of each inner loop iteration,such that the full error vectors can be written as e = [ e ( j, , · · · , e ( j,K ) ] (41) E = [ E ( j, , · · · , E ( j,K ) ] . (42)Reconstructions from two cases with data views that all have significant Ewald curvature are shownin the next two subsections. The first case is when the Ewald curvature arises due to a large incident15eam wavelength, and the second case is when a high resolution is desired. Two reconstructionsare carried out separately for each case, one uses the forward and inverse propagation operatorsas defined in equations (29) and (30), while the other reconstruction has the forward and inversepropagation operators changed to L mn = P ( z n ) (43) L − mn = P − ( z n ) , (44)i.e., with the Ewald sphere assumed flat. The reason for doing this is to allow a comparison betweenan algorithm which does not take Ewald curvature into account with the proposed algorithm whichdoes. The quality of the two different reconstructions are gauged by the Fourier shell correlation(FSC) metric Saxton & Baumeister (1982); van Heel & Schatz (2005), defined asFSC( q ) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:80) | q | = q F true ( q ) F ∗ recon ( q ) (cid:113)(cid:80) | q | = q | F true ( q ) | (cid:113)(cid:80) | q | = q | F recon ( q ) | (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (45)where F true ( q ) and F recon ( q ) are the Fourier transforms of the ground truth and the reconstructedvolumes, respectively. For this first case, the virus was calculated to d = 5 ˚A resolution and the incident electrons have awavelength λ = 0 .
34 ˚A, giving ε = 0 . . L = 29nm). The resultant virus is 68 × ×
68 voxels and the computational volume iszero-padded to 93 × ×
93 voxels. A total of 68 views were calculated from random rotations of thevirus, with randomly chosen focal distances ranging between a full length of the virus, L , on eitherside of the virus. Three example data views for this case of long wavelength are shown in Figure 8.The reconstruction from this dataset is shown in Figure 9 with the number of groups G = 3, an outerloop of K = 3 iterations, and an inner loop of J = 40 iterations. The RAAR algorithm parameterwas set to β = 0 .
7. The Fourier shell correlation for the reconstructions is shown in Figure 10.
For the second case, the virus density was calculated to d = 2 . λ = 0 . ε = 0 . L = 29nm). The resultant virus is 130 × ×
130 voxelsand the computational volume is zero-padded to 156 × ×
156 voxels. A total of 69 views werecalculated from random rotations of the virus, with randomly chosen focal distances ranging betweena full length of the virus, L , on either side of the virus. Three example data views for this case ofhigh resolution are shown in Figure 11. The reconstruction from this dataset is shown in Figure 12with the number of groups G = 2, an outer loop of K = 5 iterations, and an inner loop of J = 80iterations. The RAAR algorithm parameter was set to β = 0 .
9. The Fourier shell correlation isshown in Figure 13. 16 igure 8: Example data views for the case of long wavelength at three different focal distances z n and orientations R n with significant Ewald curvature ( ε = 0 . igure 9: Reconstruction results for the long wavelength case. Orthogonal central slices of the ground truth (Row1); the reconstructed particle with the proposed algorithm (Row 2); the reconstructed particle with the proposedalgorithm but ignoring Ewald curvature (Row 3). The two sets of reconstructions are scaled to have the same meanas the ground truth. The greyscale of the two reconstructions are also set to be the same as that for the ground truth.(Bottom row) Errors as the algorithm progresses on a logarithmic scale. The solid line is the object error, e , and thedashed line is the data error, E . Blue is the proposed algorithm, orange is the proposed algorithm but ignoring Ewaldcurvature. igure 10: Fourier shell correlations for the long wavelength case, comparing the ground truth virus density with thereconstruction from the proposed algorithm (blue open squares) and the reconstruction with the proposed algorithmbut ignoring Ewald curvature (orange open circles). The 1 / igure 11: Example data views for the case of high resolution at three different focal distances z n and orientations R n with significant Ewald curvature ( ε = 0 . igure 12: Reconstruction results for the high resolution case. Orthogonal central slices of the ground truth (Row1); the reconstructed particle with the proposed algorithm (Row 2); the reconstructed particle with the proposedalgorithm but ignoring Ewald curvature (Row 3). The two sets of reconstructions are scaled to have the same meanas the ground truth. The greyscale of the two reconstructions are also set to be the same as that for the ground truth.(Bottom row) Errors as the algorithm progresses on a logarithmic scale. The solid line is the object error, e , and thedashed line is the data error, E . Blue is the proposed algorithm, orange is the proposed algorithm but ignoring Ewaldcurvature. igure 13: Fourier shell correlations for the high resolution case, comparing the ground truth virus density with thereconstruction from the proposed algorithm (blue open squares) and the reconstruction with the proposed algorithmbut ignoring Ewald curvature (orange open circles). The 1 / . Discussion and Conclusion An algorithm for reconstructing an object from projection images affected by Ewald sphere curva-ture in cryo-electron microscopy is proposed. This algorithm was shown via simulations to be ableto reconstruct the correct 3D object from a set of 2D near-field intensities, which we have called“views,” that are affected by significant Ewald curvature. A criterion for the Ewald curvature tobecome significant is derived, relating together the three key parameters of (1) wavelength of theincoming wave, (2) desired resolution, and (3) thickness of the sample. Only the reality and positivityconstraints are applied in our reconstructions, i.e., the object is assumed to be real and positive. Asupport constraint is not used.The algorithm is based on the paradigm of iterative projection algorithms (IPAs), and a projectionoperation that makes the minimum change to a set of complex numbers such that the sum of theirimaginary parts is equal to a desired value is derived. A restructuring of the traditional IPA loopwas proposed which alleviates the memory requirement of the algorithm when the dataset containsmany views. An implication of this restructuring is that the operations no longer consist of makingthe least amount of change to the iterate to satisfy the set of all data views, because the output ofthe projection operator P M now satisfies the data views only in their respective groups. Howeverthis is necessary for practical applications when there could potentially be hundreds of thousands ofviews. Such an algorithm structure is extremely amenable to parallelization, both for the individualgroups of IPAs and also for the P M projection operation.For this proof-of-concept work, we have not parallelized the algorithm. The computational com-plexity of our algorithm in serial implementation is O ( N log N ) where N is the number of voxelsin each dimension of the 3D array. A factor N log( N ) is due to the 3D FFT, and the remainingfactor of N comes from the fact that the number of 2D views needed for a unique reconstruction isproportional to N . For unoptimized Python code and working solely with double-precision floating-point numbers, the reconstruction detailed in Section 4.1, the long wavelength case, took around 19hours. The reconstruction detailed in Section 4.2, the high resolution case, took around 8.5 days.Both reconstructions were carried out with a single core on ASU’s Agave research computing systemutilizing Intel Broadwell CPUs. Larger objects may require more iterations of the algorithm. Thecompute time can likely be reduced if the initial starting iterate of the algorithm was a preliminaryreconstruction obtained by treating the views as projections, i.e., ignoring Ewald curvature. In termsof storage requirement, our code stores the set of all ρ mn ( x, y ) as the iterate f . If we switch to storingthe set of all f n ( x, y, z ) then that would further reduce the amount of storage needed. The tradeoff istime, because ρ mn ( x, y ) will have to be generated from f n ( x, y, z ) by zero-padding and Fourier trans-forming at every iteration. As mentioned before, the greatest speed-up would come from parallelizingthe algorithm, which we are working towards.In the limit where the Ewald sphere is flat, our algorithm becomes a tomographic reconstructionmethod, capable of recovering a 3D object from conventional tomographic data where each view isjust a simple projection (sum of densities) through the object.Our method assumes the use of a conventional through-focus series for image reconstruction. Incomparison with an alternative method based on the far-out-of-focus spatial separation of spatialfrequencies Russo & Henderson (2018), while our approach lacks the benefits of masking in thatmethod, it may have particular advantages when used with a Zernike phase plate ? . Using in-focus23mages, this preserves low spatial frequencies in image formation, otherwise lost in bright-field out-of-focus images of a weak phase object. Phase plates thus avoid the need for the very large defocus(with loss of high-resolution detail) required to obtain visibility from the smallest particles, so thatour method could then give access to these smaller biomolecules in this way. It may also reducedamage by reducing the number of image recordings needed. Our approach is more computationallyintensive, but should be simpler to implement experimentally.Our algorithm assumes we know the orientation and defocus of the object in each view, both of whichshould be obtainable from existing cryo-EM software. For cryo-EM tomography, the resolution limitimposed by radiation damage means that Ewald sphere curvature is unlikely to be important, sothat we have assumed that this method will be applied to single-particle data. Additional analysisto allow for conformation variation may be considered in future developments of this method. Thecontrast transfer function (CTF) is ignored here in this work. In the event that the CTF is known,and is assumed to be the same for all views, then it can be incorporated into the current algorithmand does not require change to the overall structure of the method. If the CTF is unknown, then thelens aberration parameters may be recoverable in addition to the object potential, using a modifiedversion of the proposed algorithm, a line of inquiry we are currently investigating. Further futurework includes exploring the effect of varying some of the parameters of the algorithm, such as thenumber of groups, G , the max inner iteration, J , and the max outer iteration, K , and extending thealgorithm to deal with multiple scattering.The code for the simulations carried out in this paper can be found at: https://gitlab.com/jpchen1/em-reconstruction-with-ewald Acknowledgements : We are grateful to Prof R. M. Glaeser for useful conversations in connectionwith this work during his visit to the new Cryo-EM center at ASU. We also thank Cornelius Gatifor many helpful discussions and the two anonymous referees for their encouragement and many con-structive comments. We acknowledge the ASU Agave computing cluster and NVIDIA Corporationfor their Titan V GPU. JPJC and RAK acknowledge support from NSF STC Award DBI-1231306.JPJC, KES and RAK acknowledge support from NSF Award DBI-1565180. JCHS acknowledgesupport from ARO award AWD00035320.
Appendix A. Potential scattering
Assuming non-relativistic elastic scattering, we start with the time-independent Schr¨odinger equation − (cid:126) m ∇ ψ ( r ) − eU ( r ) ψ ( r ) = Eψ ( r ) , (A.1)where U ( r ) is the electric potential, e , m and E are the charge, mass and kinetic energy of the electron,respectively. Define k = √ mE/ (cid:126) and V ( r ) = − eU ( r ) we can write down the inhomogeneousHelmholtz’s equation ( ∇ + k ) ψ ( r ) = 2 m (cid:126) V ( r ) ψ ( r ) . (A.2)24he integral form of equation (A.2) is the Lippmann-Schwinger equation, which, assuming an in-coming plane wave e ik z , can be written as ψ ( r ) = e ik z + 2 m (cid:126) (cid:90) d r (cid:48) G ( r , r (cid:48) ) V ( r (cid:48) ) ψ ( r (cid:48) ) , (A.3)where G ( r , r (cid:48) ) is the Green’s function of the homogenous Helmholtz’s equation, G ( r , r (cid:48) ) = − e ik | r − r (cid:48) | π | r − r (cid:48) | . (A.4)Making the first Born approximation and factoring out the e ik z term in equation (A.3), gives thescattered wave ψ ( r ) = e ik z (cid:18) m (cid:126) (cid:90) d r (cid:48) G ( r , r (cid:48) ) V ( r (cid:48) ) e − ik ( z − z (cid:48) ) (cid:19) (A.5)as the solution to our scattering problem. To cast this solution into the form of freespace propagators,first write the Helmholtz Green’s function as G ( r , r (cid:48) ) = (cid:90) d k (2 π ) e i k · ( r − r (cid:48) ) k − k + iη , (A.6)where η is a positive infinitesimal. Evaluating the k z integral in equation (A.6) gives, G ( r , r (cid:48) ) = − i (cid:90) dk x π (cid:90) dk y π e ik x ( x − x (cid:48) ) e ik y ( y − y (cid:48) ) e i √ k − k x − k y | z − z (cid:48) | (cid:113) k − k x − k y , (A.7)where the square root is defined to give decaying exponentials if k x + k y > k . Substituting equa-tion (A.7) into equation (A.5), taking z > z (cid:48) , and replacing the square root in the denominator by k due to the small angle scattering geometry typical in electron microscopy experiments gives ourfinal result ψ ( r ) = e ik z (cid:20) − mi (cid:126) k (cid:90) d r (cid:48) (cid:90) dk x π (cid:90) dk y π e ik x ( x − x (cid:48) ) e ik y ( y − y (cid:48) ) e i (cid:16) √ k − k x − k y − k (cid:17) ( z − z (cid:48) ) V ( x (cid:48) , y (cid:48) , z (cid:48) ) (cid:21) . (A.8)Writing E = (cid:126) k m and λ = πk , the prefactor of the integral can be written as mi (cid:126) k = i πλE . Appendix B. Projection operator for constraining the sum of imaginary numbers
In this appendix we derive the operation that makes the minimum change to a set of complex numberssuch that the sum of their imaginary parts becomes equal to a desired value. Denote the desiredvalue by I data , and the original set of N complex numbers by z on = x on + iy on . The sum of all imaginaryparts from the N complex numbers is then I o = N (cid:88) n =1 y on . (B.1)25e would like the set of complex numbers after the projection operation, denoted by z n = x n + iy n ,to have (cid:80) Nn =1 y n = I data .The function we wish to minimize is the Euclidean distance f = N (cid:88) n =1 (cid:2) ( x n − x on ) + ( y n − y on ) (cid:3) . (B.2)The constraint equation is g = N (cid:88) n =1 y n − I data . (B.3)Applying the method of Lagrange multipliers leads us to write L = f + λg , (B.4)where λ is the Lagrange multiplier. Taking the partial derivatives of L with respect to x n and y n and setting them to zero gives 0 = ∂L∂x n = 2( x n − x on ) (B.5)0 = ∂L∂y n = 2( y n − y on ) + λ , (B.6)yielding x n = x on (B.7) y n = y on − λ . (B.8)Summing over all n in equation (B.8) and rearranging for λ gives λ = 2 N N (cid:88) n =1 (cid:0) y on − y n (cid:1) . (B.9)Substituting equation (B.9) back into (B.8) and combining the real and imaginary parts using (B.7)finally yields z n = z on + i N (cid:0) I data − I o (cid:1) . (B.10)Equation (B.10) is the projection operation that makes the minimum change to a set of complexnumbers, measured by the Euclidean distance, such that the sum of their imaginary parts is equalto I data . 26 eferences Bando, M., Morimoto, Y., Sato, T., & Tsukihara, T. (1994). Crystal structural analysis of tobacconecrosis virus at 5˚A resolution.
Acta Cryst. D , , 878–883.Bauschke, H. H., Combettes, P. L., & Luke, D. R. (2002). Phase retrieval, error reduction algorithm,and fienup variants: a view from convex optimization. JOSA A , , 1334–1345.Donatelli, J. J., & Spence, J. C. H. (2020). Inversion of many-beam bragg intensities for phasingby iterated projections: Removal of multiple scattering artifacts from diffraction data. Phys. Rev.Lett. , .Elser, V. (2003). Phase retrieval by iterated projections. J. Opt. Soc. Am. A , , 40–55.Elser, V., Rankenburg, I., & Thibault, P. (2007). Searching with iterated maps. PNAS , , 418–423.Fienup, J. R. (1982). Phase retrieval algorithms: a comparison. Appl. Opt. , , 2758–2769.Gerchberg, R. W., & Saxton, W. O. (1972). A practical algorithm for the determination of the phasefrom image and diffraction plane pictures. Optik , .Gravel, S., & Elser, V. (2008). Divide and concur: A general approach to constraint satisfaction. Phys. Rev. E , , 036706.Gureyev, T. E., Quiney, H. M., Kozlov, A., & Allen, L. J. (2019a). On electron imaging of smallmolecules, part I: relative roles of multiple scattering and fresnel diffraction, . URL: https://arxiv.org/abs/1911.11599 .Gureyev, T. E., Quiney, H. M., Kozlov, A., & Allen, L. J. (2019b). On electron imaging of smallmolecules, part II: tomographic reconstruction from defocus series, . URL: https://arxiv.org/abs/1911.11598 .van Heel, M., & Schatz, M. (2005). Fourier shell correlation threshold criteria. J. Struc. Bio. , ,250–262.Kirian, R. A. et al. (2020). reborn software. https://kirianlab.gitlab.io/reborn/index.html .[Online; accessed 9 June 2020].Lentzen, M. (2014). No surprise in the first born approximation for electron scattering. Ultrami-croscopy , , 201–210.Leong, P. A., Yu, X., Zhou, Z. H., & Jensen, G. J. (2010). Correcting for the ewald sphere inhigh-resolution single-particle reconstructions. In Methods in Enzymology (p. 369–380). Elseviervolume 482. URL: https://linkinghub.elsevier.com/retrieve/pii/S0076687910820154 .doi: .Luke, D. R. (2005). Relaxed averaged alternating reflections for diffraction imaging.
Inverse Prob-lems , , 37–50.Marchesini, S. (2007). A unified evaluation of iterative projection algorithms for phase retrieval. Rev.Sci. Instrum. , , 001301. 27illane, R. P., & Lo, V. L. (2013). Iterative projection algorithms in protein crystallography. I.theory. Acta Cryst. A , , 517–527.Ren, D., Ophus, C., Chen, M., & Waller, L. (2020). A multiple scattering algorithm for threedimensional phase contrast atomic electron tomography. Ultramicroscopy , .Russo, C. J., & Henderson, R. (2018). Ewald sphere correction using a single side-band imageprocessing algorithm. Ultramicroscopy , , 26–33. doi: .Saxton, W. O., & Baumeister, W. (1982). The correlation averaging of a regularly arranged bacterialcell envelope protein. J. Microscopy , , 127–138.Spence, J. C. (2013). High-Resolution Electron Microscopy . Oxford University Press.Unser, M., Thevenaz, P., & Yaroslavsky, L. (1995). Convolution-based interpolation for fast, high-quality rotation of images.
IEEE Transactions on Image Processing , , 1371.Wolf, M., DeRosier, D. J., & Grigorieff, N. (2006). Ewald sphere correction for single-particle electronmicroscopy. Ultramicroscopy , , 376–382. doi: .Zhu, D., Wang, X., Fang, Q., Van Etten, J. L., Rossmann, M. G., Rao, Z., & Zhang, X. (2018).Pushing the resolution limit by correcting the ewald sphere effect in single-particle cryo-em recon-structions. Nature Communications , , 1552. doi: .Zuo, J. M., & Spence, J. C. (2017). Advanced Transmission Electron Microscopy . SpringerNew York. URL: http://link.springer.com/10.1007/978-1-4939-6607-3 . doi:10.1007/978-1-4939-6607-3