PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Struct Biol. Author manuscript; available in PMC 2010 September 1.
Published in final edited form as:
PMCID: PMC2803348
NIHMSID: NIHMS119275

Ab initio maximum likelihood reconstruction from cryo electron microscopy images of an infectious virion of the tailed bacteriophage P22 and maximum likelihood versions of Fourier Shell Correlation appropriate for measuring resolution of spherical or cylindrical objects

Abstract

A maximum likelihood reconstruction method for an asymmetric reconstruction of the infectious P22 bacteriophage virion is described and demonstrated on a subset of the images used in [Lander, G.C., Tang, L., Casjens, S.R., Gilcrease, E.B., Prevelige, P., Poliakov, A., Potter, C.S., Carragher, B., Johnson, J.E., 2006. The structure of an infectious P22 virion shows the signal for headful DNA packaging. Science 312(5781), 1791–1795]. The method makes no assumptions at any stage regarding the structure of the phage tail or the relative rotational orientation of the phage tail and capsid but rather the structure and the rotation angle are determined as a part of the analysis. A statistical method for determining resolution consistent with maximum likelihood principles based on ideas for cylinders analogous to the ideas for spheres that are embedded in the Fourier Shell Correlation method is described and demonstrated on the P22 reconstruction. With a correlation threshold of .95, the resolution in the tail measured radially is greater than 0:0301 Å−1 (33.3 Å) and measured axially is greater than 0:0142 Å−1 (70.6 Å) both with probability p = 0:02.

Keywords: Bacteriophage P22, Infectious particle, 3-D reconstruction, Cryo electron microscopy

1. Introduction

Motivated by recent reconstructions from cryo electron microscopy (cryo EM) images of tailed bacteriophages epsilon15 (Jiang et al., 2006) and P22 (Lander et al., 2006), alternative maximum likelihood reconstruction and resolution calculation methodologies are described and demonstrated on the same P22 images used in Lander et al. (2006). Maximum likelihood dates back to the early 1900s (Lehmann and Casella, 1998, Section 10.1, p. 515) but continues as an important method for deriving statistical estimators in structural biology (e.g., Blanc et al. (2004) and McCoy et al. (2005) in crystallography and Scheres et al. (2007) and Singh et al. (2004) in cryo EM). In comparison with the reconstruction method described in Lander et al. (2006), the chief advantage of the reconstruction method described in the present paper is its ab initio character which manifests itself in two main differences: First, it is not necessary to have a 3-D structure of the tail machine before determining the 3-D structure of the tailed phage. Second, although the 6-fold symmetric tail machine is attached at a 5-fold symmetry axis of the capsid, it is not necessary to specify the rotational position of the tail machine relative to the symmetry axes of the capsid; Instead, all possible rotations, a range of 12 degrees (please see Supplemental material, Section L), are considered by the reconstruction method and the best is selected. Similar to the results of Lander et al. (2006), the portal end of the tail shows approximate 12-fold rotational symmetry even though no such symmetry was imposed.

The approach of this paper can be applied to any tailed bacteriophage for which an icosahedrally symmetric structure can be determined. If the tail is long and flexible, only the proximal part of the tail will be resolved in the 3-D reconstruction. More generally, the approach can probably be applied to viruses where the infection process results in a distinguished attachment site on the surface of the virus which replaces the role of the tail. Such problems are of current interest, e.g., in the case of polio virus, Bubeck et al. (2005) and Zhang et al. (2008). Finally, the methods used in this paper show how maximum likelihood approaches can be used for complicated structures by assembling the structure out of parts and estimating parameters that determine the structure of each part and the relative locations and orientations of the parts. The computations required in this approach are moderately extensive, e.g., best done on a set of tens of PCs. However, the computations are much less extensive than would be required for a statistical ab initio 3-D reconstruction of the infectious bacteriophage (i.e., capsid, tail, and genome) since such a reconstruction would be of a particle without any symmetry. Because the approach of this paper is based on combining an icosahedrally symmetric reconstruction with an ξ-fold symmetric tail reconstruction, not all possible distortions of the capsid to accommodate the tail can be represented. However, to the extent that the distortions of the capsid which allow the joining of the tail are ξ-fold symmetric, then the proximal part of the tail reconstruction will include those distortions so that the sum of the icosahedrally symmetric reconstruction and the ξ-fold symmetric tail reconstruction will accurately reflect the structure of the tailed bacteriophage.

Relative to standard Fourier Shell Correlation (FSC) calculations (please see van Heel and Schatz (2005) for new ideas on setting FSC thresholds and an extensive bibliography of FSC investigations containing 26 entries dating back to initial papers such as Frank et al. (1981) and Saxton et al. (1982)), the resolution methodology described in the present paper has several potentially attractive features: (1) It is linked to the maximum likelihood criteria used to determine the 3-D reconstruction algorithm. (2) It is possible to measure resolution independent of orientation, as is appropriate for spherical objects, or with respect to translations along and rotations around a specific axis, as may be natural for a cylindrical object such as the tail of a phage. (3) It is not necessary to perform two reconstruction calculations each with half of the entire data set. (4) It provides a probability of correctness, i.e., the answer is of the form that resolution is greater than a particular number with a certain probability.

The reconstruction method has three phases: (1) Use a standard cryo EM reconstruction algorithm to compute an icosahedrally symmetric reconstruction of the tailed phage. (2) Use the reconstruction of Phase (1) to determine origin location and projection orientation for each image by quadratic correlation. The projection orientation is only determined up to one of the 60 rotations of the icosahedral group, since the reconstruction of Phase (1) has icosahedral symmetry. (3) Use a mathematical description of the tail, the capsid reconstruction of Phase (1), the origin locations and projection orientations (up to a rotation from the icosahedral group) of Phase (2), and the maximum likelihood criteria to determine a 3-D reconstruction of the entire tailed phage. While the details of Phases (1) and (2) are described in the numerical results (Section 3), the major portion of the reconstruction part of the present paper concerns Phase (3) (Section 2).

The resolution method has two phases: (1) Compute the Hessian of the log likelihood at the maximum likelihood parameter estimates, i.e., compute the matrix of mixed second-order partial derivatives of the log likelihood with respect to the parameters evaluated at the particular vector of parameters that maximizes the likelihood. As is described in Section 4.1, the parameter estimation error, i.e., the difference between the true value of the parameter vector and the value determined by the maximum likelihood estimator, is approximately Gaussian in distribution and the negative of the inverse of this Hessian is approximately the parameter estimation error covariance matrix. (2) As is described in Section 4.4, use a Monte Carlo procedure to compute many FSC curves where the structures compared by FSC are drawn at random from the multivariate Gaussian probability density function (pdf) determined in Phase (1). From this ensemble of FSC curves, it is possible to compute a histogram which approximates the pdf for the resolution at which FSC first falls below any threshold, where the threshold might depend on the magnitude of the reciprocal space position, as is described in van Heel and Schatz (2005). From this histogram it is possible to compute the probability that the resolution exceeds a particular value. For cylindrical objects, two alternatives to FSC are described which are appropriate for measuring axial and rotational resolution, respectively.

2. Reconstruction method

For further details of the reconstruction method, please see Prust (2006).

2.1. Mathematical model for the phage capsid and tail

Real space coordinates are denoted by x with rectangular, cylindrical, and spherical coordinates denoted by (x, y, z)T, (r, ϕ, z), and (|x|, θ, ϕ), respectively. Correspondingly, reciprocal space coordinates are denoted by k, (kx, ky, kz), (kr, ϕ′, kz), and (k, θ′, ϕ′). The capsid electron scattering intensity, denoted by ρc (x), is described in a coordinate system in which the rotational symmetry axes intersect at the origin of the coordinate system, the z axis is a 5-fold symmetry axis, and the quadrant of the x–z plane that has x > 0 and z > 0 includes one of the five 2-fold symmetry axes closest to the positive z axis (Yin et al., 2003; Altmann, 1957; Laporte, 1948). The tail electron scattering intensity, denoted by ρt(x), is described in a coordinate system in which the long axis of the tail is the z axis of the coordinate system and the end of the tail that attaches to the capsid is the end at more positive z value. The entire particle is described in the same coordinate system as the capsid. The electron scattering intensity of the complete particle, denoted by ρ(x), is therefore

ρ(x)=ρc(x)+ρt(xδ)
(1)

where

δ=(0,0,zt)T.
(2)

It is not necessary to consider rotation of the tail when attaching the tail to the capsid because the mathematics to be used in the sequel can represent the tail in any rotation around the z axis. It will be convenient to assume that the tail is length-centered in the tail coordinates, i.e., ρt(x) is nonzero only in a symmetric region −z0/2 ≤ z ≤ +z0/2 (z0 > 0) in which case zt ≤ 0.

