1. Introduction

With the advent of the X-ray free electron laser, it is possible to image a single biomolecule through coherent diffraction imaging [1–4], albeit at lower wavelengths. Many promising approaches have been suggested for accomplishing this goal [5–18]. In a typical coherent X-ray imaging experiment, single biomolecules are injected into a pulsed X-ray beam. The elastically scattered photons are recorded by a detector array placed downstream. Due to the ultrashort time scale of the pulse, the scattering pattern can be recorded before the molecule is destroyed [19]. The recorded scattering data are then sorted according to the orientations of the biomolecules. Signal strength can be improved by averaging shots of the biomolecules having the same orientation, which will benefit subsequent steps of phase retrieval and image reconstruction. Due to the small size of the biomolecule and small elastic scattering cross-section of the electron, we encounter two problems for the coherent X-ray diffraction imaging. One is that the detected signal will be mostly dominated by Poisson noise. This is especially the case when the scattering angle is large, as the scattering is mostly in the forward direction for coherent illumination. The other difficulty is that the biomolecule itself is intrinsically a 3D structure and its orientation is not controlled. Thus, orientation classification is the first crucial step of data processing in coherent X-ray imaging.

Early work on classification was based on large scattering angle data where the diffraction pattern is modeled as a speckle pattern [12, 20]. In practice, this is primarily true for large-angle scattering where the detected photon level is very low with the current or near future X-ray photon budgets. Thus the classification schemes based on large-angle-scattering-data are not currently practical. Recently there are two interesting developments where Expectation-Maximization algorithms is utilized for classification [21, 22]. The signal level for successful classifications is much reduced. Another imaging scheme is to simultaneously inject many randomly oriented biomolecules [23]. The scattered photons are much more and the orientation classification is unnecessary in this setup.

In this paper we focus on scattering by a single particle and propose a classification scheme based on patterns of low scattering angles where a substantial number of photons can be detected. Further the comparison among different rotations is not directly based on the pattern itself but on compressed signals which are the radial and azimuthal distributions of the scattering pattern. As both of these distributions are signals integrated over many pixels, the variance in the radial or azimuthal distribution is much smaller than that in the individual pixels. Consequentially a classification can be relatively immune to the shot noise and the results of the classification can be trustworthy. An example of utilizing those compressed signals is given in this paper for classification of in-plane rotations (i.e., relative orientations between biomolecules are within a plane parallel to the detector plane).

2. Weak elastic scattering process and the issue of shot noise

In a typical setup for coherent diffraction imaging of single biomolecules, the samples are injected into the waist of the focused coherent X-ray beam and a detector is placed downstream to record the scattering pattern. The scattered pattern recorded by the detector array, I(x, y), can be written as

In this paper we also assume that the noise is mainly shot noise, the detected photon number in any pixel is a Poisson distributed random variable whose mean value can be calculated from Eq. (1). In Fig. 1 we show a typical scattering pattern by the single biomolecule when the incident photon number is infinite (noise free), 2 × 10¹⁴ and 2 × 10¹² per pulse. The beam is assumed to be 100 nm and the X-ray wavelength is 1.5Å. The detector array is 150 mm away from the sample and the pixel size is 0.11 mm.

 

Fig. 1 Typical simulated scattering patterns of a DNA polymerase delta (PDB ID:3IAY). The data are shown in logarithmic scale. The radius of the hole in the center of the detector array is 3.72 mm. Left: “noise free” pattern; Middle: incident photon number is 2 × 10¹⁴ and the number of collected photons is 818875; Right: incident photon number is 2 × 10¹² and the number of collected photons is 8114.

Download Full Size | PDF

From Fig. 1 we can see that for the “noise free” case, the scattering pattern looks like speckles in regions corresponding to large scattering angles. However, with a finite incident photon number such as 2 × 10¹⁴ photons per pulse, the contrast of the speckles is diminished and some even disappear. For the 2 × 10¹² photons per pulse case, which is the current available level at LCLS, the number of photons scattered with scattering angles larger than 2° is so little that basically all speckles disappear. It will be very difficult to distinguish any exposures from each other, not to mention conducting orientation classification with this incident photon level.

Consider an annular area of the scattering pattern with radius ρ_min ≤ ρ ≤ ρ_max. We divide this area into regions of radial or azimuthal segments as shown in Fig. 2(a) and 2(b). The boundaries between these radial segments are rings with radii ρ ₁, ρ ₂,..., ρ_N, which can be defined as

 

Fig. 2 Detector array is divided into radial (a) or azimuthal segments (b).

Download Full Size | PDF

In the simulation study we choose δθ = 1°. As the scattering pattern is symmetric about the origin, we limit the value of θ to be within [0, 180°]. The azimuthal distribution can be written as:

