Low-dose cryo electron ptychography via non-convex Bayesian optimization (original) (raw)

Abstract

Electron ptychography has seen a recent surge of interest for phase sensitive imaging at atomic or near-atomic resolution. However, applications are so far mainly limited to radiation-hard samples, because the required doses are too high for imaging biological samples at high resolution. We propose the use of non-convex Bayesian optimization to overcome this problem, and show via numerical simulations that the dose required for successful reconstruction can be reduced by two orders of magnitude compared to previous experiments. As an important application we suggest to use this method for imaging single biological macromolecules at cryogenic temperatures and demonstrate 2D single-particle reconstructions from simulated data with a resolution up to 5.4 Å at a dose of 20_e_ −/Å2. When averaging over only 30 low-dose datasets, a 2D resolution around 3.5 Å is possible for macromolecular complexes even below 100 kDa. With its independence from the microscope transfer function, direct recovery of phase contrast, and better scaling of signal-to-noise ratio, low-dose cryo electron ptychography may become a promising alternative to Zernike phase-contrast microscopy.

Introduction

The advent of direct electron detectors has led to a resolution revolution in the field of cryo electron microscopy in the last few years. The technique is now producing three-dimensional atomic potential maps of biological macromolecules of a few 100 kDa or lower with a resolution better than 3.5 Å1–3, such that individual amino acid side-chains can be resolved. An important role in this revolution play new image processing algorithms based on a Bayesian approach, which infer important parameters without user intervention4. Also the correction of beam-induced motion has become possible mostly due to the new generation of detectors5, 6. However, several challenges remain to be overcome in order to routinely reach 3 Å resolution also for small complexes7, 8: Firstly, beam-induced specimen charging and subsequent motion currently still render the high resolution information of the first few frames of a high repetition rate movie recorded with a direct electron detector unusable9, because the motion is too fast to efficiently correct for it. Secondly, the Detective Quantum Efficiency (DQE) of detectors is still imperfect at high spatial frequencies8, 10. Thirdly, the contrast of single images can still be improved to enable reconstructions with fewer particles and increase the throughput10.

The last of these challenges has recently been addressed with a new phase plate model11, 12, which is comparatively simple to use and provides excellent contrast at low spatial frequencies. In addition to this hardware-based approach to achieve linear phase contrast in the measured amplitudes, discovered by Zernike in the 1930s13, it is also possible to algorithmically retrieve the phase information from a set of coherent diffraction measurements. One such technique, commonly known as ptychography or scanning coherent diffractive microscopy 14, is becoming increasingly popular in the field of materials science due to experimental robustness and the possibility to obtain quantitative phase contrast over an essentially unlimited field of view15, 16. The use of ptychography for imaging radiation sensitive samples with electrons at high resolution is however precluded so far by its high dose requirements.

Here, we show how the use of non-convex Bayesian optimization to solve the ptychographic phase retrieval problem fulfills the dose requirements for imaging biological macromolecules and makes it possible to obtain 2D images from single particles with sub-nanometer resolution. After a short introduction into the technique, we will also mention how ptychography offers improvements for the other two challenges discussed above.

Despite being initially proposed as a solution to the phase problem for electrons17, 18, ptychography has seen its biggest success in X-ray imaging, due to the less stringent sample requirements and the experimental need for lensless imaging techniques. Recent developments include the introduction of iterative algorithms to enable the reconstruction of datasets collected with an out-of focus probe19, 20, which decreases the memory requirements of the method dramatically. The algorithms also have the capability for the correction of experimental difficulties such as unknown scan positions21–23, partial coherence24, probe movement during exposure25, 26, intensity fluctuations during the scan24, 27 and reconstruction of background noise15, 27.

In recent years, some of these advances have been applied in the context of electron microscopy and yielded atomic resolution reconstructions of low-atomic number materials28–30 and quantitative phase information15.

