Experimental 3D coherent diffractive imaging from photon-sparse random projections

The application of X-ray coherent diffractive imaging to structure determination of single biomolecules at free electron lasers yields signal levels down to a few hundred scattered photons per frame or less. Here we demonstrate in an analog experiment at a synchrotron source that even at a sparsity level of order 1.3 × 10−3 photons per pixel per frame the 3D structure of a nano-fabricated sample can be obtained, even if the orientation of the particle for a given individual frame is not explicitly known.

The routine atomic resolution structure determination of single particles is expected to have profound implications for probing structure-function relationships in systems ranging from energy-storage materials to biological molecules. Extremely bright ultrashort-pulse X-ray sources -X-ray freeelectron lasers (XFELs) -provide X-rays that can be used to probe ensembles of nearly identical nanoscale particles. When combined with coherent diffractive imaging, these objects can be imaged; however, as the resolution of the images approaches the atomic scale, the measured data are increasingly difficult to obtain and, during an X-ray pulse, the number of photons incident on the 2D detector is much smaller than the number of pixels. This latter concern, the signal 'sparsity', materially impedes the application of the method. An experimental analog using a conventional X-ray source is demonstrated and yields signal levels comparable with those expected from single biomolecules illuminated by focused XFEL pulses. The analog experiment provides an invaluable cross check on the fidelity of the reconstructed data that is not available during XFEL experiments. Using these experimental data, it is established that a sparsity of order 1.3 Â 10 À3 photons per pixel per frame can be overcome, lending vital insight to the solution of the atomic resolution XFEL single-particle imaging problem by experimentally demonstrating 3D coherent diffractive imaging from photon-sparse random projections. before destruction' is routinely used for serial femtosecond crystallography, in which a randomly oriented crystal is illuminated by the X-ray beam (Chapman et al., 2011;Schlichting, 2015). The crystal size has typically been large enough such that sufficient numbers of X-rays are diffracted into a recordable pattern prior to destruction to be able to determine the orientation of the crystal. This is repeated with new microcrystals until there is an adequate sampling of reciprocal space to determine the 3D structure of the unit cell of the crystal.
To obtain 3D structural information from single particles, such as virus particles or single proteins, serial diffraction data from many identical or nearly identical objects has to be measured with sufficient orientational variation. Usually, this is achieved by randomly injecting particles into the FEL beam, relying on the statistical coincidence of a single particle being hit by an FEL pulse (Barty, 2016;Spence, 2017). Due to the random nature of this process, the sample's orientation for each diffraction pattern is generally unknown and has to be recovered a posteriori in order to build up a continuous 3D diffraction volume in reciprocal space. This can then be inverted by iterative phase retrieval (Marchesini, 2007) into a real-space electron-density distribution, the last step of singleparticle CDI.
A first complete demonstration of the method was provided by the 3D structure determination of the Giant Mimivirus, approximately 450 nm in diameter, to a resolution of 125 nm (Ekeberg et al., 2015). More recently, the same type of experiment has been successfully performed using the much smaller Melbourne virus (diameter around 230 nm) to a resolution of 28 nm (Lundholm et al., 2018). At the same time, a step towards even smaller viruses -the Rice Dwarf Virus (RDV) and bacteriophage PR772, both with a diameter of around 70 nm -has been made, resulting in images at a resolution of around 17 nm  and more recently, for PR772, slighly below 10 nm (Rose et al., 2018). 1 Those data were a result of the Single-Particle-Imaging Initiative at the Linac Coherent Light Source (Aquila et al., 2015;Munke et al., 2016;Reddy et al., 2017). This also resulted in the collection of a few hundred high-resolution diffraction frames from RDV at a photon energy of 7 keV, showing that a useful diffraction signal can be collected at 5.9 Å resolution from single hits (Munke et al., 2016).
Reconstructing biological macromolecules to 3 Å resolution or better has previously been set as the ultimate goal of single-particle CDI (Aquila et al., 2015). However, when the particle is sufficiently small, only a few hundred X-rays are scattered elastically prior to destruction of the particle Fortmann-Grote et al., 2017). In this case the recorded scattering pattern contains too few X-rays to be able to determine the orientation of the particle from the single pattern alone. Such patterns are called 'sparse', i.e. the features necessary to interpret the pattern are totally dominated by Poisson noise. For example, a protein of at least ca 10 nm diameter is necessary to scatter, on average, 50 photons outside the central speckle, at a photon energy of 8 keV in a nanoscale FEL focus (see below for further details). In this case, hundreds of thousands of diffraction patterns have to be collected to build up the 3D reciprocal-space intensity, i.e. to assemble an invertible data set (Loh & Elser, 2009). To date, no such experimental data set exists and considerable method development is still required towards the realization of singleparticle CDI as an independent method of macromolecular structure determination.
However, the important case of sparse diffraction from a 3D object has only been solved experimentally in a setting different from the classical CDI problem. For example, one of the methods of orientation recovery, a statistical technique based on expectation maximization, the expand-maximizecompress (EMC) algorithm (Loh & Elser, 2009;Ayyer et al., 2016), has been applied successfully to real-space sparse radiographic data in two (Philipp et al., 2012) and three dimensions (Ayyer et al., 2014), to sparse crystallographic data limited to one (Wierman et al., 2016) and two rotation axes (Lan et al., 2017), and very recently also to synchrotron-based serial protein crystallographic data for random crystal orientations (Lan et al., 2018).
Here, we demonstrate 3D CDI from sparse random projections in the same geometry as that used for FEL-based single-particle imaging experiments. Using a synchrotron beam on a micron-scale sample, we show that, with as few as about 50 scattered photons per diffraction pattern, and without explicit knowledge of the sample's orientation for a given data frame, it is possible to robustly reconstruct the scattering distribution of the sample in reciprocal space and to invert this diffraction volume into a high-resolution 3D electron density distribution. We show that, despite strong sparsity in the data, it is possible to reconstruct a sample to a complexity of more than 40 resolution elements within the largest diameter of the sample. For ease of comparison, the term 'resolution element' is defined by analogy with the work of Loh & Elser (2009)  of the detector, we define it here in terms of the resolution as obtained after phasing.

