2.1 Problem formulation and definitions

A tomogram is described by a 3D grid of voxels that are associated with values of electron optical density. Let ={1,2,…,N} be the set of indices of all voxels in the density map, and let x_i³ be the location of voxel i in the map.

Intensity vector: then, a tomogram is represented by a N-dimensional vector of the ordered list of intensity values for all voxels.

(1)

where N is the total number of voxels and f_i is the intensity value of voxel i.

Class label vector: a pattern in a density map is defined as a collection of voxels that may have different intensities but are assigned to the same pattern class. Let us define the set of all pattern classes in a tomogram as ={1,2,…,C}, with C as the total number of classes. Each voxel in the density map can be assigned to one of the classes in . Similar to an intensity vector, we therefore define a class label vector

(2)

where l_i is the class label of voxel i. Whereas the intensity vector f is known, the class label vector l is unknown and to be inferred. The task of identifying recurrent patterns in a density map is therefore equivalent to determining the class label vector l. To achieve this goal, we apply a two-step procedure, combining a heuristic approach generating an initial tomogram classification with a subsequent refinement using a statistical inference method that maximizes the joint probability of the class assignment l and the intensity f.

2.2 Initial tomogram classification

The goal of the following section is to define an initial class label vector l, which serves in a later step as an input in the GHMRF-based refinement process.

2.2.1 Identification of segmented objects.

First, Gaussian filtering is applied to generate a blurred density map f^blur reducing the influence of noise in the segmentation. Then a seed growth approach is adapted from the watershed segmentation (Dougherty, 1993) as follows. All local intensity maxima in the filtered map are identified. Local maxima with low intensity values are likely a result of noise and are discarded. The remaining local maxima with high intensity voxels are used as seeds for an extended watershed segmentation procedure. During the seed growth segmentation, a voxel rank r is introduced that captures the order at which voxels are included to a region. A segmented region S_s is defined, where the watershed algorithm is terminated at a given voxel rank r_s^max. Let v_s,r be a specific voxel index, of segment s and voxel rank r, then the set of voxels in a region of segmented object S_s is defined as

(3)

such that for all s. Note that the filtered tomogram is only used to identify the indices of voxels in the segmented objects. Classifications are performed using the unfiltered intensity values of these voxels.

2.2.2 Rotation-invariant feature vectors.

To introduce an efficient way to compare the similarities between all the distinct objects, we introduce rotation-invariant feature vectors (Kazhdan et al., 2003; Saha et al., 2010; Xu et al., 2009).

Rotation-invariant feature vectors p(x_i) describe the intensity distribution of the tomogram in the neighborhood of a voxel located at x_i. To construct a feature vector at the voxel location x_i, we divide the neighborhood of the voxel x_i into M concentric shells (Fig. 2)

(4)

where V_j(x_i) is defined as the set of voxels that fall into the concentric shell centered at x_i and is defined by the two radii r_j−1 and r_j, where with r_j=jR/M and R is the largest chosen radius. If a concentric shell V_j(x_i) is thin (i.e. r_j+1−r_j≈ voxel length), then the voxel intensities f(x_k) with kV_j(x_i) can be approximated by a spherical function g that is defined on the surface of a sphere in spherical coordinates.

(5)

where , and θ, ϕ are the colatitude and longitude angles, respectively. g can then be approximated by a sum of its spherical harmonics (Hobson, 1931):

(6)

where L is a given bandwidth, and a_lm is a coefficient associated with the complex spherical harmonics function Y_l^m, which is independent to g.

Open in a separate window

Fig. 2.

Neighborhood volumes defined as a series of concentric shells around voxel location x_i for voxel i. Schematic view of a 2D grid with individual voxels shown as dark grey dots. Concentric shells are constructed that are centered at x_i. The largest radius is defined as R. All radii are defined as r_j=jR/M, with M as the maximal number of shells. A neighborhood volume V_j(x_i)={k:r_j−1<|x_k−x_i|≤r_j} is defined as all voxels that fall into a concentric shell defined by two radii, r_j−1 and r_j with r_j−1<r_j. As an example, the neighborhood shell V₁₀(x_i) is shown in light grey, defined as the set of voxels located between radii r₉ and r₁₀.

Based on such a decomposition, the intensity distribution of all the voxels in shell V_j(x_i) can then be described by a set of L rotation-invariant features {s_jl|l=1,…,L}(Kazhdan et al., 2003) as follows:

(7)

Such features are calculated for all the M shells. The rotation-invariant description of the density distribution around the location of a voxel x_i is then defined by the following feature vector:

