Algorithmic framework for X-ray nanocrystallographic reconstruction in the presence of the indexing ambiguity Jeffrey J. Donatellia,b and James A. Sethiana,b,1 a
Department of Mathematics and bLawrence Berkeley National Laboratory, University of California, Berkeley, CA 94720
X-ray nanocrystallography allows the structure of a macromolecule to be determined from a large ensemble of nanocrystals. However, several parameters, including crystal sizes, orientations, and incident photon flux densities, are initially unknown and images are highly corrupted with noise. Autoindexing techniques, commonly used in conventional crystallography, can determine orientations using Bragg peak patterns, but only up to crystal lattice symmetry. This limitation results in an ambiguity in the orientations, known as the indexing ambiguity, when the diffraction pattern displays less symmetry than the lattice and leads to data that appear twinned if left unresolved. Furthermore, missing phase information must be recovered to determine the imaged object’s structure. We present an algorithmic framework to determine crystal size, incident photon flux density, and orientation in the presence of the indexing ambiguity. We show that phase information can be computed from nanocrystallographic diffraction using an iterative phasing algorithm, without extra experimental requirements, atomicity assumptions, or knowledge of similar structures required by current phasing methods. The feasibility of this approach is tested on simulated data with parameters and noise levels common in current experiments.
A
lthough conventional X-ray crystallography has been used extensively to determine atomic structure, it is limited to objects that can be formed into large crystal samples ð>10 μmÞ. An appealing alternative, made possible by recent advances in light source technology, is X-ray nanocrystallography, which is able to image structures resistant to large crystallization, such as membrane proteins, by substituting a large ensemble of easier to build nanocrystals, typically 0, we use the known reference configuration B = ðb1 ; b2 ; b3 Þ to compute an approximation ~ a Þ is far from unity ~ a = BD−1 to the orientation matrix. If detðR R then the autoindexing procedure failed to compute an accurate ~ a . We then find the closest rotaorientation and we thus reject R tion matrix by computing the singular value decomposition ~ = UV T , which is used as the approxi~ a = UΣV T and setting R R mation to the image orientation up to lattice symmetry, i.e., ~ = QR where R is the full orientation and . R Crystal Size Determination To compute accurate structure factor magnitudes from measured intensities, the squared magnitude of the shape transform must be divided out of the intensity measurements in Eq. 2. Near a Bragg peak, the shape transform grows quadratically with the crystal size, which often varies up to an order of magnitude in each dimension over the nanocrystal ensemble. We determine these crystal sizes by analyzing intensities around Bragg peaks in lowangle images (Fig. 2) sampled at least twice the Nyquist rate for 1 , where W is the width the crystal, i.e., the pixel spacing is at most 2W of the crystal. A Fourier analysis of these intensities reveals the crystal sizes. Fourier Analysis of the Shape Transform. For an image I with orientation§ R, consider its restriction Ir to a small neighborhood centered at a low-angle Bragg peak with detector coordinates xo ∈ R2 corresponding to the reciprocal lattice point , where jξj is small. In , by taking a linear approximation of q and using the translation invariance of S on , the intensities are approximated, up to a constant C, by 2 [4] Ir ðxÞ ≈ CS K R ~x ;
2 where ~x = ðx; 0Þ and K = ðλDÞ−1 . If we denote GðxÞ = CSðK~xÞ , then Eq. 4 becomes the restriction of G to the rotated plane Ir ≈ GjRðR2 Þ , whose Fourier transform,{ from the Fourier projection slice theorem and the Wiener–Khinchin theorem, is approximately the X-ray projected autocorrelation of :k [5] Note that the support of this projected autocorrelation is given by the Minkowski sum of the rotated projected crystal, i.e., supp . If
we
approximate** the crystal lattice as , where N = ðN1 ; N2 ; N3 Þ are the crystal
§
Because the shape transform is symmetric with respect to rotation by elements of the use of orientations from autoindexing is sufficient here.
,
{^
Ir may be slightly smeared out due to grid alignment effects.
k
PRð3Þ is the X-ray projection operator through the detector plane normal and A is the autocorrelation operator.
**We note that alternative models of the finite lattice may be used here instead.
Donatelli and Sethian
By computing ^I r , we can deduce the crystal sizes by analyzing H as long as none of the rotated Bravais vectors RT hj is orthogonal to the detector. In general, the boundary of H consists of a series of line segments, with three normals ðn1 ; n2 ; n3 Þ along with their three antiparallel directions, which can be found by computing the convex hull of the rotated projected autocorrelated unit cell, i.e., the right-hand side of Eq. 6 with Nj = 2 for each j. In the direction ni, the extent of the convex hull bi = maxfjni · xj : x ∈ Hg to that of the unpro be equal P must jected crystal, i.e., bi = 3j=1 nTi RT hj ðNj − 1Þ. Therefore, by defining the matrix Aij = nTi RT hj , the crystal sizes N can be retrieved by solving AðN − 1Þ = b, where 1 = ð1; 1; 1Þ. If the image does not directly pass through a reciprocal lattice point, this analysis is still valid; however, the Fourier-transformed images may contain oscillations, which grow with distance to the lattice point. Image Segmentation. To retrieve the crystal sizes via the methods of the previous section, we require an estimate for the support of the projected autocorrelated crystal from the Fourier transformed images, i.e., we need to segment the support from the noisy background (Fig. 3).We begin the segmentation by initializing a set H of pixels whose ^I r value is greater than a fixed percentage of the largest†† value. Then, we traverse a sorted list of the remaining values, adding pixels to H until one reaches a point more than some threshold, typically a few pixels, away from all of the pixels currently in H, suggesting that one has reached the end of the support and has begun to see the oscillations from the background.‡‡
Structure Factor Magnitude Modeling ~ m and crystal sizes Nm are known, Once the lattice orientations R 2 we can use Eq. 2 to compute an approximation ~ F m to the
structure factor square magnitudes from the image Im, but only up to a constant factor because they are scaled by the unknown incident photon flux density Jm, which varies between images. Furthermore, due to the indexing ambiguity, one only knows the corresponding reciprocal space coordinates associated with the values of ~ F m up to the crystal lattice symmetry, i.e., the possible structure factor magnitudes for each point take the form of a multimodal distribution. Moreover, these two problems are strongly coupled together: We cannot perform the scaling correction unless we know what modes to scale to and the modes are indistinguishable in the unscaled data set. Hence, we must simultaneously determine both the scaling and the multimodal parameters. Processing the Data. We approximate the structure factor squared magnitudes at the ith reciprocal lattice point§§ ξi from the mth image by computing an average over the neighboring ball Bðξi ; rÞ with radius r:
††
We exclude the origin which picks up all of the noise within the image.
‡‡
H gives unitless coordinates which should be scaled by λD=ðNp dxÞ for a restricted image with Np × Np pixels of size dx × dx. Crystal sizes are rejected if they are outside of a set range.
§§
To keep our notation compact, we are representing the reciprocal lattice points with a single index in place of the traditional Miller indices.
PNAS | January 14, 2014 | vol. 111 | no. 2 | 595
APPLIED MATHEMATICS
j=1
frequency
frequency
Fig. 4. (Left) Histogram of the possible unscaled variance stabilized structure factor magnitudes for a reciprocal lattice point corresponding to a fourfold indexing ambiguity. (Right) Histogram of the scaled data with multimodal Gaussian model (red).
Scaling Correction. In practice, variance in incident photon flux density, noise, and errors in autoindexing and crystal size determination smear out the peaks in the histogram, making them difficult to locate via expectation maximization (Fig. 4): Data must be scaled to properly model the structure factor magnitudes. To do so, we seek scaling factors which minimize the variance in the histograms and alternate this procedure with the expectation maximization step in Eq. 8. We seek the scaling factor cm for the mth image by solving
min cm
P vi;m =
~ qðxÞ∈Bðξi ;rÞ Im ðxÞ R m
P
~ R m
2 : 2 PðxÞΔΩðxÞS ~ qðxÞ R r m e qðxÞ∈Bðξi ;rÞ m
[7]
Multimodal Analysis. Assume for now that the structure factor
magnitudes are already properly scaled. Due to the indexing ambiguity, at this point in the procedure we only know the orientations up to the lattice symmetry. Thus, for each reciprocal lattice point ξi , values of wi;m could correspond to K different structure factor . A histomagnitudes, e.g., for elastic scattering gram of fwi;m g for ξi reveals K different peaks, smeared out as noise and parameter uncertainty are increased. Our goal is to detect these peaks and model the associated multimodal distribution. To retrieve the set of possible structure factor magnitudes, we will model the computed values fwi;m g from each reciprocal lattice point ξi with a multimodal Gaussian distribution. Specifically, the associated probability density functions can be expressed in terms of multiple Gaussian distributions with means μi = ðμi;1 ; . . . ; μi;K Þ, in monotonically increasing order, and standard deviations by
,
where . M Given fwi;m gm=1 , we determine its multimodal model through an expectation maximization algorithm. In particular, given an ð0Þ ð0Þ initial guess for model parameters μi;j and σ i;j , we perform several iterations of the following:
[8]
[9]
i;j
whose solution is given by P
2 i;j wi;m μi;j Ti;j;m : 2 2 i;j wi;m Ti;j;m
cm = P
We discard any intensities below some fixed threshold from the above sum to prevent large errors from the division. Note that because we only know the orientation up to lattice symmetry, we also set{{ vt;m = vi;m for every ξt such that for some , Rξt = ξi . To simplify notation, we will assume that unmeasured values and corresponding indices are removed from the remaining sets and summations. To reduce the dependence of the standard deviation on the size of the intensities, we use two applications of variance stabilization 1=4 (12), and instead work with wi;m = vi;m .
σ i = ðσ i;1 ; . . . ; σ i;K Þ
2 X cm wi;m − μi;j Ti;j;m ;
[10]
Once the cm are computed, they are normalized so that P m cm = constant and then used to scale the images by replacing every wi;m with cm wi;m . Scaling is alternated with expectation maximization until convergence. Resolving the Indexing Ambiguity After the structure factor magnitude modeling, we know up to K possible structure factor magnitudes at each reciprocal lattice point. Resolving the indexing ambiguity amounts to correctly assigning one of these K values to each point. There are K equally valid solutions, related to each other by globally . We first use the set of multiapplying a rotation from modal model parameters to construct a graph theoretic model of the structure factor magnitude concurrency, i.e., the probability that two given structure factor magnitudes occur within the same image. Then, we resolve the indexing ambiguity by finding the maximum edge weight clique of this graph with a greedy approach. Graph Theoretic Modeling of Structure Factor Magnitude Concurrency.
Given the scaled variance stabilized structure factor magnitudes fwi;m g, means fμi;j g, and standard deviations fσ i;j g for the mth image and jth mode at the ith reciprocal lattice point, we constructkk a graph G = ðV ; EÞ with vertices V = fði; jÞg and edges E, where ðði1 ; j1 Þ; ði2 ; j2 ÞÞ ∈ E if and only if i1 ≠ i2 and , i.e., only one j can be selected at each reciprocal lattice point ξi and each j can only be selected once among its twin-related coordinates ξt = Rξi , where . Consequently, choosing a consistent set of structure factor magnitudes, where each possible value appears exactly once, is equivalent to finding a maximal clique*** in this graph.††† We define a directed weight W : E → R on G as the ratio of the sum of concurrence probabilities over the sum of the occurrence probabilities, where we sum over the sets I i1 ;i2 , consisting of all of the images which potentially measure intensities simultaneously at ξi1 and ξi2 : P m∈I i ;i Ti1 ;j1 ;m Ti2 ;j2 ;m W ðði1 ; j1 Þ; ði2 ; j2 ÞÞ = P 1 2 : [11] m∈I i ;i Ti1 ;j1 ;m 1 2
Here, Ti;j;m represents the probability that wi;m is drawn from the jth Gaussian mode. To initialize, we separate the data into K ð0Þ equal bins, define μi;j as the location in the jth bin with the ð0Þ greatest density, and set σ i;j to be less than a typical bin size. Also, we perform outlier rejection by removing any wi;m in which pðwi;m ; μi ; σ i Þ is below a given threshold.
W gives a measure of how likely the values of μi1 ;j1 and μi2 ;j2 are to occur together in one of the solutions. With no noise or error, if μi2 ;1 ; . . . ; μi2 ;K are distinct then W is exactly 1 when μi1 ;j1 and μi2 ;j2 kk
For known symmetry in the structure factor magnitudes (e.g., Friedel or Laue symmetry), we simplify G’s structure by merging corresponding symmetric vertices.
~ ***Cl⊂V is a maximal clique if for all v1 ,v2 ∈ Cl, ðv1 ,v2 Þ ∈ E with no proper superset Cl⊋Cl satisfying the same property. {{
vi,m is multivalued for an image which measures intensities at both ξ and Rξ,
596 | www.pnas.org/cgi/doi/10.1073/pnas.1321790111
.
†††
AnB = fx ∈ A : x ∉ Bg is the set difference of A and B.
Donatelli and Sethian
Table 1. Autoindexing performance Error* Jo
0.016
32,678 28,251 22,701 12,859 3,926
7,912 6,583 4,733 1,786 39
21,485 17,993 13,588 6,224 798
3,196 3,496 3,947 4,180 2,290
85 179 433 669 799
*Error in the Frobenius norm modulo
are simultaneously part of one of the solutions and is 0 otherwise. However, if there are B values of k such that μi2 ;k = μi2 ;j2 then W is 1=B. With noise present, the asymmetry of the weight function favors structure factor magnitudes with strong signal and whose histograms are highly multimodal, providing the most orientation information. We now formulate the solution to the indexing ambiguity problem. We seek the maximal clique in G with maximum edge weight: X max W ðv1 ; v2 Þ; [12] Cl∈CmðGÞ
.
Orientation Calculation. Although we can achieve a robust approximation of the complete structure factor magnitudes, accuracy can be improved by using this information to directly orient each image, and then averaging structure factor magnitudes for each corresponding reciprocal lattice point. For every image Im , we ~ m , where Qm solves compute its full orientation Rm = Qm R
v1 ;v2 ∈Cl
where CmðGÞ is the set of all maximal cliques in G. If a sufficient number of images are used, then, in the absence of noise and error, the maximizer of Eq. 12 retrieves one of the exact solutions to the indexing ambiguity problem, i.e., the maximum edge weight clique Cl assigns the correct structure factor magnitudes, up to a global rotation. For imperfect data, this maximizer will choose a solution which is most consistent with the observed structure factor magnitude concurrency.
[13] Here, A is the set of indices i where wi;m is measured by the ~ m and image at the orientation R is a set of class representatives of the quotient group , i.e., consists of the identity and twinning operators. If there are at least two orientations close to the minimum value, then we reject the computed orientation. Once the orientations are known, we obtain the structure factor magnitudes by averaging the corresponding magnitudes computed from each scaled oriented image.
Greedy Approach to the Maximum Edge Weight Clique Problem.
Even though the maximum edge weight clique problem is, in general, nondeterministic polynomial-time hard (13), when constructed from the indexing ambiguity problem via Eq. 12, we can solve it in quadratic time with a greedy approach. We initialize the clique Cl with some starting vertex‡‡‡ vs ∈ V and progressively add vertices that maximize the weight sum of the current clique. In practice, we remove any points whose associated multimodal distributions contain less than the maximum number of modes. For convenience, we use a single index for the vertices §§§ V = fvi gN i=1 and we set W ðv1 ; v2 Þ = − ∞ if ðv1 ; v2 Þ ∉ E. Algorithm 1.
The elements of the set Cl returned by Algorithm 1 are pairs of the form ði; jÞ, corresponding to choosing the jth modeled variance stabilized structure factor magnitude μi;j at the reciprocal lattice ^ → R where w ~:L ~ ðξi Þ = μi;j . With point ξi . This induces the map w a sufficient number of images and no noise or error, Algorithm 1 retrieves an exact solution to the indexing ambiguity problem, always preferring a vertex vn with nonzero weighted edges connecting to all elements of the current clique, i.e., vn is only chosen if it corresponds to the same solution as the rest of clique. This approach remains robust with imperfect data, as it considers several pairs of measured intensities over all images to choose the structure factor magnitudes at any single point.
Phase Recovery Once complete structure factor magnitudes are computed, one can determine missing phases and, thus, determine the electron density of a periodic unit with any applicable phasing method, e.g., refs. 14–17. Although these methods have been used extensively to determine structure in X-ray crystallography, each has limitations or introduces extra difficulties into the experimental setup. An appealing alternative is to directly deduce phases from Fourier magnitude data using an iterative phase retrieval technique, which only requires that Fourier magnitudes be sampled at a sufficient rate. Although infeasible in conventional crystallography, the signal from nanocrystals contains significant information between Bragg peaks, which may allow sampling at the required rates. In general, such iterative phasing is possible if P one samples the ^ i ; where Fourier magnitudes at points of the form ξ = 3i=1 n2i h ni ∈ Z (18, 19). In ref. 20 the feasibility of such an approach was demonstrated assuming that adequate signal was collected at each of the required points. However, the square magnitude of the shape transform at ξ grows quadratically in the crystal size for each dimension in which ξ · hi ∈ Z. For nanocrystallography, this typically only results in noticeable signal at reciprocal lattice vertices and edges.{{{,kkk While theory is lacking, this sampling density is sufficient in certain 3D cases (21). Alternatively, a recent approach uses Fourier magnitude information along with its gradient, assuming it can be accurately calculated, only at reciprocal lattice points (22). Here, we test the feasibility of iterative phasing using the Fourier magnitudes computed from our framework at reciprocal lattice vertices and edges. {{{
The lattice edge structure factor magnitudes may be computed via Eq. 7.
‡‡‡
§§§
vs is typically chosen from a point with a highly multimodal distribution and strong signal.
In Algorithm 1,
denotes the assignment operator.
Donatelli and Sethian
kkk
One may also need to resolve the indexing ambiguity on the used non-Bragg data if it has less symmetry than the Laue group.
****This is obtained by linearly mapping the reciprocal lattice onto a uniform grid.
PNAS | January 14, 2014 | vol. 111 | no. 2 | 597
APPLIED MATHEMATICS
Fig. 5. Electron density contours of the exact solution (Left) and the computed solution with Jo = 218 (Center) and overlay with the atomic model (Right).
21,800 2,180 218 21.8 2.18
Accepted
Table 2. Crystal size determination performance
Table 3. Orientation determination performance
Error* Jo 21,800 2,180 218 21.8 2.18
Jo
Accepted
0.4
30,506 26,260 20,195 10,639 1,851
19,652 16,080 10,941 4,071 458
9,667 8,996 7,681 4,626 458
903 992 1,166 1,388 387
98 71 126 184 112
186 121 281 370 436
*Relative error in the geometric average of the crystal sizes.
Iterative Phase Retrieval. Given a domain**** Zα = Zα1 × Zα2 × Zα3 , Fourier magnitude values a : Zα → C, and some support S⊆Zα , 2 iterative phase retrieval algorithms seek a function ρ ∈ ℓ ðZα Þ such that ^ ρ = a and ρðxÞ = 0 for all x ∉ S. Given a set Ω of points where a has a recorded value, such algorithms typically make use of the projection operators PS , where PS ρðxÞ = ρðxÞ if x ∈ S and PρðxÞ = 0 otherwise, and PM;Ω , where 8 > ρðkÞ ^ > > > aðkÞ^ρðkÞ; if ^ρ ≠ 0 and k ∈ Ω; > < Pd [14] M;Ω ρðkÞ = aðkÞ; if ^ρ = 0 and k ∈ Ω; > > > > > : ^ρðkÞ; if k ∉ Ω:
We alternate between several iterations of the error reducing algorithm (23): ρðn+1Þ = PS PM;Ω ρðnÞ , and the hybrid input-output algorithm (24): ρðn+1Þ = ðPS PM;Ω + PSc ðI − βPM;Ω ÞÞρðnÞ , β ∈ ð0; 1, to seek out the solution ρ. Furthermore, we couple these iterations with the Shrinkwrap method (25), which, starting with an initial guess such as the unit cell, updates an estimate of the true support T by convolving the current iterate with a Gaussian and then thresholding. Results We demonstrate our methodology by determining the structure of PuuE Allantoinase from simulated diffraction data using the atomic coordinates and crystal symmetry recorded in refs. 26, 27 with several different peak incident photon flux densities. Each data set consists of 33,856 diffraction images. We assume knowledge of the Bravais vector lengths and the space group, 1. Aquila A, et al. (2012) Time-resolved protein nanocrystallography using an X-ray freeelectron laser. Opt Express 20(3):2706–2716. 2. Boutet S, et al. (2012) High-resolution protein structure determination by serial femtosecond crystallography. Science 337(6092):362–364. 3. Chapman HN, et al. (2011) Femtosecond X-ray protein nanocrystallography. Nature 470(7332):73–77. 4. Johansson LC, et al. (2012) Lipidic phase membrane protein serial femtosecond crystallography. Nat Methods 9(3):263–265. 5. Kern J, et al. (2013) Simultaneous femtosecond X-ray spectroscopy and diffraction of photosystem II at room temperature. Science 340(6131):491–495. 6. Koopmann R, et al. (2012) In vivo protein crystallization opens new routes in structural biology. Nat Methods 9(3):259–262. 7. Spence JCH, Weierstall U, Chapman HN (2012) X-ray lasers for structural and dynamic biology. Rep Prog Phys 75(10):102601. 8. White TA, et al. (2012) Crystfel: A software suite for snapshot serial crystallography. J Appl Cryst 45(2):335–341. 9. Duisenberg AJM (1992) Indexing in single-crystal diffractometry with an obstinate list of reflections. J Appl Cryst 25:92–96. 10. Steller I, Bolotovsky R, Rossmann MG (1997) An algorithm for automatic indexing of oscillation images using Fourier analysis. J Appl Cryst 30:1036–1040. 11. Lovisolo L, da Silva EAB (2001) Uniform distribution of points on a hyper-sphere with applications to vector bit-plane encoding. IEEE Proc Vision Image Signal Process 148(3):187–193. 12. Guan Y (2009) Variance stabilizing transformations of Poisson, binomial and negative binomial distributions. Stat Probab Lett 14:1621–1629. 13. Alidaee B, Glover F, Kochenberger G, Wang H (2007) Solving the maximum edge weight clique problem via unconstrained quadratic programming. Eur J Oper Res 181(2):592–597. 14. Green DW, Ingram VM, Perutz MF (1954) The structure of haemoglobin. IV. Sign determination by the isomorphous replacement method. Proc. Roy. Soc. A. 225(1162):287–307.
598 | www.pnas.org/cgi/doi/10.1073/pnas.1321790111
Accepted Correct, %*
21,800
2,180
218
21.8
2.18
22,645 99.9
17,844 99.9
11,663 99.8
3,775 99.3
17 58.8
*Given the possible solution sets fR1,i g and fR2,i g, the number of correct orientations is maxj=1,2 Sj , where Sj consists of all computed orientations Ri closer to the Rj,i solution in the Frobenius norm modulo .
which, in practice, may be deduced from autoindexing information and reflection conditions. The crystal displays P4 space-group symmetry, thus diffraction data are symmetric with respect to 90° rotation about the z axis and inversion, and has a twinning operator given by 180° rotation about the x axis, or, equivalently, the y axis. Image orientations are generated by randomly sampling from a normal distribution of quaternions. Sizes were randomly generated with an average crystal width of μC = 2,948.97 Å and standard deviation σ C = 982.99 Å. For each image, we generate random incident photon flux densities J, measured in photons per square Angstrom per pulse, from a peak density Jo , via x ∼ Uð−1; 1Þ; J = Jo expð−8x2 Þ. We use experimental parameters†††† similar to ref. 3. Intensity values are computed via Eq. 2 along with shot and background noise, modeled with a Poisson distribution and a normal distribution, with a standard deviation of 1.3 photons per pixel (Figs. S1–S5). Here we present statistics‡‡‡‡ for our framework (Tables 1–3) along with a reconstruction from the processed simulated data (Fig. 5). For further details, see ref. 28. ACKNOWLEDGMENTS. We thank Stefano Marchesini for many valuable conversations. This research was supported in part by the Applied Mathematical Sciences subprogram of the Office of Energy Research, US Department of Energy (DOE) under Contract DE-AC02-05CH11231 and by the Division of Mathematical Sciences of the National Science Foundation, and used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the US DOE under Contract DE-AC02-05CH11231. J.A.S. was also supported by an Einstein Visiting Fellowship of the Einstein Foundation, Berlin. J.J.D. was also supported by a DOE Computational Science Graduate Fellowship.
††††
λ = 6:9 Å, dx = 75 μm, D = 68/141 mm for the front/rear detectors with 1,024 × 1,024 pixels, horizontal polarization, and Jo between 0.01 and 100 times current experimental levels.
‡‡‡‡
Here we use the distance in Frobenius norm modulo
.
15. Hauptman H, Karle J (1953) Solution of the Phase Problem I. The Centrosymmetric Crystal, No. 3 (American Crystallographic Association, New York). 16. Karle J (1980) Some developments in anomalous dispersion for the structural investigation of macromolecular systems in biology. Int J Quantum Chem Quantum Biol Symp 18(S7): 357–367. 17. Rossmann MG, Blow DM (1962) The detection of sub-units within the crystallographic asymmetric unit. Acta Crystallogr 15(1):24–31. 18. Hayes MH (1982) The reconstruction of a multidimensional sequence from the phase or magnitude of its Fourier transform. IEEE Trans Acoust Speech Signal Process 30(2):140–154. 19. Rosenblatt J (1984) Phase retrieval. Commun Math Phys 95:317–343. 20. Spence JCH, et al. (2011) Phasing of coherent femtosecond X-ray diffraction from sizevarying nanocrystals. Opt Express 19(4):2866–2873. 21. Millane RP (1996) Multidimensional phase problems. J Opt Soc Am A Opt Image Sci Vis 13(4):725–734. 22. Elser V (2013) Direct phasing of nanocrystal diffraction. Acta Crystallogr A 69:559–569. 23. Gerchberg RW, Saxton WO (1972) A practical algorithm for the determination of the phase from image and diffraction plane pictures. Optik (Stuttg) 35:237–246. 24. Fienup JR (1978) Reconstruction of an object from the modulus of its Fourier transform. Opt Lett 3(1):27–29. 25. Marchesini S, et al. (2003) X-ray image reconstruction from a diffraction pattern alone. Phys Rev B 68(4):140101. 26. Bernstein FC, et al. (1977) The Protein Data Bank: A computer-based archival file for macromolecular structures. J Mol Biol 112(3):535–542. 27. Ramazzina I, et al. (2008) Logical identification of an allantoinase analog (puuE) recruited from polysaccharide deacetylases. J Biol Chem 283(34):23295–23304. 28. Donatelli JJ (2013) Reconstruction algorithms for x-ray nanocrystallography via solution of the twinning problem, PhD Thesis (University of California, Berkeley).
Donatelli and Sethian
Department of Mathematics and bLawrence Berkeley National Laboratory, University of California, Berkeley, CA 94720
X-ray nanocrystallography allows the structure of a macromolecule to be determined from a large ensemble of nanocrystals. However, several parameters, including crystal sizes, orientations, and incident photon flux densities, are initially unknown and images are highly corrupted with noise. Autoindexing techniques, commonly used in conventional crystallography, can determine orientations using Bragg peak patterns, but only up to crystal lattice symmetry. This limitation results in an ambiguity in the orientations, known as the indexing ambiguity, when the diffraction pattern displays less symmetry than the lattice and leads to data that appear twinned if left unresolved. Furthermore, missing phase information must be recovered to determine the imaged object’s structure. We present an algorithmic framework to determine crystal size, incident photon flux density, and orientation in the presence of the indexing ambiguity. We show that phase information can be computed from nanocrystallographic diffraction using an iterative phasing algorithm, without extra experimental requirements, atomicity assumptions, or knowledge of similar structures required by current phasing methods. The feasibility of this approach is tested on simulated data with parameters and noise levels common in current experiments.
A
lthough conventional X-ray crystallography has been used extensively to determine atomic structure, it is limited to objects that can be formed into large crystal samples ð>10 μmÞ. An appealing alternative, made possible by recent advances in light source technology, is X-ray nanocrystallography, which is able to image structures resistant to large crystallization, such as membrane proteins, by substituting a large ensemble of easier to build nanocrystals, typically 0, we use the known reference configuration B = ðb1 ; b2 ; b3 Þ to compute an approximation ~ a Þ is far from unity ~ a = BD−1 to the orientation matrix. If detðR R then the autoindexing procedure failed to compute an accurate ~ a . We then find the closest rotaorientation and we thus reject R tion matrix by computing the singular value decomposition ~ = UV T , which is used as the approxi~ a = UΣV T and setting R R mation to the image orientation up to lattice symmetry, i.e., ~ = QR where R is the full orientation and . R Crystal Size Determination To compute accurate structure factor magnitudes from measured intensities, the squared magnitude of the shape transform must be divided out of the intensity measurements in Eq. 2. Near a Bragg peak, the shape transform grows quadratically with the crystal size, which often varies up to an order of magnitude in each dimension over the nanocrystal ensemble. We determine these crystal sizes by analyzing intensities around Bragg peaks in lowangle images (Fig. 2) sampled at least twice the Nyquist rate for 1 , where W is the width the crystal, i.e., the pixel spacing is at most 2W of the crystal. A Fourier analysis of these intensities reveals the crystal sizes. Fourier Analysis of the Shape Transform. For an image I with orientation§ R, consider its restriction Ir to a small neighborhood centered at a low-angle Bragg peak with detector coordinates xo ∈ R2 corresponding to the reciprocal lattice point , where jξj is small. In , by taking a linear approximation of q and using the translation invariance of S on , the intensities are approximated, up to a constant C, by 2 [4] Ir ðxÞ ≈ CS K R ~x ;
2 where ~x = ðx; 0Þ and K = ðλDÞ−1 . If we denote GðxÞ = CSðK~xÞ , then Eq. 4 becomes the restriction of G to the rotated plane Ir ≈ GjRðR2 Þ , whose Fourier transform,{ from the Fourier projection slice theorem and the Wiener–Khinchin theorem, is approximately the X-ray projected autocorrelation of :k [5] Note that the support of this projected autocorrelation is given by the Minkowski sum of the rotated projected crystal, i.e., supp . If
we
approximate** the crystal lattice as , where N = ðN1 ; N2 ; N3 Þ are the crystal
§
Because the shape transform is symmetric with respect to rotation by elements of the use of orientations from autoindexing is sufficient here.
,
{^
Ir may be slightly smeared out due to grid alignment effects.
k
PRð3Þ is the X-ray projection operator through the detector plane normal and A is the autocorrelation operator.
**We note that alternative models of the finite lattice may be used here instead.
Donatelli and Sethian
By computing ^I r , we can deduce the crystal sizes by analyzing H as long as none of the rotated Bravais vectors RT hj is orthogonal to the detector. In general, the boundary of H consists of a series of line segments, with three normals ðn1 ; n2 ; n3 Þ along with their three antiparallel directions, which can be found by computing the convex hull of the rotated projected autocorrelated unit cell, i.e., the right-hand side of Eq. 6 with Nj = 2 for each j. In the direction ni, the extent of the convex hull bi = maxfjni · xj : x ∈ Hg to that of the unpro be equal P must jected crystal, i.e., bi = 3j=1 nTi RT hj ðNj − 1Þ. Therefore, by defining the matrix Aij = nTi RT hj , the crystal sizes N can be retrieved by solving AðN − 1Þ = b, where 1 = ð1; 1; 1Þ. If the image does not directly pass through a reciprocal lattice point, this analysis is still valid; however, the Fourier-transformed images may contain oscillations, which grow with distance to the lattice point. Image Segmentation. To retrieve the crystal sizes via the methods of the previous section, we require an estimate for the support of the projected autocorrelated crystal from the Fourier transformed images, i.e., we need to segment the support from the noisy background (Fig. 3).We begin the segmentation by initializing a set H of pixels whose ^I r value is greater than a fixed percentage of the largest†† value. Then, we traverse a sorted list of the remaining values, adding pixels to H until one reaches a point more than some threshold, typically a few pixels, away from all of the pixels currently in H, suggesting that one has reached the end of the support and has begun to see the oscillations from the background.‡‡
Structure Factor Magnitude Modeling ~ m and crystal sizes Nm are known, Once the lattice orientations R 2 we can use Eq. 2 to compute an approximation ~ F m to the
structure factor square magnitudes from the image Im, but only up to a constant factor because they are scaled by the unknown incident photon flux density Jm, which varies between images. Furthermore, due to the indexing ambiguity, one only knows the corresponding reciprocal space coordinates associated with the values of ~ F m up to the crystal lattice symmetry, i.e., the possible structure factor magnitudes for each point take the form of a multimodal distribution. Moreover, these two problems are strongly coupled together: We cannot perform the scaling correction unless we know what modes to scale to and the modes are indistinguishable in the unscaled data set. Hence, we must simultaneously determine both the scaling and the multimodal parameters. Processing the Data. We approximate the structure factor squared magnitudes at the ith reciprocal lattice point§§ ξi from the mth image by computing an average over the neighboring ball Bðξi ; rÞ with radius r:
††
We exclude the origin which picks up all of the noise within the image.
‡‡
H gives unitless coordinates which should be scaled by λD=ðNp dxÞ for a restricted image with Np × Np pixels of size dx × dx. Crystal sizes are rejected if they are outside of a set range.
§§
To keep our notation compact, we are representing the reciprocal lattice points with a single index in place of the traditional Miller indices.
PNAS | January 14, 2014 | vol. 111 | no. 2 | 595
APPLIED MATHEMATICS
j=1
frequency
frequency
Fig. 4. (Left) Histogram of the possible unscaled variance stabilized structure factor magnitudes for a reciprocal lattice point corresponding to a fourfold indexing ambiguity. (Right) Histogram of the scaled data with multimodal Gaussian model (red).
Scaling Correction. In practice, variance in incident photon flux density, noise, and errors in autoindexing and crystal size determination smear out the peaks in the histogram, making them difficult to locate via expectation maximization (Fig. 4): Data must be scaled to properly model the structure factor magnitudes. To do so, we seek scaling factors which minimize the variance in the histograms and alternate this procedure with the expectation maximization step in Eq. 8. We seek the scaling factor cm for the mth image by solving
min cm
P vi;m =
~ qðxÞ∈Bðξi ;rÞ Im ðxÞ R m
P
~ R m
2 : 2 PðxÞΔΩðxÞS ~ qðxÞ R r m e qðxÞ∈Bðξi ;rÞ m
[7]
Multimodal Analysis. Assume for now that the structure factor
magnitudes are already properly scaled. Due to the indexing ambiguity, at this point in the procedure we only know the orientations up to the lattice symmetry. Thus, for each reciprocal lattice point ξi , values of wi;m could correspond to K different structure factor . A histomagnitudes, e.g., for elastic scattering gram of fwi;m g for ξi reveals K different peaks, smeared out as noise and parameter uncertainty are increased. Our goal is to detect these peaks and model the associated multimodal distribution. To retrieve the set of possible structure factor magnitudes, we will model the computed values fwi;m g from each reciprocal lattice point ξi with a multimodal Gaussian distribution. Specifically, the associated probability density functions can be expressed in terms of multiple Gaussian distributions with means μi = ðμi;1 ; . . . ; μi;K Þ, in monotonically increasing order, and standard deviations by
,
where . M Given fwi;m gm=1 , we determine its multimodal model through an expectation maximization algorithm. In particular, given an ð0Þ ð0Þ initial guess for model parameters μi;j and σ i;j , we perform several iterations of the following:
[8]
[9]
i;j
whose solution is given by P
2 i;j wi;m μi;j Ti;j;m : 2 2 i;j wi;m Ti;j;m
cm = P
We discard any intensities below some fixed threshold from the above sum to prevent large errors from the division. Note that because we only know the orientation up to lattice symmetry, we also set{{ vt;m = vi;m for every ξt such that for some , Rξt = ξi . To simplify notation, we will assume that unmeasured values and corresponding indices are removed from the remaining sets and summations. To reduce the dependence of the standard deviation on the size of the intensities, we use two applications of variance stabilization 1=4 (12), and instead work with wi;m = vi;m .
σ i = ðσ i;1 ; . . . ; σ i;K Þ
2 X cm wi;m − μi;j Ti;j;m ;
[10]
Once the cm are computed, they are normalized so that P m cm = constant and then used to scale the images by replacing every wi;m with cm wi;m . Scaling is alternated with expectation maximization until convergence. Resolving the Indexing Ambiguity After the structure factor magnitude modeling, we know up to K possible structure factor magnitudes at each reciprocal lattice point. Resolving the indexing ambiguity amounts to correctly assigning one of these K values to each point. There are K equally valid solutions, related to each other by globally . We first use the set of multiapplying a rotation from modal model parameters to construct a graph theoretic model of the structure factor magnitude concurrency, i.e., the probability that two given structure factor magnitudes occur within the same image. Then, we resolve the indexing ambiguity by finding the maximum edge weight clique of this graph with a greedy approach. Graph Theoretic Modeling of Structure Factor Magnitude Concurrency.
Given the scaled variance stabilized structure factor magnitudes fwi;m g, means fμi;j g, and standard deviations fσ i;j g for the mth image and jth mode at the ith reciprocal lattice point, we constructkk a graph G = ðV ; EÞ with vertices V = fði; jÞg and edges E, where ðði1 ; j1 Þ; ði2 ; j2 ÞÞ ∈ E if and only if i1 ≠ i2 and , i.e., only one j can be selected at each reciprocal lattice point ξi and each j can only be selected once among its twin-related coordinates ξt = Rξi , where . Consequently, choosing a consistent set of structure factor magnitudes, where each possible value appears exactly once, is equivalent to finding a maximal clique*** in this graph.††† We define a directed weight W : E → R on G as the ratio of the sum of concurrence probabilities over the sum of the occurrence probabilities, where we sum over the sets I i1 ;i2 , consisting of all of the images which potentially measure intensities simultaneously at ξi1 and ξi2 : P m∈I i ;i Ti1 ;j1 ;m Ti2 ;j2 ;m W ðði1 ; j1 Þ; ði2 ; j2 ÞÞ = P 1 2 : [11] m∈I i ;i Ti1 ;j1 ;m 1 2
Here, Ti;j;m represents the probability that wi;m is drawn from the jth Gaussian mode. To initialize, we separate the data into K ð0Þ equal bins, define μi;j as the location in the jth bin with the ð0Þ greatest density, and set σ i;j to be less than a typical bin size. Also, we perform outlier rejection by removing any wi;m in which pðwi;m ; μi ; σ i Þ is below a given threshold.
W gives a measure of how likely the values of μi1 ;j1 and μi2 ;j2 are to occur together in one of the solutions. With no noise or error, if μi2 ;1 ; . . . ; μi2 ;K are distinct then W is exactly 1 when μi1 ;j1 and μi2 ;j2 kk
For known symmetry in the structure factor magnitudes (e.g., Friedel or Laue symmetry), we simplify G’s structure by merging corresponding symmetric vertices.
~ ***Cl⊂V is a maximal clique if for all v1 ,v2 ∈ Cl, ðv1 ,v2 Þ ∈ E with no proper superset Cl⊋Cl satisfying the same property. {{
vi,m is multivalued for an image which measures intensities at both ξ and Rξ,
596 | www.pnas.org/cgi/doi/10.1073/pnas.1321790111
.
†††
AnB = fx ∈ A : x ∉ Bg is the set difference of A and B.
Donatelli and Sethian
Table 1. Autoindexing performance Error* Jo
0.016
32,678 28,251 22,701 12,859 3,926
7,912 6,583 4,733 1,786 39
21,485 17,993 13,588 6,224 798
3,196 3,496 3,947 4,180 2,290
85 179 433 669 799
*Error in the Frobenius norm modulo
are simultaneously part of one of the solutions and is 0 otherwise. However, if there are B values of k such that μi2 ;k = μi2 ;j2 then W is 1=B. With noise present, the asymmetry of the weight function favors structure factor magnitudes with strong signal and whose histograms are highly multimodal, providing the most orientation information. We now formulate the solution to the indexing ambiguity problem. We seek the maximal clique in G with maximum edge weight: X max W ðv1 ; v2 Þ; [12] Cl∈CmðGÞ
.
Orientation Calculation. Although we can achieve a robust approximation of the complete structure factor magnitudes, accuracy can be improved by using this information to directly orient each image, and then averaging structure factor magnitudes for each corresponding reciprocal lattice point. For every image Im , we ~ m , where Qm solves compute its full orientation Rm = Qm R
v1 ;v2 ∈Cl
where CmðGÞ is the set of all maximal cliques in G. If a sufficient number of images are used, then, in the absence of noise and error, the maximizer of Eq. 12 retrieves one of the exact solutions to the indexing ambiguity problem, i.e., the maximum edge weight clique Cl assigns the correct structure factor magnitudes, up to a global rotation. For imperfect data, this maximizer will choose a solution which is most consistent with the observed structure factor magnitude concurrency.
[13] Here, A is the set of indices i where wi;m is measured by the ~ m and image at the orientation R is a set of class representatives of the quotient group , i.e., consists of the identity and twinning operators. If there are at least two orientations close to the minimum value, then we reject the computed orientation. Once the orientations are known, we obtain the structure factor magnitudes by averaging the corresponding magnitudes computed from each scaled oriented image.
Greedy Approach to the Maximum Edge Weight Clique Problem.
Even though the maximum edge weight clique problem is, in general, nondeterministic polynomial-time hard (13), when constructed from the indexing ambiguity problem via Eq. 12, we can solve it in quadratic time with a greedy approach. We initialize the clique Cl with some starting vertex‡‡‡ vs ∈ V and progressively add vertices that maximize the weight sum of the current clique. In practice, we remove any points whose associated multimodal distributions contain less than the maximum number of modes. For convenience, we use a single index for the vertices §§§ V = fvi gN i=1 and we set W ðv1 ; v2 Þ = − ∞ if ðv1 ; v2 Þ ∉ E. Algorithm 1.
The elements of the set Cl returned by Algorithm 1 are pairs of the form ði; jÞ, corresponding to choosing the jth modeled variance stabilized structure factor magnitude μi;j at the reciprocal lattice ^ → R where w ~:L ~ ðξi Þ = μi;j . With point ξi . This induces the map w a sufficient number of images and no noise or error, Algorithm 1 retrieves an exact solution to the indexing ambiguity problem, always preferring a vertex vn with nonzero weighted edges connecting to all elements of the current clique, i.e., vn is only chosen if it corresponds to the same solution as the rest of clique. This approach remains robust with imperfect data, as it considers several pairs of measured intensities over all images to choose the structure factor magnitudes at any single point.
Phase Recovery Once complete structure factor magnitudes are computed, one can determine missing phases and, thus, determine the electron density of a periodic unit with any applicable phasing method, e.g., refs. 14–17. Although these methods have been used extensively to determine structure in X-ray crystallography, each has limitations or introduces extra difficulties into the experimental setup. An appealing alternative is to directly deduce phases from Fourier magnitude data using an iterative phase retrieval technique, which only requires that Fourier magnitudes be sampled at a sufficient rate. Although infeasible in conventional crystallography, the signal from nanocrystals contains significant information between Bragg peaks, which may allow sampling at the required rates. In general, such iterative phasing is possible if P one samples the ^ i ; where Fourier magnitudes at points of the form ξ = 3i=1 n2i h ni ∈ Z (18, 19). In ref. 20 the feasibility of such an approach was demonstrated assuming that adequate signal was collected at each of the required points. However, the square magnitude of the shape transform at ξ grows quadratically in the crystal size for each dimension in which ξ · hi ∈ Z. For nanocrystallography, this typically only results in noticeable signal at reciprocal lattice vertices and edges.{{{,kkk While theory is lacking, this sampling density is sufficient in certain 3D cases (21). Alternatively, a recent approach uses Fourier magnitude information along with its gradient, assuming it can be accurately calculated, only at reciprocal lattice points (22). Here, we test the feasibility of iterative phasing using the Fourier magnitudes computed from our framework at reciprocal lattice vertices and edges. {{{
The lattice edge structure factor magnitudes may be computed via Eq. 7.
‡‡‡
§§§
vs is typically chosen from a point with a highly multimodal distribution and strong signal.
In Algorithm 1,
denotes the assignment operator.
Donatelli and Sethian
kkk
One may also need to resolve the indexing ambiguity on the used non-Bragg data if it has less symmetry than the Laue group.
****This is obtained by linearly mapping the reciprocal lattice onto a uniform grid.
PNAS | January 14, 2014 | vol. 111 | no. 2 | 597
APPLIED MATHEMATICS
Fig. 5. Electron density contours of the exact solution (Left) and the computed solution with Jo = 218 (Center) and overlay with the atomic model (Right).
21,800 2,180 218 21.8 2.18
Accepted
Table 2. Crystal size determination performance
Table 3. Orientation determination performance
Error* Jo 21,800 2,180 218 21.8 2.18
Jo
Accepted
0.4
30,506 26,260 20,195 10,639 1,851
19,652 16,080 10,941 4,071 458
9,667 8,996 7,681 4,626 458
903 992 1,166 1,388 387
98 71 126 184 112
186 121 281 370 436
*Relative error in the geometric average of the crystal sizes.
Iterative Phase Retrieval. Given a domain**** Zα = Zα1 × Zα2 × Zα3 , Fourier magnitude values a : Zα → C, and some support S⊆Zα , 2 iterative phase retrieval algorithms seek a function ρ ∈ ℓ ðZα Þ such that ^ ρ = a and ρðxÞ = 0 for all x ∉ S. Given a set Ω of points where a has a recorded value, such algorithms typically make use of the projection operators PS , where PS ρðxÞ = ρðxÞ if x ∈ S and PρðxÞ = 0 otherwise, and PM;Ω , where 8 > ρðkÞ ^ > > > aðkÞ^ρðkÞ; if ^ρ ≠ 0 and k ∈ Ω; > < Pd [14] M;Ω ρðkÞ = aðkÞ; if ^ρ = 0 and k ∈ Ω; > > > > > : ^ρðkÞ; if k ∉ Ω:
We alternate between several iterations of the error reducing algorithm (23): ρðn+1Þ = PS PM;Ω ρðnÞ , and the hybrid input-output algorithm (24): ρðn+1Þ = ðPS PM;Ω + PSc ðI − βPM;Ω ÞÞρðnÞ , β ∈ ð0; 1, to seek out the solution ρ. Furthermore, we couple these iterations with the Shrinkwrap method (25), which, starting with an initial guess such as the unit cell, updates an estimate of the true support T by convolving the current iterate with a Gaussian and then thresholding. Results We demonstrate our methodology by determining the structure of PuuE Allantoinase from simulated diffraction data using the atomic coordinates and crystal symmetry recorded in refs. 26, 27 with several different peak incident photon flux densities. Each data set consists of 33,856 diffraction images. We assume knowledge of the Bravais vector lengths and the space group, 1. Aquila A, et al. (2012) Time-resolved protein nanocrystallography using an X-ray freeelectron laser. Opt Express 20(3):2706–2716. 2. Boutet S, et al. (2012) High-resolution protein structure determination by serial femtosecond crystallography. Science 337(6092):362–364. 3. Chapman HN, et al. (2011) Femtosecond X-ray protein nanocrystallography. Nature 470(7332):73–77. 4. Johansson LC, et al. (2012) Lipidic phase membrane protein serial femtosecond crystallography. Nat Methods 9(3):263–265. 5. Kern J, et al. (2013) Simultaneous femtosecond X-ray spectroscopy and diffraction of photosystem II at room temperature. Science 340(6131):491–495. 6. Koopmann R, et al. (2012) In vivo protein crystallization opens new routes in structural biology. Nat Methods 9(3):259–262. 7. Spence JCH, Weierstall U, Chapman HN (2012) X-ray lasers for structural and dynamic biology. Rep Prog Phys 75(10):102601. 8. White TA, et al. (2012) Crystfel: A software suite for snapshot serial crystallography. J Appl Cryst 45(2):335–341. 9. Duisenberg AJM (1992) Indexing in single-crystal diffractometry with an obstinate list of reflections. J Appl Cryst 25:92–96. 10. Steller I, Bolotovsky R, Rossmann MG (1997) An algorithm for automatic indexing of oscillation images using Fourier analysis. J Appl Cryst 30:1036–1040. 11. Lovisolo L, da Silva EAB (2001) Uniform distribution of points on a hyper-sphere with applications to vector bit-plane encoding. IEEE Proc Vision Image Signal Process 148(3):187–193. 12. Guan Y (2009) Variance stabilizing transformations of Poisson, binomial and negative binomial distributions. Stat Probab Lett 14:1621–1629. 13. Alidaee B, Glover F, Kochenberger G, Wang H (2007) Solving the maximum edge weight clique problem via unconstrained quadratic programming. Eur J Oper Res 181(2):592–597. 14. Green DW, Ingram VM, Perutz MF (1954) The structure of haemoglobin. IV. Sign determination by the isomorphous replacement method. Proc. Roy. Soc. A. 225(1162):287–307.
598 | www.pnas.org/cgi/doi/10.1073/pnas.1321790111
Accepted Correct, %*
21,800
2,180
218
21.8
2.18
22,645 99.9
17,844 99.9
11,663 99.8
3,775 99.3
17 58.8
*Given the possible solution sets fR1,i g and fR2,i g, the number of correct orientations is maxj=1,2 Sj , where Sj consists of all computed orientations Ri closer to the Rj,i solution in the Frobenius norm modulo .
which, in practice, may be deduced from autoindexing information and reflection conditions. The crystal displays P4 space-group symmetry, thus diffraction data are symmetric with respect to 90° rotation about the z axis and inversion, and has a twinning operator given by 180° rotation about the x axis, or, equivalently, the y axis. Image orientations are generated by randomly sampling from a normal distribution of quaternions. Sizes were randomly generated with an average crystal width of μC = 2,948.97 Å and standard deviation σ C = 982.99 Å. For each image, we generate random incident photon flux densities J, measured in photons per square Angstrom per pulse, from a peak density Jo , via x ∼ Uð−1; 1Þ; J = Jo expð−8x2 Þ. We use experimental parameters†††† similar to ref. 3. Intensity values are computed via Eq. 2 along with shot and background noise, modeled with a Poisson distribution and a normal distribution, with a standard deviation of 1.3 photons per pixel (Figs. S1–S5). Here we present statistics‡‡‡‡ for our framework (Tables 1–3) along with a reconstruction from the processed simulated data (Fig. 5). For further details, see ref. 28. ACKNOWLEDGMENTS. We thank Stefano Marchesini for many valuable conversations. This research was supported in part by the Applied Mathematical Sciences subprogram of the Office of Energy Research, US Department of Energy (DOE) under Contract DE-AC02-05CH11231 and by the Division of Mathematical Sciences of the National Science Foundation, and used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the US DOE under Contract DE-AC02-05CH11231. J.A.S. was also supported by an Einstein Visiting Fellowship of the Einstein Foundation, Berlin. J.J.D. was also supported by a DOE Computational Science Graduate Fellowship.
††††
λ = 6:9 Å, dx = 75 μm, D = 68/141 mm for the front/rear detectors with 1,024 × 1,024 pixels, horizontal polarization, and Jo between 0.01 and 100 times current experimental levels.
‡‡‡‡
Here we use the distance in Frobenius norm modulo
.
15. Hauptman H, Karle J (1953) Solution of the Phase Problem I. The Centrosymmetric Crystal, No. 3 (American Crystallographic Association, New York). 16. Karle J (1980) Some developments in anomalous dispersion for the structural investigation of macromolecular systems in biology. Int J Quantum Chem Quantum Biol Symp 18(S7): 357–367. 17. Rossmann MG, Blow DM (1962) The detection of sub-units within the crystallographic asymmetric unit. Acta Crystallogr 15(1):24–31. 18. Hayes MH (1982) The reconstruction of a multidimensional sequence from the phase or magnitude of its Fourier transform. IEEE Trans Acoust Speech Signal Process 30(2):140–154. 19. Rosenblatt J (1984) Phase retrieval. Commun Math Phys 95:317–343. 20. Spence JCH, et al. (2011) Phasing of coherent femtosecond X-ray diffraction from sizevarying nanocrystals. Opt Express 19(4):2866–2873. 21. Millane RP (1996) Multidimensional phase problems. J Opt Soc Am A Opt Image Sci Vis 13(4):725–734. 22. Elser V (2013) Direct phasing of nanocrystal diffraction. Acta Crystallogr A 69:559–569. 23. Gerchberg RW, Saxton WO (1972) A practical algorithm for the determination of the phase from image and diffraction plane pictures. Optik (Stuttg) 35:237–246. 24. Fienup JR (1978) Reconstruction of an object from the modulus of its Fourier transform. Opt Lett 3(1):27–29. 25. Marchesini S, et al. (2003) X-ray image reconstruction from a diffraction pattern alone. Phys Rev B 68(4):140101. 26. Bernstein FC, et al. (1977) The Protein Data Bank: A computer-based archival file for macromolecular structures. J Mol Biol 112(3):535–542. 27. Ramazzina I, et al. (2008) Logical identification of an allantoinase analog (puuE) recruited from polysaccharide deacetylases. J Biol Chem 283(34):23295–23304. 28. Donatelli JJ (2013) Reconstruction algorithms for x-ray nanocrystallography via solution of the twinning problem, PhD Thesis (University of California, Berkeley).
Donatelli and Sethian