Figure 1 shows the experimental set-up for an out-of focus ptychography experiment. A ptychographic dataset is collected by scanning a spatially confined, coherent beam, subsequently called ‘probe’, over the specimen and recording far-field diffraction patterns at a series of positions such that the illuminated regions of neighboring positions overlap. The diffraction-limited resolution r d (half-period) of the final image is given by rd=λ⋅ΔzNpix⋅dpix, where N pix is the number of detector pixels, d pix is the detector pixel size and λ is the de-Broglie-wavelength of the electrons. Given the set of positions r⃗1…K and a realistic forward model for the formation of the corresponding diffraction patterns I(1…K), the complex-valued transmission function T(r⃗), which describes the atomic properties of the specimen31, can then be retrieved by solving a non-convex inverse problem. The electron dose used for successful reconstruction has exceeded 1 × 103 e −/Å2 so far, limiting the usability of ptychography to radiation-hard specimens. Table 1 lists recently published electron ptychography experiments and the used average electron doses. The lowest dose was reported in ref. 32, which used an estimated 3.33 × 103 e −/Å2 at a resolution of dr = 58.4 pm, resulting in a dose of 1.1 × 103 e −/pixel and achieving a line resolution of 2.3 Å, demonstrated by resolving the lattice spacing of gold nanoparticles. We will demonstrate via simulations that it is possible to reduce this dose by a factor of 100, thus reaching the dose range allowed for imaging biological macromolecules.

Figure 1.

Figure 1

Experimental geometry in ptychography. The coherent electron wave function ψ(r⃗) illuminates several regions (centered at r⃗1…K) across the sample, which is characterized by the transmission function T(r⃗). For each position, a 2D diffraction pattern I(1…K) is recorded in the far field at distance Δ_z_. The sample thickness t can be neglected for biological macromolecules in the reconstruction at the resolutions presented in this paper.

Table 1.

List of previously published electron ptychography experiments and used electron dose.

Reference resolution e −/Å2
D’Alfonso et al.30 ~1.5 Å 1.77 × 104
Yang et al.29 atomic 1.3 × 104
Putkunz et al.28 ~1 Å 9.2 × 106
Humphry et al.32 ~2.3 Å 3.33 × 103

The problem of beam-induced sample movement has already been addressed before the development of fast direct detectors. Scanning with small spots of several 10 nm in size over a vitrified sample has shown to reduce beam induced specimen movement33–35 in real-space imaging and it has been noted that the remaining movement may be due to radiation damage, not sample charging33. Ptychography naturally operates with a confined beam, thus minimizing the area where charge can build up, such that the movement should be reduced compared to the illumination of large areas in cryo-EM. The sampling requirements given by the experimental setup allow to operate the detector with large effective pixel sizes, such that the DQE and MTF are at near-perfect values and can be neglected in the reconstruction. The overhead resulting from the need to take multiple exposures can therefore be greatly reduced by using fast detectors with few pixels.

Results

Image formation of cryo-TEM and ptychography

We perform multislice simulations of three different biological macromolecules with molecular weights ranging from 64 kDa to 4 MDa. We choose the 64 kDa hemoglobin36, the 706 kDa 20S proteasome from yeast2, and the 4 MDa human ribosome37. We create atomic potential maps using the Matlab code InSilicoTem38, with a thickness of 50 nm and at an electron energy of 300 keV. We use the isolated atom superposition approximation, without solving the Poisson-Boltzmann equations for the interaction between the molecule and the ions. We also do not model the amorphousness of the solvent, which was performed in ref. 38 using molecular dynamics simulations, but was seen to have a negligible effect at very low doses. As described in ref. 38, we model the imaginary part of the potential via the inelastic mean free path, creating a minimal transmission contrast between the vitreous ice and the protein. Using these potential maps, we simulate a ptychography experiment by cropping three-dimensional slices from the potential at several positions and propagate a coherent incoming wave through the slices using the methods described in ref. 39 in custom code. The final model for the formation of the intensity on the detector is then

for the diffraction pattern and

| I0(q⃗)=|ℱ−1[ℱ[ψexit(r⃗)]⋅CTF(q⃗)]|2 | 2 | | ----------------------------------------- | - |

for the cryo-EM image, where ℱ and ℱ−1 denote the forward and inverse Fourier transform, and CTF the contrast transfer function, respectively.

The detector properties are modeled as in ref. 38, by multiplying the Fourier transform of the exit-wave intensity with the square root of the detective quantum efficiency DQE(q⃗) 40, before applying shot noise by sampling from a Poissonian distribution. We finally convolve the noisy signal with the noise transfer function NTF(q⃗) to yield the measured intensity.