Experimental
The experiment was performed on the undulator beamline ID10, end station EH2, of the European Synchrotron Radiation Facility (ESRF) (Chushkin et al., 2014). The photon energy was set to 8.1 keV using a water-cooled Si(111) pseudochannel-cut monochromator with an intrinsic energy resolution of ÁE/E ' 1.4 Â 10 À4 . The sample was a solid gold object with a largest diagonal length of about 1.1 m, fabricated by electroplating and supported by a silicon nitride membrane. It was placed in the beam on a high-precision tomographic stage at a distance of 4.0 m from the detector. Defined by several sets of slits (Chushkin et al., 2014), the lateral beam size at the sample plane was approximately 10 Â 10 m. A schematic of the setup is shown in Fig. 1. For further details of all the major aspects of our data treatment and the algorithms employed in the reciprocal-and real-space reconstruction of our data, we refer the reader to the supporting information.
For data collection, the mixed-mode pixel array detector (MM-PAD) was used. It is a wide dynamic range integrating detector developed at Cornell University that is capable of collecting data at a kilohertz frame rate with a high signal-tonoise ratio (SNR), ranging from one X-ray photon per pixel to a maximum rate exceeding 10 8 photons per pixel per second at 8 keV photon energy (Tate et al., 2013). The beam was attenuated using a polished single-crystal Si attenuator with a transmission of ca 0.3. As a result, the overall flux reaching the detector was approximately 10 8 photons s À1 .
The data set analyzed here comprises diffraction patterns from 227 unique 3D orientations. With 2000 collected frames per orientation, this amounts to a total of M data = 454 000 data frames. Each frame was collected with an illumination time of 25 ms. The orientations were obtained from two independent tomographic series (see Fig. 1), during each of which only , the angle about the tomographic rotation axis y, was varied.