(8)

whose elements consist of the ordered sequence of s_jl for all consecutive shells {V_j(x_i)|:j=1,…,M}. These feature vectors are rotation invariant, which means that they are independent of the relative orientation of the local density distribution. Therefore, feature vector-based similarities between density regions can be calculated even if these regions are at different relative orientations to each other. Previously, we have shown that such rotational invariant features are also robust against noise in maps (Xu et al., 2009).

2.2.3 Comparing the similarity of segmented objects.

The goal is to cluster segmented objects based on their feature vector dissimilarity. The dissimilarity o_a,b between objects a and b is defined as one minus the fraction of voxels with similar feature vectors in both objects S_a and S_b.

(9)

where is the number of similar feature vectors that appear both in objects S_a and S_b.

(10)

_i^fea is the set of feature vectors that are most similar to feature vector i, but are located spatially apart from i in the tomogram. and are the total number of feature vectors in S_a and S_b, respectively. is the sum of the total number of feature vectors in object a for which a similar feature vector exists in object b and the number of feature vectors in b with a similar feature vector in a.

To determine _i^fea, the voxels in the tomogram are clustered based on the Euclidian distance between their feature vectors. Cryo electron tomograms are generally so large that not all feature vectors can be stored in computer memory at the same time. We, therefore, employ a large-scale clustering technique, (BIRCH) (Harrington and Salibián-Barrera, 2008; Zhang et al., 1997). The clustering cutoff is chosen heuristically by sampling the similarities between randomly selected feature vectors. Good performances have been achieved when the cutoff comprises the top 5% of sampled feature vector similarities. Then all voxels are detected that are in the same cluster but are apart in grid space so that they are not direct neighbors in the tomogram. For each voxel i in each cluster A a voxel set is defined as U_i={k:|x_k−x_i|>r_d,kA}, where r_d is a predefined grid distance.

It is computationally beneficial to restrict ^fea_i to only a fixed number of the closest voxel neighbors. Here, ^fea_i is defined as the subset of U_i with the 20 voxels whose feature vectors are most similar to the feature vector of i. Slight variations of this number do not affect the outcome of our calculations.

2.3 Classification refinement by Gaussian Hidden Markov Random Field Model

In this section, the initial classification is refined using an GHMRF model. In the GHMRF framework, class labels are not directly observable variables and are modeled as a hidden random field. The class labels are determined by maximizing the joint probability P(l,f) of class labels l and the intensity vector f. GHMRF has been used for image analysis and segmentation (Li et al., 2009; Zhang et al., 2001). Here, we apply and extend the method to 3D density maps for the classification of frequently occurring macromolecular complexes in cryo electron tomograms.

2.3.1 Gaussian Hidden Markov Random Field Model.

In the GHMRF framework, the voxel labels are modeled as a hidden random field. Their conditional independence follows the Markov property and can be described by an undirected graph where each node corresponds to a voxel.

(11)

where l_i={l_j|∀j_i} are the labels of neighbor voxels of i.

In standard GHMRF models, the voxel neighborhood is defined by an undirected graph with voxels as vertices and edges as the cubic grid connecting the vertices (Fig. 3). In our method, this graph is augmented by additional edges between those voxels that share similar feature vectors but are at far distance in the tomogram grid. In other words, for a given voxel i, its neighborhood list _i includes all direct grid neighbors in the tomogram and all voxels that have a similar density environment even if they are located at far distance (Fig. 3). By augmenting the list of grid neighbors, we are able to connect those voxels in the graph that have the same class label, even if these voxels are part of different copies of the same complex located at different regions of the tomogram. The neighborhood _i of voxel i is, therefore, defined as

(12)

where ^map_i is the set of all voxels that are adjacent to i in the tomogram grid and ^fea_i is the set of voxels that have similar feature vectors.

Open in a separate window

Fig. 3.

Gaussian Hidden Markov Random Field. GHMRF with observable intensity random field (above in red) and hidden class random field (below in blue). In the hidden field, the Markov property graph is defined by the direct neighbors of voxels in the grid (grey connections, for simplicity only a 2D grid is shown) and also by voxels with similar feature vectors (green dotted connections) (i.e. a green connection is formed if two voxels are defined as neighbors in the feature space).

Voxel intensities f_i form an observable random field. Given the class label l_i=l, the voxel intensity f_i is assumed to follow a conditional probability distribution

(13)

where Φ(f;μ_l,σ_l) denotes Gaussian distribution with mean μ_l and SD σ_l that are specific to class label l. In addition, given one instance of l, the intensities f are assumed to be conditionally independent.

