4. Best-alignment algorithms

Taking Eq. (2) as the signal model and further assuming that the noise is Gaussian and white, the likelihood and the log-likelihood of the set of observed images f gIi given the structure and the nuisance variables are:

Likelihood:

$$L(I\_i \mid \mathbf{S}, \mathcal{N}, T) = \prod\_{i=1}^{N} \frac{1}{(2\pi)^{A/2} \sigma^4} \exp\left(-\frac{||R\_{\theta, t\_i} I\_i - C\_i P\_{\mathfrak{n}\_i}(\mathbf{S})||^2}{2\sigma^2}\right) \tag{4}$$

$$L(I\_i \mid \mathbb{S}, \mathcal{N}, T) = \prod\_{i=1}^{N} L(I\_i \mid \mathbb{S}, \mathfrak{n}\_i, \theta\_i, t\_i) \tag{5}$$

Having aligned every image, line 5 of the algorithm updates the structure. Because the log-likelihood is quadratic in S, this step has a closed-form solution:

> X N

!

: (8)

–10<sup>4</sup>

i¼1 P∗ nk i CiRθ<sup>k</sup> i ,t k i Ii

is the back-projection operator defined earlier, and we use the fact that

Thus, best-alignment reconstruction iteratively aligns images to the projections of the current estimate of the structure and updates the estimate of the structure

The algorithm, as I have described above, is computationally far too expensive to implement directly. To see why, consider some numbers: a typical high-resolution

number of projection directions (vertices of the spherical grid), and the order of 10<sup>2</sup> rotations and translations per image. A naïve implementation would require Num.

cryo-EM reconstruction uses the order of 104–10<sup>5</sup> images, the order of 10<sup>3</sup>

<sup>S</sup><sup>k</sup>þ<sup>1</sup> <sup>¼</sup> <sup>X</sup>

Algorithm 1: Simple Best-Alignment Reconstruction

DOI: http://dx.doi.org/10.5772/intechopen.90099

the CTF operator Ci is self-adjoint.

Here, P<sup>∗</sup> nk i

2: <sup>k</sup> 1, <sup>S</sup><sup>k</sup> <sup>S</sup>

4: <sup>N</sup> <sup>k</sup>

6: k k þ 1 7: end while 8: return S<sup>k</sup> 9: end procedure

3: while S<sup>k</sup> not converged do

5: Skþ<sup>1</sup> arg max SI I f gj <sup>i</sup> <sup>S</sup>, <sup>N</sup> <sup>k</sup>

Figure 5.

49

Best-alignment reconstruction.

N

, <sup>T</sup><sup>k</sup> arg max <sup>N</sup> ,TI I f gj <sup>i</sup> Sk, <sup>N</sup> , <sup>T</sup> � �

!�<sup>1</sup>

1: procedure SIMPLE-BA f g Ii ð Þ , S f g Ii = images, S = initial structure, k = iteration index

A Gentle Introduction to Cryo-EM Single-Particle Reconstruction Algorithms

, T<sup>k</sup> � �

i¼1 P∗ nk i C2 <sup>i</sup> Pn<sup>k</sup> i

using the aligned images. Figure 5 illustrates the algorithm.

Log-likelihood:

$$L(\{I\_{i}\} \mid \mathcal{S}, \mathcal{N}, T) = -\sum\_{i=1}^{N} \frac{\left\| R\_{\theta\_{i}, t\_{i}} I\_{i} - C\_{i} P\_{\mathfrak{n}\_{i}}(\mathcal{S}) \right\|^{2}}{2\sigma^{2}} + const. \tag{6}$$

where σ<sup>2</sup> is the variance of the noise, which (for simplicity) I will assume is known and the same for every image; A is the number of pixels in an image; and Const: is a constant independent of S, N , and T. Also, by slight abuse of notation, I have expressed each term in the product on the right-hand side of Eq. (4) as L Ið <sup>i</sup>j S, ni, θi,tiÞ in Eq. (5).

The maximum likelihood estimates of S, N, and T are available as:

$$\hat{S}, \hat{\mathcal{N}}, \hat{T} = \arg\max\_{\mathcal{S}, \mathcal{N}, T} L(\{I\_i\} \mid \mathcal{S}, \mathcal{N}, T). \tag{7}$$

