Journal Of Structural Biology - Yale University
Journal of Structural Biology 171 (2010) 256–265Contents lists available at ScienceDirectJournal of Structural Biologyjournal homepage: www.elsevier.com/locate/yjsbiAn adaptive Expectation–Maximization algorithm with GPU implementationfor electron cryomicroscopyHemant D. Tagare a,b, Andrew Barthel b, Fred J. Sigworth c,*aDepartment of Diagnostic Radiology, Yale University, New Haven CT 06520, United StatesDepartment of Biomedical Engineering,Yale University, New Haven CT 06520, United StatescDepartment of Cellular and Molecular Physiology, Yale University, CT 06520, United Statesba r t i c l ei n f oArticle history:Received 18 February 2009Received in revised form 30 May 2010Accepted 2 June 2010Available online 9 June 2010Keywords:Cryo-EMSingle-particle a b s t r a c tMaximum-likelihood (ML) estimation has very desirable properties for reconstructing 3D volumes fromnoisy cryo-EM images of single macromolecular particles. Current implementations of ML estimationmake use of the Expectation–Maximization (EM) algorithm or its variants. However, the EM algorithmis notoriously computation-intensive, as it involves integrals over all orientations and positions for eachparticle image. We present a strategy to speedup the EM algorithm using domain reduction. Domainreduction uses a coarse grid to evaluate regions in the integration domain that contribute most to theintegral. The integral is evaluated with a fine grid in these regions. In the simulations reported in thispaper, domain reduction gives speedups which exceed a factor of 10 in early iterations and which exceeda factor of 60 in terminal iterations.Ó 2010 Elsevier Inc. All rights reserved.1. IntroductionSingle-particle reconstruction is the process by which noisytwo-dimensional images of individual, randomly-oriented macromolecular ‘‘particles” are used to determine one or more threedimensional electron-scattering-density ‘‘maps” of the underlyingmacromolecules. All of the algorithms that perform such reconstructions are iterative. Each iteration begins with a guess of theparticle structure, then aligns the cryo-EM images with the structure, and averages the aligned images to update the structure.The alignment step is especially critical in overcoming noise andimproving resolution.For the alignment step, traditional algorithms choose the bestalignment to the current guess of the structure. This ‘‘best alignment” strategy is simple to implement, but can be troublesomewhen used with noisy images because noisy images can match atwrong alignments. The problem is that the ‘‘best alignment” strategy ignores alignments that are almost as good as – but are not –the best alignment.One class of reconstruction algorithms which does not have thislimitation is similar to the maximum-likelihood (ML) principle.These algorithms do not use the best alignment for structure update, but instead use the Expectation–Maximization (EM) algorithm or variants thereof, which form a weighted average overall alignments for the structure update. The weighted averaging al* Corresponding author. Fax: 1 203 785 4951.E-mail address: fred.sigworth@yale.edu (F.J. Sigworth).1047-8477/ - see front matter Ó 2010 Elsevier Inc. All rights reserved.doi:10.1016/j.jsb.2010.06.004lows the non-best-match alignments to contribute. In cryo-EM, theML principle was proposed for 2D image restoration (Sigworth,1988) and subsequently generalized to the estimation of multipleimage classes (Scheres et al., 2005b) and for the optimization ofunit-cell alignment in images of crystals (Zeng et al., 2007). TheML principle has also been extended to the problem of 3D reconstruction from 2D cryo-EM images by Yin et al. (2001, 2003),Scheres et al. (2005b), Lee et al. (2007), Scheres et al. (2007a,b).Even though it has appealing properties, the EM algorithm has alimitation: it is computationally slow. Calculating the expectationover all alignments (i.e. averaging over all alignments) is expensiveand the EM algorithm can easily take CPU-months to converge. Inthis paper, we report two strategies for speeding up the EM algorithm for cryo-EM. The first strategy is similar to the work of Sanderet al. (2003) and speeds up the EM algorithm by computing theexpectation only over those alignments that contribute significantlyto the weighted average. This is the ‘‘adaptive” part of our algorithm.Our experiments show that it significantly speeds up the EM algorithm without losing accuracy. The second strategy is the use ofgraphics processing units (GPUs) to accelerate the most computation-intensive parts of the calculations (Castano-Diez et al., 2008).To put previous attempts to speed up the EM algorithm in context, we first explain the basic EM iteration. The EM iteration hastwo steps, in the first step the latent data probability is calculatedand in the second step this probability is used for calculating theweighted averages. Both calculations involve integration over allpossible alignments and contribute equally to the computationalcomplexity.
257H.D. Tagare et al. / Journal of Structural Biology 171 (2010) 256–265Previous attempts to speed up the EM algorithm can be classified into two groups. The first group of algorithms uses sphericalharmonics as basis functions for the structure (Yin et al., 2001;Yin et al., 2003; Lee et al., 2007). These algorithms have the computational advantage that integrations with respect to rotationscan be calculated in closed form. However, integration with respectto translations still need to calculated numerically. The secondgroup contains the algorithm of Scheres et al. (2005b) and is similar in its approach to our algorithm. This algorithm decreases computational cost by calculating the EM integrals over a smallerintegration domain. Scheres et al.’s strategy is quite complex. Thefirst EM iteration is carried out in its entirety over alignments.For all subsequent iterations, every image is compared with everyclass mean using the most significant translations from the previous iteration for this pair. Using these translations, the latent probability that the image comes from that class is calculated for everyrotation. All rotations for which the calculated probability exceedsa fraction of the maximum value of all calculated probabilities areretained. The EM integration is carried out over alignments formedby the surviving rotations and all translations.One problem with this strategy is that it is not entirely clearwhether thresholding the probability to reduce the domain givesgood approximations to the integral. If the domain of integrationis to be reduced, then the reduction strategy ought to be basedon how changing the domain affects the integral rather than onthresholding a function. Our algorithm includes such a strategy,and is inspired by adaptive integration techniques in numericalanalysis (Atkinson, 1978). We first use a coarse grid to estimatethe contribution to the integral at every alignment. Then, we retainthe smallest set of alignments over which the integral contributes afraction (e.g. 0.999) of the net integral.2. The ML formulationSuppose that l is a projection of a structure along a specificdirection. An observed cryo-EM image I is l with additive noiseand a random in-plane rotation and shift. Letting T s representthe in-plane image rotation and translation operator, wheres ¼ ðh; tx ; ty Þ is the rotation and translation, the observed image isI ¼ T s ðl þ nÞ, or, T s ðIÞ ¼ l þ n, where, n is additive white Gaussian noise. Thus the probability density function (pdf) of observingan image I ispðIjl; rÞ ¼Zpg ðT s ðIÞjl; rÞpðsÞds;ð1ÞXwhere, pðsÞ is the density of s; X is the support of pðsÞ, andpg ðT s ðIÞjl; rÞ ¼1 expð2pr2 ÞP 2! kT s ðIÞ lk2;2r2ð2Þwhere, kk is the usual Euclidean norm, P is the number of pixels inthe image, and r2 is the noise variance of each pixel. The supportX ¼ ½0; 2pÞ ½ t max ; tmax ½ t max ; t max for some maximum translation t max .A small aside to emphasize an important point – typically in Eq.(2) the value of kT s ðIÞ lk2 is several orders of magnitude largerthan the value of r2 (the former is the pixel-wise squared difference summed over the entire image while the latter is the noisevariance of a single pixel). Because small changes in kT s ðIÞ lk2are vastly amplified by the exponentiation, small changes inkT s ðIÞ lk2 can cause large changes in the value of pg .Cryo-EM obtains images from unknown random projectiondirections. It is common to model this phenomenon as follows:Letlj ; j ¼ 1; . . . ; M be the projection of the particle along the jth direction. Then, assuming that images are random draws from one ofthe projections, the pdf of an image I ispðIjl1 ; . . . ; lM ; r; a1 ; . . . ; aM Þ ¼MXaj pðIjlj ; rÞ;ð3Þj¼1where, the coefficients aj are non-negative and sum to 1, and thedensities pðIjlj ; rÞ are given by Eq. (1) with l ¼ lj . This probabilitydensity model is popularly called a mixture model (McLachlan andPeel, 2000). The densities pðIjlj ; rÞ are called class densities andthe coefficients aj are called mixture coefficients. The means lj arecalled class means.Suppose that N images Ik ; k ¼ 1; . . . ; N are obtained in this way.Then, the joint density of the images ispðI1 ; . . . ; IN jl1 ; . . . ; lM ; r; a1 ; . . . ; aM Þ¼NYpðIk jl1 ; . . . ; lM ; r; a1 ; . . . ; aM Þ:ð4Þk¼1A maximum-likelihood estimate ofgiven by 1 ; . . . ; a M g ¼ arg 1; . . . ; l M; r ; afll1 ; . . . ; lM ; r; a1 ; . . . ; aM ismaxl1 ;.;lM ;r;a1 ;.;aM log pðI1 ; . . . ; IN jl1 ; . . . ; lM ; r; a1 ; . . . ; aM Þ:ð5Þ2.1. The EM algorithmThe EM algorithm iteratively converges to the maximum-likelihood estimate of Eq. (5). The EM algorithm and its application tocryo-EM is not new. But we present some its details in order to motivate the adaptive-EM algorithm. The EM algorithm works by introducing additional random variables, called latent variables. For themaximum-likelihood problem of Eq. (5) the EM algorithm introduces two latent variables per image. Together they denote the component of and the transformation for the kth image. For the kth imagethe latent random variables are denoted yk and sk . The variable yk is adiscrete valued random variable taking values yk 2 f1; . . . ; Mg. Theevent yk ¼ j indicates that the kth image comes from the jth component. The random variable sk ¼ ðhk ; tx;k ; ty;k Þ takes values sk 2 X.Because yk takes values in yk 2 f1; . . . ; Mg and sk in X, the jointrandom variable ðyk ; sk Þ takes values in f1; . . . ; Mg X. In Fig. 1a,we show f1; . . . ; Mg X as a column of M X domains. There areN such columns corresponding to N images. Below we will need to refer to the M N copies of the domain X in Fig. 1a individually and collectively. We will refer to the domain in the jth row and kth column asXjk . We will collectively refer to all domains as f, i.e. f ¼ [j;k Xjk .The EM iterations proceed as given below. The superscripts½n 1 and ½n refer to the values of the variables in the n 1stand nth iteration:2.1.1. The EM algorithm½0 ½0 1. Initialize: Set n ¼ 0 and initialize aj ; lj ; r½0 for j ¼ 1; . . . ; M.2. Start Iteration: Set n ¼ n þ 1.3. Calculate Latent Probabilities: For all points in Xjk calculatepðyk ¼ j; sk jIk ; h½n 1 Þ½n 1 a½n 1 pg ðT sk ðIk Þjlj ; r½n 1 Þpðsk Þj:R½n 1 ½n 1 ; r½n 1 Þpðsk Þdski¼1 aiXik pg ðT sk ðI k Þjlið6Þpðy ¼ j; sk jIk ; h½n 1 Þ ðyk ¼ j; sk jIk ; h½n 1 Þ ¼ PN R k:p½n 1 Þdsii¼1 Xji pðyi ¼ j; si jI i ; hð7Þ¼ PMandThe variable i in the both denominators of the above formulae is asummation variable. In Eq. (6) the variable i sums over rows, andin Eq. (7) the variable i sums over columns.
258H.D. Tagare et al. / Journal of Structural Biology 171 (2010) 256–265N images (k)RownormalizationVertex vat center ofthe cube.M means (j)11Cube C(v)jktxtyMNColumnnormalizationFig. 1. Structures for the EM algorithm.½n 4. Update Parameters: The parameter ajprobability density p of Eq. (6)½n jais updated using the3.1. Speeding up the EM algorithmN Z1X¼pðyk ¼ j; sk jIk ; h½n 1 Þdsk :N k¼1 Xjkl½n j ¼ð8Þ½n of Eq.l½n are updated using the density pj ;rThe parameters(7):N ZXXjkk¼1 ðyk ¼ j; sk jIk ; h½n 1 Þdsk ;T sk ðIk Þpð9Þandðr2 Þ½n ¼M XN Z1 X½n 1 ðykkT sk ðIk Þ lj k2 pMP j¼1 k¼1 Xjk¼ j; sk jIk ; h½n 1 Þdsk :Note that the updates ofð10Þl and r are weighted averages of T sk ðIk Þ½n 1 2and kT sk ðIk Þ lj as the weight.k with p½n ½n 5. Loop: Stop if aj ; lj ; r½n have converged. Else go to step 2.A useful visualization of the EM iterations is presented inFig. 1a:½n 1 21. Begin by calculating kT sk ðIk Þ ljk in Xjk .½n 1 3. Adaptive EM; r½n 1 Þpðsk Þ in every Xjk .2. From this calculate ak pg ðT sk ðIk ÞjljNormalize this function column wise in f so that the net integral of the function in each column is 1. This is illustrated inFig. 1a and referred to as column normalization. The normalizedfunction is the latent data probability of Eq. (6). The denominator in Eq. (6) achieves the normalization.3. Next, normalize the function obtained in the above step so thatits net integral in each row of f is 1. This is also illustrated inFig. 1a and called row normalization. The normalization isachieved by the denominator in Eq. (7). The result of row nor .malization is the function p4. Update lj and r2 according to Eqs. (9) and (10) using weighted as the weight.averages with pThe EM iteration described above is computationally expensive.Eqs. (6)–(10) require integration over the domains Xjk . These integrals are not available in closed form and have to be evaluatednumerically. To do this, we introduce a grid in every Xjk (seeFig. 1b) and approximate the integrals with a Riemann sum of theintegrand over the vertices of the grid. The grid consists of cubesof size Dh Dt x Dty . The center of each cube is a vertex of the grid.We use v to refer to a vertex, and Cðv Þ to refer to its cube.The computationally expensive part of the algorithm is the cal½n 1 2culation of the kT sk ðIk Þ ljk term at each vertex of the grid.We employ two strategies to speed up the EM iteration. Bothstrategies require various formulae in Eqs. (6)–(10) to be approximated. The first strategy is called domain reduction and it replacesthe integration over Xjk with integration over smaller subsets. Thesecond strategy is called grid interpolation and uses two grids – acoarse grid for estimating the reduced domain, and a fine grid inthe reduced domain for calculating the formulae. Analogous strategies have been used by Sander et al. (2003) for angular search in aconventional reconstruction algorithm. The difference is that ourstrategy is for approximating an integral whereas Sander et. al.’sstrategy is for approximating a maximum-seeking search.Both strategies are based on observations which are illustratedin Fig. 2. The (a) part of the figure shows an image which was obtained by projecting a ribosome structure in two different directions, summing the projections, adding noise, and blurring. Thisimage is similar to an estimated class meanl½n 1 in early iterajtions. The (b) part shows a single projection with white noiseadded to it. This image is similar to the observed images. The (c)½n 1 2part shows kT sk ðIk Þ ljk as a function of rotation (transla-tions are held fixed for simplicity) with the image in part (a) as ,l½n 1 and the image in part (b) as Ik . The (d) part shows how pj½n 1 2which is calculated from kT sk ðIk Þ lj has strong spikes.The function pFig. 2 illustrates several effects:k , behaves with rotation. has more than one peak. A peak-seeking align1. The function pment method would align Ik at the strongest of these peaks.But it is not clear that just using the strongest peak is the best
H.D. Tagare et al. / Journal of Structural Biology 171 (2010) 256–265259Fig. 2. Domain reduction and grid interpolation.alignment decision, especially since it is the best alignment to½n 1 an unconverged mean lj.On the other hand, recall from Eqs. (9) and (10) that the EM . Thus, all peaks contribalgorithm works by using all values of pute in the EM algorithm, and it avoids the premature decision ofusing a single peak. This shows why the EM algorithm is preferable to a peak-seeking alignment algorithm. is essentially2. In spite of multiple peaks, Fig. 2d suggests that pzero over a significant part of Xjk . If the numerical integration concan be restricted to only that part of the domain where ptributes significantly, then considerable computational gaincan be made without losing too much accuracy. This is themotivation behind domain reduction.½n 1 23. Note that kT sk ðIk Þ ljk is a much smoother function than p. This suggests a way to estimate the reduced domains: introduce a coarse grid in the domains Xkj , calculate½n 1 2kT sk ðIk Þ ljk at the vertices of the coarse grid, interpolateusing the vertex values and use the interpolated function to and the reduced domains. This is grid interpolation.estimate p½n 1 2The solid line in Fig. 2c shows kT sk ðIk Þ ljk evaluated on a
260H.D. Tagare et al. / Journal of Structural Biology 171 (2010) 256–265fine grid (spacing 1 ). The dashed line in Fig. 2c shows½n 1 2kT sk ðIk Þ ljk evaluated on a coarse grid (spacing 12 ) ’s are shown inand interpolated by a B-spline. The resulting p calculated from the B-spline gives anFig. 2d. Note that the p and the figure suggests that the reducedapproximation of pdomain can be calculated from this approximation.The procedure that obtains the reduced domain is describedbelow in detail. The result of using that procedure gives thereduced domain shown in Fig. 2d (the parameter f is explainedbelow).We now describe both strategies in detail:3.2. Domain reductionThe idea in domain reduction is to replace integration over Xjkwith integration over a smaller subset X jk . The X jk s are estimated .to contain a significant fraction f of the probability mass of pThe fraction of the probability mass is equal to a parameter fwhose value is user-chosen but constrained to be 0 f 1. Because we want to retain the averaging properties of the EM algorithm, we set the value of f close to 1, e.g. f ¼ 0:999. Eventhough the f is very close to 1, the domain is reduced considerably is spiky.because p is obTo describe more precisely how X jk are found, recall that ptained by column and row normalization as shown in Fig. 1. Thus sums to 1 in each row of f. Let Xj ¼ [k Xjk bethe integral of pthe union of all Xjk in a row, then row normalization implies thatR Xj pds ¼ 1. Let O Xj be any subset of Xj , thenZ dspO in O. Let O be the subset of Xjmeasures the probability ‘‘mass” of pwith the smallest volume that has a probability mass of f:O ¼ argmin volðOÞ subject to the constraints O Xj ;andZO ds ¼ f:pOð11Þ jk Then, the reduced domain is X ¼ O \ Xjk .The reduced domain O is easy to find: Let T be a threshold and are greater than orXTof Xj in which the values of pj be the subset R ds is a monotonically decreasing function ofequal to T. Then XT pjT and its range of values is ½0; 1 . Thus there is a unique T forR ds ¼ f. It is straight forward to show that for this Twhich XT pj the set XTj equals the set O of Eq. (11).The algorithm for domain reduction is a binary-search algorithm for the appropriate T is:3.2.1. The domain reduction algorithm1. Initialize: Set the upper and lower limit of T to 0 and 1respectively.2. Binary search: Carry out a binary search in within ½0; 1 for theT that solves:Z ds ¼ f:pXTj3. Calculate theX jk ¼ XTj \ Xjk .domains:Returnthereduceddomains3.3. Grid interpolationGrid interpolation uses two grids, a coarse grid Gc and a fine gridGf . The two grids are chosen such that any cube of the coarse gridcontains an integer number of the cubes of the fine grid. The twogrids are used to estimate the reduced domain as follows:3.3.1. Domain reduction with grid interpolation1. Coarse calculation: Calculate kT sk ðIk Þ ljk2 at the vertices ofthe coarse grid Gc .2. Interpolate:
Likelihood Expectation–Maximization abstract Maximum-likelihood (ML) estimation has very desirable properties for reconstructing 3D volumes from noisy cryo-EM images of single macromolecular particles. Current implementations of ML estimation make use of the Expectation–Maximization
Yale Yalelift 360 Yale VS /// CM Series 622 Yale Yalelift ITP/ITG Yale Yalelift LHP/LHG Coffing LHH Table of Contents 1-2 3-4 5 6-7 8-9 10-11 Hand Chain Hoist Ratchet Lever Hoist Trolleys Beam Clamps Specialities 12 13 14-15 16 17 Yale Yalehandy Yale UNOplus Yale C 85 and D 85 Yale D 95 Yale AL 18 19-20 CM CB
animation, biology articles, biology ask your doubts, biology at a glance, biology basics, biology books, biology books for pmt, biology botany, biology branches, biology by campbell, biology class 11th, biology coaching, biology coaching in delhi, biology concepts, biology diagrams, biology
the Yale Access App The Yale Access App is needed to control your Yale products from your mobile device. The Yale Access App is available for iPhone and Android. Download the Yale Access App from the App Store or Google Play depending on your device. Once you have downloaded the Yale Access App, log in to
1 Yale University Press, “Books from Yale at Christmas” (print advertisement), Yale Daily News, 10 December 1958, 9. 2 Gideon Gordon, “From the Couch: Study of Yale Mind Published,” Yale Daily News, 22 October 1958, 1. 3 Bryant M. Wedge, preface to Psychosocial Problems of College M
Yale Law School 2019-2020 Yale Law School 2019-2020. BULLETIN OF YALE UNIVERSITY Series 115 Number 11 August 10, 2019 (USPS 078-500) is published seventeen times a year (one time in May and October; three times in June and September; four times in July; five times in August) by Yale University, 2 Whitney
Yale Smart Lift can reduce cycle times by up to 25%. Exclusive Yale Smart Slow Down To further ensure that every load remains stable, the Yale MPB-VG truck features optional Yale Smart Slow Down technology. When the operator turns the truck to change direction, the Yale Smart Slow Down feature intuitively reduces the truck's speed,
DAT Study Tips* Biology Materials: DAT Destroyer, Feralis Biology Notes, Cliff's AP Bio 3rd Edition, DAT Bootcamp (Both Cliff’s AP Bio and Feralis Notes are free online) Biology is one of the most time consuming sections to study for, given that the scope of the material covered in DAT biology is so randomly big. Cliff's AP Bio 3rdFile Size: 527KBPage Count: 9Explore furtherDAT Bootcamp Biology Flashcards Quizletquizlet.comHow to Study for the DAT Biology Section the Right Way .datbootcamp.comFeralis Biology Notes DAT Study Tips Free Downloadferalisnotes.comFeralis Biology Notes? Student Doctor Network Communitiesforums.studentdoctor.netBiology Cumulative Exam Flashcards Quizletquizlet.comRecommended to you b
4 n Yale Nursing THE PATH STRATEGIC VISION FOR A RESILIENT FUTURE FORWARD Education Health Systems Partnerships Science 4 Download the full Strategic Blueprint: nursing.yale.edu/blueprint Like other graduate and professional schools at Yale University, Yale School of Nursing (YSN) maintain
