The emergence of super-bright light sources - X-ray free electron lasers(XFELs) combined with Single Particle Imaging(SPI) method, makes it possible to obtain nanometer resolution 3D structure of biological particles such as proteins or viruses without needing to freeze them. SPI relies on the \di raction before destruction" principle, meaning that each sample only produces a single di raction image before being destroyed by an X-ray pulse. The orientation of the particle in the beam is random for each shot. This gives rise to the problem of orientation recovery, in which an array of 2D di raction images has to be combined into a single 3D image, necessary for the reconstruction of 3D structure of the studied particle. The orientation recovery problem is most commonly solved by the EMC algorithm, which is the most computationally expensive part of data analysis for SPI experiments. In this work we introduce several modi cations to the EMC algorithm aimed at improving the quality of reconstruction and/or increasing the algorithm's speed of convergence. We analyse the e ectiveness of these modi cations using simulated di raction data.
The emergence of fourth generation light sources - X-ray Free Electron Lasers (XFEL), opens up new possibilities in many elds of science. Compared to third generation synchrotron light sources, XFELs produce extremely bright (more than 1012 photons) and extremely short (tens of femtoseconds) x-ray pulses [6]. These properties provide a new way to study the 3D structure of bioparticles such as proteins and viruses - Single Particle Imaging (SPI) [8]. Compared to other methods such as X-Ray crystallography or cryogenic electron microscopy, SPI has the advantage of being able to study non-crystalline bioparticles in their natural state (suspended in water).
SPI relies on "di raction before destruction"[
This gives rise to orientation recovery problem: SPI produces a number of at di raction images of the particle in an unknown orientation. In order to reconstruct the 3D structure of said particle these images need to be combined to produce a single 3D di raction density. The established way to solve this problem is by using the EMC algorithm [5].
In this work we introduce three modi cations of orientation recovery algorithm EMC, which aim to improve quality of reconstruction and to increase the speed of convergence. We examine their e ectiveness using simulated di raction data. 2 2.1
In orientation recovery problem we have an array of Mdata di raction patterns scattered by identical particles in unknown orientations. Each such pattern is a spherical slice of a single 3D density in reciprocal space, and the goal of an orientation recovery algorithm is to reconstruct this 3D density W by nding to which orientation each di raction pattern belongs to. This problem is most commonly solved by the EMC algorithm.
EMC algorithm starts with input data (di raction images K) and an initial approximation of 3D di raction density W 0. Then this initial approximation is iteratively improved until it converges to the nal value of W T . Each iteration of the algorithm consists of three steps which give the name to the algorithm: Expand, Maximize and Compress.
During Expansion step the current approximation W, which is typically represented by its values on regular 3D grid, is converted to tomographic representation Wij - values of W in points corresponding to i-th pixel of the di raction image with scattering particle in j-th orientation. In order to do this we sample the 3D rotation group SO(3) by a nite number of "evenly spaced" rotations Rj (j = 1 : : : Mrot). And for each pixel of the detector we calculate a corresponding point in reciprocal space qi (i = 1 : : : Mpix). For typical at detectors all these points will be lying on Ewald's sphere. After calculating rotations Rj and points qi we can de ne Wij as W (Rj qi), which is calculated via linear interpolation.
The Maximization step updates current approximation of 3D di raction
density Wij ! Wi0j based on maximizing log-likelihood function Q(W 0). This step is
equivalent to one iteration of EM algorithm [
Mdata
P k=1 Mdata
P k=1
PjkDik
Pjk :
Finally, Compression step converts Wi0j back onto regular 3D grid by reversing the interpolation procedure used in the expansion step. In essence, EMC algorithm is equivalent to EM where after each iteration we perform combination of compression and expansion steps, which enforce extra constraints on current model Wij . Whereas for EM all Wij are treated as independent variables, in reality they are derived from values on 3D intensity grid, and when the distance between points Rj qi and Rj0 qi0 is small, values Wij and Wi0j0 are not independent.
In this work we propose three modi cations to the EMC algorithm: Incremental EMC modi es the maximization procedure, instead of performing a single update of 3D intensity W using data from all images, only one image is used on each iteration. For randomly selected image Dk probabilities Pjk are calculated and then the current values of Wij are immediately updated: Wi0j =
Mdata ( P k=1
Mdata ( P k=1
: Such an update can be performed in O(1) time if the numerator and denominator of Wij are saved separately, as well as all the values of Pjokld [7]. This way Mdata of such updates can be performed in the same time as the single maximization step of EMC algorithm. This incremental M-step is followed by compression and expansion steps same as in regular EMC algorithm. However, expansion steps outputs only values of Wij , and in order to use out incremental maximization step again we rst need to perform a normal M-step to obtain separate values of numerators PMdata PjkDik, denominators PMdata Pjk and probabilities Pjokld.
k=1 k=1 E ectively, this modi cation of EMC alternates its M-steps between normal and incremental, and only performs compression and expansion steps after every other iteration.
Stepwise EMC uses batches of images in its maximization step. First a new 3D di raction density Wij is calculated using only a subset of images D 2 D. Then we take weighted average between current desity Wij and Wij as the result of maximization:
Wi0j = (1 (2 + t) )Wij + (2 + t) where t is the number of iteration and 0:5 < 1 is a coe cient that ensures that the weight before newly calculated density W exponentially decreases, thus guaranteeing convergence of the algorithm [10]. This maximization step is then performed several more times with di erent subsets of images, until all Mdata images have been used. Then follow normal compression and expansion steps. Adaptive EMC is the nal proposed modi cation which does not modify any of the three steps of EMC, and instead it only takes e ect after the M-step, when the new values of Wi0j are calculated. Instead of taking the point W 0 itself, which maximizes the expected log-likelihood function Q(W 0), we try to go further in the same direction:
^ Wij = Wij + (Wi0j
Wij ); where coe cient 1 is adaptively changed based on log-likelihood Q(W^ ). If Q(W^ ) Q(W ) then is increased and the W^ is accepted as the result of current iteration, otherwise is reset to 1 and W 0 is taken as a result instead. This approach guarantees that log-likelihood never decreases and the algorithm thus converges [9]. 2.2
To evaluate the performance of our modi cations and to compare them to regular EMC we tested them using simulated di raction data.
In the rst test we used a binary contrast torus in reciprocal space as our object. Using such a simple model allows us to successfully perform orientation recovery in a short time and without tuning of any reconstruction parameters. This model allows us to easily compare di erent algorithms, but it has a drawback of being very far removed from actual experimental data. For that reason we performed a second more life-like test.
In the second test we used di raction patterns of keyhole limpet hemocyanin
[4] simulated by Dragon y [
In both tests for each algorithm we performed the reconstruction 5 times, using a random initial approximation of 3D di raction density W 0. Then we compared the output of each iteration W t with the initial density Wtrue, that was used to generate the di raction images. We used root mean square di erence (RMS di erence) between these two 3D di raction densities as our metric of quality of reconstruction. 3
In the rst test we used 1000 images with 4900 pixels in each image. Number of possible rotations Mrot was increasing every 10 iterations. Starting from 420 0.24 0.22 possible orientations on the rst ten iterations and up to 10860 orientations for iterations 41 through 50. The results of this test are presented on Fig. 1.
One can see that out of three proposed modi cations, only incremental EMC demonstrated better results than unmodi ed algorithm. Stepwise EMC proved to be too dependent on the random selection of the rst batch of images used in reconstruction. Adaptive EMC didn't provide any boost to the speed of convergence, due to the fact that the coe cient never became greater than 1 without decreasing the expected log-likelihood.
For these reasons stepwise and adaptive modi cations were ruled out as ine ective, and for the second test only Incremental EMC was evaluated. In this test we used 5000 images with 10000 pixels and 10860 possible orientations. The results of this test are presented on Fig. 2.
Both algorithms demonstrated similar results, however our chosen metric of quality doesn't perform very well on this test. The di erence between 0-th iteration which is just random noise and the nal output of the algorithm is quite small, when looking only at the RMS di erence between the reconstructed and initial 3D di raction density. Due to this we cannot conclusively evaluate the performance of incremental modi cation of EMC algorithm. 4
In this work we proposed, developed, and tested three modi cations to the EMC algorithm, used to solve orientation recovery problem in SPI experiments. After the rst preliminary test, adaptive and stepwise modi cations have shown worse results than unmodi ed algorithm and were thus ruled out as non-viable. For the more promising incremental modi cation we performed a second set of tests with more life-like input data. In these tests both incremental and regular EMC 0 5 10
15 Iteration 20 25 30 demonstrated similar results, but due to instability of reconstruction process and di culty of the evaluation of reconstruction quality, we can not de nitively say that one algorithm is better than the other. Additional tests may be required to establish that. 5
This research was supported by Russian Science Foundation (project No. 18-4106001) and Helmholtz Associations Initiative Networking Fund (HRSF-0002).