In best-alignment algorithms, the maximization in Eq. (7) is carried out numerically, by first initializing S and then alternately maximizing over (N , T) and S in several iterations. See Algorithm 1 (the superscript k is the iteration index).

The maximization with respect to N , T is carried out by some form of exhaustive search over a discrete grid. That is, a spherical grid is created over the unit sphere, and the maximization with respect to each ni is restricted to a search over the vertices of this grid. Similarly, the angle interval [0, 2π) is covered by a 1D grid, and the translation square is covered by a 2D grid, and the maximization with respect to each θ<sup>i</sup> and ti are restricted to the vertices of these grids.

The maximization in line 4 of the algorithm can be decomposed into independent maximizations with respect to (n1, θ1, t1), (n2, θ2, t2), … because the righthand side of Eq. (6) is a sum of terms, each term depending only on one tuple (ni, θi, ti). The maximization of a single term is carried out by minimizing R<sup>θ</sup>i,ti Ii � CiPni <sup>S</sup><sup>k</sup> � � � � � � 2 by exhaustive search with (ni, θi, ti) restricted to the vertices of the grids mentioned above. By minimizing (ni, θi, ti), the image Ii is best aligned to a projection of S<sup>k</sup> , in the sense that R<sup>θ</sup>i,ti Ii � CiPni Sk � � � � � � 2 is minimized (hence the name of this algorithm).

A Gentle Introduction to Cryo-EM Single-Particle Reconstruction Algorithms DOI: http://dx.doi.org/10.5772/intechopen.90099


Having aligned every image, line 5 of the algorithm updates the structure. Because the log-likelihood is quadratic in S, this step has a closed-form solution:

$$\mathbf{S}^{k+1} = \left(\sum\_{i=1}^{N} P\_{\mathfrak{n}\_i^k}^\* \mathbf{C}\_i^2 P\_{\mathfrak{n}\_i^k}\right)^{-1} \left(\sum\_{i=1}^{N} P\_{\mathfrak{n}\_i^k}^\* \mathbf{C}\_i \mathbf{R}\_{\theta\_i^k, \mathfrak{k}\_i^k} I\_i\right). \tag{8}$$

Here, P<sup>∗</sup> nk i is the back-projection operator defined earlier, and we use the fact that the CTF operator Ci is self-adjoint.

Thus, best-alignment reconstruction iteratively aligns images to the projections of the current estimate of the structure and updates the estimate of the structure using the aligned images. Figure 5 illustrates the algorithm.

The algorithm, as I have described above, is computationally far too expensive to implement directly. To see why, consider some numbers: a typical high-resolution cryo-EM reconstruction uses the order of 104–10<sup>5</sup> images, the order of 10<sup>3</sup> –10<sup>4</sup> number of projection directions (vertices of the spherical grid), and the order of 10<sup>2</sup> rotations and translations per image. A naïve implementation would require Num.

Figure 5. Best-alignment reconstruction.

estimate the structure as well as the nuisance variables. I will call these algorithms (for reasons that will become clear below) best-alignment algorithms. On the other hand, there are algorithms that only determine the structure. These are based on the

Taking Eq. (2) as the signal model and further assuming that the noise is Gaussian and white, the likelihood and the log-likelihood of the set of observed images

N

i¼1

Rθi,ti

In best-alignment algorithms, the maximization in Eq. (7) is carried out numerically, by first initializing S and then alternately maximizing over (N , T) and S in several iterations. See Algorithm 1 (the superscript k is the iteration index).

The maximization with respect to N , T is carried out by some form of exhaustive search over a discrete grid. That is, a spherical grid is created over the unit sphere, and the maximization with respect to each ni is restricted to a search over the vertices of this grid. Similarly, the angle interval [0, 2π) is covered by a 1D grid, and the translation square is covered by a 2D grid, and the maximization with

The maximization in line 4 of the algorithm can be decomposed into independent maximizations with respect to (n1, θ1, t1), (n2, θ2, t2), … because the righthand side of Eq. (6) is a sum of terms, each term depending only on one tuple (ni, θi, ti). The maximization of a single term is carried out by minimizing

of the grids mentioned above. By minimizing (ni, θi, ti), the image Ii is best aligned