Detector calibration and further data treatment
As the collected data were very sparse, reliable identification of single-and few-photon events becomes crucial for the analysis which explicitly takes into account the Poissonian nature of the noise in the data (see below). Calibration of the detector data, i.e. the transformation of the raw detector A schematic diagram of the experiment. The sample, a gold nanostructure supported on a silicon nitride membrane, was rotated about the y axis by an angle to obtain diffraction patterns at different orientations with respect to the optical axis z. A first rotation series about the y axis was followed by an in-plane rotation of the sample about z (angle ) and a subsequent second rotation series about y. The beam attenuation and illumination time were adjusted so that each data frame contains only about 50 scattered photons. An example of a single diffraction pattern is shown in the inset on the upper left-hand side (see also animation in the supporting information). The area of a single pixel in comparison with the field-of-view area has been enlarged for visualization. output into photon counts, is therefore an essential step of the analysis.
The applied calibration procedure consisted of several steps. A binary mask was used to reject all pixels that were inactive or were within gaps between detector modules. Second, the dark signal was subtracted for each pixel. The gain of the detector was determined from the histogram of pixel values over many frames within a region of interest where the maximum intensity per pixel is only a few photons per frame. The resulting histogram shows discrete peaks corresponding to zero, one, two etc. photons per pixel. In terms of analog-todigital units (ADU), a gain of 11.1 ADU photon À1 at 8.1 keV was determined from the spacing of these peaks. This gain, together with the widths of these peaks, yields an SNR of 5.2 for single photons. For the transformation of ADUs into single photons, a threshold energy of E t = E À HWHM n = 0.77E was used, with the single-photon energy E and the noise peak half-width at half-maximum HWHM n .
With around 40 000 active detector pixels in a region of interest (ROI) of 255 Â 255 pixels around the beam center, this leads to a false-positive probability P(1|0) ' 3 Â 10 À5 or between 1 and 2 events per frame. At the sparsity level of the data and with more than 400 000 frames comprising the complete data set, it was found that cosmic radiation made up a considerable contribution to the total recorded intensity. We therefore applied a threshold-based removal procedure. The sparsity level of the single frames is reflected by an average of 49.3 scattered photons detected per frame (median 48 and standard deviation 9.6), after masking out pixels dominated by empty-beam scattering. Due to the low level of background scattering from the instrument, e.g. optical components or apertures, no background or 'empty-beam' subtraction was performed on the data.

Orientation determination
The goal of orientation determination is to obtain the 3D reciprocal-space intensity W(q) that is proportional to the modulus-squared Fourier transform of the 3D electron density of the sample. Here, q denotes the 3D Cartesian reciprocalspace coordinate. To reconstruct W(q), the EMC algorithm correlates each data frame, K d (d = 1, . . . , M data ), with tomographic slices W j (j = 1, . . . , M rot ) of W, corresponding to M rot possible sample orientations, based on the current iterate of our model of W(q). Each iteration comprises expanding the current model W into slices W j , an update W j ! W j 0 by maximizing a log-likelihood function Q(W 0 ), and compressing slices W j 0 into a new 3D model W 0 (q). The update itself consists of forming the weighted sum The index i specifies a pixel and P jd (W) is the probability of a frame K d having been collected at an orientation j, based on the slice W j . For very noisy diffraction data from weak scatterers this probability can be well described using Poisson statistics (Loh & Elser, 2009).
Each orientation j of the sample is represented by a unit quaternion, q j . For the present experiment, the set fq j g of possible orientations was derived from the experimental setup and procedure (for details, see supporting information). As a result, the optimized set fq j g of orientations could be used to generate a Fourier intensity W ref (q), assembled using full knowledge of the orientations, as a reference for the EMCbased intensity reconstruction. EMC received the same set of orientations as an input, together with all data frames, but without explicit knowledge of any frame's orientation. As a further input, a geometry file was included that contained the reciprocal-space coordinate of each detector pixel and a 3D binary mask, S, identifying those voxels in the cubic domain of W(q) that are to be excluded from the analysis process, e.g. because they are never reached by any Ewald sphere slice. As S is non-symmetric with respect to q = 0, the Friedel symmetrization step included in EMC (Loh & Elser, 2009) was modified accordingly. To account for the fact that some pixels in each frame K d contain a large fraction of sample scattering, but still have a significant contribution of 'parasitic' or beamline scatter, we defined a binary mask on the detector ROI to identify those pixels to be included in the update rule W ! W 0 , but not into the calculation of probabilities P jd (W) .
As EMC was always initiated with a random intensity distribution, independent runs of the algorithm show some statistical variation. To reduce the associated uncertainty, the EMC algorithm was run 20 times for 500 iterations, followed by an averaging procedure similar to that described by . A small fraction of the ensemble, 2 out of 20 reconstructions, exhibited artifacts due to localized overweighting of certain orientations. These could be automatically discarded by rejecting highly non-homogeneous distributions of orientations. Within the remaining results, two main classes could be observed which are related by an overall rotation of about 180 around an axis close to one of the coordinate axes. This is in accordance with a previous study for an isotropic orientational distribution . After manual attribution to one of the two classes, the results were averaged and their relation was verified by orientational registration. The final averaged 3D reciprocal-space volume hW(q)i was obtained as an average of 13 individual results in the same orientation as W ref (q).