In order to study the effect of the noise, the mean value of the radial and azimuthal distributions are first calculated from the “noise free” diffraction patterns. The noisy versions of the radial or azimuthal distribution are calculated from the noisy diffractions when the incident photon number is 2 × 10¹⁴ and 2 × 10¹² respectively. A typical radial and azimuthal distributions and their mean values are shown in Fig. 3(a) and 3(b) respectively. We see that when the incident photon number is 2 × 10¹⁴, the noisy version of the radial and azimuthal distribution is similar to its mean value. However when the incident photon number is 2 × 10¹², the noisy version of the radial and azimuthal distribution deviates substantially from their corresponding mean values, which will affect the reliability of the classification adversely.

 

Fig. 3 Typical radical (a) and azimuthal (b) distributions. Three cases are studied: the mean value case (dotted line); the 2 × 10¹⁴ photon number case (solid line) and the 2 × 10¹² photon number case (dashed line).

Download Full Size | PDF

In Figs. 4 we show surface plots when the biomolecule rotates about the x, y, and z axes. We see that with rotations about x or y, both the radial and the azimuthal distributions change smoothly with rotation angles. For rotations about z, the radial distribution stays the same while the azimuthal distribution shifts continuously with the rotation angle. Thus rotation about z (also called an in-plane rotation) can be recognized by studying the evolution of the radial distributions, and the rotation angles can be identified by finding the amount of the shift of the azimuthal distribution. For rotations other than along the z-axis, both the radial and azimuthal distributions are going to change smoothly with the rotation angles. Thus a distance between two distributions, e.g. Euclidean distance, can be a measure describing the global changes during continuous rotation. A small distance between distributions means that the two distributions are quite similar to each other and the relative rotation angle between them is small. In a separate simulation we modified a large biomolecule (PDB ID: 3I55, size ∼ 200Å) by removing all atoms greater than 90Å from the centroid and showed that the same process can be easily applied to spherical molecules.

 

Fig. 4 Surface plots of the radial (left column) and azimuthal (right column) distributions of the scattering patterns when the biomolecule rotates about the x-axis (a,b), y-axis (c,d), and z-axis (e–f). The incident photon number is 2 × 10¹⁴.

Download Full Size | PDF

3. Euler angles and the possibility of in-plane rotation

In the previous sections we have devised a method to preprocess the scattering data to improve the data consistency for the classification of orientations. In this section, we will use introduce Euler angles, which have been widely used in the cryo-EM field, to represent arbitrary orientations [24–26].

An arbitrary rotation of an angle ϕ about an arbitrary axis can be described by three consecutive rotations, i.e.,

As the low-angle scattering pattern is proportional to the square of the amplitude of the Fourier transform of the scattering potential, which is a real and even function, we can limit the ranges of α, β and γ to [0, π].

The last γ-rotation is the so-called in-plane rotation because the z-axis of the intrinsic coordinate and the z-axis of the laboratory coordinates are co-axial. If we can classify the in-plane rotation out of all possible rotations and know the angle of the in-plane rotation of a selected pattern with respect to a reference pattern, we can rotate that pattern to overlap with the reference pattern. This will give us an opportunity to increase the SNR by averaging over all the in-plane rotations. The extent of the SNR improvement depends on the number of in-plane rotations.

Consider the random orientation of a single biomolecule when injected into the passage of the X-ray beam, the probability of the orientation which can be described by the Euler angles (α, β, γ) and can be written as [27]

Thus the probability density function is uniform within a rectangular volume of 2π ² with Euler rotations defined by α, γ and cosβ. With a total number N_total of random rotations, the number density for α and γ is $N_{\mathit{\text{total}}}^{1/3}/\mathit{\text{pi}}$ and $N_{\mathit{\text{total}}}^{1/3}/2$ for cos β. Thus the mean number of in-plane rotations is expected to be:

4. In-plane rotation recognition through the radial distribution

We can also see from the surface plots shown in Fig. 4(e) that the radial distribution stays unchanged with continuous rotation about z. Thus, we may identify the in-plane rotation by studying the radial distribution.

Consider a distance defined as

 

Fig. 5 The radial distance as a function of the rotation angle when the protein rotates about y, and z respectively. The curve for rotation along z is shown in pink; in blue for rotation along y. Two incident photon levels are studied here: (a) 10¹⁴; (b) 10¹²; The mean values of the distances are shown in the solid lines; The color-shaded areas represent ±σ around the mean values for both the rotation about y and z respectively. σ is the variance. (c)–(d) are zoomed-in versions of (a) and (b) respectively.

Download Full Size | PDF

5. Finding the in-plane rotation angles through the azimuthal distribution