where σ<sup>2</sup> is the variance of the noise, which (for simplicity) I will assume is known and the same for every image; A is the number of pixels in an image; and Const: is a constant independent of S, N , and T. Also, by slight abuse of notation, I have expressed each term in the product on the right-hand side of Eq. (4) as

Ii � CiPni k k ð Þ<sup>S</sup> <sup>2</sup>

<sup>S</sup>, <sup>N</sup>b, <sup>T</sup>^ <sup>¼</sup> arg max <sup>S</sup>, <sup>N</sup> ,TL I ðf g<sup>i</sup> <sup>j</sup> <sup>S</sup>, <sup>N</sup> , <sup>T</sup>Þ: (7)

by exhaustive search with (ni, θi, ti) restricted to the vertices

� 2

is minimized (hence the

Ii � CiPni Sk � � � � �

<sup>σ</sup><sup>A</sup> exp � <sup>R</sup>θi,ti k k Ii � CiPnið Þ<sup>S</sup> <sup>2</sup>

2σ<sup>2</sup> !

L Ið <sup>i</sup>j S, ni, θi, tiÞ (5)

<sup>2</sup>σ<sup>2</sup> <sup>þ</sup> Const: (6)

(4)

expectation-maximization algorithm.

4. Best-alignment algorithms

Likelihood:

Log-likelihood:

L Ið <sup>i</sup>j S, ni, θi,tiÞ in Eq. (5).

R<sup>θ</sup>i,ti

48

Ii � CiPni <sup>S</sup><sup>k</sup> � � � � �

to a projection of S<sup>k</sup>

name of this algorithm).

� 2

f gIi given the structure and the nuisance variables are:

Technology, Science and Culture - A Global Vision, Volume II

N

1 ð Þ <sup>2</sup><sup>π</sup> <sup>A</sup>=<sup>2</sup>

N

i¼1

The maximum likelihood estimates of S, N, and T are available as:

respect to each θ<sup>i</sup> and ti are restricted to the vertices of these grids.

, in the sense that R<sup>θ</sup>i,ti

L I<sup>ð</sup> <sup>i</sup><sup>j</sup> <sup>S</sup>, <sup>N</sup> , <sup>T</sup>Þ ¼ <sup>Y</sup>

i¼1

L I ðf g<sup>i</sup> <sup>j</sup> <sup>S</sup>, <sup>N</sup> , <sup>T</sup>޼�<sup>X</sup>

^

L I<sup>ð</sup> <sup>i</sup><sup>j</sup> <sup>S</sup>, <sup>N</sup> , <sup>T</sup>Þ ¼ <sup>Y</sup>

images Num. projections Num. rotations and translations (>104 103 <sup>10</sup><sup>2</sup> ) evaluations of R<sup>θ</sup>i,ti Ii CiPni <sup>S</sup><sup>k</sup> 2 for every execution of line 4 of the algorithm. Moreover, calculating R<sup>θ</sup>i,ti Ii CiPni <sup>S</sup><sup>k</sup> 2 requires one image rotation plus translation. This is computationally expensive. A number of "tricks" have been developed to keep the computational cost manageable. I will discuss these below.

particularly sophisticated form of gridding incorporating a variant of the Pipe

A Gentle Introduction to Cryo-EM Single-Particle Reconstruction Algorithms

different images and also to account for nonwhite noise.

� 2

mizing the correlation coefficient between R<sup>θ</sup>i,ti

DOI: http://dx.doi.org/10.5772/intechopen.90099

5. Expectation-maximization algorithms

Ii � CiPni <sup>S</sup><sup>k</sup> � � � � �

estimated. Only S is estimated.

included for S, if necessary.

Algorithm 2: Simple EM Reconstruction

3: while S<sup>k</sup> not converged do

<sup>Γ</sup><sup>k</sup>þ<sup>1</sup> L Ið j<sup>i</sup> Sk, <sup>n</sup>i, <sup>θ</sup>i, ti<sup>Þ</sup> <sup>p</sup> <sup>n</sup>i, <sup>θ</sup>i, tijSk ð Þ

L I<sup>ð</sup> <sup>i</sup>jSk, <sup>n</sup>i, <sup>θ</sup>i, ti<sup>Þ</sup> <sup>p</sup> <sup>n</sup>i, <sup>θ</sup>i, tijSk ð Þ <sup>d</sup>n<sup>i</sup> <sup>d</sup>θidti