Reciprocal space (intensity)
A comparison of hW(q)i with W ref is shown in Fig. 2. Visually, the orthogonal slices through the reconstructed and reference intensities are in very good agreement. This observation is reflected by an overall Pearson correlation coefficient

Validation (intensity reconstruction)
To assess the reliability of the reconstructed 3D intensity in reciprocal space, we randomly assigned the frames of the data set to two independent half-data sets and reconstructed two independent 3D reciprocal-space volumes as described before. A Fourier shell correlation (FSC) curve ( van Heel & Schatz, 2005), obtained here directly from the two reciprocal-space volumes, is shown in Fig. 3 (red line). It intersects the half-bit threshold curve ( van Heel & Schatz, 2005) at a value beyond q = 120, indicating a self-consistent reconstruction of the reciprocal-space volume close to the Nyquist limit. Here, q is measured in units of ÁXk/D, where ÁX = 150 m is the detector pixel pitch, k is the wave number and D is the sampleto-detector distance. (Note that, in the limit of a flat Ewald sphere patch, which is well fulfilled here, q then corresponds to the distance to the central beam in units of detector pixels. For further details, see the supporting information.) This analysis shows the effect of EMC alone rather than merging the effects of orientation determination and phasing, as would result for a traditional FSC analysis on phased results in real space. For comparison, the FSC curve resulting from correlating hW(q)i with W ref , as obtained from known orientations, is also shown (blue line). The high degree of similarity between these two curves indicates strong agreement between the result obtained by EMC and the reference intensity distribution, assembled using full knowledge of the orientations.

Real space (density)
To obtain the real-space electron-density distribution, the missing phases of the 3D Fourier intensity need to be determined. To this end, we applied standard iterative phase retrieval to the 3D diffraction data, i.e. a combination of the hybrid input-output (HIO) algorithm and the error reduction (ER) algorithm (Fienup & Wackerman, 1986;Marchesini et al., 2003;Xiong et al., 2014). In total, 600 iterations were applied, i.e. 420 iterations of HIO with a feedback parameter = 0.9, followed by 180 iterations of ER. For further details see the supporting information.
To ensure the reproducibility of the obtained result, 60 reconstructions were performed in total. The results were filtered in a two-step selection process. In the first, manual, step, we discarded images which visually deviated from the most abundant reconstruction result. In the second step, the Orthogonal slices through (a)-(c) the EMC-reconstructed and (d)-(f) the reference 3D diffraction volumes. The EMC-reconstructed diffraction volume hW(q)i results from averaging the results of 13 independent EMC runs, each starting with a random intensity distribution. The reference diffraction volume W ref (q) was constructed based on the known orientations of the sample for each frame during the measurement. The dashed circles in panels (c) and (f) indicate a radius of 127 voxels, whereas the solid circles indicate radii of 20 and 50, respectively. All slices are drawn on the same scale with dimensionless lateral coordinates in units of ÁXk/D.