A notable difference both in simulation and practice is the fact that for cryo-EM, usually no pixel binning is applied to maximize the imaged area and increase throughput. Therefore, also high spatial-frequency regions with low values of DQE and NTF are used for image formation41. For ptychography, on the other hand, the detector can be heavily binned, as long as the real-space patch given by λΔz/dpix≡rd⋅Npix still encompasses the probe beam on the sample. For typical detectors this condition is fulfilled at bin sizes equivalent to a few percent of the Nyquist frequency. This leads to a near-constant DQE and a near-unity NTF, such that they can be omitted in the ptychography reconstructions, whereas we still include them in the simulation of the diffraction data. We note, however, that a convolution with a detector transfer function can be modeled with a partially coherent beam if necessary, as demonstrated in refs 42 and 43. We choose the Gatan K2 Summit as the detector for our simulations because it has the highest published DQE and MTF values at low spatial frequencies at 300 keV41. We note that direct detection cameras with frame rates of 1 kHz and above may be more suitable for high-throughput scanning experiments44–46, but characteristics for these cameras at 300 keV are either not published or inferior to the K2 Summit. Assuming the K2 Summit for both ptychography and phase-contrast TEM simulations also simplifies a direct comparison between the two methods. The intensity after detection is modeled as38:

I(q⃗)=ℱ−1[ℱ[Poisson(ℱ−1[ℱ[I0(q⃗)]⋅DQE(q⃗)])]⋅NTF(q⃗)], 3

where NTF and DQE are properties of the detector40, 45 and Poisson(x) samples from a Poisson distribution with mean x.

Single-particle reconstruction

Figure 2 shows a comparison of low-dose ptychography reconstructions with currently used methods for single-particle imaging with electrons: defocus-based cryo-EM, and Zernike phase contrast cryo-EM with a Volta phase-plate. We choose exemplary doses of 5 e −/Å2 as the typical threshold where the highest resolution details are destroyed47 and 20 e −/Å2 as a typical dose at which experiments are performed. We have reversed the contrast in the cryo-EM images to simplify the visual comparison with the ptychography reconstructions. To quantitatively assess the image quality, we have computed the 2D Fourier Ring Correlation (FRC)48 with the ground truth for the both ptychographic reconstruction and simulated cryo-EM images of the macromolecules, as shown in Fig. 3. As ground truth for the images we use the electron counts in a noiseless, aberration-free phase-plate image. Using the 1-bit criterion as a resolution threshold48, the achieved resolutions at 5 e −/Å2 and 20 e −/Å2, respectively, are 12 Å and 8.9 Å for hemoglobin; 10.9 Å and 9.1 Å for 20S proteasome; and 10.3 Å and 5.4 Å for human ribosome. In the case 20S proteasome, these values are identical to the FRC threshold for the phase plate image; for hemoglobin and human ribosome, the phase plate image yields a slightly better resolution of 8.7 Å and 5.1 Å respectively at a dose of 20 e −/Å2 and 10 Å at a dose of 5 e −/Å2. As the FRC is insensitive to very small and very large values of signal-to-noise ratio (SNR), we also show the spatial-frequency resolved SNR in Fig. 3(d)–(f). We define the SNR as

| SNR(q)=10⋅log10(|ℱ[T(r⃗)]|2|ℱ[T(r⃗)]−ℱ[Tmodel(r⃗)]|2)dB. | 4 | | ---------------------------------------------------------------- | - |

Figure 2.

Figure 2

Cryo-electron ptychography reconstructions from simulated data and simulated cryo-EM images for different doses and 3 macromolecules with growing molecular weights in columns 1–3. Row (a) Phase of the transmission function, the ground truth for the ptychography reconstructions. The scale bar next to the figures is in rad. Rows (b) and (e) ptychography reconstruction at doses of 5 e −/Å2 and 20 e −/Å2. Rows (c) and (f) Simulated cryo-EM image with a defocus of 1.6 μm at a dose of 5 e −/Å2 and 20 e −/Å2. Rows (d) and (g) Simulated cryo-EM image with a Zernike phase plate and a defocus of 50 nm at doses of 5 e −/Å2 and 20 e −/Å2. Column (1) hemoglobin, column (2) 20S proteasome, column (3) human ribosome.

Figure 3.

Figure 3

FRC (a)–(c) and SNR (d)–(f) as a function of spatial frequency for the cryo-electron ptychography reconstructions and simulated cryo-EM images in Fig. 2.