2: <sup>k</sup> 1, Sk <sup>S</sup>

Ð

51

Algorithm 2.

The L I<sup>ð</sup> <sup>i</sup>jS<sup>k</sup>

The simplifications I made to describe best-alignment reconstruction are easy to discard. It is straightforward to allow for unknown and unequal noise variances for

Some packages (e.g., SPIDER, EMAN) carry out the alignment step by maxi-

In classical statistics, the estimate of any set of parameters improves as more data are added, provided that the number of parameters stays fixed. The number of parameters in best-alignment reconstructions does not stay fixed; the number of parameters in ð Þ N , T increases linearly with the number of images. This can potentially limit the asymptotic accuracy of best-alignment reconstructions. Expectationmaximization algorithms attempt to overcome this limitation by treating ð Þ N , T as latent variables, i.e., as variables that influence the likelihood, but which are not

The theory of the expectation-maximization algorithm can be found in many textbooks, e.g., [16]. The EM algorithm works iteratively, where in each step the conditional mean of the latent variables is used to update the parameter estimate. For the cryo-EM reconstruction problem, variables N , T are taken to be the latent variables, and S is taken to be the parameter to be estimated. A prior can also be

RELION [8, 9] uses the EM algorithm with a Gaussian prior on the amplitude of the Fourier coefficients of S. The resulting algorithm is rather complex, and instead of discussing all of the details of the algorithm, I will discuss a simplified version in

, <sup>n</sup>i, <sup>θ</sup>i, tiÞterm in line 4 comes from Eq. (5). The term <sup>p</sup> <sup>n</sup>i, <sup>θ</sup>i, tijS<sup>k</sup> � � in line

Line 4 in Algorithm 2 calculates the conditional probability of the alignment parameters ni, θi, ti given the image Ii and the current estimate of the structure S<sup>k</sup>

be set to a uniform probability density. The denominator on the right-hand side of the assignment in line 4 is the normalizing constant, which makes <sup>Γ</sup><sup>k</sup>þ<sup>1</sup> <sup>n</sup>i,θi,ti a probability density. Also, note that <sup>Γ</sup><sup>k</sup>þ<sup>1</sup> <sup>n</sup>i,θi,ti is a function of <sup>n</sup>i, <sup>θ</sup>i, ti and is calculated for all values of ni, θi, ti on the vertices of the spherical, angular, and translation grids.

4: For each image Ii calculate the following conditional probability for all values of ni, θi, ti

4 is the prior probability of the alignment parameters given S<sup>k</sup>

1: procedure SIMPLE-EM f g Ii ð Þ , S f g Ii = images, S = initial structure, Initialization

5: Use the conditional probabilities to update the structure as:

Ii and CiPni Sk � � rather than mini-

.

, and typically this can

. This corresponds to using the signal model of Eq. (3)

and Menon method [15].

4.3 Comments

mizing R<sup>θ</sup>i,ti

in the alignment step.

Another problem with the simple algorithm is the structure update step of Eq. (8). The matrix representations of the operators in this equation are far too large to compute with, and, in practice, tricks are also used to simplify this calculation. I will discuss these below as well.

Different cryo-EM reconstruction packages differ in the tricks that the packages use to speed up computation. The combination of various tricks can occasionally be overwhelming to understand, but it helps to keep in mind that the underlying algorithm is just Algorithm 1 or a minor variation of it. SPIDER, EMAN, FREALIGN, and cryoSPARC implement versions of best-alignment reconstruction.

### 4.1 Speeding up alignment

Approaches to speed up the alignment step are:


#### 4.2 Speeding up structure update

The structure update of Eq. (8) has a simpler form in the Fourier domain. The Fourier slice theorem [4, 12] suggests that the projection and back-projection operators simplify to 2D slice extraction and slice insertion operators in the Fourier domain. The CTF filter operator also reduces to a point-wise multiplication by the CTF. With these simplifications, Eq. (8) becomes tractable.

Inserting or extracting a 2D slice from a 3D volume requires careful numerical interpolation. A method called gridding is used for this [13, 14]. All cryo-EM reconstruction packages use some form of gridding. RELION uses a

particularly sophisticated form of gridding incorporating a variant of the Pipe and Menon method [15].