Figure 3
The red line illustrates the Fourier shell correlation (FSC) between two EMC-retrieved reciprocal-space volumes resulting from splitting the data set into two equal halves and performing the same analysis on them as on the whole data set. The blue line indicates the FSC between hW(q)i and W ref , i.e. the reciprocal-space volume resulting from analyzing the whole data set using EMC and the reference intensity assembled using known orientations. The green line denotes the half-bit threshold curve, used as a common criterion for resolution determination in analysis of FSC curves. remaining reconstructions (43 for the reference and 37 for EMC data set) were aligned with sub-pixel precision (Guizar-Sicairos et al., 2008) and averaged. Then the 20 reconstructions showing the highest correlation with this average were selected for the final average. Note that the procedure applied in the second step could also be used to avoid any manual intervention. However, in such a case, several iterations would likely be required in order to avoid bias by strong outliers in the average. The resulting real-space reconstructions, i.e. the real part of the final average, from both reciprocal-space intensities are shown in Fig. 4.
All details of the reference reconstruction are reproduced in the EMC-based reconstruction down to a resolution level of very few pixels. A comparison of a scanning electron microscopy image of the sample with an iso-surface rendering of the EMC-based reconstruction shows that height variations due to imperfections in the fabrication process are well reproduced by the reconstruction. This identifies the sample as a true 3D structure with features in all coordinate directions.