The SNR of the ptychographic reconstruction is significantly lower than the SNR of the phase-plate image for all three particles at spatial frequencies below 0.1 Å−1. It does however scale better to high spatial frequencies, and does not drop below −15 dB for all particles and all resolution up to nearly 0.5 Å−1, where the SNR is close to two orders of magnitude better than the phase-plate image SNR. This helps ptychography perform better when multiple reconstructions are averaged, because a positive single-digit SNR can be reached with fewer particles.

Effect of averaging

In single-particle cryo-EM, a three-dimensional reconstruction can be obtained by collecting a large ensemble of 2D images, before orienting and averaging them in three dimensions, such that a satisfactory SNR is achieved. A similar 3D reconstruction from ptychographic data is out of the scope of this paper. A straightforward approach would be to use the reconstructed 2D phase images as an input to existing 3D reconstruction algorithms4, 49, as is done routinely in ptychographic X-ray tomography50; however, it is likely that a dedicated algorithm which reconstructs the 3D structure directly from the diffraction patterns will achieve the best results. Coarse orientation alignment for such an algorithm could be done in real space from 2D reconstructions shown here. We leave the evaluation of such approaches to future studies, and for now concentrate on the achievable 2D resolution when averaging a larger ensemble of particles.

To give a rough estimate how the resolution and SNR achieved by our algorithm scales with averaging over multiple datasets, we simulate 30 independent datasets with identically oriented particles, and average the reconstruction results. We found that the resolution corresponding to the diffraction limit as defined by the probe-forming aperture is achieved with roughly 40 averaged reconstructions (see the Supplementary Information). Superresolution beyond this point is in principle possible, but due to the low electron counts and the low contrast transfer for higher angles, more images are needed for further improvements. To limit the amount of necessary computation we use an average of 30 images, where a resolution of roughly 3 Å is reached, close to the resolution corresponding to the probe-forming aperture of 2.7 Å. Figure 4(a)–(f) show images of the averaged reconstructions of the three samples, at doses of 5 e −/Å2 and 20 e −/Å2 respectively. We also compare the FRC between respectively 30 averaged reconstructions of 60 independently created ptychographic data sets, to give a resolution estimate. We use here the 1/2-bit resolution threshold discussed in ref. 48, which gives a slightly more conservative estimate than the 0.143-criterion commonly used for the 3D Fourier Shell correlation in averaged reconstructions for 3D cryo-EM. With averaging, a resolution of 3.4 Å is achieved for hemoglobin, 3.1 Å for 20S proteasome and 2.9 Å for human ribosome, all at a dose of 20 e −/Å2. This shows that cryo-electron ptychography can recover atomic resolution information in 2D from only tens of averaged images, facilitated by the favorable scaling of the SNR with spatial frequency, as discussed above. It may therefore be possible to drastically reduce the number of particles that is required for a successful 3D reconstruction at atomic resolution. A rough tomographic estimate51 suggests that the ribosome could be reconstructed to 3 Å resolution in 3D with less than 10000 particles.

Figure 4.

Figure 4

Average over 30 ptychographic reconstructions from independent data sets for (a) hemoglobin with 5 e −/Å2, (b) hemoglobin with 20 e −/Å2, (c) proteasome 20S with 5 e −/Å2, (d) proteasome 20S with 20 e −/Å2, (e) human ribosome with 5 e −/Å2, (f) human ribosome with 20 e −/Å2. FRC of averaged reconstructions from independent data sets for (g) hemoglobin, (h) proteasome 20S (i) human ribosome.

Probe and dose dependence

It is well-known that the phase profile of the ptychographic probe can heavily influence the reconstruction quality52–56. Here, we numerically test three different probes, depicted in Fig. 5, and their influence on the reconstruction SNR at low and high doses: (1) a standard defocused probe with defocus aberration of 400 nm, (2) a defocused Fresnel Zone Plate (FZP), and (3) a randomized probe generated by a holographic phase plate and a conventional lens. Figure 5 depicts these probe in real and Fourier space, and typical diffraction patterns at infinite dose and low dose. The FZP was recently suggested as a phase modulator for bright-field STEM57, because its simple phase modulation allows analytical retrieval of linear phase contrast. However, diffractive optics typically have imperfections due to the manufacturing process, which introduce errors and dose inefficiency if the phase modulation is obtained by a simple fitting procedure. Iterative ptychography algorithms allow for the simultaneous retrieval of the probe wave function19, 20, and therefore offer full flexibility in the design of the phase profile. Empirically, probes with a diffuse phase profile result in better reconstructions; therefore, we test as a third probe a random illumination generated by a holographic phase plate and a focusing lens.