Because of the cylindrical shape of the tail and the rotational symmetry of the tail around its long axis, a cylindrical coordinate system is used. Because the tail is by definition periodic with period 2π in ϕ, ρt(x) can be expressed as a Fourier series with radially- and axially-dependent weights denoted by cl(r, z). Assuming that the tail has a known maximum radius, denoted by R+, it follows that the radial dependence of the weights can be expanded in a sum of weighted Bessel functions in r. The axial dependence of the weights can be expanded as a sum of weighted complex exponentials in z. Assume that the tail has rotational symmetry of order ξ (ξ = 1 is permitted). It follows from Eqs. 185, 174, 180, and 92 in the Supplemental material that

ρt(r,φ,z)=l=+p=1n=+cl,p,nq(z)fn(z)hl,p(r)exp(ilξφ)
(3)

where cl,p, n are unknown complex-valued coefficients,

q(z)={1,|z|z0/20,otherwise,
(4)
fn(z)=exp(i(2π/z0)nz),
(5)
hl,p(r)={0,rR+J|lξ|(γ|lξ|,pr/R+),r<R+,
(6)

Jl(x) is the lth Bessel function of the first type, and γl,p is the pth zero of Jl(x). Since ρ(x) is real, it follows by Hermitian symmetry (see Eq. 179 in the Supplemental material) that cl,p,n has the symmetry

cl,p,n=cl,p,n.
(7)

Note that Eq. (7) implies that

[c0,p,0]=0
(8)

where I indicates the imaginary part. The reciprocal-space representation of the electron scattering intensity ρt(x) is denoted by Pt(k) and is the 3-D Fourier transform of ρt(x) which is defined by

Pt(k)=ρt(x)exp(i2πkTx)dx.
(9)

It then follows from Eqs. 202, 203, 198, and 118 in the Supplemental material that

Pt(k)=l=+p=1n=+Lt(kr,φ,kz),(l,p,n)cl,p,n
(10)

where

Lt(kr,φ,kz),(l,p,n)=Q(kzn/z0)exp(ilξ(φπ/2))Hl,p(k),
(11)
Q(kz)=z0sinc(kzz0),
(12)
Hl,p(kr)=R+2γ|lξ|,pJlξ1(γ|lξ|,p)J|lξ|(2πkrR+)(2πkrR+)2γ|lξ|,p2,
(13)

and sinc(z) = sin(πz)/(πz).

The reciprocal-space representation of the electron scattering intensity of the capsid (complete particle), i.e., of ρc(x) [ρ(x)], is denoted by Pc(k) [P(k)] and is the 3-D Fourier transform of ρc(x) [ρ(x)] where the 3-D Fourier transform is defined by Eq. (9). Eq. (1) implies that

P(k)=Pc(k)+exp(i2πkTδ)Pt(k).
(14)

Results related to the icosahedral average of the complete particle (capsid plus tail) are necessary in the approach of this paper. The icosahedral group has Ng = 60 operations each of which is a rotation that can be expressed (for a specific coordinate system) as a 3 × 3 matrix. For β [set membership] {0, …, Ng − 1}, let Sβ [set membership] R3×3 be the matrices which, since they are rotation matrices, satisfy Sβ1=SβT and det Sβ = +1. If a function f is rotated to yield a function f′ and the rotation is described by the rotation matrix R then the definition used in this paper is that f′(x) = f(R−1x). With these preliminary results, the icosahedral average of the complete particle, denoted by [rho with macron](x), is

ρ¯(x)=˙1Ngβ=0Ng1ρ(Sβ1x)
(15)
=1Ngβ=0Ng1[ρc(Sβ1x)+ρt(Sβ1xδ)]
(16)
=ρc(x)+1Ngβ=0Ng1ρt(Sβ1xδ)
(17)

since ρc(x) has icosahedral symmetry, i.e., ρc(Sβ1x)=ρc(x) for β [set membership] {0, …, Ng − 1}. The icosahedral average, [rho with macron](x), also has icosahedral symmetry, i.e., ρ¯(Sβ1x)=ρ¯(x) for β [set membership] {0, …, Ng − 1}. Let P(k) be the reciprocal space representation of [rho with macron](x). Eq. (17) implies that

P¯(k)=Pc(k)+1Ngβ=0Ng1exp(i2πkTSβδ)Pt(Sβ1k).
(18)

As is described in Supplemental material Section B.6, the mathematics used in this paper does not uniquely represent the tail since if cl,p,n represents the tail ρt(x) then exp(−ilϕ0)cl,p,n represents the same tail rotated around the z axis by the angle ϕ0 where ϕ0 [set membership] [0, 2π) is arbitrary. However, when the tail is attached to the capsid, only 5 of these rotations are equivalent for the combination of capsid and tail since only under the 5 rotations ϕ0 [set membership] {n2π/5 : n [set membership] {0, …, 4}} is the capsid unaltered since the z axis is a 5-fold symmetry axis of the capsid.

2.2. Mathematical model for the image formation process and the difference image

A standard image formation equation is used. Let σi(χ) be the ith real-space image and Σi(κ) be the corresponding reciprocal space image which is its 2-D Fourier transform defined analogously to Eq. (9) where χ [set membership] R2 and κ [set membership] R2 are the 2-D coordinate vectors in real and reciprocal space, respectively, where κ=̇|κ| and χ=̇|χ|. Let χ0,i be the offset between the location of the particle's center in the ith image and the center of the ith image. Let G(κ) be the contrast transfer function (CTF) (Baker et al., 1999, p. 873; Scherzer, 1949). Let Ri be the rotation matrix that describes the orientation of the particle in the microscope or, equivalently, the projection orientation. Then (Erickson, 1973, Eq. 11c; Yin et al., 2003, Eq. 10),

i(κ)=exp(iκTχ0,i)G(κ)P(Ri1(κT,0)T).
(19)

As is described in Section 1, the approach of this paper is based on difference images where the difference is between the experimental image and the predicted image where the prediction is based on an icosahedrally symmetric reconstruction of the complete particle. As a part of the reconstruction, estimates are made of the origin offset for each image, the particle orientation for each image, and the icosahedrally symmetric electron scattering intensity. In this paper, the following assumptions are made:

  1. The origin offset estimate, denoted by [chi]0, i, is exact, i.e., [chi]0, i = χ0, i.
  2. The estimate of the rotation matrix describing the particle orientation, denoted by Ri, is exact up to a rotation from the icosahedral group, i.e.,
    R^i=RiSβi.
    (20)
  3. The estimate of the icosahedrally averaged electron scattering intensity of the complete particle, denoted by [rho with circumflex](x), is exact, i.e., [rho with circumflex](x) = [rho with macron](x).

A predicted image is needed in order to form the difference image and the natural predicted image, denoted by [Sigma]i(k), is

^i(κ)=exp(iκTχ^0,i)G(κ)P^(R^i1(κT,0)T)
(21)
=exp(iκΤχ0,i)G(κ)P¯((RiSβi)1(κT,0)T)
(22)
=exp(iκΤχ0,i)G(κ)P¯(Sβi1(Ri1(κT,0)T))
(23)
=exp(iκΤχ0,i)G(κ)P¯(Ri1(κT,0)T)
(24)

since P(k) has icosahedral symmetry.

The ith difference image, denoted by Δi(κ;Ri), is

Δi(κ;Ri)=i(κ)^i(κ)
(25)
=exp(i2πκTχ0,i)G(κ)[P(Ri1(κT,0)T)P^(Ri1(κT,0)T)]
(26)
=exp(i2πκTχ0,i)G(κ)[exp(i2π(κT,0)Riδ)Pt(Ri1(κT,0)T)1Ngβ=0Ng1exp(i2π(κT,0)RiSβδ)Pt((RiSβ)1(κT,0)T)].
(27)

Since the origin offset is known by Assumption 1, it is natural to assume that the boxed images are shifted such that the origin offset in the shifted image is zero. In that case the factor exp(−i2πκTχ0,i) has value 1. The icosahedral averaging that creates P(k) causes the single tail to be replicated Ng = 60 times in 12 groups of 5 with one group at each 5-fold symmetry axis where each replication is at 1/60 the scattering intensity of the true tail. We refer to these replicated low-scattering-intensity tails as “ghosts”. Let Σg(κ;Ri) denote the ghost tail reciprocal-space image with definition

g(κ;Ri)=1Ngβ=0Ng1exp(i2π(κT,0)RiSβδ)Pt((RiSβ)1(κT,0)T).
(28)

Note that

g(κ;RiSβ)=g(κ;Ri)
(29)

for all β [set membership] {0, …, Ng − 1} because the {Sβ} form a multiplicative group. Using this definition, the difference image can be written as

Δi(κ;Ri)=exp(i2πκTχ0,i)G(κ)[exp(i2π(κT,0)Riδ)Pt(Ri1(κT,0)T)g(κ;Ri)]
(30)

which describes the difference image as the superposition of the true tail and the ghost tails.

2.3. Statistical model for the noisy images

The statistical model falls within the general class of models described in Doerschuk and Johnson (2000) and Yin et al. (2003) and only the main characteristics are briefly summarized here. The central feature, and the key difference from the numerical examples described in Doerschuk and Johnson (2000) and Yin et al. (2003), is that the orientations of the particles are still independent random variables but the probability density functions (pdfs) for these random variables are not identical. In particular, each particle has its pdf concentrated on the Ng = 60 icosahedrally related orientations that are the outcome of orienting the image of a nonsymmetrical particle with predicted images of an icosahedrally symmetric particle. As is described in Section 2.2, it is assumed that the difference boxed images are shifted if necessary so that the center of the capsid is projected to the center of the image. Therefore, there is no uncertainty in the location of the center of the capsid in the image and so, in the notation of Doerschuk and Johnson (2000) and Yin et al. (2003), χi,j = (0,0)T. In the calculations described in this paper, we assume that there is only one class of capsid and one class of tail. Therefore, in the notation of Doerschuk and Johnson (2000) and Yin et al. (2003), Nη = 1. This restriction could be removed. The reciprocal space image is assumed to be corrupted by additive zero-mean white Gaussian noise with known variance σ2. The variance is, in fact, estimated from the images in a preliminary calculation.

Removing the one class restriction would require a precise definition of how multiple classes occur. For instance, if the capsid has only one class but the tail has multiple classes, then the following generalization would be natural: (1) Merge all of the data to compute an icosahedrally symmetric reconstruction. (2) Use the icosahedrally symmetric reconstruction to compute difference images. (3) Use a multiclass generalization of the algorithm described in this paper, exactly following the multiclass algorithm of Doerschuk and Johnson (2000), to reconstruct multiple tail structures. Alternatively, if the capsid has multiple classes but only one possible tail exists, then the following generalization would be natural: (1′) Use the multiclass algorithm of Doerschuk and Johnson (2000), to reconstruct multiple icosahedrally symmetric structures and label each image with its estimated class. (2′) Based on the estimated class label, compute difference images using the appropriate icosahedrally symmetric reconstruction. (3′) Apply the algorithm described in this paper to the difference images to determine a reconstruction of the single class of tail. Finally, if there are multiple classes of capsid and multiple classes of tail and all possible mixtures of capsid and tail are present in the data, then a combination of these two generalizations would be necessary, in particular, (1′), (2′), and (3).

Let the ith difference image be arrayed in a vector denoted by yi. Let the unknown coefficients, cl,p,n, be arrayed in a vector denoted by c. Let the additive pixel noise for the ith difference image be arrayed in a vector denoted by wi which is, therefore, Gaussian with mean 0 and covariance Qi = σ2INy where Ny is the number of pixels in the image. From Eq. (30), the ith difference image depends linearly on Pt(k). From Eq. (10), Pt(k) depends linearly on the unknown coefficients cl,p,n (which are the elements of the vector c). Therefore, there is a matrix, which is denoted by LΔ, that relates the ith difference image to the unknown coefficients cl,p,n (which are the elements of the vector c). (In Eq. (11) the elements of the simpler matrix relating Pt(k) to c, where c has elements cl,p,n, is shown explicitly). The matrix depends on the identity of the image, i.e., on i, because it depends on the orientation of the particle, i.e., on Ri. The matrix also depends on random variables, such as which of the Ng = 60 icosahedrally related orientations is present, and the collection of such random variables for the ith image is denoted by zi. In the calculations of this paper, the only random variables on which the matrix depends are the orientations and therefore, the integrals in Eqs. (35)(37) are actually discrete sums though in more general problems the matrix could depend on additional variables such as translations in which case the integrals would not reduce to discrete sums. The conclusion of this paragraph is that there is an equation,

yi=LΔ(i,zi)c+wi,
(31)

analogous to Doerschuk and Johnson (2000, p. 1718) and Yin et al. (2003, Eq. 19, p. 31), which describes the entire imaging system.

2.4. Maximum likelihood reconstruction method and expectation-maximization algorithm

Use of the maximum likelihood estimation ideas and formulas of Doerschuk and Johnson (2000) and Yin et al. (2003) allows a reconstruction of the tail by estimating values for the unknown cl,p,n coefficients which in turn specify the tail through Eq. (3). With the value of the matrix LΔ and the pdfs for the random variables on which LΔ depends, both new in this paper for this application, the algorithm is as follows. First, pre-compute the following quantities:

ai(yi;zi)=j=1NT(i){ln[(2π)Ny/2detQi,j(zi)]+12yi,jTQi,j1(zi)yi,j}
(32)
bi(yi;zi)=j=1NT(i)LΔT(i,j,zi)Qi,j1(zi)yi,j
(33)
Di(zi)=j=1NT(i)LΔT(i,j,zi)Qi,j1(zi)LΔ(i,j,zi)
(34)

where, for the ith virion, these equations allow NT(i) tilt images, denoted by yi,j, to be processed and allow the pixel noise covariance Qi,j and LΔ to depend on both the virion index i and the tilt series index j. These quantities allow rapid evaluation of a Gaussian pdf which is the product of NT(i) independent Gaussian pdfs with means LΔc the covariances Qi,j which is the pdf needed in the expectation of the expectation-maximization algorithm (zi are the so-called nuisance parameters of the algorithm). Second, determine an initial condition for the values of the cl,p,n coefficients. In the numerical results presented in Section 3, that initial condition is random as is described in Section 3.1. Third, starting at this initial condition, iterate the following actions until the values of the cl,p,n coefficients have converged.

  1. Compute
    γi(c0,yi)=zip(yi|zi,c0)p(zi)dzi
    (35)
    βi(c0,yi)=zibi(yi;zi)p(yi|zi,c0)p(zi)dzi
    (36)
    Δi(c0,yi)=ziDi(zi)p(yi|zi,c0)p(zi)dzi
    (37)
    where c0 is the current value of the vector constructed from the cl,p,n coefficients.
  2. Combine these results to compute
    g=i=1Nv1γi(c0,yi)βi(c0,yi)
    (38)
    F=i=1Nv1γi(c0,yi)Δi(c0,yi).
    (39)
    where Nv is the number of virions, i.e., the number of images if each virion has a tilt series containing only one image.
  3. Solve the linear system
    Fc=g
    (40)
    for the vector c which is the new value of the vector constructed from the cl,p,n coefficients.

Action 1 evaluates the expectations of the expectation-maximization algorithm while Actions 2 and 3 evaluate the maximization of the expectation-maximization algorithm where the function that must be maximized turns out to be a quadratic form in c so the location of the maximum can be computed by solving a linear system, i.e., Eq. (40). These calculations generalize immediately to the case where each virion is from one of a finite number of different classes and the class label for the virion shown in a particular image is not known (Doerschuk and Johnson, 2000; Yin et al., 2003). Because the approach of this paper is based on difference images (Eq. 30), multiclass algorithms require care that the image and the prediction of the image used to compute the difference image come from the same class and three multiclass situations and algorithms are described in the second paragraph of Section 2.3.

2.5. Parallel computation methods

The calculations implied by the algorithm described in this paper are large and so efficiency and parallel computing on a cluster of commodity PCs has been critical. Prust (2006, Section 2.2.6, pp. 19–21) provides methods based on Eq. (7) that allow fast computation of Lt and therefore of LΔ. The parallel software is based on modifications of the software described in Zheng (2002). The modifications are the new LΔ and the new pdfs for the random variables zi on which LΔ depends. (In the calculations reported in this paper, the zi are the orientation parameters for each image and take 1 of 60 values for each image but the possible values differ from image to image. For each image, the pdf on the zi is uniform over the possible values for that image). With regard to the pdfs, the key modification, a generalization, is to provide a different pdf for each difference image. This has major implications for the storage footprint of the software, specifically the storage required for the D matrices (Eq. 34), which will be returned to in Section 5.

3. Numerical results 1: The reconstruction of P22

For further details concerning the reconstruction of P22, please see Prust (2006).

3.1. Practical issues

The boxed images were hand selected based on absence of adjacent particles, broken particles, or junk in the image. No preference was given to images based on the visibility of the tail. Fig. 14 shows samples of both accepted images and discarded images. No masking of the images was performed because of the difficulty of designing a procedure that did not mask side-pointing tails.

Fig. 1
Examples of accepted and rejected P22 images. (a and b) Show a pair of P22 images that were included in the 3-D reconstruction while (c and d) show images that were rejected.

The selected images were oriented, modulo a rotation from the icosahedral group, and centered by quadratic correlation. A library of 5000 reference images with projection directions uniformly distributed through the asymmetric unit of the icosahedral group was computed by Spider (Frank et al. (1996)) using command PJ3Q from a high-resolution icosahedrally symmetric reconstruction of the capsid of P22 (Lander et al. (2006)). The rest of the processing was performed in Matlab.5 For each boxed image, estimates of the orientation (modulo a rotation from the icosahedral group) and origin offset are computed by maximizing the normalized quadratic correlation between the boxed image at a variety of shifts and the reference images each with a different orientation. Let σiref(χ) denote the ith reference image and hence the ith orientation and let χ0 denote the origin offset. Let σ(χ) denote one of the boxed images. The estimates for that boxed image are

i^,χ^0=arg maxi{1,...,5000},χ0{(m1Δ,m2Δ)T:m1,m2{mmax,...,+mmax}}J1(σ,σiref,χ0)
(41)

where

J1(σ,σiref,χ0)=χσ(χχ0)σiref(χ)[χσ2(χ)][χ[σiref(χ)]2],
(42)

Δ is the image sampling interval, and mmax = 2.

Difference images were computed by subtracting the reference image from the shifted boxed image after computing, by least squares, an optimal gain to apply to the reference image in order to compensate for the unknown scaling of the boxed image. The optimal gain for a particular boxed image, denoted by ĝ, is the gain that minimizes a cost function, denoted by J2:

g^=arg mingJ2(g,σ,σi^ref,χ^0)
(43)

where

J2(g,σ,σi^ref,χ^0)=χ[σ(χχ^0)gσi^ref(χ)]2.
(44)

The calculation of ĝ can be done explicitly in terms of σ, σi^ref, and [chi]0. Fig. 2 shows sample difference images for the P22 reconstruction.

Fig. 2
Generation of difference images for the P22 reconstruction. (a and d) Show a pair of raw P22 images. (b and e) Show the corresponding reference image from the 5000 image library. (c and f) Show the resulting difference images using the optimal gains.

Expectation-maximization is an iterative algorithm that requires an initial condition. As is described in more detail in Supplemental material Section C, the calculations described in this paper used random initial conditions computed in two steps. First, compute pseudo random variables cl,p,n from a pdf which is uniform on the interval [−ω, + ω] subject to the energy constraint that

ω2/4l=+p=1n=+[cl,p,n]23ω2/2.
(45)

Second, set cl,p,n equal to cl,p,n for those values of (l,p,n) where cl,p,n [set membership] R (please see Eq. 8) and set cl,p,n equal to cl,p,nexp(iξl,p,n) for the remaining values of (l,p,n) where ξl,p,n are pseudo random variables from a pdf which is uniform on the interval [0, 2π]. The value of ω was set to 0.00002 which yielded satisfactory performance. Because expectation-maximization is guaranteed only to converge to a local maximum of the likelihood function, a multi-start optimization was performed in which 99 initial conditions were tested and the best answer, i.e., the answer with highest likelihood, was retained.

Similar to most cryo EM and X-ray crystallography reconstruction algorithms, the resolution of the reconstruction is increased in a series of steps. Resolution of the model is controlled by truncating the l, p, and n sums in Eq. (3) or equivalently Eq. (10) to −lmaxllmax, 1 ≤ ppmax, and −nmaxnnmax. Table 1 lists the values of lmax, pmax, and nmax for each step and the corresponding number of cl,n,p coefficients, denoted by Nc. The multi-start optimization described in the previous paragraph was used for Step 1 and the answer from Step 1 was the best of the multi-start answers. In Steps 2–6, the single initial condition was always the answer from the previous step augmented with additional cl,p,n coefficients with value 0.

Table 1
Parameters at each step of the reconstruction algorithm. ±lmax, pmax, and ±nmax are the truncation limits for the infinite sums in the tail model (Eq. 3 or equivalently Eq. 10). Nc is the total number of cl,p,n coefficients.

In Steps 1 and 2 the difference image is predicted by Eq. (30) which includes both the true tail and the ghost tails. However, in Steps 3–6, in order to save computation, the ghost tails are omitted, i.e., Σg(κ; Ri) is deleted from Eq. (30).

The parameters for the reconstruction were as follows: A total of Nv = 276 images were used. Each image had a sampling interval of Δ = 4.04Å/pixel and measured 288 × 288 pixels. The CTF was unity. The cylinder containing the tail had radius R+ = 130Å and length z0 = 380Å. The symmetry of the tail was ξ = 6 fold rotational symmetry. The tail coordinate system was displaced by zt = −380Å from the capsid coordinate system which is the same as the whole-particle coordinate system. (Recall that the tail is nonzero in the tail coordinate system for the region −z0/2 ≤ z ≤ +z0/2 and that the capsid, tail, and whole-particle coordinate systems are described in Eqs. (1) and (2)).

3.2. The 3-D cube

Fig. 3 shows the final answer (i.e., Step 6) of the reconstruction at a single viewing angle but at various contour levels. Fig. 4 shows cross sections through the final answer (i.e., Step 6). Fig. 5 shows two views of the assembled P22 virus structure. Supplemental material Section D contains additional figures showing the 3-D real-space reconstruction of the tail. Supplement Fig. 12 shows the resulting structure at each step of the reconstruction process (see Table 1). Supplement Fig. 13 shows the final answer (i.e., Step 6) from various viewing angles. Note that there is no ambiguity in the relative location and orientation of the tail and capsid structures since the tail coordinate system is locked to the capsid coordinate system and because the tail reconstruction algorithm can reconstruct the tail in any rotation as implied by the data. The fact that the tail is not rotationally blurred demonstrates that the tail structure attaches to the capsid in a specific way (i.e., the 6-fold symmetric tail specifically connects to the 5-fold symmetric capsid). The resolution of the reconstruction is presented in Section 5.

Fig. 3
Final P22 tail reconstruction (i.e., Step 6) rendered at increasing contour levels from left to right.
Fig. 4
Cross sectional plots in the x–y plane of the final P22 tail reconstruction (i.e., Step 6). The cross sections are at distances 240, 380, and 530 Å from the center of the capsid in (a–c), respectively, which implies that (a–c) ...
Fig. 5
3-D reconstruction of the complete P22 virus structure. Side and end-on views are shown in (a and b), respectively. Note that the rotational uncertainty in assembling the capsid and tail structures was uniquely determined by the tail reconstruction algorithm. ...

3.3. The angle between a tail spike and the icosahedral symmetry axes

The mathematics used in this paper can represent a ξ-symmetric tail independent of the rotation of the tail around the z axis relative to the capsid. Therefore, the rotational relationship between the protein molecules in the tail and the icosahedral symmetry of the capsid can be determined free of any initial assumptions.

One method to visualize the rotational relationship is to select a region of z values and integrate ρ(x) over this region to create a 2-D averaged cross sectional real-space image of the tail. Assume that the region of z is specified in the tail coordinate system (which is displaced by zt from the complete particle coordinate system as is shown in Eqs. (1) and (2)) by −z0/2 ≤ zzz+ ≤ z0/2. The integral can be done analytically because the only factor of Eq. (3) that must be integrated is fn(z) and that integral has the value (please see Supplemental material Section G, Eq. 209)

z=zz+fn(z)dz=z0πnexp(iπn(z++z)/z0)sin(πn(z+z)/z0).
(46)

Alternatively, a 3-D real-space cube, such as is visualized in Section 3.2, can be summed in the z direction over the appropriate subset of planes to create an averaged 2-D cross sectional real-space image of the tail.

Taking the second approach, Fig. 6 shows averaged cross sections of the real-space 3-D cube visualized in Fig. 3 and Fig. 12(f) (Supplementary material) and Fig. 13 in which the sampling interval is 2 Å in all three coordinates. Specifically, Fig. 6(a) shows the sum of planes between z = 120Å and z+ = 150Å in the tail coordinate system (−260Å and −230Å in the total particle coordinate) and Fig. 6(b) shows the sum of planes between z = −100Å and z+ = 50Å in the tail coordinate system (−480Å and −330Å in the total particle coordinate system) where the center of the capsid is at the origin in the total-particle coordinate system. Fig. 6(a and b) correspond to the portal-end and mid-tail portion of the tail structure, respectively. The image shown in Fig. 6(b) shows the 6-fold symmetric locations of the protein molecules that make up the tail as present in the mid-tail portion of the structure. The center of mass of the molecule closest to the x axis in the positive rotational direction is 33.74° from the x axis. Fig. 6(a) shows similar information for the portal end of the structure. Here, in addition to the exact 6-fold symmetry, there is an approximate 12-fold symmetry, also seen in Lander et al. (2006), that was not imposed on the structure. The center of mass of the molecule closest to the x axis in the positive rotational direction is 0.83° from the x axis. As discussed in the final paragraph of Section 2.1, the structure described here is in one of 5 equivalent coordinate systems where the 5 systems are related by rotations by 2πn/5 (n [set membership] {0, …,4}) around the z axis. If the coordinate system is chosen in order to make these angles as small as possible, which is a way to make the angles unique, then 33.74° becomes 9.74° and 0.83° is unaltered.

Fig. 6
Averaged cross sections of the tail (a and b) and the relationships between the tail molecules and the icosahedral capsid symmetries (c and d). (a): Averaged cross section near the portal end of the tail, specifically, averaged over distances of 230–260 ...

The diagrams shown in Fig. 6(c and d) display the relationships between the molecules at the portal-end or the mid-tail portion of the tail structure and the icosahedral symmetries of the capsid structure. Specifically, the diagrams show the x and y axes, the centers of mass of the six (Fig. 6c) or 12 (Fig. 6d) molecules, and the x–y projection of the five 2-fold rotational symmetries of the capsid icosahedral symmetry that lie closest to the negative z axis, i.e., lie closest to the attachment point of the tail. Essentially the same information is provided in Fig. 4B of Lander et al. (2006). Given the small range of angles that are relevant (as is shown in Supplemental material Section L, the range of angles is 6° for 12-fold and 12° for 6-fold) and the inaccuracies of the structures, these angles appear to support those of Lander et al. (2006). Unlike the method of Lander et al. (2006), in which a 3-D structure for the tail machine is attached at a specific arbitrary rotation angle relative to the capsid symmetry axes and then the combined capsid-tail structure is refined, in the approach described in this paper no assumption about this angle is ever made at any stage of the algorithm and so this is an ab initio determination of the value of the angle. The reason that no value for this angle is ever assumed is that the mathematics used to represent the tail can represent the tail in any rotational position as is described in the final paragraph of Section 2.1.

At any cross sectional level, by using the maximum likelihood estimation error ideas of Section 4.1 and Monte Carlo ideas analogous to those of Section 4.4 to compute many structures and therefore many averaged cross sections of the tail, it would be possible to determine statistical information about the location of the molecules in the cross section of the tail. For instance, such statistical information might be the 2 × 2 covariance matrix describing the uncertainty in the center of mass location or a histogram describing the uncertainty in the angle. Given the level of precision with which this angle can be determined, these calculations would not provide substantial additional insight.

4. Resolution methods

A standard definition of resolution is based on the Fourier Shell Correlation (FSC) function (van Heel, 1987, Eq. 2; Harauz and van Heel, 1986, Eq. 17; Baker et al., 1999, p. 879) which compares the reciprocal-space scattering densities of two structures. The purpose of FSC, as described for instance in Saxton et al. (1982) (starting on p. 131 line 17), is to determine how the image noise effects the 3-D reconstruction where the determination is done in 3-D reciprocal space as a function of the magnitude of the spatial frequency vector. In the standard approach, this determination is made by comparing two 3-D reconstructions made from nonoverlapping subsets of the available images. In the approach described here, this determination is made based on the shape of the likelihood function at the maximum where the shape is measured as the matrix of mixed second-order partial derivatives of the log likelihood function with respect to the parameters evaluated at the parameter values that maximize the likelihood function. In the approach described here, the goal is to compute the probability density functions (pdfs) of certain random variables. We are not able to perform these calculations symbolically. Therefore, we use Monte Carlo methods to compute histograms which approximate the pdfs. So, while the Monte Carlo methods are important to implement our approach, they are not intrinsic to our approach.

Let Pa(k) and Pb(k) be the two reciprocal-space scattering densities to be compared. The FSC function [denoted by pFSC(k)] is a function of the magnitude of the reciprocal-space frequency vector (k) and is defined by

pFSC(k)=˙Pa(k)[Pb(k)]dΩ|Pa(k)|2dΩ|Pb(k)|2dΩ
(47)
=SPa,PbΩ(k)SPa,PaΩ(k)SPb,PbΩ(k)
(48)

where dΩ′ is integration over the angles of spherical coordinates (i.e., dΩ=φ=02πθ=0πsin(θ)dθdφ where θ′ and ϕ′ are the angles of spherical coordinates in reciprocal space) and

SPa,PbΩ(k)=˙Pa(k)[Pb(k)]dΩ
(49)

for any pair of reciprocal-space functions Pa(k) and Pb(k). Note that pFSC(k) is real valued because ρ(x) is real valued and that |pFSC(k)| ≤ 1 by the Cauchy–Schwarz inequality. The two structures, Pa(k) and Pb(k), are often the reconstructions based on even and odd numbered images, respectively. Once the FSC has been computed, the resolution is defined as the smallest value of k such that PFSC(k) is less than a threshold which may depend on k (van Heel and Schatz, 2005).

Our interest is a resolution measure with the following properties:

  1. The resolution measure distinguishes between axial and radial directions because that distinction is intrinsic in the mathematics used in this paper, e.g., the cutoff for the n sum versus the l and p sums in Eq. (3) influence axial versus radial directions.
  2. The resolution measure is related to the maximum likelihood criteria used to compute the reconstruction.
  3. The resolution measure attaches a probability to its results, analogous to the probability included in tests like t tests.
  4. The resolution measure provides the resolution of the reconstruction based on the full set of images rather than by comparing two probably lower resolution reconstructions based on even and odd numbered images, respectively.

In the remainder of this section a statistical version of FSC is described, suitable for separate axial and radial resolution measures, that is based on the general theory of errors in maximum likelihood estimation. As described above for standard FSC, it is still necessary to have a threshold and the ideas described do not include new ideas about the specification of the threshold. However, all current threshold ideas of which we are aware, including k dependent thresholds such as those advocated by van Heel and Schatz (2005), can be used.

4.1. General theory of estimation error covariance for maximum likelihood estimators

The standard theory (Efron and Hinkley, 1978) for the estimator error covariance of a maximum likelihood estimator is described in this section. Let y be the vector of data and c be the vector of unknown parameters. Let the estimate of c, which is a function of y, be denoted by ĉ(y). Let the Hessian of the log likelihood function, the matrix of mixed second-order partial derivatives of the log likelihood function, be denoted by H(c) with i, jth element defined by [partial differential]2 ln p(y|c)/[partial differential]ci[partial differential]cj where p(y|c) is the conditional probability density function on the data y given the unknown parameters c. Let c* be the true value of the parameters. The key result (Efron and Hinkley, 1978) is that the estimation error, ĉ(y) − c*, is approximately Gaussian distributed with mean vector 0 and covariance matrix − [H(ĉ(y))]−1.

4.2. A simpler example

In this section the general theory of Section 4.1 is demonstrated on the simpler example of Yin et al. (2003, Section 2.1) both to provide intuition and to demonstrate the accuracy of the general theory. Demonstrating the accuracy in almost any nonlinear problem requires Monte Carlo simulation and hence can only be done on simpler problems.

The simple example is to assume that each experiment produces a data point, denoted by yν, which is the sum of a random variable, denoted by Lν, times the unknown quantity, denoted by c, plus a noise denoted by nν:

yν=Lνc+nν.
(50)

The integer ν [set membership] {1, …, Nν} is an index describing independent realizations of the experiment. The quantities yν, Lν, c, and nν are all real numbers. The goal is to estimate the value of c from all of the yν data. Note the similarity with the cryo EM situation (Eq. 31). Assume that the sets of random variables {Lν : ν [set membership] {1, …, Nv}} and {nν : ν [set membership] {1, …, Nv}} are independent of each other, that the sequence of random variables Lν is independent and identically distributed according to a Gaussian pdf with mean mz and variance σz2, and that the sequence of random variables nν is independent and identically distributed according to a Gaussian pdf with mean 0 and variance σ2. Then, by direct computation (Yin et al., 2003, Eq. C14), the log likelihood function is

lnp(y|c)=Nv2ln(c2σz2+σ2)Nv2ln(2π)121c2σz2+σ2×ν=1Nv(yνcmz)2
(51)

and it is possible to show (Yin et al., 2003, Eqs. 57) that the maximum likelihood estimate of c, denoted by ĉ, is one of the three roots of the polynomial

0=c3σz4c2mzσz2y¯+c[σz2ryσ2(σz2+mz2)]+mzσ2y¯
(52)

where

y¯=1Nvν=1Nvyν
(53)
ry=1Nvν=1Nvyν2.
(54)

The Hessian required by the general theory is just the negative of the inverse of the second derivative with respect to the unknown c because c is a single real number. Starting from Yin et al. (2003, Eq. C15), the second derivative can be computed directly for arbitrary value of c with the result that

2lnp(y|c)c2=Nv(c2σz2+σ2)3[σz6c4+σz2σ4+3σz4ryc2σz2ryσ22σz4c3mzy¯+6σz2cmzy¯σ23σz2c2mz2σ2+mz2σ4]
(55)

Each iteration of Monte Carlo (indexed by the integer t) includes the following computations:

  1. Compute Nv pairs of pseudo-random variables Lν and nν drawn from the pdfs described previously.
  2. Compute Nv measurements yν from Eq. (50).
  3. Compute y and ry from Eqs. (53) and (54), respectively.
  4. Compute the coefficients and then the three roots of the polynomial described by Eq. (52).
  5. Among the real roots of Step 4, the estimate, denoted by ĉ(t), is the root that maximizes the log likelihood given by Eq. (51). Since the polynomial is of third order, at least one real root is guaranteed.
  6. Using Eq. (55), compute the Hessian at the estimate ĉ(t) and denote the result by h(t). The value −1/h(t) is the estimate of the estimation error covariance.
  7. Compute the true error, denoted by δ(t) and defined by δ(t) = cĉ(t).

The Monte Carlo estimates of the mean and variance of the error after T Monte Carlo iterations are the sample mean and variance of the δ(t) values, specifically,

δ¯=˙1Tt=1Tδ(t)
(56)
s2=˙1Tt=1T[δ(t)δ¯]2,
(57)

respectively.

Consider a case where the true value of c is 5.0, Nv = 100, mz = 2.0, σz = 0.1, and σ = 0.3. Based on T = 106 Monte Carlo iterations, the Monte Carlo estimate for the mean and variance of the estimation error are [delta with macron] = 0.000120 and s2 = 0.000849. In addition, the histogram of the T different estimation error variance results from the general theory of Section 4.1, i.e., the T different values of the covariance estimate −1/h(t) computed from the Hessian h(t), are shown in Fig. 7. The fact that the histogram is narrow and is centered around s2 = 0.000849 demonstrates the accuracy of the general theory of Section 4.1 on a problem that resembles a scalar version of the cryo EM problem.

Fig. 7
The histogram of results from the general theory of Section 4.1, i.e., the histogram of different values of −1/h(t) where h(t) is the Hessian. The Monte Carlo estimate for the mean and variance of the estimation error are [delta with macron] = ...

The calculation described in the previous paragraph and Fig. 7 does not make clear how the measure of performance proposed in this paper depends on the SNR of the original data. Therefore, in Fig. 8 are plotted the Monte Carlo variance of the estimation error, i.e., s2, (based on 103 Monte Carlo trials) and the result of the general theory for maximum likelihood estimators, i.e., −1/h(t), as a function of the SNR of the original data. The parameters are the same as in the previous paragraph except that σ varies in order to vary the SNR of the original data. The fact that −1/h(t) tracks s2 accurately as the SNR changes and the fact that both change dramatically as the SNR improves illustrate the high quality of the general theory for this specific example and the fact that the general theory is really about how SNR of the data influences SNR of the estimates.

Fig. 8
The variance of the estimation error, i.e., s2, computed by Monte Carlo using 103 trials and the approximate variance based on the general theory for maximum likelihood estimators, i.e., −1/h(t), both as a function of standard deviation of the ...

4.3. Fourier Axial Correlation (FAC) and Fourier Radial Correlation (FRaC) for cylindrical objects

As shown in Eq. (47), the standard FSC quantity is an integration over the two angles of spherical coordinates. For a cylindrical object, especially when the resolution can be independently controlled in the axial and radial directions, it is natural to consider integrations over various combinations of cylindrical coordinates. If the integration is taken over cylindrical shells of reciprocal space, i.e., ϕ′ and kz, then the criteria is called Fourier Radial Correlation (FRaC)6 and is denoted by pFRaC(kr) and the criteria measures radial resolution. Alternatively, if the integration is taken over cross sectional planes of reciprocal space, i.e., ϕ′ and kr, then the criteria is called Fourier Axial Correlation (FAC) and is denoted by pFAC(kz) and the criteria measures axial resolution. Define the integrals over cylindrical shells and over cross sectional planes, both in reciprocal space, by

SPa,Pbφ,kz(kr)=˙kz=+φ=02πPa(k)[Pb(k)]dφdkz
(58)
SPa,Pbφ,kr(kz)=˙kr=0φ=02πPa(k)[Pb(k)]dφkrdkr
(59)

for any pair of reciprocal-space functions Pa(k) and Pb(k). Then FRaC and FAC are defined by

pFRaC(kr)=˙SPa,Pbφ,kz(kr)SPa,Paφ,kz(k)SPb,Pbφ,kz(kr)
(60)
pFAC(kz)=˙SPa,Pbφ,kr(kz)+SPa,Pbφ,kr(kz)[SPa,Paφ,kr(kz)+SPa,Paφ,kr(kr)][SPb,Pbφ,kr(kz)+SPb,Pbφ,kr(kz)]
(61)
=[SPa,Pbφ,kr(kz)]SPa,Paφ,kr(kz)SPb,Pbφ,kr(kz)
(62)

where R indicates the real part. As is shown in Supplemental material Section K, both pFRaC(kr) and pFAC(kz) are real valued because ρ(x) is real valued. In addition, by the Cauchy-Schwarz inequality, |p FRaC(k)| ≤ 1. Finally, starting with Eq. (62) and using |[SPa,Pbφ,kr(kz)]||SPa,Pbφ,kr(kz)| followed by the Cauchy-Schwarz inequality, it follows that |pFAC(k)| ≤ 1. Defining pFAC(kz) as the sum of terms at kz and −kz allows it to combine the positive and negative frequencies which have the same interpretation with respect to resolution. In addition, summing the terms makes pFAC(kz) a function of |kz| which, like kr, ranges from 0 to ∞ while, without the sum of terms, it would be a function of kz, which ranges from −∞ to +∞. Finally, it is only with the sum of terms that ρ(x) real implies that pFAC(kz) is also real (Supplemental material Section K).

For the tail model of Eq. (3), SPa,Pbφ,kz(kr) and SPa,Pbφ,kr(kz) can be computed symbolically in terms of the cl,p,na and cl,p,nb for the two structures Pa(k) and Pb(k) under the assumption that both structures share the same values of z0 (Eqs. 4 and 5) and R+ (Eq. 6). The reason that the calculations can be done symbolically is the orthogonality of the factors of Lt(kr, ϕ′, kz), (l,p,n) (Eq. 11):

φ=02πexp(ilξφ)exp(ilξφ)dφ=2πδl,l
(63)
kr=0Hl,p(kr)Hl,p(kr)krdkr=R+22[J|l|+1(γ|l|,p)]2δp,p
(64)
kz=Q(kzn/z0)Q(kzn/z0)dkz=z0δn,n
(65)

where Eq. (63) is an elementary integral, Eq. (64) follows from Supplemental material Eq. 127 since Hl,p(kr) [set membership] R, and Eq. (65) is Supplemental material Eq. 225. Using these results,

Spa,pbφ,kz(kr)=2πz0l=+p=1n=+p=1cl,p,na(cl,p,nb)Hl,p(kr)Hl,p(kr)
(66)
Spa,pbφ,kz(kz)=2πl=+p=1n=+n=cl,p,na(cl,p,nb)×R+22[J|l|+1(γ|l|,p)]2Q(kzn/z0)Q(kzn/z0).
(67)

As is shown in Supplemental material Section J, Eq. (67) implies that Eq. (62) can be simplified and that the simplified equation is rational in kzz0. Taking the limit as kzz0 grows large leads to the result (Eq. 225) that

limkzpFAC(kz)=lpjl,pnn[cl,p,nacl,p,nb](1)n+n[lpjl,pnn|cl,p,na|2(1)n+n][lpjl,pnn|cl,p,nb|2(1)n+n]
(68)

where

jl,p=˙R+22[J|l|+1(γ|l|,p)]2.
(69)

4.4. The connection between estimation error covariance and FAC, FRaC, and FSC

Let c* denote the vector containing the true cl,p,n coefficients, let ĉ(y) denote the vector containing the maximum likelihood estimates of the cl,p,n coefficients from the image data y, and let H(ĉ(y)) be the Hessian of the log likelihood evaluated at the maximum likelihood estimate of the cl,p,n coefficients. The Hessian can be computed starting from Doerschuk and Johnson (2000, Eq. 15, p. 1721) by taking partial derivatives of the log likelihood with respect to two components of the vector containing the cl,p,n coefficients and using the fact that the first derivative of the log likelihood when evaluated at the maximum likelihood estimate is zero because the estimate is the maximum. The resulting formula is Prust (2006, Section 3.2, pp. 35–36)

H(c^(y))=i=1Nυ(1γi(c^,yi))Δi(c^,yi,k)
(70)

where it is important to note that H(ĉ(y)) can be computed in terms of the γ and Δ variables of Eqs. (35) and (37) so its computation does not add to the memory footprint or computational cost of the algorithm. Similar to the comment at the end of Section 2.4, these calculations can be generalized to the case where each virion is from one of a finite number of different classes and the class label for the virion shown in a particular image is not known (Doerschuk and Johnson, 2000; Yin et al., 2003) with the interesting result that the Hessian matrix has a block diagonal structure where each class corresponds to a different block. From the general theory of Section 4.1, the estimation error ĉ(y) − c* is Gaussian distributed with mean 0 and covariance −[H(ĉ(y))]−1. Therefore, if the entire experiment and computation were repeated many times, the estimates ĉ(y(n)) would be Gaussian distributed with mean c* and covariance − H[(ĉ(y))]−1 where n indexes the repetitions. If the true c* is approximated by the maximum likelihood estimate ĉ(y), then the distribution of ĉ(y(n)) is fully specified and sampling from this distribution is a generalization of the two reconstructions traditionally computed by using even versus odd numbered images.

Suppose there are two structures that are either computed from different images or are different samples from the distribution of the previous paragraph. Denote the estimated cl,p,n coefficients by ĉa and ĉb. From Eqs. (60) and (66) (Eqs. 62 and 67), FRaC (FAC) can be computed for any value of kr (kz) from ĉa and ĉb. As described in Section 4, resolution in FSC is measured as the smallest value of k for which pFSC(k) is below the threshold. Adopting the same definition for pFRaC(kr) and pFAC(kz) implies that the resolution, denoted by kr* and kz*, respectively, is defined to be

kr=minkr0{kr:pFRaC(kr)<tFRaC(kr)}
(71)
kz=minkz0{kz:pFAC(kz)<tFAC(kz)}
(72)

where tFRaC(kr) [tFAC(kz)] is the radial (axial) threshold which is possibly kr (kz) dependent. Though it is not shown in the notation, kr* (kz*) is a function of pFRaC(kr) [pFAC(kz)] which is a function of ĉa and ĉb. Therefore, kr* and kz* are derived random variables, derived from ĉa and ĉb.

Because kr* and kz* are derived random variables, their pdfs can be computed by Monte Carlo. Before starting the iterative Monte Carlo calculation, it is most efficient to compute the Cholesky factorization of H(ĉ(y)), which is denoted by H1/2 and which has the property that H(ĉ(y)) = H1/2(H1/2)T. Each iteration of Monte Carlo (indexed by the integer t) includes the following computations:

  1. Compute υa,υb [set membership] RNc whose components are independent and identically distributed Gaussian pseudo random variables with mean 0 and variance 1 where Nc is the number of cl,p,n coefficients.
  2. Compute ca = H1/2 υa + ĉ(y) and cb = H1/2 υb + ĉ(y) which are samples from the Gaussian distribution for the estimates. ca and cb play the role of the structures computed from even versus odd numbered images.
  3. From Eqs. (60, 66, and 71), compute kr*(t).
  4. From Eqs. (62, 67, and 72), compute kz*(t).

After T Monte Carlo iterations, an estimate of the pdf for kr* (kz*), denoted by pk*(kr) [pkz*(kz)], can be determined by computing the histogram for the set kr*(t) [kz*(t)] for t [set membership] {1,…,T}. Once a probability p is chosen, e.g., p = .01, the estimates for kr* and kz* denoted by kr and kz, respectively, can be determined as the values such that

kr=0k^rpkr(kr)dkr=p
(73)
kz=0k^zpkz(kz)dkz=p.
(74)

Alternatively, the same pdfs can be used to compute confidence intervals. Let kr (kz) be the sample mean of kr*(t) [kz*(t)], i.e.,

k¯r=1Tt=0Tkr(t)
(75)
k¯z=1Tt=0Tkz(t).
(76)

Then the symmetric 100q% confidence intervals are [krδr, kr + δr] and [kzδz, kz + δz] where δr and δz are defined by

kr=k¯rδrk¯r+δrpkr(kr)dkr=q
(77)
kz=k¯zδzk¯z+δzpkz(kz)dkz=q,
(78)

respectively.

While the equations are not shown here, these ideas can also be applied to the FSC based on Yin et al. (2003, Eqs. 23 and 25).

5. Numerical results 2: The resolution of the reconstruction of P22

As described in Section 2.5, the storage footprint of the current software system is large. Each D matrix requires Nc(Nc + 1)/2 locations (where Nc is the number of cl,p,n coefficients used) and the number of D matrices is the number of abscissas, which is Ng = 60 since the orientation uncertainty is due to the uncertainty in which icosahedrally related orientation is present, times the number of images (denoted by Nν) since each image has a different set of Ng = 60 possible orientations. Fitting the D matrices into memory constrains the number of cl,p,n coefficients Nc and/or the number of images Nν and the largest calculations reported here have Nc = 385 (implied by lmax = 2, pmax = 7, and nmax = 5) and Nν = 276. With so few images, the traditional method of performing reconstructions with even and with odd numbered images and then comparing the two 3-D cubes by FSC is not attractive.

Independent of image quality and image number, Nc sets an upper bound on the resolution that can be achieved. The Nc cl,p,n coefficients describe a function that is nonzero in a cylinder of length z0 and radius R+. Since the function is constrained to have ξ-fold rotational symmetry, the function is uniquely defined by its values on 1/ξ of the cylinder's volume. Alternatively, the same volume might be represented by Nc voxels where each voxel measures T × T × T. Equating the volume measured in voxels and the volume of 1/ξ of the cylinder gives NcT3=πR+2z0/ξwhich implies that

T=[πR+2z0ξNc]1/3.
(79)

The resolution that a model with Nc coefficients is able to represent when averaged in all directions, independent of the number and quality of the images, is unlikely to exceed about 2T unless the basis functions used in Eq. (3) much more efficiently represent the electron scattering intensity function in comparison with voxel basis functions (i.e., basis functions which are 1 in a voxel and 0 outside of the voxel). For ξ = 6, Nc = 385, R+ = 130Å, and z0 = 380Å, the value of T is T = 20.59Å. Therefore the achieved resolution may be limited by the value of Nc rather than the number and quality of the images.

A second measure of the upper bound on resolution that is set by Nc independent of the image quality and image number can be determined by considering special 3-D structures composed of just those cl,p,n and corresponding basis functions with the highest spatial frequencies. In the truncated system used for numerical calculations, that special structure has cl,p,n coefficients defined by

cl,p,nHF={1,l{±lmax},p=pmax,n{±nmax}0,otherwise
(80)

where “HF” stands for “high frequency”. Let PHF(k) be the corresponding reciprocal space function. The averaged distribution in reciprocal space of the energy in the special structure is determined by SPHF,PHFφ,kz(kr) (averaged over cylindrical shells) and SPHF,PHFφ,kr(kz) (averaged over cross sectional planes) and these functions are plotted in Fig. 9. While the curves oscillate, the curve has permanently decreased to below 10% of its maximum value by 0.0495Å−1 for SPHF,PHFφ,kz(kr) and by 0.0150Å−1 for SPHF,PHFφ,kr(kz) and to below 1% of its maximum value by 0.0568Å−1 for SPHF,PHFφ,kz(kr) and by 0.0226Å−1 for SPHF,PHFφ,kr(kz). With Nc = 385 (implied by lmax = 2, pmax = 7, and nmax = 5) the mathematical model cannot achieve higher spatial resolution than somewhere between the 10% and 1% values of kr and kz independent of the number and quality of the images used to determine the values of the cl,p,n coefficients.

Fig. 9
The distribution of energy in reciprocal space (0 ≤ kr, kz ≤ 0.08Å−1) for the highest frequency cl,p,n and corresponding basis function. Both SPHF,PHFφ,kz(kr), showing energy averaged over cylindrical ...

Due to the considerations of the previous two paragraphs, the resolution cannot be above kFRaC=̇0.06Å−1 or kFAC=̇0.023Å−1 for FRaC or FAC, respectively, and plots are stopped at k*=̇ max(kFRaC, kFAC) = 0.06Å−1. Fig. 10 shows 4 realizations of pFRaC(kr) (Eq. 60) and pFAC(kz) (Eq. 62). The plots of pFAC(kz) clearly show the approach of pFAC(kz) to the asymptotic value as kzz0 (z0 = 380Å) grows large as is expected from Eq. (68). It is apparent that resolution is not limited by the number or quality of the images, since the curves remain high over the entire range of 0 to kFRaC or 0 to kFAC for FRaC or FAC, respectively, but rather is limited by the number of cl,p,n coefficients (i.e., by Nc) that our current software can accommodate. Therefore, in order to compute resolutions less than kFRaC and kFAC, it is necessary to choose a strict threshold in Eqs. (71) and (72). Using T = 1000 Monte Carlo iterations and threshold functions tFRaC(kr) = 0.95 (Eq. 71) and tFAC(kz) = 0.95 (Eq. 72) the histograms of kr* (Eq. 71) and kz* (Eq. 72) values are shown in Fig. 11. The reason that the kr* histogram of Fig. 11(a) is bimodal is that the kr value at which pFRaC(kr) first drops below the threshold typically occurs in one of two successive oscillations of pFRaC(kr) where the oscillations are due to the Hl,p(kr) functions (defined in Eq. 13) via Eq. (66) in Eq. (60). Choosing p = 0.02 in Eqs. (73) and (74) leads to a radial resolution of kr = 0.0301Å−1 (33.3 Å) and an axial resolution of kz = 0.0142Å−1 (70.6 Å) where both are with probability p = 0.02, that is, the probability that the resolution is actually less than kr = 0.0301Å−1 (33.3 Å) and kz = 0.0142Å−1 (70.6 Å) is p = 0.02.

Fig. 10
Four example realizations of pFAC(kz) (a) and pFRaC(kr) (b) defined in Eqs. (60) and (62), respectively, and plotted on the range from 0 to k* = 0:06Å−1. A different color is used for each realization so that the curve of an individual ...
Fig. 11
Histograms of kz* (a, defined in Eq. (72), range 0:01 ≤ kz* ≤ 0:06Å−1 ) and kr* (b, defined in Eq. (71), range 0:024 ≤ kr* ≤ 0:044Å−1 ).

6. Discussion

Two connected methodological contributions are presented in this paper. The first is a maximum likelihood reconstruction method for an asymmetric reconstruction of an object of the form “sphere plus cylinder” where the sphere part has icosahedral symmetry, the cylinder part has ξ-fold rotation symmetry around the axis of the cylinder, and the axis of the cylinder and a 5-fold axis of the icosahedron are coincidental. While an icosahedrally symmetrized reconstruction of the object is used, in particular, the method described in this paper is really based on images that are the difference of the experimental image and the predicted icosahedrally symmetric image, no previous structure for the cylinder is required. In addition, the rotation angle between the cylindrical and spherical objects is determined without any assumptions directly from the data. The second methodological contribution is a statistical method of measuring resolution that combines standard Fourier Shell Correlation (FSC) ideas with standard statistical maximum likelihood ideas and can measure resolution axially and radially in a cylindrical object as well as radially in a spherical object as is done by FSC.

The reconstruction and the resolution methods are used to study the infectious P22 bacteriophage virion using a subset of the images used in Lander et al. (2006). Without making any assumptions at any stage of the reconstruction algorithm, the rotational positioning of the components of the tail are determined relative to the icosahedral symmetry axes of the capsid as is described in Fig. 6. The combination of the reconstruction and resolution calculations is important because limitations of the reconstruction software make a resolution calculation based on reconstructions from even versus from odd numbered images unattractive and the methods described here can be applied to reconstructions from the complete set of images. With a correlation threshold of 0.95, the resolution in the tail measured radially is greater than 0:0301Å−1 (33.3Å) and measured axially is greater than 0:0142Å−1 (70.6Å) both with probability p = 0.02.

Supplementary Material

01

Acknowledgments

C.J.P. and P.C.D. gratefully acknowledge the support provided by the National Institutes of Health under Grant 1R01EB000432-01 and the National Science Foundation under Grants CCR-0098156 and CCF-0735297 and the computations were carried out by C.J.P. and P.C.D. at Purdue University using facilities supported by the Army Research Office under Contract DAAD19-99-1-0015 and the National Science Foundation under Grant EIA-0130538. Electron microscopic imaging was conducted at the National Resource for Automated Molecular Microscopy, which is supported by the NIH through the National Center for Research Resources' P41 program (RR17573).

Footnotes

4All 2-D images (boxed images, cross sections, etc.) in this paper were made with Matlab (http://www.mathworks.com/) and all 3-D surface plots were made with the Spider/Web system (Frank et al., 1996).

5http://www.mathworks.com/.

6“FRaC” is used in order not to conflict with Fourier Ring Correlation which is abbreviated FRC.

Appendix A. Supplementary data: Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.jsb.2009.04.013.

References

  • Altmann SL. On the symmetries of spherical harmonics. Proc Camb Phil Soc. 1957;53:343–367.
  • Baker TS, Olson NH, Fuller SD. Adding the third dimension to virus life cycles: three-dimensional reconstruction of icosahedral viruses from cryo-electron micrographs. Microbiol Mol Biol Rev. 1999;63(4):862–922. [PMC free article] [PubMed]
  • Blanc E, Roversi P, Vonrhein C, Flensburg C, Lea SM, Bricogne G. Refinement of severely incomplete structures with maximum likelihood in BUSTER-TNT. Acta Crystallogr. 2004;D60(Part 121):2210–2221. [PubMed]
  • Bubeck D, Filman DJ, Cheng N, Steven AC, Hogle JM, Belnap DM. The structure of the poliovirus 135S cell entry intermediate at 10-Angstrom resolution reveals the location of an externalized polypeptide that binds to membranes. J Virol. 2005;79(2):7745–7755. [PMC free article] [PubMed]
  • Doerschuk PC, Johnson JE. Ab initio reconstruction and experimental design for cryo electron microscopy. IEEE Trans Info Theory. 2000;46(5):1714–1729.
  • Efron B, Hinkley DV. Assessing the accuracy of the maximum likelihood estimator: observed versus expected Fisher information. Biometrika. 1978;65(3):457–487.
  • Erickson HP. The Fourier transform of an electron micrograph—first order and second order theory of image formation. In: Barer R, Cosslett VE, editors. Advances in Optical and Electron Microscopy. Vol. 5. Academic Press; London and New York: 1973. pp. 163–199.
  • Frank J, Radermacher M, Penczek P, Zhu J, Li Y, Ladjadj M, Leith A. SPIDER and WEB: processing and visualization of images in 3D electron microscopy and related fields. J Struct Biol. 1996. pp. 190–199. Available from: < http://www.wadsworth.org/resnres/frank.htm>. [PubMed]
  • Frank J, Verschoor A, Boublik M. Computer averaging of electron micrographs of 40S ribosomal subunits. Science. 1981;214(4527):1353–1355. [PubMed]
  • Harauz G, van Heel M. Exact filters for general geometry three dimensional reconstruction. Optik. 1986;73(4):146–156.
  • Jiang W, Chang J, Jakana J, Weigele P, King J, Chiu W. Structure of epsilon15 bacteriophage reveals genome organization and DNA packaging/injection apparatus. Nature. 2006;439(7076):612–616. [PMC free article] [PubMed]
  • Lander GC, Tang L, Casjens SR, Gilcrease EB, Prevelige P, Poliakov A, Potter CS, Carragher B, Johnson JE. The structure of an infectious P22 virion shows the signal for headful DNA packaging. Science. 2006;312(5781):1791–1795. [PubMed]
  • Laporte O. Polyhedral harmonics. Z Naturforschg. 1948;3a:447–456.
  • Lehmann EL, Casella G. Theory of Point Estimation. 2nd. Springer-Verlag; New York: 1998.
  • McCoy AJ, Grosse-Kunstleve RW, Storoni LC, Read RJ. Likelihood-enhanced fast translation functions. Acta Crystallogr. 2005;D61(Part 4):458–464. [PubMed]
  • Prust CJ. Ph D thesis. School of Electrical and Computer Engineering, Purdue University; West Lafayette, Indiana, USA.: 2006. Model-based statistical inference problems concerning non-linear 3-D tomography with applications to the structural biology of asymmetric virus particles.
  • Saxton WO, Baumeister W. The correlation averaging of a regularly arranged bacterial cell envelope protein. J Microsc. 1982;127(2):127–138. [PubMed]
  • Scheres SHW, Gao H, Valle M, Herman GT, Eggermont PPB, Frank J, Carazo JM. Disentangling conformational states of macromolecules in 3D-EM through likelihood optimization. Nat Methods. 2007;4(1):27–29. [PubMed]
  • Scherzer O. The theoretical resolution limit of the electron microscope. J Appl Phys. 1949;20:20–29.
  • Singh V, Marinescu DC, Baker TS. Image segmentation for automatic particle identification in electron micrographs based on hidden markov random field models and expectation maximization. J Struct Biol. 2004;145(1–2):123–141. [PMC free article] [PubMed]
  • van Heel M. Similarity measures between images. Ultramicroscopy. 1987;21:95–100.
  • van Heel M, Schatz M. Fourier shell correlation threshold criteria. J Struct Biol. 2005;151:250–262. [PubMed]
  • Yin Z, Zheng Y, Doerschuk PC, Natarajan P, Johnson JE. A statistical approach to computer processing of cryo electron microscope images: virion classification and 3-D reconstruction. J Struct Biol. 2003;144(12):24–50. [PubMed]
  • Zhang P, Mueller S, Morais MC, Bator CM, Bowman S, nad Hafenstein Valorie D, Wimmer E, Rossmann MG. Crystal structure of CD155 and electron microscopic studies of its complexes with polioviruses. PNAS. 2008;105(47):18284–18289. [PubMed]
  • Zheng Y. Master's thesis. School of Electrical and Computer Engineering, Purdue University; West Lafayette, Indiana, USA.: 2002. Parallel implementations of 3-D reconstruction algorithms for cryo electron microscopy: a comparative study.