(14)

2.3.2 Iterative optimization method for voxel classification.

The voxel re-classification is formulated as a maximization problem:

(15)

is the best estimation for the true class label vector. P(f,l) is the joint probability of the intensities f and class labels l, where θ={μ_l,σ_l} is the collection of mean and SD parameters for the intensity values of the voxels in each of the pattern classes. To perform the maximization, we use an iterative approach (Zhang et al., 2001) that consists of two steps. First, for a given estimation of the parameter set at iteration t, a classification is performed to obtain updated class labels . Then expectation maximization determines the optimal parameter set given obtained from the first step.

Step 1: classification.

Given , we want to obtain

(16)

According to (Zhang et al. (2001), the minimization problem can be formulated as

(17)

with the energy functions

(18)

and

(19)

where κ is the set of all cliques on the graph. Following Zhang et al. (2001), we choose V_c=−δ_lilj, where δ is the Kronecker delta function. Calculating the global minimum of Equation (17) is computationally infeasible. However, because both U(f|l) and U(l) can be expanded into sums over all voxels i, the optimal solution can be approximated using the iterated conditional modes (ICM) algorithm (Besag, 1986), which is popular in solving such optimization problems. This method optimizes U(l_i|l_{i},f) for each i, while keeping all other labels l_{i} fixed.

Step 2: Expectation maximization for model fitting.

The expectation maximization method updates parameters θ by maximizing the conditional expectation

(20)

Following Zhang et al. (2001), new parameters are estimated as

(21)

and

(22)

where

(23)

Following Li (2009) and given the class labels l_i obtained from the classification step, the probability for a voxel i to be assigned to class label l is calculated as

(24)

The iterations between expectation and maximization step are repeated until convergence is reached, leading to the newly refined class label vector l.

Finally, a class label can be defined for each object as the label that is assigned to the majority of voxels in the object. The classification performance is increased, if all the voxels with labels different to the corresponding object label are reassigned to the background class.

2.4 Generating simulated cryo electron tomograms

For a reliable assessment of the method, tomograms must be generated by simulating the actual tomographic image reconstruction process, allowing the inclusion of noise, tomographic distortions due to missing wedge and electron optical factors such as contrast transfer function (CTF) and Modulation Transfer Function (MTF) (Fig. 4).

Open in a separate window

Fig. 4.

Simulated electron tomograms including missing wedge effects, CTF and MTF for a tomogram at different SNR levels. (A) Contour volume representation of the tomogram and (B) a slice through x–y plane of the tomogram (top panel) and a slice through the x–z plane of the tomogram (bottom panel).

We follow a previously applied methodology for simulating the tomographic image formation mechanism as realistically as possible (Beck et al., 2009; Förster et al., 2008; Nickell et al., 2005). The electron optical density of a macromolecule is proportional to its electrostatic potential and the density map can be calculated from the atomic structure by applying a low pass filter at a given resolution. Here, density maps are generated at 4 nm resolution using the PDB2VOL program of the Situs 2.0 package (Wriggers et al., 1999) with voxel length of 1 nm. These initial density maps are then used as samples for simulating electron micrograph images at different tilt angles. In cryoET, the sample is tilted in small increments around a single-axis. At each tilt angle, a simulated micrograph is generated from the sample. We set the tilt angle rotating from −70 to 70^○ with steps of 2^○, which is a typical procedure for experimental tomograms. As a result, our data contains a wedge-shaped region in Fourier space for which no data have been measured (missing wedge effects), similar to experimental measurements. The missing wedge effect leads to distortions of the density maps along the tilt-axis. To generate realistic micrographs, noise is added to the images and the resulting image map is convoluted with a CTF, which describes the imaging in the transmission electron microscope in a linear approximation. Any negative contrast values beyond the 1st zero of the CTF are eliminated. We also consider the MTF of a typical detector used in whole cell tomography, and convolute the density map with the corresponding MTF. The CTF and MTF describe distortions from interactions between electrons and the specimen and distortions due to the image detector (Frank, 2006; Nickell et al., 2005). Typical acquisition parameters that were also used during actual experimental measurements of whole cell tomograms (Beck et al., 2009) were used: voxel grid length=1 nm, the spherical aberration=2×10⁻³ m, the defocus value=−4×10⁻⁶m, the voltage=200 kV, the MTF corresponded to a realistic electron detector (McMullan et al., 2009), defined as sinc(πω/2) where ω is the fraction of the Nyquist frequency.

Finally, we use a backprojection algorithm to generate a tomogram from the individual 2D micrographs that were generated at the various tilt angles (Beck et al., 2009). To test the influence of increasing noise, we add different amount of noise to the images, so that the SNR levels range between 5 and 0.5, respectively (Fig. 4).

2.5 Benchmark set

The classification is tested on simulated tomograms containing 40 complexes of four different classes (Fig. 4). To assess bias arising from specific combinations of complexes, a set of 24 types of complexes of variable sizes were selected from the PDB databank (Fig. 5). From this pool, 52 benchmark sets were generated that each contained four randomly chosen types of complexes. To generate a tomogram for each of these combinations, 10 instances of each complex were randomly oriented and placed into a grid of size 150 nm×150 nm×100 nm (Fig. 4). For two of the 52 benchmark sets (set 1 and set 2), 50 independent density maps were generated each with randomly placed complexes. For each density map, a tomogram was simulated at four different signal-to-noise ratios (SNR=0.5, 0.1, 5, ∞) (Fig. 4) leading to a total of 200 simulated tomograms per set. Set 1 contains GroEL [1KP8], heat-shock protein ACR1 [2BYU], carboxipeptidase [2BO9] and propionyl-CoA carboxylase [1VRG] (PDB ID in squared brackets). Set 2 contains ornithine carbamoyltransferase [1A1S], octomeric enolase[1W6T], RNA polymerase [2GHO] and ClpP [1YG6]). For the remaining 50 benchmark sets of complexes, tomograms were generated at SNR=0.5.