After having identified scattering patterns corresponding to in-plane rotations, we can find the angle of the in-plane rotation by studying the azimuthal distribution I_θ. This angle information is needed if we hope to average over in-plane rotations. In this section we describe a simple method to find this angle.

If $I_{\theta }^{\mathit{\text{ref}}}$ is the azimuthal distribution of the reference pattern and $I_{\theta }^{\mathit{\text{test}}}$ is the azimuthal distribution of the test pattern with only an in-plane rotation γ with respect to the reference pattern, then these two distributions should be shifted with respect to each other. This can be clearly seen in the surface plot shown in Fig. 4(f). The amount of the shift should be the angle of the in-plane rotation of the biomolecule with respect to its original position.

Consider a new function, $I_{\theta }^{\prime }$, which is a circularly-shifted version of $I_{\theta }^{\mathit{\text{test}}}$, i.e.,

 

Fig. 6 Azimuthal distance reaches minimum at the in-plane rotation angle. Two incident photon number levels are studied: 2 × 10¹⁴(solid line) and 2 × 10¹²(dashed line).

Download Full Size | PDF

With in-plane rotations classified (section 5) and the angle of in-plane rotation determined as described in this section, we can rotate those in-plane rotated patterns to overlap the reference patterns, and average those patterns to improve the signal-to-noise ratio. In Fig. 7 we show a typical single shot of the scattering pattern and the averaged pattern from 19 shots with different in-plane rotations. Comparing these two images we can clearly see that the signal-to-noise ratio is improved after the averaging, especially in regions away from the center.

 

Fig. 7 (a) Scattering pattern from a single measurement; (b) Scattering pattern after averaging over 19 patterns with different in-plane rotations; (c) “Noise free” pattern.

Download Full Size | PDF

6. Accuracy of the classification vs photon budget

As we saw in the previous sections, the simple distances can be used to compare the radial or azimuthal distributions to distinguish between in-plane rotations and out-of-plane rotations, and to find the degree of in-plane rotation. In this section we are going to study the performance of the classification scheme for different incident photon numbers. To be specific, we want to know the ratio of the correct classification results from the results of the proposed classification method.

In order to visualize the success rate of the proposed algorithm, we have generated a number of random orientations and their corresponding diffraction patterns, then used the method described in Section 5 to find the in-plane rotations. The Euler angles for those orientations are limited to {0°, 180°}. Twenty random values for each α, cosβ, γ are generated independently. Thus a total number of orientations is 20 × 20 × 20 = 8000, ignoring the inconsequential fact that the ranges for those three random variables (α, cosβ, γ) are not all equal to each other. The diffraction patterns for the protein with those orientations are generated using Eq. (1). Then we select an arbitrary pattern as a reference frame and assign the rest of the patterns to be the test patterns. The corresponding radial distribution for the reference frame and the test patterns are computed according to Eq. (3). Then the distance is calculated using Eq. (9) to compare the reference and test frames.

Consider a total number N_total of patterns from randomly oriented biomolecules, the total number of pairwise distances, N_d, will be N_total × (N_total – 1) if each pattern is compared to other patterns in the data set. Note that the comparison between any two patterns are performed twice since each one has the opportunity to be a reference pattern. We also assume that the number of random values for each Euler angle is $N_{\mathit{\text{total}}}^{1/3}$. Thus the percentage for in-plane rotations among all possible rotations shall be $N_{\mathit{\text{total}}}^{1/3}/N_{\mathit{\text{total}}}$. The total number of distances that describe in-plane rotations, $N_{\mathit{\text{in}}}^d$, is expected to be

As we have generated the patterns and know the orientations of the test frames, we can verify the results of the classification with the known information. Thus we can define the accuracy of the classification as:

Table 1. Accuracy of the In-Plane Rotation Classification^a

View Table

We can see that even with the current incident level, the success of the classification can be as high as 79% with a class size of 1°, which is required for atomic resolution. If we can relax the requirement on the sensitivity of the classification, we can achieve 92% accuracy for the class size of 3° and 97% accuracy for 5°. In this way we can break the data set into smaller groups and use more sophisticated algorithms such as described in [21] or [22] to find additional information. Research in this direction is underway. With higher photon levels, almost 99% accuracy can be achieved for class size of 1°. The accuracy from a mere guess is expected to be N_inplane/N_total, which is only 0.25% in this case.