Figure 5.

Figure 5

Different probes evaluated in this paper and corresponding diffraction patterns. Row 1: defocused beam with defocus aberration of 400 nm, convergence half-angle 9.2 mrad; row 2: defocused beam created by a Fresnel zone plate, 600 nm from focus; row 3: randomized beam, generated by a holographic phase plate and focused by a conventional lens. Column (a) beam in real space, at the sample position, scale bar is 8.5 nm; column (b) beam at the probe forming aperture, scale bar is 4.5 mrad; column (c) diffraction pattern of human ribosome at unlimited dose, normalized to the maximum intensity; column (d) diffraction pattern for a scan with an electron dose of 20 e −/Å2. The inset in 1a shows the color wheel that is used to represent amplitude and phase in columns (a) and (b).

Figure 6 shows the SNR of the three proposed probes as a function of spatial frequency for doses of 20 e −/Å2 and 80 e −/Å2. It can be seen that the simple defocused probe has almost 2 orders of magnitude worse SNR than the FZP and the random probe at the lowest spatial frequencies. At low spatial frequency the FZP achieves the best SNR, while at high spatial frequencies the random probe does slightly better. We have therefore used the random probe for the reconstructions shown in Figs 2 and 4. We give a qualitative explanation of this fact, but emphasize that a theory for optimal measurement design in ptychography and a practically feasible implementation of it is still outstanding, and may drastically improve upon the results presented here. From Fig. 5 column c) it can be seen, that the randomized probe induces the strongest intensity fluctuations (speckle) in the diffraction pattern. These very localized fluctuations vary strongly when scanning the sample, while the diffraction pattern from the defocused probe has slowly-varying features, which correlate strongly with the real space transmission function and vary only weakly when scanning the sample, which leads to a less diverse dataset and can cause problems in the reconstruction.

Figure 6.

Figure 6

SNR of reconstructions of the human ribosome at different radiation doses using the defocused probe, the Fresnel zone plate and the randomized probe.

Methods

Mathematical framework of Ptychography

We define the two-dimensional grid with size n1×n2∈ℕ×ℕ and length scale r > 0 as Drn1×n2:=(rα,rβ)α,β=0n1,n2⊂ℝ2. The two-dimensional complex transmission function of the object is discretized as a n 1 × n 2 matrix and denoted as T:Drdn1×n2→ℂ, where r d > 0 is the diffraction-limited length scale as introduced above. The object is illuminated by a small beam with known distribution, and discretized as a m 1 × m 2 matrix, denoted as ψ:Drdm1×m2→ℂ. For simplicity, in this paper we only consider the case n 1 = n 2 and m 1 = m 2, i.e. a uniform discretization in both axes. In the experiment, the beam is moved over the sample to positions r⃗i, and illuminates K > 1 subregions to obtain K diffraction images. The intensity measured for position i is then

| Ii(q⃗)=|ℱ[ψ(r⃗+r⃗i)⋅T(r⃗)]|2,i∈{0,…,K}, | 5 | | ------------------------------------------- | - |

where the real-space coordinates are discretized in steps of r d, and reciprocal-space coordinates in steps of (m {1,2} r d)−1. Mathematically, ptychographic reconstruction can be understood as a special case of the generalized phase retrieval problem: given a phase-less vector of measurements y∈ℝ+m find a complex vector z∈ℂn such that

where A:ℂn→ℂm is an arbitrary linear operator. We follow the notations in ref. 52 to write the ptychography problem in this form. First, we vectorize the transmission function as TV∈ℂN with N=n1⋅n2∈ℕ. We introduce the matrix R(i)∈ℝM×N, which extracts an m 1 × m 2 sized area centered at position r⃗i from T V. With these notations in place, the relation between the noise-free diffraction measurements collected in a ptychography experiment and T V can be represented compactly as

where P is constructed by cropping K regions from T V and multiplying by the incoming beam in Q, and applying a 2D discrete Fourier transform F, i.e. P = FQ:

