Cryo electron microscopy: 3-D reconstruction of stochastic objects having statistics with geometric symmetry from 2-D noisy conditionally-linear imaging

2.6 MB·

by  Nan Xu

1
Cryo electron microscopy: 3-D reconstruction of
stochastic objects having statistics with geometric
symmetry from 2-D noisy conditionally-linear
imaging
Nan Xu, Student Member and Peter C. Doerschuk, Senior Member
7
Abstract—Theimagesobtainedfromcryoelectronmicroscopy following a single instance of the particle through the 4-D
1
aretheprojectionofmanyheterogeneousinstancesoftheobject trajectory that it would experience in its natural environment.
0
under study (e.g., a virus). When the object is heterogeneous
2 Single-particlecryoelectronmicroscopy(cryoEM)istrans-
but has an overall geometric symmetry, it is natural to describe
forming structural biology [4] because, due to better image
n the object as stochastic with symmetrical statistics. A maximum
a likelihood reconstruction approach is described in which all detectors [5], phase plates [6] for improved image contrast,
J instances of the object lack symmetry but the statistics of the andmoresophisticatedstatisticalcomputing[7]–[9],cryoEM
5 ensemble of objects are constrained to have symmetry. This is beginning to provide atomic-resolution structures similar
algorithm is demonstrated on bacteriophage HK97 and is con- to those of x-ray crystallography without the need to create
] trasted with existing algorithms in which each object, while still a crystal. This is a very important advance because many
P heterogeneous, has the symmetry. Reconstruction results show
A that the proposed algorithm eliminates long-standing distortions particles simply do not form crystals for a variety of biophys-
in previous heterogeneity calculations associated with symmetry ical reasons. However, more important for understanding the
.
t axes, and provides estimates that make more biological sense functioning of the particles, cryo EM has the capability of
a
than the estimates of existing algorithms.
t displaying heterogeneity and dynamics among the particles.
s
[ Index Terms—cryo electron microscopy, virus structural bi- In a cryo EM experiment, typically 104–106 noisy 2-D
ology, 3-D signal reconstruction from projections, maximum projectionimagesofthe3-Delectronscatteringintensityofthe
1 likelihood estimation, random signal with symmetrical statistics
flash-frozenparticlesarerecorded,oneimageforeachinstance
v
of the particle. Typical particle weights are in the range of
6
0 I. INTRODUCTION 105–108Daltons.Theexperimentaldesignisdominatedbythe
2 THE goal of structural biology is to determine the 3-D desiretoavoiddamagetotheparticlebytheelectronbeam.In
1 particular,damageminimizationisachievedbylowdose(lead-
geometry of a biological macromolecular complex. Such
0
ing to noisy images), by recording only one projection image
. complexes,forexample,avirus,functionasmachinesandare
1 per particle (unlike medical x-ray computed tomography),
typicallycalled“particles”.The3-Dgeometry,typicallycalled
0
and by performing the experiment on flash-frozen particles at
7 the “structure”, can be as detailed as the 3-D location of each
cryogenictemperatures.Becausetheorientationoftheparticle
1 atom in the particle or, if only lower resolution is achievable,
: the electron density as a function of position in 3-D (in x- duringtheflashfreezingisnotcontrolled,theorientationofthe
v
projectionimageisrandomandbecausetheimageisnoisy,the
i ray crystallography) or the electron scattering intensity as a
X projection direction cannot be determined from an individual
function of position in 3-D (in electron microscopy). Under-
r standing the 3-D geometry of synthetic nanoparticles is also image alone. Standard computational methods assume that all
a instances of the particle are identical or that all instances of
becoming important in materials science [1]–[3]. Structural
the particle belong to one of a few discrete classes (“discrete
biology has long been dominated by x-ray crystallography
heterogeneity”, the class label for each image is not known,
which has very limited capability for dealing with hetero-
typicallytherearelessthantenclasses)andthateveryinstance
geneity and dynamics among the particles. In particular, the
within a class is identical. Then a structure is computed for
particles making up the crystal must be sufficiently uniform
each class, possibly at atomic resolution. More sophisticated
to allow the formation of a crystal and then the additional
approaches characterize the continuous heterogeneity of the
non-biological interactions between neighboring particles in
particles within each class.
the crystal further restrict the individual particles. Finally the
The geometric symmetry of the particle is an important
data is a function of the entire crystal, i.e., a function of all
aspect of structural studies of both biological and synthetic
of the particles. Therefore x-ray crystallography provides a
particles. Many so-called spherical viruses, which include
3-D snapshot of a particular state of the particle rather than
poliovirusandotherhumanpathogens,haveaso-calledcapsid
N.XuiswiththeSchoolofElectricalandComputerEngineering,Cornell that is constructed out of many copies of a peptide building
University,Ithaca,NY,14853USAe-mail:[email protected] block which are arranged in pattern that has icosahedral
P.C.DoerschukiswiththeSchoolofElectricalandComputerEngineering
symmetrymuchlikeahuman-madegeodesicdome.Symmetry
and the Nancy E. and Peter C. Meinig School of Biomedical Engineering,
CornellUniversity,Ithaca,NY,14853USAe-mail:[email protected] is described mathematically as a rotational point group. This
2
paper describes how to impose symmetry on the statistics of and near symmetry axes of the particle [31, p. 173]. Since the
the particle rather than on each realization of the particle symmetry axes can be the location of important biology, e.g.,
(Sections II, III, IV, and VII) and how to integrate this inFlockHouseVirus(FHV)[32],[33],thebiologyconcerning
mathematicalmodeloftheparticleintoamaximumlikelihood the pentamer of γ peptides near the 5-fold symmetry axis of
reconstructionmethod(SectionV).Theresultisdemonstrated the icosahedrally symmetric particle [34], [35], so anomalous
on images of a Virus Like Particle of bacteriophage HK97 results near symmetry axes is not a satisfactory situation.
(Section VI). (3) Assume that each particle in a class is different, that
The approach of this paper is focused on the estimation no particle has the symmetry, but the statistics of the i.i.d.
of first and second order statistics of the electron scattering ensemble of particles has the symmetry as is described in the
intensity of the particles in order to characterize the het- following paragraph. The authors are unaware of prior work
erogeneity of the particles and related work is discussed in of this type.
Section II. However there are other approaches to character- Oneapproachtointroducesymmetricstatisticsistorequire
izing the differences among the particles visualized by cryo symmetry in all of the finite-dimensional probability mea-
EM. Prominent approaches include direct reconstruction of sures that together, via Kolmogorov’s extension theorem [36,
principal components [10], low-dimensional manifold embed- Theorem 2.1.5 p. 11], define the electron scattering intensity
ding of multiple reconstructions with interpolation to create stochasticprocess.Iftheprobabilitymeasurescanbedescribed
trajectories [11], [12], and normal modes of a references byprobabilitydensityfunctions(pdfs),thenthesymmetrycon-
model [13]. These papers do not concern symmetry (the dition is p (ρ ,...,ρ ) = p (ρ ,...,ρ )
numerical examples are all ribosome problems) and so are for all β ∈Rβ−{11x,1.,.....,R,Nβ−1x}k, k1∈ {1,2k,...},x1x,...,,.x.k.,x1 ∈ IRk3,
g 1 k
outside of the scope of this paper.
and ρ ,...,ρ ∈ IR. Define the mean function and covari-
1 k
Notation: E denotes expectation, i.i.d. stands for “indepen-
ance function of the electron scattering intensity, denoted by
dent and identically distributed”, x∈IR3, and x=(cid:107)x(cid:107) . The
2 ρ¯(x) and C(x ,x ), respectively, by ρ¯(x) = E[ρ(x)] and
quantity x/x is a unit vector and, since x∈IR3, is shorthand 1 2
C(x ,x )=E[(ρ(x )−ρ¯(x ))(ρ(x )−ρ¯(x ))],respectively.
1 2 1 1 2 2
for the two angles of spherical coordinates. I is the n×n
n Fromthek =1andk =2instancesofthesymmetrycondition
(cid:82)
identity matrix and 0 is the m × n zero matrix. dΩ
m,n for the pdfs, it is straightforward to compute that
is integration over the surface of the sphere in 3-D which,
in spherical coordinates, is (cid:82)π (cid:82)2π sin(θ)dθdφ. Define the ρ¯(R−1x)=ρ¯(x) (1)
θ=0 φ=0 β
Kroneckerdeltafunctionδm,n byδm,n =1ifm=nand=0 for any x∈IR3 and β ∈{1,...,N },
g
otherwise where m and n are integers or vectors of integers.
C(R−1x ,R−1x )=C(x ,x ) (2)
β 1 β 2 1 2
for any x ,x ∈IR3 and β ∈{1,...,N }.
II. SYMMETRYANDSYMMETRICALSTATISTICS 1 2 g
ThesymmetriesthatarerelevanttothecryoEMapplication If the electron scattering intensity is assumed to be Gaussian,
are the 3-D rotational point group symmetries which are the the benefits of which are described in Ref. [25] and outlined
icosahedral I, octahedral O, tetrahedral T, cyclic C (n ∈ inSectionV,thenEqs.1–2implythatthesymmetrycondition
n
{1,2,...}), and dihedral D (n ∈ {1,2,...}) groups. Once for the pdfs is satisfied for all values of k.
n
a coordinate system is chosen, each of these groups can be
described as a collection of matrices denoted by R ∈IR3×3
β III. BASISFUNCTIONS
(β ∈ {1,...,N }) which are rotation matrices, i.e., R−1 =
g β
RT and detR =+1. The value of N is 60, 24, 12, n, and The particle is described as a linear combination of basis
β β g
2n for I, O, T, C , and D , respectively. functions, specifically,
n n
The symmetry can be applied in three different ways. (cid:88)
ρ(x)= c F (x), (3)
(1) Assume that the particles in a class are identical and ω ω
that the common electron scattering intensity, which is de- ω
noted by ρ(x), has the symmetry, i.e., ρ(R−1x) = ρ(x) where c are random variables and F (x) are known func-
β ω ω
for all β ∈ {1,...,N } and x ∈ IR3. This is available in tions. The goal of the reconstruction calculation (Section V)
g
essentially all software systems, e.g., Refs. [7]–[9], [14]–[16]. is to estimate the joint probability distribution of the c . To
ω
(2) Assume that each particle in a class is different but that achieve the desired constraints on the first and second order
each particle has the symmetry so that the electron scattering statisticsofρ(x)(Eqs.1and2)requiresconstraintsonthefirst
intensityoftheithparticle,whichisdenotedbyρ (x),satisfies andsecondorderstatisticsofc .Determiningtheseconstraints
i ω
ρ (R−1x)=ρ (x) for all β ∈{1,...,N } and x∈IR3. This is simplified if F (x) have specific properties under rotations
i β i g ω
is available in some software systems by resampling [17]– from the symmetry point group of the particle.
[21] (implemented, for instance, in Ref. [7]), moment esti- ThenumericalexampleofSectionVIconcernsaviruswith
mators [22]–[24], or maximum likelihood estimation of the icosahedral (I) symmetry. For I (and also for O and T), the
parameters of a probabilistic model for the Fourier series rotationalsymmetryaxesareorientedinavarietyofdirections
coefficients of the electron scattering intensity [25]–[30]. It leading to particles that tend to fill their minimum bounding
is known that the resulting spatial variance functions (i.e., the sphere.ThereforethechoiceofF (x)thatisconsideredinde-
ω
variance as a function of position) have anomalous peaks on tailusessphericalcoordinatesandistheproductofaso-called
3
angularfunction(whichdependsontheanglesofsphericalco- where T is transpose not Hermitian transpose and R is
g
ordinatesandisalinearcombinationofsphericalharmonics1) the3×3rotationmatrixcorrespondingtothegthelement
and a so-called radial function (which depends on the radius of the group.
of spherical coordinates), i.e., F (x)=I (x/x)ψ (x) where Examples of I functions for l = 15 are visualized in
ω ω ω p,j;l,n
x = (cid:107)x(cid:107) and ω is an index. For C and possibly D , it is Figure 1. Note that p=1 exhibits all of the symmetries of an
2 n n
possible to use the same framework but, especially for C , icosahedron.
n
it might be more natural to use cylindrical coordinates and a
product of a so-called axial function (in z), a polar function
(in φ), and a radial function (in r). Because rotations of x do
not effect x=(cid:107)x(cid:107), in this section we focus on I (x/x).
ω
Ref.[38]derivesabasisI (x/x)thathaspreciseproperties
ω
under the operations of a rotational point group. In order
to describe the properties of the basis, it is necessary to
describethetheideaofrepresentations.Arepresentation(rep) p=1 p=2
of a rotational point group, which is a finite group, is a
set of matrices which are homomorphic under matrix mul-
tiplication to the group elements [39, p. 61]. The irreducible
representations (irreps) of a finite group are a set of unitary
reps such that all unitary reps can be decomposed into a
direct sum of the irreps [39, Theorem 3.3 p. 69]. For the
icosahedral and octahedral cases, all irreps can be chosen to p=3 p=4 p=5
bereal-valuedorthonormalratherthancomplex-valuedunitary
and, because of the biological importance of the icosahedral Fig. 1: An icosahedron with a 5-fold axis pointing out of
case, the remainder of this section focuses on the real-valued the page and example angular basis functions in the same
orthonormal icosahedral case. Let Nrep be the number of orientation with l=15 and p∈{1,...Nrep}. The surfaces of
irreps of the Ng-order group G. Let the set of real-valued 3-D objects defined by ξ(x)=1 for x≤κ1+κ2Ip;l,n(x/x)
matrices in the pth irrep be denoted by Γp(g) ∈ Rdp×dp for and 0 otherwise, where κ1 and κ2 are chosen so that κ1 +
all g ∈ G where p ∈ {1,...,Nrep}. For the icosahedral κ2Ip;l,n(x/x) varies between 0.5 and 1.0, are visualized by
group, N = 5, N = 60, and d = 1,3,3,4,5 for p = UCSF Chimera [41] where the color indicates the distance
rep g p
1,2,3,4,5, respectively. The basis functions I (x/x) were from the center of the object.
p,ω
derived [38] by applying real-valued generalized projection
operators constructed from real-valued irreps to real-valued
IV. CONSTRAINTSONTHEMOMENTSOFTHEWEIGHTS
sphericalharmonicsresultingintheindexω beingtheordered
CAUSEDBYSYMMETRICSTATISTICS
pair (l,n) where l indexes the subspace defined by spherical
In this section we suppress the index η which describes the
harmonics of order l which contains the basis function and n
class of the particle. The electron scattering intensity, denoted
indexes the different basis functions within the subspace. The
byρ(x),isdescribedasalinearcombinationofbasisfunctions
basis functions I (x/x) have the properties
p;l,n
wherethebasisfunctionsaredenotedbyI (x/x)(pindexes
p,ζ
1) Each I is a d -dimensional real-valued vector func-
p;l,n p theirrep,ζ isashorthandforallotherindices,andx=(cid:107)x(cid:107) )
tion, i.e., Ip;l,n ∈Rdp. and take values in IRdp (d is the dimension of the pth irrep2)
p
2) TheI functionsareorthonormalonthesurfaceofthe
p;l,n and the weights in the linear combination are denoted by
sphere, i.e., c (x), also take values in IRdp, and depend on x=(cid:107)x(cid:107) . In
p,ζ 2
(cid:90) summary,
I (x/x)IT (x/x)dΩ=δ δ δ I (4)
p,l,n p(cid:48),l(cid:48),n(cid:48) p,p(cid:48) l,l(cid:48) n,n(cid:48) dp Ng
(cid:88)(cid:88)
ρ(x)= cT (x)I (x/x). (6)
p,ζ p,ζ
3) The subspace of square integrable functions on the sur-
p=1 ζ
faceofthespheredefinedbysphericalharmonicsofindex
The electron scattering intensity ρ(x) is real-valued. In this
l, contains a set of I functions with a total of 2l+1
p;l,n section we deal with the situation where there are a complete
components.
set of real-valued orthonormal irreps (including the important
4) EachI functionhasaspecifictransformationproperty
p;l,n icosahedral case) so that I (x/x) can be taken as real and
p,ζ
underrotationsfromtherotationalpointgroup[40,p.20],
thereforec canalsobetakenasreal.InSectionVIIthecase
p,ζ
in particular,
where I (x/x) and c are complex-valued is described.
p,ζ p,ζ
The constraint on the mean function of c (x) is derived
I (R−1x/x)=(Γp(g))TI (x/x) (5) p,ζ
p;l,n g p;l,n inthisparagraph.Defineρ¯(x)=E[ρ(x)]∈IR1 andc¯ (x)=
p,ζ
E[c (x)]∈IRdp. From Eq. 6,
p,ζ
l ∈1SNph,ethriecaolrhdaerrmmoniscastiasrfieedsemnot∈ed{b−ylY,l.,.m.(,θl},φa)nwdh(eθr,eφth)eadreegtrheeealnsgatliessfieosf (cid:88)Ng (cid:88)
sphericalcoordinateswith0≤θ≤π and0≤φ≤2π [37,Section14.30, ρ¯(x)= c¯Tp,ζ(x)Ip,ζ(x/x) (7)
pp.378–379].
p=1 ζ
4
and from Eqs. 6 and 5 (g ∈{1,...,N }) From Eq. 15 it follows that
g
Ng N(cid:88)rep (cid:88)N(cid:88)rep (cid:88)
ρ¯(Rg−1x)=(cid:88)(cid:88)c¯Tp,ζ(x)(Γp(g))TIp,ζ(x/x). (8) C(x1,x2)= IpT1,ζ1(x1/x1) (16)
p=1 ζ p1=1 ζ1 p2=1 ζ2
×C (x ,x )I (x /x ).
From Eq. 1 it follows that the left hand sides of Eqs. 7 and 8 p1,ζ1;p2,ζ2 1 2 p2,ζ2 2 2
are equal so that for all g ∈{1,...,Ng}, and from Eqs. 15 and 5 (g ∈{1,...,Ng})
(cid:88)pN=g1(cid:88)ζ c¯Tp,ζ(x)(cid:2)Idp −(Γp(g))T(cid:3)Ip,ζ(x/x)=0. (9) C(Rg−1x1,Rg−1x2)=pN(cid:88)1r=ep1(cid:88)ζ1 pN(cid:88)2r=ep1(cid:88)ζ2 IpT1,ζ1(x1/x1) (17)
×Γp1(g)C (x ,x )[Γp2(g)]TI (x /x ).
p1,ζ1;p2,ζ2 1 2 p2,ζ2 2 2
Multiplying by IT (x/x) and integrating over the surface
p(cid:48),ζ(cid:48) FromEq.1itfollowsthatthelefthandsidesofEqs.16and17
of the sphere implies [via the orthogonality of the I (x/x)
p,ζ are equal so that for all g ∈{1,...,N },
functions (Eq. 4)] that for all g ∈{1,...,N }, g
g
(cid:104) (cid:105) N(cid:88)rep(cid:88)N(cid:88)rep (cid:88) (cid:20)
c¯Tp(cid:48),ζ(cid:48)(x) Idp(cid:48) −(Γp(cid:48)(g))T =01,dp(cid:48). (10) IpT1,ζ1(x1/x1) Cp1,ζ1;p2,ζ2(x1,x2) (18)
p1=1 ζ1 p2=1 ζ2
Taking transposes and renaming p(cid:48) to p and ζ(cid:48) to ζ gives that (cid:21)
for all g ∈{1,...,N }, −Γp1(g)Cp1,ζ1;p2,ζ2(x1,x2)[Γp2(g)]T Ip2,ζ2(x2/x2)=0.
g
Multiplying on the left by I∗ (x /x ) and on the right by
(cid:2)Idp −Γp(g)(cid:3)c¯p,ζ(x)=0dp,1. (11) IH (x /x )andintegratingp(cid:48)1o,vζ1(cid:48)erth1esu1rfaceofthespherefor
p(cid:48),ζ(cid:48) 2 2
2 2
Theidentityirrepisthep=1caseandforthiscase,d =1 bothx1 andx2 implies[viatheorthogonalityoftheIp,ζ(x/x)
p=1
and Γp=1(g) = 1 for all g ∈ {1...,Ng} so that c¯p=1,ζ(x) functions(Eq.4withIpT(cid:48),l(cid:48),n(cid:48)(x/x)replacedbyIpH(cid:48),l(cid:48),n(cid:48)(x/x))]
is unconstrained by Eq. 11. Eq. 11 implies that for all g ∈ that for all g ∈{1,...,Ng},
{1,...,Ng}, Cp(cid:48),ζ(cid:48);p(cid:48),ζ(cid:48)(x1,x2)=Γp(cid:48)1(g)Cp(cid:48),ζ(cid:48);p(cid:48),ζ(cid:48)(x1,x2)[Γp(cid:48)2(g)]T.
Γp(g)c¯p,ζ(x)=c¯p,ζ(x). (12) 1 1 2 2 1 1 2 2 (19)
Consider p ∈ {2,...,Nrep} such that dp > 1. If c¯p,ζ(x) (cid:54)= Renamingp(cid:48)1,ζ1(cid:48),p(cid:48)2,andζ2(cid:48) top1,ζ1,p2,andζ2,respectively,
0 then there is at least a one-dimensional subspace which multiplying on the right by [Γp2(g)]∗, taking advantage of the
isdpin,1variant under the action of Γp(g) for all g ∈{1,...,N } unitarity of the irrep matrices Γp(g), and moving the second
g
and therefore, by definition [39, p. 67], Γp(g) is not an irrep. term to the right hand side of the equation gives
But this contradicts the assumption that Γp(g) is an irrep.
C (x ,x )[Γp2(g)]∗ =Γp1(g)C (x ,x )
Therefore,thestatementc¯ (x)(cid:54)=0 mustbefalse.Finally, p1,ζ1;p2,ζ2 1 2 p1,ζ1;p2,ζ2 1 2
p,ζ dp,1 (20)
consider p ∈ {2,...,N } such that d = 1. If c¯ (x) (cid:54)= 0
rep p p,ζ which must be true for all g ∈ {1,...,N }, p ,p ∈
thenEq.12impliesthatΓp(g)=1forallg ∈{1,...,N }but g 1 2
g {1,...,N }, ζ ,ζ in their unspecified range, and x ,x ∈
this is the identity irrep which is the p = 1 irrep. Therefore, rep 1 2 1 2
IR+ ∪ {0}. Finally, use the fact that the irreps Γp(g) (p ∈
c¯ (x)=0. In summary,
p,ζ {1,...,N }, g ∈ {1,...,N }) are real-valued orthonormal
rep g
(cid:26) rather than complex-valued unitary to get
0 , p∈{2,...,N }
c¯ (x)= dp,1 g . (13)
p,ζ arbitrary, p∈{1} C (x ,x )Γp2(g)=Γp1(g)C (x ,x )
p1,ζ1;p2,ζ2 1 2 p1,ζ1;p2,ζ2 1 2
(21)
The constraint on the covariance function of c (x) is
p,ζ which must be true for all g ∈ {1,...,N }, p ,p ∈
derived in this paragraph acting as if c (x), I (x/x), and g 1 2
p,ζ p,ζ {1,...,N }, ζ ,ζ in their unspecified range, and x ,x ∈
Γp(g) are complex valued which simplifies the presentation rep 1 2 1 2
IR+ ∪{0}. This set of matrix equations each of dimension
in Section VII. Define ρ˜(x) = ρ(x)−ρ¯(x) and c˜ (x) =
p,ζ d ×d hassubstantialstructurebecausetheΓp(g)(p ,p ∈
c (x)−c¯ .DefineC(x ,x )=E[(ρ(x )−ρ¯(x ))(ρ(x )− p1 p2 1 2
p,ζ p,ζ 1 2 1 1 2 {1,...,N }, g ∈ {1,...,N }) matrices are irreps that are
ρ¯(x ))] = E[ρ˜(x )ρ˜(x )] ∈ CC1 (as in Section II) and rep g
2 1 2 distinctforp (cid:54)=p .Schur’sLemma[39,Corollary3.2,p.70]
C (x ,x )=E[(c (x )−c¯ (x ))(c (x )− 1 2
p1,ζ1;p2,ζ2 1 2 p1,ζ1 1 p1,ζ1 1 p2,ζ2 2 implies that
c¯p2,ζ2(x2))T] = E[c˜p1,ζ1(x1)c˜Tp2,ζ2(x2)] ∈ CCdp1×dp2. Eq. 6 (cid:26)
implies that c (ζ ,x ;ζ ,x )I , p =p
Cp1,ζ1;p2,ζ2(x1,x2)= 0p1 1 ,1 2 2 dp1 ot1herwi2se
Ng dp1,dp2 (22)
(cid:88)(cid:88)
ρ˜(x)= c˜Tp,ζ(x)Ip,ζ(x/x) (14) where cp1(ζ1,x1;ζ2,x2)∈CC.
p=1 ζ The results for the mean function and covariance function
(cid:88)Ng (cid:88) of cp,ζ(x) ∈ IRdp given in Eqs. 13 and 22 are infinite-
= IpT,ζ(x/x)c˜p,ζ(x). (15) dimensional results. To reduce these to finite-dimensional
p=1 ζ resultsforcomputation,weassumethatc (x)isdescribedas
p,ζ
5
a linear combination of scalar real-valued orthonormal basis V. HETEROGENEOUSRECONSTRUCTIONWITH
functions with vector real-valued weights. In particular, SYMMETRICALSTATISTICS
Nq This section describes the maximum likelihood estimator
(cid:88)
c (x)= c ψ (x) (23) that computes the first and second order statistics of the
p,ζ p,ζ,q p,ζ,q
q=1 weightsintheorthonormalexpansionoftheelectronscattering
intensity ρ(x) from the image data under the constraint of
where c ∈ IRdp, ψ (x) ∈ IR, and
(cid:82)∞ ψ p,ζ(,xq)ψ (x)x2dx =p,ζ,q δ . The goal is symmetrical statistics for ρ(x) (Eqs. 1 and 2).
x=0 p,ζ,q1 p,ζ,q2 q1,q2 The standard imaging model is that the 2-D image is
to determine the constraints on the mean vector c¯
p,ζ,q the projection along the axis of the microscope of the 3-D
and covariance matrix V of the weights
p1,ζ1,q1;p2,ζ2,q2 electron scattering intensity of the particle linearly modified
c where c¯ = E[c ] and V =
p,ζ,q p,ζ,q p,ζ,q p1,ζ1,q1;p2,ζ2,q2 by the so-called contrast transfer function of the microscope
E[(c − c¯ )(c − c¯ )T]. In this
p1,ζ1,q1 p1,ζ1,q1 p2,ζ2,q2 p2,ζ2,q2 and observed with an unknown translation in additive noise.
notation,
Let N be the number of 2-D projection images. For the
Nq v
(cid:88) ith particle and ith image, let c be the random variable
c¯ (x)= c¯ ψ (x). (24) i
p,ζ p,ζ,q p,ζ,q
coefficients in the orthonormal expansion for ρ (x). Since
q=1 i
ρ (x)isalinearfunctionofc andthestandardimagingmodel
i i
To determine the mean vector c¯ , substitute Eq. 24 into
p,ζ,q is a linear function of ρ (x), there is a linear relationship
Eq. 13, multiply by ψ (x)x2, integrate over x ∈ [0,∞), i
p,ζ,q(cid:48) between the ith image and c , which can be written in the
and rename q(cid:48) to q to find that i
spatial frequency domain (so-called reciprocal space) by the
c¯ =0 for all p∈{2,...,N }. (25) projection slice theorem. Specifically [25],
p,ζ,q dp,1 g
TodeterminethecovariancematrixVp1,ζ1,q1;p2,ζ2,q2,firstnote yi =L(zi)ci+wi (30)
from the definition C (x ,x ) = E[(c (x ) −
p1,ζ1;p2,ζ2 1 2 p1,ζ1 1 wherey istheith2-Dprojectionimageinreciprocalspace,z
c¯ (x ))(c (x )−c¯ (x ))T]that(viaEqs.23and24) i i
p1,ζ1 1 p2,ζ2 2 p2,ζ2 2 isavectorofnuisanceparameterswhichincludetheprojection
Nq Nq orientation of the ith image and the projected location of the
(cid:88) (cid:88)
Cp1,ζ1;p2,ζ2(x1,x2)= ψp1,ζ1,q1(x1) (26) particle in the ith image, ci ∈RNc is the vector of weights of
q1=1q2=1 theithimage,L(zi)isthetransformationmatrixfromweights
×V ψ (x ), to the ith reciprocal-space image, and w is the noise which
p1,ζ1,q1;p2,ζ2,q2 p2,ζ2,q2 2 i
corrupts the ith image. We assume that the noise w and the
so that i
weights c are jointly Gaussian and independent with pdfs
(cid:90) ∞ (cid:90) ∞ i
Vp1,ζ1,q1(cid:48);p2,ζ2,q2(cid:48) = x1=0 x2=0ψp1,ζ1,q1(cid:48)(x1) (27) wiA∼mNax(i0m,Qum) alnikdelcihio∼odNe(sc¯ti,mVa)t.or is used to estimate c¯ and
×C (x ,x )ψ (x )x2x2dx dx . V from the image data, as is described in Ref. [25]. The
p1,ζ1;p2,ζ2 1 2 p2,ζ2,q2(cid:48) 2 1 2 1 2
maximumlikelihoodestimatoriscomputedbyanexpectation-
Define vp1(ζ1,q1;ζ2,q2)∈IR by maximization algorithm with nuisance parameters {zi}Ni=v1.
(cid:90) ∞ (cid:90) ∞ ThejointlyGaussianassumptiononc andw isquiteimpor-
i i
v (ζ ,q ;ζ ,q )= ψ (x ) (28)
p1 1 1 2 2 x1=0 x2=0 p1,ζ1,q1 1 wtaintthbwec,ausspeectihfiecnaltlhy,e stochastic part of ci can be combined
×c (ζ ,x ;ζ ,x )ψ (x )x2x2dx dx . i
p1 1 1 2 2 p2,ζ2,q2 2 1 2 1 2
Substitute Eq. 22 into the right hand side of Eq. 27 and then yi =L(zi)c¯i+wi(cid:48) (31)
replace q(cid:48) and q(cid:48) by q and q , respectively, to find that
1 2 1 2 where E[w(cid:48)|z ] = 0 and E[w wT|z ] = L(z )VLT(z )+Q.
i i i i i i i
Vp1,ζ1,q1;p2,ζ2,q2 =(cid:26) v0p1(ζ1,q,1;ζ2,q2)Idp1, opt1h=erwpi2se . Tpahriasmtreatnesrsfoartmtahteiocnorsetmoofvheasv{incgi}aNi=nv1oifsreomw(cid:48)thweistehtaofstnruucistaunrecde
dp1,dp2 (29) covariance matrix. Alternating updates of c¯i and V are used.
Thegoalofthispaperistoestimatec¯ andV Additional algorithmic results are described in Ref. [25] in-
p,ζ,q p1,ζ1,q1;p2,ζ2,q2
or equivalently v (ζ ,q ;ζ ,q ) from image data. cluding results for multiple discrete classes of particle where
p1 1 1 2 2
Eq. 29 implies structure for the complete covariance matrix the a priori distribution is to be estimated.
V. Suppose that the basis functions are enumerated with Our previous calculations [25]–[27], [29] were for icosahe-
indices changing in the order p (slowest), ζ, q, and j (fastest) dral symmetry and used only the p = 1 subset of the basis
where j is the index that enumerates the elements of the functions described in Section III. By direct calculation or
vector I (x/x). Then the V matrix is a block matrix with basedontheresultsinSectionIV,noconstraintsonc¯orV are
p,ζ
N ×N blocks(determinedbythevaluesofthepindices) requiredinordertoachieveEqs.1and2.Theoptimizationof
rep rep
andonlythediagonalblocksarenonzero(duetothe“p =p ” V uses Matlab’s fmincon (option “trust-region-reflective”)
1 2
condition in Eq. 29). Furthermore, each block is constructed with symbolic cost, gradient of the cost, and Hessian of the
of N ×N subblocks (where N is the total number of cost.Inordertoincorporatethecompletesetofbasisfunctions
ζ,q ζ,q ζ,q
(ζ,q) pairs) and each subblock is proportional to the d ×d (i.e., all values of p) into this software [42], we think of the
p p
identitymatrixwithproportionalityconstantv (ζ ,q ;ζ ,q ). mean vector and covariance matrix as functions of a further
p1 1 1 2 2
6
levelofparameterization.Thefurtherlevelofparameterization
isquitesimple.Specifically,themeaniseitherarbitraryorzero
(Eq.25)sothefurtherparametersarethearbitraryvalueswhile
thecovariancehascertainelementsthatmuchbeequalbutare
otherwise arbitrary (Eq. 29) so the further parameters are the
arbitraryvalues.Usingtheadditionallevelofparameterization,
the gradient and Hessian can be computed via the chain rule.
These equations compute the necessary gradient and Hessian
intermsoflargervectorsandmatriceswhicharethenreduced
in size. While this approach fits the software easily, more
efficient approaches may be possible.
VI. RECONSTRUCTIONRESULTS
Numerical results on both synthetic and experimental im-
ages of a Virus Like Particle (VLP) derived from bacterio-
phage HK97 [43] are described in this section. The VLP is
essentiallythebacteriophageminusthebacteriophage’sportal,
tail, and genome leaving only the icosahedrally symmetric Fig. 2: Four simulated real-space 2-D images of HK97. The
capsid. HK97 has a complicated lifecycle wherein it first self- color map includes only the lower 98% of the pixel values
assemblesinnearequilibriumconditionsandthenundergoesa and saturates on the upper 2%.
sequenceofessentiallyirreversiblematurationtransformations
which result in a robust particle capable of surviving outside
[47, p. 879]) between two structures (the two estimates
of the host cell. One of the first steps of maturation is the
of mean electron scattering intensity) computed from dis-
digestion by a virally-encoded protease of the so-called δ
joint sets of images. This measure is applied in all cal-
domainofthe60×7copiesofthecapsidpeptidethattogether
culations, evaluated by the formulas of Ref. [48, Eqs. 22–
make up the capsid of the bacteriophage. The experimental (cid:112)
images [44] come from the particle denoted by PhIPro+ [29] S25]: F(SkC)(=k)(cid:82) P=(k)SPP1∗,(Pk2)(dkΩ)/(cid:48) anSdP1th,Pe1(iknt)eSgPra2,lPc2a(kn)bewdhoenree
which contains a protease that is defective so that the particle P1,P2 1 2
symbolically [48, Eqs. 23].
is trapped in the Prohead I step of maturation. The average
outer radius of the capsid is 254A˚.
A. Simulated images
Both the synthetic and experimental image calculations are
performed in two steps. First a homogeneous reconstruction, Simulated2-Dprojectionimagesdirectlyinreciprocalspace
i.e., a reconstruction in which all particles are assumed to were generated from Eq. 30 which is the mathematical model
be identical and have icosahedral symmetry, is performed. In usedinthereconstructionalgorithm.Therearemultiplepossi-
terms of Section V, this is the case where p = 1, c = c¯, blesourcesforc¯,includingfromatomicresolutioncoordinates
i
and V = 0 so this algorithm is denoted by Hetero fromx-raycrystallographyexperiments,butfewsourcesforV.
(zero covar)N. cS,Neccond, a heterogeneous reconstruction, i.e., a Thevaluesofc¯andV usedinthesimulation(“groundtruth”)
reconstructioninwhichnoparticlesarethesame,isperformed comefromtheHetero(SymStat)reconstructionresults.Each
usingthehomogeneousreconstructionc¯asaninitialcondition. of the 1200 images measures 100×100 pixels with a pixel
Two types of heterogeneous reconstruction are computed and sampling interval of 5.52A˚ and a SNR of 0.5 measured on
contrasted, one uses only the p=1 basis functions (Hetero variance not standard deviation. Examples of the simulated
(SymPart)) so that each particle has the icosahedral symmetry real-space 2-D images are shown in Fig. 2.
The sampling interval is 5.52A˚ so the Nyquist spatial
while the other uses the full set of p ∈ {1,...,N = 5}
basis functions with the constraints of Section IV (rHepetero frequency for the images is 0.5/5.52 = 0.0906A˚−1. Over
(SymStat)) so that no particle has the icosahedral symme- this range the FSC curve based on reconstructions from the
try but the first and second order statistics of the electron non-overlapping subsets of images is very close to 1 so that
scattering intensity ρ(x) do have the icosahedral symmetry. the structures are essentially identical. Due to computation
The Hetero (zero covar) followed by Hetero (SymPart) limits imposed by our computer system, the number of basis
reconstructions used basis functions with p = 1, l ≤ 55, and functionsthatcanbeusedislimitedwhichrestrictsthespatial
q ≤ 20 (1060 total basis functions) while the Hetero (zero resolutionthatcanbeachieved.Onemeasureofthisrestriction
covar) followed by Hetero (SymStat) reconstructions used is that the spherically averaged energy (the denominator of
p ∈ {1,...,N }, l ≤ 10, and q ≤ 20 (2020 total basis FSC) is quite small by k = 0.0906A˚−1, approximately 10−6
rep
functions)whereonlysmallerl valuescanbeusedbecauseof times it peak value.
the substantially larger number of basis functions that occur
B. Experimental images
when p ∈ {1,...,N } versus p = 1. Only diagonal V
rep
matrices are considered. Limitations of our computer system make it impossible to
A standard measure of performance in cryo EM is the process the entire set of images in Ref. [44] so two non-
Fourier shell correlation (FSC) ( [45, Eq. 2], [46, Eq. 17], overlaping sets of 1200 images each showing one instance of
7
theparticlewererandomlyselectedforprocessing.Eachimage
measures 200×200 pixels with a pixel sampling interval of
2.76A˚.
From the mean vector c¯ and the covariance matrix V, it
is possible to compute the mean function ρ¯(x) = E[ρ(x)]
(Eqs. 7 and 24 or Ref. [25, Eq. 16]), and covariance function
C(x ,x ) = E[(ρ(x ) − ρ¯(x ))(ρ(x ) − ρ¯(x ))] (Eqs. 16
1 2 1 1 2 2
and26orRef.[25,Eq.18])oftheelectronscatteringintensity
(cid:112)
ρ(x) of the particle. We often visualize s (x) = C(x,x)
ρ
which has the same units as ρ¯(x) and is only 3-D versus 6-
D. Two pairs of mean function ρ¯(x) and standard deviation
functions (x)werecomputedbasedontheestimatesprovided
ρ
byHetero(SymPart)andHetero(SymStat).Thetwopairs
were computed twice based on the different non-overlapping
subsets of the images.
FSC was computed for the Hetero (SymPart) recon-
structions from the non-overlapping subsets of images with
the result that resolution is 0.061A˚−1 (equivalently 16.4A˚)
(a) symmetrical particles (b) symmetrical statistics
based on the frequency at which the FSC curve crosses the
(Hetero (SymPart)) (Hetero (SymStat))
0.5 level. The same calculation for the Hetero (SymStat)
reconstruction leads to an FSC curve that is almost 1 for the
Fig. 4: 3-D reconstructions of ρ¯ for PhIPro+ from Hetero
entire range of frequencies out to the Nyquist frequency for
the images which is 0.5/2.76 = 0.181A˚−1 since the pixel (SymPart)andHetero(SymStat)displayedwithdifferentcol-
ormaps.Theshapeisasurfaceofconstantintensity(5×10−4)
sampling interval is 2.76A˚. Therefore these two structures are
of ρ¯(x) colored by the standard deviation s (x), which is
essentially identical. ρ
visualized by UCSF Chimera [41]. All markings are scaled
by 10−3.
row. Results are shown for both a Hetero (SymPart) and a
Hetero (SymStat) reconstruction.
The primary conclusion from Figs. 3–5 is that the assump-
tionofsymmetricalparticles(i.e.,Hetero(SymPart))leadsto
unphysical peaks in the variance s (x) near and on symmetry
ρ
(a) 222 coordinate (b) symmetric (c) symmetric axeswhilethesepeaksarenotpresentwhenusingtheassump-
system particles statistics tion of symmetrical statistics (i.e., Hetero (SymStat)). Using
a different approach, specifically, a resampling approach, but
Fig. 3: 3-D reconstructions of the standard deviation s (x) still assuming that each particle is symmetric, this same
ρ
for PhIPro+. The shape is a surface of constant intensity type of artifact is well-known to be a common result [31,
(0.0038) of the standard deviation s (x), which is visualized p. 173] from standard software systems, e.g., EMAN2 [7].
ρ
by UCSF Chimera [41]. The panel “symmetric particles” The radially-directed artifacts in Figs. 3–5 are inconsistent
(“symmetric statistics”) is computed by Hetero (SymPart) withtheannularlystructuredPhIPro+particle,whichhas,from
(Hetero (SymStat)). outsidetoinside,annuliofwellorderedcapsidprotein,poorly
ordercapsidδdomainandproteaseprotein,andpoorlyordered
Figs. 3(b) and (c) show the standard deviation s (x) com- DNAgenome.Theannularcharacteristicofthestandarddevia-
ρ
puted by Hetero (SymPart) and Hetero (SymStat), respec- tionfunctions (x)fromtheHetero(SymStat)reconstruction
ρ
tively. In Fig. 4, ρ¯(x) computed by Hetero (SymPart) and isparticularlyclearinFig.3(b)andevenmoresoinFig.4(a).
the two heterogeneous solutions s (x) computed by Hetero Theannularly-directedartifactsaretheoriginofthelightblue
ρ
(SymPart) and Hetero (SymStat) are jointly visualized by radial lines in Fig. 5(a) while the actual physical structure of
computing a 3-D surface of constant value of ρ¯(x) and color- the particle is the origin of the annular character of Fig. 5(b).
ing the surface by s (x). The ρ¯and s functions in Figs. 3–4 Insummary,theresultsbasedonassumingsymmetricstatistics
ρ ρ
have the same orientation as is shown in Fig. 3(a), wherein (Hetero(SymStat))ratherthansymmetricparticles(Hetero
the x-, y-, and z-axes are all 2-fold rotational symmetry axes (SymPart)) make more biological sense. Previous analyses
of the icosahedron which is one of two possible so-called 222 basedonHetero(SymPart)primarilyusedthesquarerootof
coordinatesystems.Fig.5showscrosssectionsofthestandard thesphericalaverageofthevarianceC(x,x)andthereforethe
deviation function s (x) with contours of the mean function peaks along radially-directed lines were averaged out leading
ρ
ρ¯(x). Each cross section is either perpendicular to one of the tointerpretableinformation.However,aswetrytoincreasethe
three types of symmetry axis (5-, 3-, and 2-fold) or includes spatialdetailofourinterpretation,suchbehaviorisdifficultto
thesymmetryaxisasisdescribedinthediagramsintheupper understand.
8
morelimited[39,p.70],e.g.,amatrixmaynothaveanyeigen-
values since the roots of the characteristic polynomial could
potentially all be complex. In this section, the constraints on
themomentsoftheweightscausedbysymmetricstatisticsare
investigated when the irreps and therefore the basis functions
and the weights are complex.
The results for the mean function of c (x) are unaltered
p,ζ
from Eqs. 7–13.
The results for the covariance function of c (x) are
p,ζ
changed from Eqs. 15–20. Because the variables are now
complex, it is necessary to compute and to apply symmetry
constraints to both Cρ,ρ(x ,x ) and Cρ,ρ∗(x ,x ) which are
1 2 1 2
definedbyCρ,ρ(x ,x )=E[(ρ(x )−ρ¯(x ))(ρ(x )−ρ¯(x ))]
1 2 1 1 2 2
(the same as C(x ,x ) in Section IV) and Cρ,ρ∗(x ,x ) =
1 2 1 2
E[(ρ(x ) − ρ¯(x ))(ρ(x ) − ρ¯(x ))∗], respectively. The cal-
1 1 2 2
culations for Cρ,ρ are the same as Eqs. 15–20. In a more
detailednotationwhereCc,c (x ,x )=E[(c (x )−
p1,ζ1;p2,ζ2 1 2 p1,ζ1 1
c¯ (x ))(c (x )−c¯ (x ))T], the result is
p1,ζ1 1 p2,ζ2 2 p2,ζ2 2
Cc,c (x ,x )[Γp2(g)]∗ =Γp1(g)Cc,c (x ,x ).
p1,ζ1;p2,ζ2 1 2 p1,ζ1;p2,ζ2 1 2
(32)
(a) symmetrical individuals (Hetero (SymPart))
ItisnotpossibletocontinueontoEq.21becauseΓp(g)isno
longer real and a new solution is described in the following
paragraph. The calculations for Cρ,ρ∗ follow the same plan
with the result that [Cc,c∗ (x ,x ) = E[(c (x ) −
p1,ζ1;p2,ζ2 1 2 p1,ζ1 1
c¯ (x ))(c∗ (x )−c¯∗ (x ))T]]
p1,ζ1 1 p2,ζ2 2 p2,ζ2 2
Cc,c∗ (x ,x )Γp2(g)=Γp1(g)Cc,c∗ (x ,x ).
p1,ζ1;p2,ζ2 1 2 p1,ζ1;p2,ζ2 1 2
(33)
Eq. 33 is of the same form as Eq. 21 and so has the same
form of solution which is
Cc,c∗ (x ,x )
p1,ζ1;p2,ζ2 1 2
(cid:26) cc,c∗(ζ ,x ;ζ ,x )I , p =p
= p1 1 1 2 2 dp1 1 2 . (34)
0 , otherwise
dp1,dp2
(b) symmetrical statistics (Hetero (SymStat))
The solution to Eq. 32 is described in this paragraph. Two
irreps are the same when there exists an invertible matrix
Fig.5: Crosssectionsofthestandarddeviationfunctions (x)
for PhIPro+ using both Hetero (SymPart) and Hetρero S such that Γ1(g) = S−1Γ2(g)S for all g. Since Γp2(g)
is an irrep, it follows that [Γp2(g)]∗ is a rep. By Ref. [39,
(SymStat)anddisplayedwithacommoncolormap.Thewhite
Theorem 3.5 p. 70], [Γp2(g)]∗ is an irrep. There is a fixed set
lines depict contour levels of the mean ρ¯(x). The geometrical
of irreps, so there must be a function τ from {1,...,N }
orientation of the s (x) function is shown in the top row rep
and for each orientaρtion, cross sections in the x-y and x- to {1,...,Nrep} such that [Γp2(g)]∗ = S−1Γτ(p2)(g)S for
some invertible matrix S. The function τ must be a permu-
z planes are displayed so that, for instance, a cross section
perpendiculartoa5-foldaxisandacrosssectionthatcontains tΓapti∗o(ng)sianncde iSf(cid:48)[τΓ(pp(cid:48)22()g)=]∗[τS((cid:48)p]−(cid:48)21) ==pΓ∗p∗th(gen).SE[qΓupa2t(ign)g]∗tSh−e1le=ft
a 5-fold axis are both shown. Panel (a) shows the recon-
hand sides of these two equations gives S[Γp2(g)]∗S−1 =
struction results based on symmetrical particles, i.e., Hetero S(cid:48)[Γp(cid:48)2(g)]∗[S(cid:48)]−1 which implies T−1Γp2(g)T = Γp(cid:48)2(g)
(SymPart), whereas Panel (b) shows the reconstruction results
(where T = [S−1S(cid:48)]∗) which implies that p and p(cid:48) are the
based on symmetrical statistics, i.e., Hetero (SymStat). 2 2
same irrep. Therefore, Eq. 32 is equivalent to
VII. CASESWHEREREAL-VALUEDIRREPSARENOT Cc,c (x ,x )Γτ(p2)(g)=Γp1(g)Cc,c (x ,x ).
AVAILABLEANDNON-SPHERICALCOORDINATESYSTEMS p1,ζ1;p2,ζ2 1 2 p1,ζ1;p2,ζ2 1 2(35)
Eq. 35 is of the same form as Eq. 21 and so has the same
Standard references for rotational point group irreps give
form of solution which is
the irrep matrices in complex-valued unitary form. In some
instances, there is no similarity transformation to a real- Cc,c (x ,x )
valuedorthonormalform[38].Whileitispossibletoconsider p1,ζ1;p2,ζ2 1 2
(cid:26) cc,c(ζ ,x ;ζ ,x )I , p =τ(p )
representation theory over a real-valued vector space rather = p1 1 1 2 2 dp1 1 2 . (36)
0 , otherwise
than a complex-valued vector space, the real-valued theory is dp1,dτ(p2)
9
Rather than Cc,c (x ,x ) and Cc,c∗ (x ,x ) with random-variable weights c and there is no mean or
p1,ζ1;p2,ζ2 1 2 p1,ζ1;p2,ζ2 1 2 p,ζ
it is probably more natural to work with the real- covariance function intermediate on the path to mean vector
valued functions C(cid:60)c,(cid:60)c (x ,x ) = E[(cid:60)[c (x ) − andcovariancematrixconstraintsanalogoustoEqs.25and29.
p1,ζ1;p2,ζ2 1 2 p1,ζ1 1
c¯ (x )](cid:60)[c (x )−c¯ (x )]T], C(cid:61)c,(cid:61)c (x ,x )=
p1,ζ1 1 p2,ζ2 2 p2,ζ2 2 p1,ζ1;p2,ζ2 1 2
E[(cid:61)[c (x ) − c¯ (x )](cid:61)[c (x ) − c¯ (x )]T], VIII. SYMMETRIZINGNON-SYMMETRICFUNCTIONS
p1,ζ1 1 p1,ζ1 1 p2,ζ2 2 p2,ζ2 2
and Cp(cid:60)1c,,ζ(cid:61)1;cp2,ζ2(x1,x2) = E[(cid:60)[cp1,ζ1(x1) − Suppose s is a function from IR3 to IR. Suppose Rβ
c¯p1,ζ1(x1)](cid:61)[cp2,ζ2(x2) − c¯p2,ζ2(x2)]T] which can easily (β ∈{1,...,Ng}) are the rotation matrices for a point-group
be computed via expectations of products of the terms symmetry. An important result in theory and software is that
(cid:60)[c (x ) − c¯ (x )] = [(c (x ) − c¯ (x )) +
(c p1,ζ(1x )1−c¯ p1(,ζx1 ))1∗]/2and(cid:61)p[c1,ζ1 (1x )−c¯p1,ζ1 (x1)]= . (cid:88)Ng
p1,ζ1 1 p1,ζ1 1 p1,ζ1 1 p1,ζ1 1 ρ(x)= s(R−1x) (37)
[(c (x )−c¯ (x ))−(c (x )−c¯ (x ))∗]/(2i). β
Eitph1e,rζ1 1Cc,c p1,ζ1 (x1 ,x ) p1,ζa1nd1 Cpc1,c,∗ζ1 1 (x ,x ) β=1
or C(cid:60)c,p(cid:60)1c,ζ1;p2(,xζ2,x1),2 C(cid:61)c,(cid:61)c (px1,ζ,1x;p2),ζ,2 1 an2d has the symmetry described by Eq. 1. The proof is via the
p1,ζ1;p2,ζ2 1 2 p1,ζ1;p2,ζ2 1 2 Rearrangement Theorem [40, Theorem II, p. 24] which states
Cthpe(cid:60)1c,s,ζa(cid:61)1m;cp2e,ζm2(exth1o,xds2)usecdanin SbeectimonadIeV. fiFninitael-lyd,imafetnesriochnoalosinbgy thatasumoveralltheelementsofagroup,i.e.,(cid:80)Nβ=g1f(Rβ),
isunaltered(althoughthetermsinthesumoccurinapermuted
anorderingoftheindices(p,ζ)andtherealversusimaginary
order)ifR isreplacedbyR R whereβ isafixedelement
parts, a real-valued positive semi-definite matrix V can be β β0 β 0
of the group. The purpose of this section is to describe a
assembled which has block structure similar to the situation
similarsymmetrizationprocedureforstochasticprocessessuch
in Section IV.
that the mean (covariance) of the result has the symmetry
The calculations in this section allow ρ(x) to be complex.
described by Eq. 1 (Eq. 2).
However, ρ(x) is real in the cryo EM application. A “penalty
One construction that leads to a stochastic process with the
function” approach would be to create synthetic measure- desired statistics starts with a stochastic process s(x)∈IRNg
ments of (cid:61)ρ(x) with value 0 and variance σ2 and process
(cid:61) with components si(x) that is exchangeable relative to first
the resulting complex-valued ρ(x), possibly allowing σ(cid:61)2 to and second order statistics, i.e., if s(cid:48)(x) ∈ IRNg is a permu-
decrease as the iterative algorithm converges. A constraint tation of the elements of s(x), then s(cid:48)(·) and s(·) have the
approach would be to determine additional constraints on
same first and second order statistics. This implies that there
cp,ζ(x) such that ρ(x) is guaranteed to be real. We believe exist functions s¯(x), C=(x ,x ), and C(cid:54)=(x ,x ) such that
that such constraints involve the behavior of basis functions . 1 2 1 2
s¯ (x) = E[s (x)] = s¯(x) (for all β ∈ {1,...,N }) and
underconjugation(i.e.,forFourierserieswithperiodT in1-D, β β . g
C (x ,x ) = E[(s (x )−s¯(x ))(s (x )−s¯(x ))] =
ψarn∗e(nt)ot=cl[oesxeply(ir2eπlantetd/Tto)]t∗he=roetxapti(o−nai2lπponitn/tT-g)ro=upψs−ynm(mt)e)traineds δββ11,β,β22C=1(x12,x2)+(1β−1 δβ11,β2)C(cid:54)=(1x1,xβ22) (f2or all β12,β2 ∈
{1,...,N }) .
g
discussed in this paper.
Define ρ(·) by
Except in this paragraph, results in this paper are stated
in the spherical coordinate system (e.g., Eq. 6). Spheres, (cid:88)Ng
ρ(x)= s (R−1x). (38)
which are contours of constant radius in spherical coordi- β β
nates, enclose typical so-called spherical viruses reasonably β=1
closely and therefore lead to parsimonious parameterizations. We will directly compute the first and second order statistics
This is typical for particles obeying the symmetries of I, ofρ(·)andfindthattheysatisfyEq.1andEq.2,respectively,
.
O, and T because the symmetry axes have widely varying as is desired. For the first order statistics, ρ¯(x) = E[ρ(x)] =
directions. However, for particles obeying the symmetries (cid:80)Ng s¯(R−1x) which is in the form of Eq. 37 so that ρ¯(·)
β=1 β
of C (and possibly D ), e.g., the tail of bacteriophage has the symmetry described by Eq. 1. For the second order
n n
P22 [49], cylindrical coordinates (r,φ,z) can be more natural statistics,
with the single symmetry axis of Cn on the z axis of the .
C(x ,x )=E[(ρ(x )−ρ¯(x ))(ρ(x )−ρ¯(x ))] (39)
coordinate system. For the case of C , Eq. 6 is replaced 1 2 1 1 2 2
n
by ρ(x) = (cid:80)Ng (cid:80) cT (r,z)I (φ) and the constraint Ng Ng
p=1 ζ p,ζ p,ζ = (cid:88) (cid:88) C (R−1x ,R−1x ) (40)
of symmetrical statistics specifies the mean and covariance β1,β2 β1 1 β2 2
functions of cp,ζ(r,z) rather than the mean (Eq. 13) and β1=1β2=1
covariance (Eqs. 22, 34, and 36) functions of cp,ζ(x). The (cid:88)Ng (cid:88)Ng (cid:104)
case of Dn is somewhat more complicated because there are = δβ1,β2C=(Rβ−11x1,Rβ−21x2) (41)
additional rotational symmetry axes beyond z axis. In either β1=1β2=1
(cid:105)
case,theeffectofthesymmetricalstatisticsconstraintisacon- + (1−δ )C(cid:54)=(R−1x ,R−1x ) .
straint on the covariance function of the coordinate-dependent β1,β2 β1 1 β2 2
weights where the weights depend on the coordinates that Independent of the properties of C= and C(cid:54)=, application of
are not involved in the symmetry operation. Finally, if basis the Rearrangement Theorem [40, Theorem II, p. 24] demon-
functions F (x) are used which directly have symmetry strates the C(x ,x ) satisfies the symmetry of Eq. 2. For
p,ζ 1 2
properties that are the obvious generalization of Eq. 5, then Gaussianstochasticprocesses,meanandcovariancesatisfying
the orthonormal expansion is ρ(x) = (cid:80)Ng (cid:80) cT F (x) the constraints of Section IV or VII provides a simulation
p=1 ζ p,ζ p,ζ
10
method based on the basis functions with special symmetry [7] G.Tang,L.Peng,P.R.Baldwin,M.D.S.,W.Jiang,I.Rees,andS.J.
properties but for non-Gaussian processes, problems with Ludtke, “EMAN2: An extensible image processing suite for electron
microscopy,”J.Struct.Biol.,vol.157,no.1,pp.38–46,2007.
manybasisfunctionshavingsymmetryproperties,orproblems
[8] S.H.W.Scheres,“RELION:ImplementationofaBayesianapproachto
with basis functions that lack symmetry properties, Eq. 38 cryo-EMstructuredetermination,”J.Struct.Biol.,vol.180,pp.519–530,
provides an alternative method. 2012.
[9] T.R.Shaikh,H.Gao,W.T.Baxter,F.J.Asturias,N.Boisset,A.Leith,
and J. Frank, “Spider image processing for single-particle reconstruc-
IX. CONCLUSION tion of biological macromolecules from electron micrographs,” Nature
Protocols,vol.3,no.12,pp.1941–1974,2008.
In this paper we formulate and solve 3-D modeling and [10] H.D.Tagare,A.Kucukelbir,F.J.Sigworth,W.Hongwei,andM.Rao,
reconstructionproblemsforstochasticsignalswherethestatis- “Directlyreconstructingprincipalcomponentsofheterogeneousparticles
fromcryo-EMimages,”J.Struct.Biol.,vol.191,pp.245–262,2015.
tics of the stochastic signal are required to have a geometric
[11] J. Frank and A. Ourmazd, “Continuous changes in structure
symmetry.Theworkismotivatedbycryoelectronmicroscopy, mapped by manifold embedding of single-particle data in cryo-
which is of great importance in structural biology and of EM,” Methods, vol. 100, pp. 61–67, 2016, single Particle Cryo-
EM, from sample to reconstruction. [Online]. Available: http:
growing importance in the material science of nanoscale
//www.sciencedirect.com/science/article/pii/S1046202316300251
particles. [12] A.Dashti,P.Schwander,R.Langlois,R.Fung,W.Li,A.Hosseinizadeh,
Previous algorithms assumed that the realizations of the H. Y. Liao, J. Pallesen, G. Sharma, V. A. Stupina, A. E. Simon, J. D.
Dinman, J. Frank, and A. Ourmazd, “Trajectories of the ribosome
stochastic signal have the geometric symmetry. While this
as a brownian nanomachine,” Proceedings of the National Academy
implies that the statistics have the symmetry, it is a sufficient of Sciences, vol. 111, no. 49, pp. 17492–17497, 2014. [Online].
condition rather than a necessary condition and it has lead to Available:http://www.pnas.org/content/111/49/17492.abstract
[13] Q. Jin, C. O. S. Sorzano, J. M. de la Rosa-Trev´ın, J. R. Bilbao-
inaccurateresultsforthesecondorderstatisticsoftheelectron
Castro, R. Nu´n˜ez-Ram´ırez, O. Llorca, F. Tama, and S. Jonic´,
scattering intensity on and near symmetry axes [31, p. 173]. “Iterative elastic 3D-to-2D alignment method using normal modes
Sincetheseareoftenlocationsofimportantbiology,inaccurate for studying structural dynamics of large macromolecular complexes,”
Structure, vol. 22, no. 3, pp. 496–506, 2014. [Online]. Available:
results greatly reduce the value of the entire experiment.
http://www.sciencedirect.com/science/article/pii/S0969212614000136
Beyond correcting the problem in existing computations, [14] S. J. Ludtke, P. R. Baldwin, and W. Chiu, “EMAN: semiautomated
this is the first publication of which the authors are aware in software for high-resolution single-particle reconstructions,” J. Struct.
Biol.,vol.128,no.1,pp.82–97,Dec.11999.
whichafull6-Dcovariancefunctioniscomputed.Interpreting
[15] S. H. W. Scheres, “A Bayesian view on cryo-EM structure determina-
the 6-D information is challenging and we are working with tion,”J.Mol.Biol.,vol.415,no.2,pp.406–418,13January2012.
our biological collaborators on this problem. [16] J. Frank, M. Radermacher, P. Penczek, J. Zhu, Y. Li, M. Ladjadj, and
A. Leith, “SPIDER and WEB: Processing and visualization of images
in3Delectronmicroscopyandrelatedfields,”J.Struct.Biol.,vol.116,
ACKNOWLEDGMENTS pp.190–199,1996.
[17] P.A.Penczek,M.Kimmel,andC.M.T.Spahn,“Identifyingconforma-
WearegratefultoProfessorsJohnE.Johnson(TheScripps tionalstatesofmacromoleculesbyeigen-analysisofresampledcryo-EM
ResearchInstitute)andDavidVeesler(UniversityofWashing- images,”Structure,vol.19,pp.1582–1590,9November2011.
[18] C.M.T.SpahnandP.A.Penczek,“Exploringconformationalmodesof
ton) for the HK97 data and helpful discussions; to Drs. Yili
macromolecularassembliesbymultiparticlecryo-EM,”CurrentOpinion
Zheng and Qiu Wang and to Ms. Yunye Gong for their effort inStructuralBiology,vol.19,pp.623–631,2009.
on the Hetero software system; and for support from NSF [19] W. Zhang, M. Kimmel, C. M. T. Spahn, and P. A. Penczek, “Hetero-
1217867. genneity of large macromolecular complexes revealed by 3D cryo-EM
varianceanalysis,”Structure,vol.16,pp.1770–1776,2008.
[20] P.A.Penczek,C.Yang,J.Frank,andC.M.T.Spahn,“Estimationof
REFERENCES varianceinsingle-particlereconstructionusingthebootstraptechnique,”
J.Struct.Biol.,vol.154,no.2,pp.168–183,2006.
[1] S.E.Kim,L.Zhang,K.Ma,M.Riegman,F.Chen,I.Ingold,M.Conrad, [21] A. Simonetti, S. Marzi, A. G. Myasnikov, A. Fabbretti, M. Yusupov,
M. Z. Turker, M. Gao, X. Jiang, S. Monette, M. Pauliah, M. Gonen, C. O. Gualerzi, and B. P. Klaholz, “Structure of the 30S translation
P.Zanzonico,T.Quinn,U.Wiesner,M.S.Bradbury,andM.Overholtzer, initiation complex,” Nature, vol. 455, no. 7211, pp. 416–420, 18 Sep.
“Ultrasmallnanoparticlesinduceferroptosisinnutrient-deprivedcancer 2008.[Online].Available:http://dx.doi.org/10.1038/nature07192
cells and suppress tumour growth,” Nature Nanotechnology, vol. 11, [22] H. Y. Liao, Y. Hashem, and J. Frank, “Efficient estimation of three-
no.11,pp.977–985,Nov.2016. dimensionalcovarianceanditsapplicationintheanalysisofheteroge-
[2] K. Ma, H. Sai, and U. Wiesner, “Ultrasmall sub-10nm near-infrared neous samples in cryo-electron microscopy,” Structure, vol. 23, no. 6,
fluorescent mesoporous silica nanoparticles,” J. American Chemical pp.1129–1137,2June2015.
Society,vol.134,no.32,pp.13180–13183,15Aug2012. [23] H. Y. Liao and J. Frank, “Classification by bootstrapping in single
[3] Y. He, T. Ye, M. Su, C. Zhang, A. E. Ribbe, W. Jiang, and C. Mao, particle methods,” in Proceedings of the 2010 IEEE International
“Hierarchical self-assembly of DNA into symmetric supramolecular SymposiumonBiomedicalImaging,14Apr.2010,pp.169–172.
polyhedra,” Nature, vol. 452, pp. 198–201, 13 March 2008. [Online]. [24] E.Katsevich,A.Katsevich,andA.Singer,“Covariancematrixestima-
Available:doi:10.1038/nature06597 tion for the cryo-EM heterogeneity problem,” SIAM Imaging Science,
[4] S. Subramaniam, W. Ku¨hlbrandt, and R. Henderson, “CryoEM at vol.8,no.1,pp.126–185,2015.
IUCrJ: a new era,” IUCrJ, vol. 3, no. 1, pp. 3–7, Jan 2016. [Online]. [25] Y. Zheng, Q. Wang, and P. C. Doerschuk, “3-D reconstruction of the
Available:https://doi.org/10.1107/S2052252515023738 statistics of heterogeneous objects from a collection of one projection
[5] G.McMullan,A.Faruqi,D.Clare,andR.Henderson,“Comparisonof image of each object,” Journal of the Optical Society of America A,
optimal performance at 300 keV of three direct electron detectors for vol.29,no.6,pp.959–970,Jun.2012.
use in low dose electron microscopy,” Ultramicroscopy, vol. 147, pp. [26] Q. Wang, T. Matsui, T. Domitrovic, Y. Zheng, P. C. Doerschuk, and
156 – 163, 2014. [Online]. Available: http://www.sciencedirect.com/ J. E. Johnson, “Dynamics in cryo EM reconstructions visualized with
science/article/pii/S030439911400151X maximum-likelihood derived variance maps,” Journal of Structural
[6] M.Khoshouei,M.Radjainia,A.J.Phillips,J.A.Gerrard,A.K.Mitra, Biology,vol.181,no.3,pp.195–206,Mar.2013.
J.M.Plitzko,W.Baumeister,andR.Danev,“Voltaphaseplatecryo-EM [27] J. Tang, B. M. Kearney, Q. Wang, P. C. Doerschuk, T. S. Baker, and
of the small protein complex Prx3,” Nature Communications, vol. 7, J.E.Johnson,“DynamicandgeometricanalysesofNudaureliacapensis
2016.[Online].Available:http://dx.doi.org/10.1038/ncomms10534 ω virusmaturationrevealtheenergylandscapeofparticletransitions,”