Open Access
Add to BibSonomyAbstract
In this paper, we propose the SPR (sparse phase retrieval) method, which is a new phase retrieval method for coherent x-ray diffraction imaging (CXDI). Conventional phase retrieval methods effectively solve the problem for high signal-to-noise ratio measurements, but would not be sufficient for single biomolecular imaging which is expected to be realized with femto-second x-ray free electron laser pulses. The SPR method is based on the Bayesian statistics. It does not need to set the object boundary constraint that is required by the commonly used hybrid input-output (HIO) method, instead a prior distribution is defined with an exponential distribution and used for the estimation. Simulation results demonstrate that the proposed method reconstructs the electron density under a noisy condition even some central pixels are masked.
© 2011 Optical Society of America
1. Introduction
X-ray free electron lasers (XFELs) can potentially provide us a novel mean to determine the three-dimensional (3D) structure of biomolecules from the diffraction images of single molecules [1,2]. Conventionally, x-ray crystallography has been the principal tool to determine high-resolution 3D structures of proteins, nucleic acids, and their complexes. However, a critical process is the crystallization, where a sufficiently large crystal must be prepared. It is known that about 40% of biomolecules, particularly membrane proteins, do not crystallize. The coherent x-ray diffraction imaging (CXDI) by XFELs does not require a crystal and will change the situation.
In order to realize the single molecule imaging, a straightforward scenario has been proposed [3, 4] and demonstrated with single mimivirus particles [5]. Firstly, a series of single molecules with the same structure are injected into vacuum and exposed to a very short (less than 5 fs) pulse x-ray beam. The diffraction image of the molecule with unknown orientation is recorded before radiation-induced structural degradation. This process is repeated until a sufficient number of images are recorded. The diffraction data are then classified into groups according to the orientations of the molecule and averaged to improve the signal-to-noise ratio [6,7]. Next, the relative orientations of the classes are determined to assemble a 3D coherent diffraction pattern in reciprocal space. Finally, the 3D electron density of the molecule is obtained by a phase retrieval method.
The above scenario was proved to be effective for the 3D reconstruction of mimivirus particles [5]. However, the condition for a single molecule CXDI will be more severe since the number of the photons received by each pixel will be very limited. Some methods of reconstructing 3D diffraction intensity have been proposed for this problem [8–10]. They are likely to be computationally intensive and the feasibility remains to be seen.
Another scenario might be possible. After collecting a sufficient number of diffraction images with unknown orientations, 2D electron density images are first reconstructed from the 2D diffraction data, followed by the 3D reconstruction method used in transmission electron microscopy [11]. In either scenario, the difficulty is due to the small number of photons scattered by a single molecule even though x-ray beams that linac coherent light source (LCLS) has started to commission give us illumination which is around 109 times brighter (∼1019 photons/pulse/mm2 in the unfocused beam) than that of current synchrotrons [12].
We propose a new method for retrieving phases of a diffraction image which will be obtained by XFELs. The method is based on the maximum a posteriori (MAP) estimation of the Bayesian statistics. The posterior distribution is computed from the likelihood and the prior. For the present problem, the Poisson distribution is appropriate for the likelihood. The prior reflects the assumption of the electron density. Since the electron density is known to have a lot of zeros, we employ the exponential distribution which encourages the MAP estimate to have a lot of zeros (this is called a sparse solution [13, 14]. A similar idea is employed in [15]). Therefore, we call our method, sparse phase retrieval (SPR) method.
This combination is effective for the problem and is distinct from other Bayesian approaches proposed for different phase retrieval problems [16, 17]. Compared to the widely used hybrid input-output (HIO) method [18–20], the SPR method does not require the support region and is robust against the limited photon counts and the missing central pixels. In this paper, we focus on 2D phase retrieval problem. The resulting SPR method can be used directly for the second scenario. If 3D diffraction data are available, it is also possible to apply the SPR method. But reconstruction of 3D diffraction data from 2D images with unknown orientations is beyond the scope of the present paper.
2. Proposed method
2.1. Diffraction image and noise
We first formulate the phase retrieval problem of a single diffraction image. Let f (x,y) ≥ 0 be the electron density of a molecule projected onto a two-dimensional (2D) plane. The coordinate is discretized, that is, x,y = 1, ⋯ , M. We introduce the following notation
where, F is the Fourier transform of f defined as follows, The coefficient of the Fourier transform differs from the commonly used 1/M2. We prefer the above definition because the power of F and f is the same. We assume that the frame size is sufficiently larger than the molecule in the image and many components of f are zero.In the ideal situation, each pixel of a diffraction image receives the number of photons proportional to the power spectrum |Fuv|2, where (u,v) is the 2D index of the pixel. However, the XFEL measurement of a single molecule diffraction will be different because the flux of the x-ray laser is not sufficiently large. The number of the photons detected by each sensor will be small and behave stochastically.
Let Nuv be the number of photons detected by the sensor at (u,v). We study the distribution of Nuv. The total number Nall = Σuv Nuv is a stochastic variable which follows the Poisson distribution. The expected value of Nall is proportional to the intensity IX of x-ray. Assuming each Nuv is independent and the distribution belongs to the same family, Nuv also follows a Poisson distribution. Let the expected value of Nuv be Suv and the distribution of Nuv becomes
We only consider this counting process and do not consider other observation noise. Approximately, Suv is denoted as Suv = α|Fuv|2 cos3θ, where α is a positive constant and θ is the scattering angle (see Fig. 1). Because θ is a function of u and v, we further rewrite it as α|Fuv|2cuv, where cuv = c(u,v) = cos3θ.The Fourier transform, F, is a linear function of f, and the probability of N, where N = {Nuv} = {N11, ⋯ , NMM}, is denoted as follows
Note that the scaling of |Fuv|2 cannot be determined only from the observed scattering pattern and α is set to 1. This scaling indeterminacy is ignored here. It would be recovered from other knowledge, such as the total energy.In the following, we consider the density reconstruction as an estimation of f. Commonly used estimate is the maximum likelihood estimate (MLE) which maximizes p(N|f). The likelihood in Eq. (2) is maximized by setting |Fuv|2 = Nuv/cuv, but this is not sufficient for the estimation of f because the phase is lost. Some additional information is needed.
2.2. SPR method
Instead of the MLE, the MAP estimate of the Bayesian statistics is used for the SPR method. The MAP estimate maximizes the posterior distribution, where the posterior distribution is defined as the product of the likelihood function and the prior distribution. Similar framework has been proposed for different phase retrieval problems [16, 17].
The likelihood function is defined as a Poisson distribution as in Eq. (2) and we discuss the prior distribution in this subsection. For a 2D image of a single molecule diffraction, it is natural to assume that the electron density is high in the central part and becomes zero in the peripheral part. But we do not know which pixel is zero. The prior should reflect this knowledge.
It is known that the prior distribution with an exponential distribution encourages the MAP estimate to have a lot of zeros [13, 14] and we use this for the SPR method
where each fxy is assumed to be independent. The hyperparameter ρxy reflects the prior belief of fxy = 0. Larger the value of ρxy, stronger the belief. In order to reflect the prior knowledge, we define the following one parameter function ρ(μ)xy The two parameters a and b were adjusted to make wxy equal to 0 at the center and 1 at the corners (see Fig. 2) and μ adjusts the overall strength of the prior. We have employed a parabolic function, but if there is any prior knowledge, ρ(μ)xy should reflect it.The posterior distribution of f observing N has the following relation,
It is convenient to take the logarithm of p(f|N) in order to compute the MAP estimate. By taking the logarithm and collecting terms related to f, the following function is obtained The first summation corresponds to the likelihood with the Poisson distribution (similar idea is used in [16]), the second summation represents the prior with the exponential distribution, and μ controls the balance between two terms. This combination is the key idea of the proposed method. The SPR method is summarized as follows. We start with a large μ and decrease it. The MAP estimate is computed for each μ, and the best value of μ is selected according to some criterion. This issue will be discussed later.We note that our formulation is related to the compressed sensing (CS) [21], and the idea of CS has been applied to the phase retrieval problem [15]. Our approach is different from the work since the Poisson likelihood is included.
2.3. Algorithm
A gradient based algorithm is used to compute the MAP estimate. The algorithm is shown below starting with the following relation,
By defining the inverse Fourier transform as the derivative of ℓμ(f|N) becomes where g(F;N) = {guv(Fuv; Nuv)} and In order to compute the derivative of ℓμ(f|N), ℱ is applied to f, then ℱ−1 is applied to g(F; N). Thus, the computational cost is similar to the HIO method [19].We employed a naïve iterative updating rule for the MAP estimate of fxy, x, y = 1, ⋯ , M.
where t is the iteration number. The “max” operation keeps fxy nonnegative. The positive variable ηt controls the step size and a simple line search of ηt accelerated the convergence largely.The initial value is generally important for gradient based algorithms. But for the SPR method, the influence is effectively reduced by controlling μ. The MAP estimate is first computed for a large μ and then for smaller μ’s successively. For a very large μ, most components of the optimal f are 0 and the influence of the initial value is negligible. The MAP estimate is then used for the initial value of a slightly smaller μ, which is a good starting point in general.
2.4. HIO method and Bayesian framework
In this subsection, we discuss the relation between the widely used HIO method and the SPR method.
The HIO method [18–20] exploits the knowledge that the positive density is localized in a small support region. Let us define the support region γ, and the method is written as follows,
This problem can be formulated as the MAP estimation problem by setting the Gaussian likelihood function and the uniform prior distribution, that is, where C is a positive number. If fxy is unbounded, above p̃u(fxy) is improper. The MAP estimate is computed by solving Eq. (7).Formally, the difference between the HIO method and the SPR method is the likelihood and the prior. But more importantly, the proposed SPR method does not require the support region. In the SPR method, a lot of components of the estimated electron density f become 0 automatically.
3. Numerical Experiment
3.1. Simulated data
The SPR method was tested with simulated diffraction data. The electron density of lysozyme, a 14 kD protein, was projected onto a 2D plane. The projected density was then Fourier transformed into reciprocal space (u,v) to simulate a diffraction image. The density and the ideal diffraction image are shown in Figs. 3(a) and 3(b), respectively.
Fig. 3 Simulated 2D electron density (a), the diffraction image without noise (b), and simulated diffraction patterns (c,d) of a lysozyme protein are shown. Each image is 308×308 pixels. The spatial frequency of the diffraction patterns is 0.0812 nm−1 and the unit of the x-ray fluxes I is photons/pulse/mm2.
Experimental diffraction images were simulated by converting the diffraction intensity into the number of photons per effective pixel (λ/σL)2 according to the Poisson distribution with the mean, Suv, defined as follows,
Here, I is the incident x-ray flux (photons/pulse/mm2), rc is the classical electron radius (2.82×10−12 mm), λ is the x-ray wave length (0.1 nm), σ is the oversampling ratio per one dimension (2), and L is the molecule diameter (6.16 nm) (see Fig. 1). The Ewald sphere curvature was ignored but this is a good approximation for small angles.Two examples for scattered photons of x-ray fluxes with I = 5.0×1021 and 1.0×1021 photons/pulse/mm2 are shown in Figs. 3(c) and 3(d), respectively. Although the fluxes in the above simulation are around 100 times larger than that of the current LCLS, the sizes of the molecules of interest are generally larger than lysozyme. When the size is around 1000 kD, the total number of photons arriving at the 2D plane will be in the similar range as Figs. 3(c) and 3(d).
3.2. Results
The SPR method was applied to the data to reconstruct the density. Figure 4 shows the results of the SPR method applied to the ideal diffraction image. As μ decreases, more components become positive. The reconstructed electron density becomes similar to the true density when μ = 1. Thus, it is important to select the best μ. A simple strategy is to rank the results with a criterion. Many criteria have been proposed in statistics, such as, Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) [22], but here we use a different one. Let f̂(μ) be the maximizer of ℓμ(f|N) and the following error is used as the criterion
where F̂(μ)uv = ℱ(f̂(μ))uv. This error function has been widely used for the phase retrieval problem and is employed here as well in order to rank the results. Roughly speaking, the error decreases as μ shrinks since it controls the balance between the prior and the likelihood. But the error does not necessarily decrease monotonically because the log posterior in Eq. (4) is different from the error. Thus, we choose the μ which minimizes ErrorF(μ).Fig. 4 Reconstruction of the electron density from the ideal diffraction image by the SPR method with different values of μ.
Figure 5 shows the ErrorF (μ) for diffraction images in Fig. 3. The value of μ was varied from 1×104 to 1×10−4. Note that although the horizontal axis of Fig. 5 increases from left to right, the value of μ was decreased in the experiments. The best μ was selected for each diffraction image by ErrorF(μ) and the reconstructed electron densities are shown in Figs. 6(a), 6(b), and 6(c).
Fig. 5 The values of ErrorF (μ) for simulated diffraction images in Fig. 3. The black circle and triangles on the curves correspond to the reconstructed density images of Figs. 6.
We also applied the HIO method for the data. The central square (154×154 pixels) was set as the support region. The method has a problem with convergence, and we stopped the iteration loop when the incremental difference of the estimated density becomes small enough. Figure 7 shows the reconstructions obtained by the HIO method for the diffraction images in Fig. 3. The reconstructed electron densities are similar to Fig. 6, but the SPR method obtained better ErrorF (μ) for all of the diffraction images.
We note the computational time of the SPR method. The method was implemented with C and run on a desktop machine with an Intel Xeon X5570 (2.93 GHz 4 cores). The optimal electron density was computed for all the values of μ’s (33 different values of μ) within 10 minutes.
3.3. Missing centers
One of the problems of the phase retrieval from diffraction images is missing centers. The intensities around the center cannot be measured experimentally because the scattered beam is superimposed by the direct beam around the center and a beam stop is used to block the direct beam.
The SPR method can be applied to the case easily by modifying Eq. (4) as follows,
where α is the active pixel set where the central pixels are excluded. The error function is also modified asFirstly, we masked the central pixel from each diffraction image in Fig. 3 and applied the SPR method. By masking the central pixel, Fig. 3(b) lost around 5.7% of the photons, Fig. 3(c) lost 294 out of 4630 photons and Fig. 3(d) lost 56 out of 957 photons. Figure 8 shows the errors. Note that the error values are better than those in Fig. 5, but it does not mean the reconstructed densities in Fig. 9 are better than those in Fig. 6, because the errors are computed only for the active pixel set. Figure 10 shows the results of the HIO method. Since the central pixel corresponds to a scaling factor of the electron density of the protein, the reconstructed images are quite similar to those reconstructed from the complete diffraction images. We see that the error values of the SPR method are better than those of the HIO method.
Fig. 8 The values of ErrorF (μ) for simulated diffraction images in Fig. 3 without the central pixel. The black circle and triangles on the curves correspond to the reconstructed density images of Fig. 9.
Next, we further masked the central 3 × 3 pixels from each diffraction image and applied the SPR method. Note that Fig. 3(b) lost around 35% of the photons, Fig. 3(c) lost 1595 out of 4630 photons and Fig. 3(d) lost 312 out of 957 photons. More than 30% of the photons are lost and the results are expected to be degraded largely. Figure 11 shows the errors. The reconstructed densities are shown in Fig. 12. Figure 13 shows the results of HIO method. We see that the error values of the SPR method are better and the reconstructed results are more stable than the HIO results. It is worth noting that the SPR method worked well even if the central 5 × 5 pixels, a half area of the centrospeckle pattern, were masked. In that condition, 54.4 % of the total photons are lost.
Fig. 11 The values of ErrorF (μ) for simulated diffraction images in Fig. 3 without the central 3 × 3 pixels. The black circle and triangles on the curves correspond to the reconstructed density images of Figs. 12.
4. Discussion
We have proposed a new phase retrieval method. The SPR method is based on the MAP estimation of the Bayesian statistics, it will be effective for single molecule diffraction images which will be obtained by XFELs.
From a Bayesian viewpoint, the difference between other phase retrieval methods including the HIO method and the SPR method is summarized in the likelihood and the prior. Many methods use the squared loss, which corresponds to a Gaussian likelihood, while the present method employs a Poisson distribution. This is suitable for CXDI data because the number of the scattered photons detected at each pixel is a counting process. Similar idea has been found in a related work [16]. The prior in HIO method corresponds to a uniform distribution with a bounded support, while the proposed method uses the exponential distribution. Similar prior has been widely used in statistics [13]. This combination shows a new promising direction for phase retrieval.
The proposed method has been tested with simulated data. The electron densities are reconstructed with a reasonable computational cost, and errors of the results are smaller than the HIO method. We further tested the method with the diffraction data where some central pixels are masked. The proposed method is applied without any major modification, and the results are quite stable even some pixels are masked.
The SPR method reconstructs the density from an incomplete observation utilizing the sparsity. This idea is similar to CS [21]. In [15], CS is applied to the phase retrieval problem, however, in our problem, the number counts are very limited, and we needed to extend the model.
Finally, we list some of our future works to improve the SPR method. One is the hyperparameter. The hyperparameter ρ(μ)xy can be modified to reflect the knowledge of the true density. This may result in a better estimate. Another issue is the algorithm. The densities for the whole sequence of μ in Fig. 5 are computed with a simple line search algorithm but advanced optimization methods may speed it up.
Acknowledgments
We would like to thank Nobuhiro Go for fruitful discussion and Atsushi Tokuhisa for providing the electron density data. We are also grateful to Satoshi Ito for helpful discussion on the optimization method. This work was partly supported by X-ray Free Electron Laser Utilization Research Project of the Ministry of Education, Culture, Sports, Science and Technology of Japan.
References and links
1. D. Sayre, “Prospects for long-wavelength X-ray microscopy and diffraction,” in “Imaging Processes and Coherence in Physics,”, vol. 112 of Springer Lecture Notes in Physics, M. Schlenker, M. Fink, J. P. Goedgebuer, C. Malgrange, J. C. Viénot, and R. H. Wade, eds. (Springer, 1980), pp. 229–235. [CrossRef]
2. K. Gaffney and H. Chapman, “Imaging atomic structure and dynamics with ultrafast x-ray scattering,” Science 316, 1444–1448 (2007). [CrossRef] [PubMed]
3. R. Neutze, R. Wouts, D. van der Spoel, E. Weckert, and J. Hajdu, “Potential for biomolecular imaging with femtosecond X-ray pulses,” Nature (London) 406, 752–757 (2000). [CrossRef]
4. G. Huldt, A. Szőke, and J. Hajdu, “Diffraction imaging of single particles and biomolecules,” J. Struct. Biol. 144, 219–227 (2003). [CrossRef] [PubMed]
5. M. M. Seibert, T. Ekeberg, F. R. N. C. Maia, M. Svenda, J. Andreasson, O. Jönsson, D. Odić, B. Iwan, A. Rocker, D. Westphal, M. Hantke, D. P. DePonte, A. Barty, J. Schulz, L. Gumprecht, N. Coppola, A. Aquila, M. Liang, T. A. White, A. Martin, C. Caleman, S. Stern, C. Abergel, V. Seltzer, J.-M. Claverie, C. Bostedt, J. D. Bozek, S. Boutet, A. A. Miahnahri, M. Messerschmidt, J. Krzywinski, G. Williams, K. O. Hodgson, M. J. Bogan, C. Y. Hampton, R. G. Sierra, D. Starodub, I. Andersson, S. Bajt, M. Barthelmess, J. C. H. Spence, P. Fromme, U. Weierstall, R. Kirian, M. Hunter, R. B. Doak, S. Marchesini, S. P. Hau-Riege, M. Frank, R. L. Shoeman, L. Lomb, S. W. Epp, R. Hartmann, D. Rolles, A. Rudenko, C. Schmidt, L. Foucar, N. Kimmel, P. Holl, B. Rudek, B. Erk, A. Hömke, C. Reich, D. Pietschner, G. Weidenspointner, L. Strüder, G. Hauser, H. Gorke, J. Ullrich, I. Schlichting, S. Herrmann, G. Schaller, F. Schopper, H. Soltau, K.-U. Kühnel, R. Andritschke, C.-D. Schröter, F. Krasniqi, M. Bott, S. Schorb, D. Rupp, M. Adolph, T. Gorkhover, H. Hirsemann, G. Potdevin, H. Graafsma, B. Nilsson, H. N. Chapman, and J. Hajdu, “Single mimivirus particles intercepted and imaged with an X-ray laser,” Nature (London) 470, 78–81 (2011). [CrossRef]
6. G. Bortel and G. Faigel, “Classification of continuous diffraction patterns: A numerical study,” J. Struct. Biol. 158, 10–18 (2007). [CrossRef]
7. G. Bortel, G. Faigel, and M. Tegze, “Classification and averaging of random orientation single macromolecular diffraction patterns at atomic resolution,” J. Struct. Biol. 166, 226–233 (2009). [CrossRef] [PubMed]
8. R. Fung, V. Shneerson, D. K. Saldin, and A. Ourmazd, “Structure from fleeting illumination of faint spinning objects in flight,” Nat. Phys. 5, 64–67 (2009). [CrossRef]
9. N.-T. D. Loh and V. Elser, “Reconstruction algorithm for single-particle diffraction imaging experiments,” Phys. Rev. E 80, 026705 (2009). [CrossRef]
10. D. K. Saldin, V. L. Shneerson, R. Fung, and A. Ourmazd, “Structure of isolated biomolecules obtained from ultrashort x-ray pulses: exploiting the symmetry of random orientations,” J. Phys.: Condens. Matter 21, 134014 (2009). [CrossRef]
11. J. Frank, Three-Dimensional Electron Microscopy Of Macromolecular Assemblies: Visualization Of Biological Molecules In Their Native State (Oxford University Press, 2006). [CrossRef] [PubMed]
12. L. Young, E. P. Kanter, B. Krässig, Y. Li, A. M. March, S. T. Pratt, R. Santra, S. H. Southworth, N. Rohringer, L. F. DiMauro, G. Doumy, C. A. Roedig, N. Berrah, L. Fang, M. Hoener, P. H. Bucksbaum, J. P. Cryan, S. Ghimire, J. M. Glownia, D. A. Reis, J. D. Bozek, C. Bostedt, and M. Messerschmidt, “Femtosecond electronic response of atoms to ultra-intense X-rays,” Nature (London) 466, 56–61 (2010). [CrossRef]
13. R. Tibshirani, “Regression shrinkage and selection via the lasso,” J. R. Stat. Soc. Ser. B 58, 267–288 (1996).
14. P. Zhao and B. Yu, “On model selection consistency of lasso,” J. Machine Learning Res. 7, 2541–2563 (2006).
15. S. Marchesini, “Ab initio compressive phase retrieval,” ArXiv Physics Optics e-prints, arXiv:0809.2006, (2008).
16. R. Irwan and R. G. Lane, “Phase retrieval with prior information,” J. Opt. Soc. Am. A 15, 2302–2311 (1998). [CrossRef]
17. S. Baskaran and R. P. Millane, “Bayesian image reconstruction from partial image and aliased spectral intensity data,” IEEE Trans. Image Process. 8, 1420–1434 (1999). [CrossRef]
18. J. Fienup, “Reconstruction of an object from the modulus of its Fourier transform,” Opt. Lett. 3, 27–29 (1978). [CrossRef] [PubMed]
19. J. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21, 2758–2769 (1982). [CrossRef] [PubMed]
20. J. Miao, D. Sayre, and H. Chapman, “Phase retrieval from the magnitude of the Fourier transforms of nonperiodic objects,” J. Opt. Soc. Am. A 15, 1662–1669 (1998). [CrossRef]
21. D. Donoho, “Compressed sensing,” IEEE Trans. Inform. Theory 52, 1289–1306 (2006). [CrossRef]
22. H. Zou, T. Hastie, and R. Tibshirani, “On the “degrees of freedom” of the lasso,” Ann. Stat. 35, 2173–2192 (2007). [CrossRef]