Validation (density reconstruction)
The resolution of the final image was estimated via the phase retrieval transfer function (PRTF) according to a procedure similar to the one described by Chapman et al. (2006). More specifically, before summation of the complexvalued reconstructions, their constant phases were adjusted so that the real part of each reconstruction was maximized. The PRTF curves for the results shown in Fig. 4 are shown in Fig. 5. The full-period resolution, as determined by the spatial frequency corresponding to a PRTF value of 1/e, amounts to a value between 40 and 45 nm. This corresponds to 24 to 26 (full-period) resolution elements within the largest linear extension of the particle, as given by the smallest sphere completely containing the particle (diameter ' 1.1 m).

Significance of the observed level of sparsity
To assess the significance of the sparsity level for a single detector frame in the present experiment, we have calculated the expected average total number of scattered photons outside the central speckle for a selection of 35 000 human protein structures from the RCSB Protein Data Bank (PDB; http://www.rcsb.org) (Berman et al., 2000). Here, a focal-spot diameter of 300 nm was assumed, at a photon energy of 8 keV and a pulse energy of 1 mJ (with 20% beamline transmission). As a result, it could be shown that, under these realistic conditions, the minimum diameter for a protein to scatter 50 photons outside the central speckle amounts to 10.6 nm. For further details, see the supporting information.
This clearly shows that the signal level in the present experiment, obtained at a synchrotron source from a nanofabricated gold structure, is comparable with what can be expected under realistic conditions from a relevant protein structure at an FEL source.

Particle complexity, rotation group sampling, SNR and sparsity
Another parameter to be discussed is the particle complexity R, as measured in half-resolution units per particle radius (Loh & Elser, 2009). Reconstructing a particle with 10 nm diameter down to a resolution of 3 Å results in a complexity of R ' 33, far beyond the current state of the art for FEL-based SPI (Lundholm et al., 2018) which is R ' 8 for a globular virus particle. The present structure reaches a Reconstruction of the 3D electron density. (a) Reconstruction from the result derived by EMC. The electron density projected along an axis perpendicular to the drawing plane is shown here. (b) Reconstruction from the reference Fourier volume. Again, the projected electron density is shown. (c) 3D iso-surface rendering of the reconstructed electron density shown in panel (a). The threshold of the iso-surface has been set to 0.2, given a normalized density with values between 0 and 1 (see also animation in the supporting information). (d) Scanning electron micrograph from the original sample.

Figure 5
Phase retrieval transfer functions for the reconstructions from the EMCgenerated Fourier space intensity hW(q)i and the reference intensity W ref .
The curves decay to a value of 1/e between q = 90 and q = 100, corresponding to a half-period resolution between 20 and 23 nm. complexity >20 in two dimensions, being constrained in the height direction to a value between 2 and 3. A comparison between a scanning electron microscopy (SEM) image and an isosurface rendering of the reconstructed particle density shows that all features in the height direction are very well reproduced. Despite its flat shape, this clearly identifies the particle as a true 3D structure and underlines the significance of the present result as a step forward towards the complexity level required for real protein structures.
Furthermore, compared with a serial imaging experiment at an FEL, where thousands of particles in random 3D orientations contribute to a full data set, the number of unique 3D orientations contributing to the present data set seems relatively low (M rot = 227). It can be shown, however, that at the given resolution and complexity this does not restrict the relevance of the result. The required minimum angular separation between adjacent orientations for a sufficient sampling of the 3D rotation group is linked to the complexity R of the particle (Loh et al., 2010): = 1/R. In case of a nonglobular shape, the maximum complexity in a given coordinate direction allows for a conservative estimate, leading here to = 1/R max = 2.2 . This shows that a finer sampling for the tomographic series contributing to the present data set would not have added more information, at the obtained resolution.
Evidently, the present experiment profits from a high SNR (see the supporting information) which would have been impossible without a setup well optimized for forward-scattering CDI (Chushkin et al., 2014). Most importantly, this consists of a set of accurately placed apertures upstream of the sample which are adjusted to define the beam incident on the sample and, at the same time, to suppress scattering arising from upstream apertures by those further downstream. Similar schemes can be applied at FEL sources to make them compatible with CDI experiments (Munke et al., 2016). Even though in the FEL case the aerosol jet in which sample particles are injected through a stream of carrier gas causes an additional source of background scatter , the present data set gives an experimental benchmark for a signal-to-noise level which would likely allow FEL-based single-particle imaging.
To provide some understanding of how orientation information can be retrieved despite the sparsity of photon measurements per pattern, consider the two key attributes of a scheme that succeeds in this retrieval. First, this scheme must recognize when it is presented with the correct 3D Fourier intensity distribution, and second, it must contain instructions to improve the compatibility between a candidate intensity distribution and the set of measured diffraction patterns. To satisfy the first attribute of recognizability, the EMC algorithm has to model the detection noise and the randomness of the orientations as accurately as possible. In the present case, this is achieved by a Poissonian noise model and an equal weight for each possible orientation. If a compatible 3D intensity distribution is recognized, then the most likely orientations of each photon pattern within this distribution can be estimated naturally. To fulfill the second attribute, EMC uses this noise model to modify any candidate 3D intensity distribution to increase its likelihood of generating the set of measured diffraction patterns. A derivation of the corresponding maximization rule is given by Loh & Elser (2009). If there are a sufficient number of photons per pattern, rotation group sampling, signal-to-background ratio and total number of measurements, then EMC can seek and recognize a family of compatible 3D intensity distributions even if the initial guess were completely random (Loh & Elser, 2009). A more quantitative definition of sufficiency for this orientation problem, beyond what is already in the literature (e.g. Loh & Elser, 2009;Elser, 2009;Philipp et al., 2012;Ayyer et al., 2014;Loh, 2014), is the subject of a future publication by some of the co-authors here.

Summary and conclusions
In summary, we have demonstrated experimentally, in the same geometry as used for FEL-based single-particle imaging, the reconstruction of a complex three-dimensional object using coherent diffractive imaging (CDI) from photon-sparse random projections at a sparsity level to be expected for a typical protein at an FEL source. To this end, we have collected 454 000 data frames with about 50 scattered photons per frame, evenly distributed over 227 unique orientations, and reconstructed a consistent 3D reciprocal-space volume without explicit knowledge of the orientation of a given frame.
It has been shown that, by application of the expandmaximize-compress (EMC) algorithm, both the reconstructed reciprocal-space intensity and the real-space density of the sample agree to a high level with reconstructions obtained using complete knowledge of the frame orientations.
The data set is made freely available in the CXI data bank (Maia, 2012; http://cxidb.org/), entry 84, to be used as a testbed for algorithm development for CDI-based single-particle imaging, e.g. by alternative methods of orientation determination. In addition, the data set can serve as a target for a signal-to-noise level enabling FEL-based single-particle imaging in the future. ments were performed on beamline ID10 at the European Synchrotron Radiation Facility (ESRF), Grenoble, France. The work was done partially while K. Giewekemeyer was visiting the Institute for Mathematical Sciences, National University of Singapore, in 2018. The visit was supported by the Institute. Duane Loh and Colin Teo would like to acknowledge funding support from the Singapore National Research Foundation's Competitive Research Program funding (NRF-CRP16-2015-05). This research used resources from the DOE Early Career Program and Stanford Synchrotron Radiation Lightsource, SLAC National Accelerator Laboratory, supported by the US Department of Energy, Office of Science, Office of Basic Energy Sciences, under contract number DE-AC02-76SF00515. Part of this work was performed at the Stanford Nano Shared Facilities (SNSF), supported by the National Science Foundation under award ECCS-1542152.