Note that with the limited number of generated random orientations due to available computing speed, the orientations which are close to each other in the R ³ space may not have been fully simulated in the above example. Thus we perform another round of simulations where the Euler angles are limited to α ∈ (0, 10°), cosβ ∈ (0.9, 1). As the radial distributions are invariant to in-plane rotations, we choose the limit on the in-plane rotation angle still to be (0, 180°). For each Euler angle we generate 10 random numbers resulting in 1000 patterns. With a much smaller range of Euler angles, the percentage for two orientations whose Euler angles (α_i, β_i, i = 1, 2) are within 1° becomes much higher. In our simulation, the percentage is about 48%, i.e., almost half of the orientations have a similar orientation. The accuracy in this small data set can be lower compared to the entire data set. In Table 1 we list the accuracy (numbers inside the parenthesis) for different photon numbers and orientation class sizes.

We see that even in this smaller data set where almost half of the orientations have a nearby orientation, the success of the classification is still almost 50%. Note that with a smaller data set, the expected success rate by a mere guess is bigger. In the case we have described above, it is about 10/1000 = 1% and so the proposed classification scheme still gives a much better result.

7. Discussion and summary

In this paper we have shown a method to compress the 2-dimensional diffraction patterns into 1-dimensional distributions, namely, radial and azimuthal distributions. With the compressed signals, the variance is greatly reduced. Those signals can be utilized as a first step in processing the data from coherent X-ray experiments. As an example, we have used a simple distance measure computed from the radial and azimuthal distributions to compare the test frame and the reference frames. As the radial distributions are invariant to in-plane rotations, we can use the radial distribution to find the in-plane rotations. The data is divided into classes of (α, β). With the proposed classification scheme we can achieve 79% accuracy even with current incident photon number. With a protein that is bigger than used in this paper or possesses additional symmetries, a good classification may be achieved with an even lower incident photon level. We have performed two simulations to test the performance of the proposed classification scheme. The performance for a real situation should be between these two cases.

With the degree of in-plane rotation identified, we can rotate one pattern to overlap with the corresponding reference pattern. This provides us an opportunity to be able to average over more patterns improving the signal-to-noise level further compared to averaging over only patterns with the same orientations. The method is non-iterative in nature, which makes it a good candidate for data preprocessing.

Another interesting aspect of the suggested compression method is that the radial distribution is invariant to in-plane rotations. This property could be used in conjunction with the manifold method described in [22]. A hyper-space constructed from the radial distribution shall exhibit a 2-manifold depending only on α and β. With fewer parameters to estimate, searching the manifold will be easier, which could potentially lower the required photon level even more. Further as the variance in the radial distributions is much smaller than those in the individual pixels, the searching of the 2-manifold will be more robust to noise, thus offering us a possibility to work with even fewer photons per pixel.

As far as atomic resolution is concerned, it is associated with large-angle scattering data, which lies on the Ewald sphere. As the signal-to-noise ratio is extremely low in these regions, averaging diffraction patterns within the same orientation class is necessary to obtain meaningful large-angle-scattering data. If the size of the orientation class is large, information will be lost during averaging. For a biomolecule of size a, the maximum displacement of the same atom within one orientation class of size δθ can be estimated as

In this paper, the biomolecules are assumed to be rigid during the injection and measuring time. Other than their orientations and relative positions in the X-ray beam, the biomolecules are also assumed to be identical, i.e. no shot-to-shot variation of the structure or conformation is considered. Different conformations of the molecules will result in different diffraction patterns. The small-angle scattering data and the resulting radial and azimuthal distributions will encode information on different confirmations [28]. The classification scheme proposed in this paper will hence distinguish different classes to different conformations. Whether a certain class belongs to an out-of-plane orientation of the reference molecule in the same conformation or an orientation class of the molecule in a different conformation can be addressed in the next step of the classification. One possible method for such purpose is the manifold method [22]. A general heterogeneous sample may not be solely classified by means of the low-angle scattering data, large-angle-scattering may have to be additionally included in the classification algorithm. This is a research direction we are currently pursuing.

A common technique for solving the orientation problem in cryo-EM is the common line approach [24–26]. Because of the extremely low signal-to-noise level, this method is not applicable in current x-ray diffraction experiments. Current classification work on X-ray diffraction data [21, 22] is primarily focused on the direct 2D diffraction patterns. We propose to compress the diffraction data by radial and azimuthal integrals, so that the orientation information can still be extracted, but the variance is greatly reduced. As we have shown in Section 2, we can use the simple difference measure as described in Section 4 and 5 for classification of in-plane rotations. We can also use more elaborate methods for a complete orientation classification, such as [21, 22], with the advantage of improved SNR.

In summary, the radial and azimuthal distributions of scattering patterns from single biomolecules has shown a great success in classifying in-plane rotations by using a simple distance measure. With more sophisticated algorithms such as the manifold method, we may be able to perform a complete orientation classification and produce the expected diffraction patterns with enough sensitivity to obtain atomic resolution.