| [I1⋮IKM]⏞I∈RKM=|[F⋯0⋮⋱⋮0⋯F]⏞F∈CKM×KM[diag(ψ)R(1)⋮diag(ψ)R(K)]⏞Q∈CKM×NTV|2. | 8 | | ---------------------------------------------------------------------------------- | - |

The matrix P∈ℂKM×N is sometimes called design matrix, because its entries determine the measurement outcome and reflect the experimental design. In the last decades many algorithms to solve this problem have been devised, only a few of which we will review with regards to low-dose reconstruction in the following section. For the subsequent analysis, we denote the KM row vectors of P as p i.

Bayesian optimization with truncated gradients

The most prominent iterative algorithms to solve the ptychographic phase retrieval problems are the difference map (DM) algorithm14, and the extended ptychographic iterative engine (ePIE)20. The difference map belongs to the family of algorithms which use projections onto non-convex sets to reach a fix-point, i.e., the solution lying at the intersection of the two sets. It can be shown that the standard algorithm of alternating projections is equivalent to steepest-descent optimization with a Gaussian likelihood, and is not suited for low-dose reconstructions52, because the for this case Poisson distribution arising from discretized count events differs too strongly from a Gaussian. While this argument does not hold for the more elaborate projection algorithms like DM and relaxed averaged alternating reflections (RAAR)58, they also fail in practice at low doses25, 59, and statistical reconstruction methods have to be used. Thibault and Guizar-Sicairos24 have analyzed maximum likelihood methods in conjunction with a conjugate gradient update rule as a refinement step, after the DM algorithm has converged. They demonstrate improved SNR by two orders of magnitude compared to the DM algorithm alone. They note, however, that starting directly with maximum likelihood optimization often poses convergence problems.

Due to the lack of algorithms with convergence guarantees, the mathematical community has recently picked up the problem, and a host of new algorithms with provable convergence has been developed. While we do not elaborate on them here we point the interested reader to the summary articles60, 61 and the article62, which refers to the most recent developments.

Here, we focus on developments which specifically target low-dose applications. Notable in this area is the work by Katkovnik et al.63, which in addition to the maximum likelihood estimate introduces a transform-domain sparsity constraint on the object and optimizes two objective functions in an alternating fashion: one for the maximizing the likelihood, and one for obtaining a sparse representation of the transmission function. However, instead of including the Poissonian likelihood directly, an observation filtering step is performed with a Gaussian likelihood. To obtain a sparse representation of the object, the popular BM3D denoising filter is used64.

During the writing of this paper, Yang et al. suggested using the Wigner Distribution Deconvolution technique for low-dose ptychography65, however no statistical treatment of the measurement process is included so far.

In this work, we formulate ptychographic phase retrieval as a Bayesian inference problem, by rewriting the probability of the transmission function T V given a set of measurements y=(y1,y2,…,yKM)T∈ℝ+KM according to Bayes’ rule:

| P(TV|y)=P(y|TV)P(TV)P(y). | 9 | | --------------------------- | - |

Since the measurements y i follow the Poissonian distribution

the likelihood function is given by

P(y|TV)=∏i=0KMIi(TV)yiyi!e−Ii(TV). 11

The prior distribution P(T V) is usually chosen such that it favors realistic solutions, so that noise is suppressed in the reconstructed image. Here we evaluate two different models. A simple prior, suggested in ref. 42, penalizes large gradients in the image with a Gaussian distribution on the gradient of the transmission function, which is also known as Tikhonov regularization:

PTikhonov(T)=exp(−μ0κ‖∇T(r⃗)‖2)=exp(−μ0κ∑i=1N(DxTV)i2+(DyTV)i2) 12

with κ=8Npix2Nm‖I‖1 chosen as in ref. 42. N m is the total number of valid measurements, N pix * K in the case when the detector has no hot pixels. This scales the numerical value of the prior to be close to the likelihood, such that the weight μ 0 can take values between 1 × 10−1 and 1 × 10−2. D_x_ and D_y_ are the discrete forward difference operators. The second prior we evaluate is based on the work by Katkovnik et al.63 and uses sparse modeling to denoise the transmission function:

Psparse(TV)=exp(−μ‖TV−TsparseV‖2) 13