Open in a separate window

Fig. 5.

Average minimum feature vector distance maps. Contour level plot of the average minumum feature vector distance maps for 24 complexes. The value assigned to each grid voxel in a map is the average minimum distance between its feature vector and the feature vectors in the density maps of all the other complexes. The contour plane contains the maximal value. Colors are based on a rainbow scheme with red as the maximum and blue as the minimum values. The PDB ID of each complex structure is also shown.

3.1 Analysis of feature vectors

We now determine the locations of the most unique feature vectors in a complex. For each of the 24 benchmark complexes, a simulated density map is generated at 4 nm resolution (Wriggers et al., 1999) and its rotation-invariant feature vectors are calculated. For each voxel in a complex, the distance is calculated between its feature vector and all the feature vectors in the other complexes. For each of the 24 complexes, the minimum distance is determined. For the given voxel, the average of all minimum distances is a measure of the feature vector's uniqueness with respect to the feature vectors in all the other complexes. The average minimum distance correlates strongly with the position of the voxel from the mass center of the complex (average Pearson's correlation −0.68±0.05). Voxels at the complex center have feature vectors that are most different from those in other complexes and are therefore most discriminative for classification purposes (Fig. 5). Moreover, larger complexes have generally a greater number of unique feature vectors (Fig. 5). These findings provide an important guide for increasing accuracy and computational efficiency of the method. Instead of processing all feature vectors in a tomogram, one should focus only on those vectors that are located at the central regions of a complex. Considering only the most informative feature vectors will not only increase the accuracy of the method but also greatly increase the computational efficiency and reduce the necessary computational memory consumption. Noise and distortions in a map will reduce the number of unique feature vectors, however, also then the central regions are most discriminative for classification. Accordingly, the similarity between the segmented objects in a tomogram is calculated by using only its central core regions.

3.2 Parameter definition

The classification accuracy is increased, if the similarity between segmented objects is calculated from their core regions. Here, the core region is set to be 500 nm³ and is defined by the termination of the adapted watershed algorithm after 500 voxels are selected: S_s^core={v_s,r|r=1,…,r_s^max=500}. The good performance for using a fixed rather than a relative sized core region is presumably a result of the limited size range of soluble complexes in a cell, which ranges typically up to ≈2700 nm³ (ribosome) with a mean of ≈656 nm³ in our sample. A slight variation of the core region size does not change the outcome of our analysis. However, using more extended regions will not improve but often decrease the classification accuracy while reducing computational efficiency significantly.

While the object similarities are determined from the core regions (S^core), the determined class labels are assigned to all voxels in an extended region S_s^ext={v_s,r|r=1,…,r_s^max=2000}.

The authors thank Dr. Shihua Zhang for his valuable suggestions and discussions.

Funding: Human Frontier Science Program (RGY0079/2009-C to F.A.); National Institute of Health (NIH) (1R01GM096089; and 2U54RR022220-06 to F.A.), Alfred P. Sloan Research foundation (to F.A.); F.A. is a Pew Scholar in Biomedical Sciences, supported by the Pew Charitable Trusts.

Conflict of Interest: none declared.