Here, TsparseV is built up by applying the BM3D collaborative filtering algorithm64, 66, which we briefly describe in the Supplementary Material. As input for the BM3D algorithm we transform T V into hue-saturation-value format using domain coloring. The prior P sparse(T) reduces the difference between the denoised version of the current transmission function and the transmission function itself. We do not take into account the marginal likelihood P(y) due to the high dimensionality of the problem. Given the likelihood function P(y|T V) and the prior distribution P(T V), we can now write the objective function for the maximum-a-posteriori (MAP) estimate:

TMAPV:=argminTVℒMAP(TV)=argminTV(−log(P(y|TV)P(TV)P(y))). 14

The log-likelihood is given as

| ℒ(TV)=∑i=1KM[|piTV|2−yilog(|piTV|2)], | 15 | | ------------------------------------------ | -- |

and the MAP objective functions are

ℒTikhonov−MAP(TV)=ℒ(TV)+μ0κ‖∇T(r⃗)‖2 16

and

ℒBM3D−MAP(TV)=ℒ(TV)+μ1‖TV−TsparseV‖2, 17

for the two prior models, respectively. We calculate the gradients of both expressions:

| ∇ℒTikhonov−MAP(T)=∑i=1KM2piTV(1−yi|piTV|2)pi†+2μ0κ∑i=1N(DxTV)i+(DyTV)i, | 18 | | --------------------------------------------------------------------------- | -- |

| ∇ℒBM3D−MAP(T)=∑i=1KM2piTV(1−yi|piTV|2)pi†+μ1(TV−TsparseV) | 19 | | ------------------------------------------------------------ | -- |

Since equations (16) and (17) are non-convex functions, there is no guarantee that standard gradient descent converges to a global minimum. Recently, a non-convex algorithm for the generalized phase retrieval problem with Poisson noise was presented67, that provably converges to a global minimum with suitable initialization. It introduces a iteration-dependent regularization on the gradients of the likelihood to remove terms which have a negative effect on the search direction. Namely, it introduces a truncation criterion

| ℰi(TV)={|yi−|piTV|2|≤αhKM‖y−I‖1|piTV|‖TV‖2}, | 20 | | ---------------------------------------------- | -- |

that acts on the gradient of the likelihood and suppresses the influence of measurements that are too incompatible with the reconstruction. The truncation parameter α h ≥ 5 is described in ref. 67. The regularized likelihood gradient is then

| ∇ℒℰi(TV)=∑i∈ℰi(TV)KM[|piTV|2−yilog(|piTV|2)]. | 21 | | ------------------------------------------------- | -- |

We compute the next step using conjugate gradient descent68, 69, since this lead to much faster convergence compared to the update procedure described in ref. 67.

Initialization

Truncated spectral initialization for ptychography was first proposed by Marchesini et al.52, based on the notion that the highest intensities in the diffraction pattern carry the strongest phase information. They compute the phase of the largest eigenvector of the following hermitian operator:

| 1|yi|>εFQ(Q†Q)−1Q†F†1|yi|>ε, | 22 | | ------------------------------- | -- |

where ε is chosen such that the largest 20 percent of the intensities are allowed to contribute and F and Q are defined as above. The largest eigenvalue of a sparse hermitian matrix can be efficiently computed either with power iterations70, or with the Arnoldi method71. In ref. 67, truncated spectral initialization with a truncation rule with 1|yi|<α02λ02 is used, with λ0=∑i=1KMyi and α 0 a free parameter. It is also important to note that the truncated spectral initialization only produces visually correct initial phase to a dose of roughly 100 e −/Å2. Figure 7(b) shows an example initialization for a dose of 100 e −/Å2. For doses below this value, we initialized the transmission function with unity transmission and normal-distributed phase with mean 0.1 and variance of 0.1. With this random initialization we found no problem of convergence for all algorithms tested in this paper.

Figure 7.

Figure 7

(a) Convergence behavior of different gradient update rules. The normalized root mean square error (NRMSE) is defined in the Supplementary material. MAP refers to a constant prior. (b) Example for the transmission function initialization T 0 after 70 power iterations, for an electron dose of 100 e −/Å2, intensities were truncated at the 80th percentile. (c) TsparseV for human ribosome after 60 iterations of BM3D-MAP. Scale bar is 10 nm.

Reconstruction parameters

All ptychography reconstructions were performed with a probe area overlap of 75% in real-space, where the probe area is defined by all pixels contributing more than 1% of the maximum intensity. This corresponds to a step size of roughly 3 nm, depending on the probe used. At a dose of 20 e −/Å2, this corresponds to 8540 electrons per diffraction pattern, i.e. an average electron count of 0.52 electrons per detector pixel, and 58 electrons per image pixel for a pixel size of 1.7 Å. For the regularization parameters we performed a grid search evaluating the final normalized root mean square error and found the best values to be μ 0 = 1 × 10−2, μ 1 = 8 × 10−2. We choose the biorthogonal spline wavelet transform as the linear transform for BM3D as it achieves the best PSNR for high noise72. For the figures shown in the paper we use the BM3D denoiser. A typical intermediate denoised image for ribosome after 60 iterations can be seen in Fig. 7(c).

Implementation Details

The algorithms presented in this paper were implemented with the Torch scientific computing framework73. The gradient update routines were adapted from the optim package for Torch69. For efficient computing on the graphics processing unit (GPU) with complex numbers, the zcutorch library for CUDA was developed74. Hyperparameter optimization was done with the hypero75 package for Torch. For BM3D denoising we use the C++ implementation76. The code was run on an Intel i7-6700 processor with 32 GB RAM and a NVidia Titan X GPU with 12 GB RAM. The run time for optimization with ℒMAP was 26 s, and for optimization with ℒBM3D−MAP 35 s. This is expected, because the BM3D algorithm used here is not yet implemented on the GPU, and the BM3D denoising is computationally more intensive.

Data Availability

The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.

Conclusion

In this paper we have, via numerical experiments, demonstrated the feasibility of retrieving high-resolution electron transmission phase information of biological macromolecules using ptychography and Bayesian optimization. With the methods presented in this paper, it should be possible to achieve a resolution better than 1 nm for true single-particle imaging of molecular complexes with molecular weights ranging from below 100 kDa to a few MDa, and a resolution around 3 Å with simple averaging of 30 datasets. We have given a detailed explanation of the optimization and initialization procedures used, and have emphasized the importance of choosing an appropriate illumination function. We note that, while the high data redundancy in a ptychographic dataset empirically makes it experimentally very robust, there is much room for improvement in terms of measurement complexity. For the results presented here, the measurement dimension KM is larger than the problem dimension N by a factor of at least 30, while the theoretical limit for successful phase retrieval is at KM = 4_N_ 77. By reducing the number of measurements, the SNR of each individual diffraction measurement could be increased, yielding an improved image SNR in the reconstruction. Therefore, the development of an optimized experimental scheme, including design of the illumination function and scanning scheme is a promising direction of research and may enable significant improvement to the results presented here.

We would like to point out two obstacles that one may have to overcome in the experimental implementation of our method. Firstly, the best results are to be expected when recording zero-loss diffraction patterns with the use of an energy filter. The energy filter may introduce phase distortions into the diffraction patterns, which may need to be accounted for in the reconstruction algorithm. Secondly, although beam-induced movements are expected to be reduced by a large amount due to spot-scanning, the remaining movement may cause problems in the reconstruction. Statistically stationary sample movements can be accounted for in the reconstruction algorithm25, 78, but beam-induced motions are likely to be non-stationary, and dedicated algorithms may need to be developed to account for it.

Cryogenic ptychographic imaging of biological samples is also being developed in the X-ray sciences79, and our results could equally be implemented there to improve the dose-effectiveness. Finally, the methods presented here may find application in electron phase imaging of radiation-sensitive samples under non-cryogenic conditions, and the incorporation of Bayesian methods into in-focus ptychographic reconstruction procedures18, 65, may provide similar gains in SNR as the ones discussed here, while also keeping the analytical capabilities of traditional scanning TEM imaging.

Electronic supplementary material

Acknowledgements

This work was funded by the Max Planck Society through institutional support. P.M.P. acknowledges support from the International Max Planck Research School for Ultrafast Imaging & Structural Dynamics.

Author Contributions

P.M.P., R.J.D.M. and R.B. conceived the idea. P.M.P. and W.Q. wrote the simulation code. P.M.P. wrote the ptychographic reconstruction code. P.M.P. wrote the paper with input and discussions from R.B., G.K., and R.J.D.M.

Competing Interests

The authors declare that they have no competing interests.

Footnotes

Electronic supplementary material

Supplementary information accompanies this paper at doi:10.1038/s41598-017-07488-y

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

References

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

Data Availability Statement

The datasets generated and analysed during the current study are available from the corresponding author on reasonable request.