Comparing cryo-EM methods and molecular dynamics simulation to investigate heterogeneity in ligand-bound TRPV1

Domain: physics
Published: 2025-12-29
Editorial outcome: major_revision

Abstract
Cryogenicelectronmicroscopy(cryo-EM)hasemergedasapowerfulmethod
forresolvingthestructureofbiologicalmacromolecules.Recently,several
computationalmethodshavebeendevelopedtostudytheheterogeneityof
moleculesinsingle-particlecryo-EM.Inthisstudy,weanalyzeapubliclyavailable
datasetofTRPV1usingfivesuchmethods:3DFlex,3DVA,cryoDRGN,Mani-
foldEM,andBayesianensemblereweighting.Wefindsignificantheterogeneity,
buteachmethodproducesdifferentresults,withsomedetectingonlycompo-
sitionalorconformationalheterogeneity.Tocomparethesediverseresults,we
developAnaVoxtoquantitativelydetermineagreementbetweenheterogeneity
methods.Furthermore,applyingBayesianensemblereweightingcombinedwith
moleculardynamicssimulationssupportsthepresenceoftheserarerstates
withinthesample.Thisstudyshowsthatalthoughcurrentmethodsrevealthe
presenceofheterogeneity,theirstochasticityandpotentialbiaspresentchal-
lengesfortheirroutineuse.However,withfuturedevelopment,thesetoolswill
enabletheuseofcryo-EMdataforquantitativebiophysicalinvestigations.

Introduction
The field of structural biology has been dominated by the drive to moreclassicalBayesianinferenceapproaches12–14,wepointthereader
achievehigh-resolution,atomicstructuresofbiomolecules.Cryogenic elsewhereforamorecomprehensivereview5.Theincreasingnumber
electronmicroscopy(cryo-EM)hasemergedastheleadingtechnique ofheterogeneitytoolsintheliteraturemotivatedustoperformthe
to determine these high-resolution structures for larger mostrigorousstudyofheterogeneityinabiologicallyrelevantsample
biomolecules1,2.Recently,theabilityofcryo-EMtocapturenotjusta to date and determine if these tools provide complementary or
single structure but a continuous distribution of states has been orthogonalinformation.Whilerecentworkmotivatedbysimilargoals
attracting growing interest. Access to this information hidden in comparedresultsfromheterogeneitymethodswhenappliedtosyn-
existing cryo-EM data has implications for understanding biomole- theticdata15,thestudypresentedhereprovidescomparisonmetrics
cularmechanismsandtreatingdisease.Heterogeneityhasbeenknown and validation tools when studying real data. Data collected in
toexistinvitrifiedsamplessincetheinceptionofcryo-EM3–5.However, experimentstypicallyexhibitssignificantlygreatercomplicationsthan
itisanopenquestionwhethercontinuouslandscapesofconforma- idealizedsyntheticdatacreatedforbenchmarks,suchasamixtureof
tionalheterogeneityandthepopulationsoftheirrepresentativestates bothcompositionalandconformationalheterogeneity,anon-uniform
canbereliablyrecoveredfromthesesamples.Recentanalysistechni- distributionofposes,andunknowngroundtruthstructures.Forthis
questakeadvantageofthelatestadvancesinmachinelearningand more practical context, we have developed a software framework,
statisticalinferenceandrangefromgeometricmachinelearning6to AnaVox, to quantitatively determine whether two heterogeneity
variational autoencoders7–9 to flow field generators10,11 to linear or methodsagreewitheachother.

Results
3DFlexrevealsTMDtwistingandcytosolicdomainlengthening
Forallanalysesinthisstudy,weusedthe218,717particleimagesof
TRPV1boundtoDkTxandRTXavailableinEMPIAR-1005921.Following
signal subtraction of the nanodisc (see Methods for details), we
reconstructed this stack to a resolution of 3.0 Å using cryoSPARC.
Recently, 3D Flexible Refinement (3DFlex) was added to the cryoS-
PARC software suite. It was designed to increase the resolution of
flexibleportionsofabiomoleculebyaccountingfortheirconforma-
tional heterogeneity10. By embedding the molecule’s consensus
reconstructioninaflexiblemeshgrid,3DFlexintroducescontinuous
deformationstoavolume.Thenatureofthesedeformationsisdic-
tated by a generative neural network that uses a flow generator to
convertalatentcoordinatevectorintoadeformationfieldthatcon-
vectsthecanonicalmapoftheconsensusreconstruction.Itistrained
toincreasetheoverlapbetweenadeformedvolumeandtheelectron
densityinparticleimages,therebyincreasingthesignal-to-noiseratio
inflexibledomainsofthemolecule.3DFlexmodelsareverysensitive
toThechoicesofhyperparameters.Forexample,morerigidmodels
exhibitednoticeableasymmetriesintheirreconstructionofdifferent
subunitsofthecytosolicdomains,despitethesedomainsbeingsym-
metric(SupplementaryFig.S1A).Theseasymmetrieswouldpersiston
symmetry-expandedparticlestacks,andwithmodelstrainedwitha
lessrigidmodel(Rigidity=1),buttheywerelesspronounced.Further,
wefoundthatbyindependentlytrainingtwomodelswiththesame
hyperparameters, we could obtain noticeably different reconstruc-
tions and learned deformations (Supplementary Fig. S1B, Supple-
mentaryMovie2).Despitethesenoticeabledifferences,thesemodels
yieldthesametrainingscoresandglobalresolutions,raisingquestions
surroundingwhichsetof3DFlex-derivedmotionstobelieve.
Theoriginal3DFlexpaperanalyzedtheTRPV1datasetwefocuson
here,andshowedincreasedlocalresolutionofregionsofthecytosolic
domainsthathadpreviouslybeenpoorlyresolved.Inourhands,we
achievedsimilarimprovementinthelocalresolutionofthesecytosolic
domains,tobelow3Åinsomehelices(Fig.2A,B).Bytrainingatwo
dimensional latent space, we find two dominant conformational
changesinthesample(Fig.2C).ThefirstisadeformationoftheTMDs
andcytosolicdomainstwistinginoppositedirections,whilethelatter
domainssimultaneouslylengthen,similartowhatonewouldexpect
foratransitionbetweentheligand-boundopenandtheapoclosed
states (Fig. 2D, SupplementaryMovie 3).During thesemotions, the
height of the density increases by approximately 5 Å. The second
motioncharacterizedbyour3DFlexmodelcorrespondstoanasym-
metricsqueezing motionof the cytosolicdomains,with the helices
movingapproximately5Å(Fig.2E).Overall,themotionswehaveseen
aresimilartothoseseenintheinitial3DFlexanalysisofthisdataset,
though we found the optimization of hyperparameters nontrivial.
Sincetheprimaryobjectiveofthismethodistoincreasetheresolution
of the consensus refinement, and because the extra density of the
toxinmaycauseartifactsinthevolume-preservingflowfield,wewere
motivated to useother methods toinvestigate the existenceofthe
motions seen in this ligand-bound TRPV1 dataset by the 3DFlex
analysis.

3DVAdetectsTMDtwisting,cytosolicdomainlengthening,and
variabledensityinDkTx
3DVariabilityAnalysis(3DVA)isasimpleranalysistoolthan3DFlex,
andisalsoavailableincryoSPARC.Itperformsprincipalcomponent
analysis (PCA) on particle images in Fourier space14, introducing
deviations from the average volume along linear eigenvectors to
minimize the difference between volumes deformed along these
eigenvectorsandtheensembleofimages.WhenappliedtothisTRPV1
dataset,3DVAdetectsconformationalchangescomparabletothose
seenin3DFlex:thefirstmotiondetectedby3DVAistheasymmetric
squeezingofthecytosolicdomains(Fig.3A),whilethesecondmotion
isthetwistingofthetransmembraneandcytosolicdomains(Fig.3B,
SupplementaryMovie4).Notethedifferentordercomparedto3DFlex
(Fig.3C).Additionally,theresultsfrom3DVAappeartoberobustto
changes in hyperparameters. Changing the filter resolution of the
3DVAorchoosingtoexpandthesymmetryoftheparticlestackdidnot
noticeablychangeourresults.However,inspectionsofhigherorder
principalcomponentsrevealedcomparativelysmallerconformational
changes,withhigheramountsofasymmetrycomparedtothefirsttwo
principalcomponentspresentedinFig.3.Thesemotionsappearasa
‘wobbling’ofthecytosolicdomainsrelativetotheTMDs,similarto
motionsfoundinthehigherorderprincipalcomponentsoftheMD
simulations, which we will analyze later (Supplementary Movie 5,
SupplementaryFig.S7).
In addition to the conformational changes detected by 3DVA,
explorations of the latent space also produce volumes that exhibit
compositional heterogeneity (Fig. 3A, B). Both the first and second
componentsofthelatentspacerevealfadingofthedensityofDkTx,
indicatingthattheremaybecompositionalheterogeneityinthispar-
ticlestack—insomeparticles,thetoxinispartiallyorfullyunbound
fromtheTRPV1protein.However,wecannotbesurethatthisistrue
compositional heterogeneity, as experiments with synthetic data
demonstrate that under some conditions, a non-uniform pose dis-
tributioninadatasetmaycause3DVAtodetectcompositionalhet-
erogeneitywherenoneexists(SupplementaryFig.S16).Givenitsfast
runtime(SupplementaryTableS1)anditsabilitytodetectbothcom-
positional and conformational heterogeneity, we found that 3DVA
couldbeusedtoquicklyidentifythemainsourcesofheterogeneityin
adataset.However,thelinearbasisofprincipalcomponentscalculated
inaFourierbasismayinfluencewhichkindsofmolecularchangesare
detected24.

TheCryoDRGNlatentspaceencodeschangesinDkTxdensity
While3DFlexand3DVAhavetheconvenienceofbeingimplementedin
cryoSPARCanotherpopulartoolforanalyzingheterogeneityinsingle
particle cryo-EM is cryoDRGN (Deep Reconstructing Generative
Networks)7.CryoDRGNisavariationalautoencoder(VAE)thatconsists
oftwoneuralnetworksbuiltwithanimage-encoder-volume-decoder
architecture.Thefirstnetworktakesrawparticleimagesasinput,and
outputsalow-dimensionalvector,embeddingtheimagesinalatent
space.Thesecondneuralnetworktakesapointinthislatentspaceas
input and outputs a 3D volume. These neural networks are trained
suchthatreconstructionsfromimagesinaspecificregionofthislatent
spacewillproducedensitiesthataresimilartothevolumespredicted
bythedecodernetwork.
Applying cryoDRGN to this TRPV1 dataset resulted in a latent
spacethatcapturedsignificantfadingofthedensityofDkTx,similarto
whatwasseenin3DVA(Fig.4A,SupplementaryMovie6).However,
there are comparatively few conformational changes to the TRPV1
proteinitself,unlike3DVAand3DFlex.Byidentifyingregionsofthe
latentspacethatproducevolumeswithweakDkTxdensity,wecan
identifyclustersofparticlesthatmaynotexhibitgoodviewsofDkTx
(Fig.4B).Directreconstructionsofimagesintheseclustersrevealed
volumeswithtoxinoccupancyconsistentwiththepredictionsofthe
cryoDRGN VAE (Fig. 4C, Supplementary Movie 7).Theseresults are
consistentwithwhatwouldbeexpectedifsomeparticlesinthesample
didnothaveDkTxboundtoTRPV1,orifsomeparticleshaveDkTx
bound to only certain subunits, as has been observed in other
samples22.Givenourotherresults,thelackofconformationalhetero-
geneity in the cryoDRGN result was unexpected, and suggests that
compositional heterogeneity can overwhelm small conformational
changes with thisVAE architecture. However, whether these results
representtruecompositionalheterogeneityorahallucinationdueto
thedataset’snon-uniformposedistributionisunclear(Supplementary
Figs.S9,S16).

ManifoldEMresolvessmallmotions
Manifold embedding, or ManifoldEM, represents one of the first
approaches to the study of conformational heterogeneity in cryo-
EM6,25.Thismethodhasbeenusedtoresolveconformationaltransi-
tionsintheSARS-CoV-2spikeprotein26,27,theribosome6,andtherya-
nodine receptor28. However, despite previous versions of the code
beingmadeavailable28–30,themethodsadoptionwasnotwidespread
duetoitsslowspeedandlackofstability.Inthiswork,wehaveuseda
modernPythonversionofManifoldEM,whichistwoordersofmag-
nitudefasterandactivelymaintained23.Unlikeothermethods,Mani-
foldEMfirstanalyzesimagesfromthesameprojectiondirection(PrD)
beforestitchingtheresultstogetheroverthefullorientationsphere.
PrDs with a number of images belowa user-specified threshold are
ignored. After an initial Euclidean distance calculation followed by
diffusion mapping, non-linear Laplacian spectral analysis (NLSA) is
then performed to project the diffusion map results back into the
interpretablespaceof2Dimages.InitialversionsofManifoldEMthen
requiredtheusertomanuallyannotatetheresulting2Deigenvectors
foreveryPrD,butamorerecentdevelopmentallowsthisprocesstobe
semi-automatedwithonlyafew‘anchornodes’definedbytheuserand
beliefpropagationofthemotionsdeterminedbythecorrespondence
of optical flow vectors30. Together, these algorithms determine the
globalconformationalcoordinateforthemotionofinterestinthree
dimensions.
Viewingthe2DNLSAeigenvectorsderivedbyManifoldEMreveals
heterogeneousmotionsconsistentwiththosefoundbyothertools.
WhenviewingPrDsderivedfromsideviewsofTRPV1,changestothe
conformationsofthecytosolicdomainsarevisible(Fig.5A).Whereas
thePrDsfromtopandbottomviewsofTRPV1allowonetoseetwistsin
the transmembrane domains (Fig. 5B). These results are consistent
with the conformational changes observed in 3DVA and 3DFlex.
Additionally,thefocusonprojectiondirectionsallowsManifoldEMto
detectmotionsnotvisibletoothermethods.Forexample,theNLSA
eigenvectorsfromthetopandbottomviewsofTRPV1showchangesin
intensitywithinthepore.However,whilethelargermotionsdepicted
in the 2D NLSA movies appear consistent with other methods, the
volumesproducedinthefinalstagesofManifoldEMarenotofhigh
enoughresolutiontocapturechangeswithinthe pore,meaning we
cannotdiscountthepossibilitythatthesechangesmaybeanartifact
leftoverfromsignalsubtraction.Itislikelythatthereasonforthelower
qualityvolumesproducedbythistoolisthat,duringpreprocessingfor
thismethod,anyPrDswithfewerthan399imagesarediscarded.By
increasingtheminimumthresholdto400asopposedto200,itwas
foundthatmorenoisyPrDswereexcluded,leadingtoanincreasein
thequalityofthereconstructedvolumes.However,thisthresholding
procedureresultsinthelossofroughlyhalfoftheparticlestackfrom
theanalysis.Eventhoughtheyarelowerresolutionthanothermeth-
ods,thetrajectoryofvolumesatdifferentlocationsoftheManifoldEM
latent space displays changes in the density of DkTx,similarto the
resultsofcryoDRGN(Fig.5C,D,SupplementaryMovie8).Intriguingly,
thesecompositionalchangeswerenotobviousintheindividualNLSA
eigenvectorsofindividualPrDs.

MDsimulationsrevealmotionssimilartothosefoundin
cryo-EM
Themethodswehavestudiedinprevioussectionshavereliedsolely
ontheanalysisofcryo-EMimages.Theymadenoassumptionsabout
the underlying physics of the molecule under study. By contrast,
numerical simulation techniques like molecular dynamics (MD) use
physicalpriorstogenerateadistributionofstatesforabiomolecule.
Such simulations are often used to investigate the dynamics of ion
channels31. Hence, we ran atomistic MD simulations of the TRPV1
homotetramer bound to DkTx and RTX in a POPC bilayer as an
orthogonalmethodtocomparethemotionsfoundbythepreviously-
describedimage-basedmethods(Fig.6A).Thesesimulationsrevealed
significant motions that offer interesting insights into the allosteric
nature of the TRP channel. Allostery is known to be crucial for the
functionofTRPchannels32,33.SinceMDoffersdirectaccesstoatomistic
details,severalalgorithmshavebeendevelopedtodetectanddescribe
allostericnetworksfromsimulations.Here,weemployedCorrelation
ofAllRotamericandDynamicalStates(CARDS)forthispurpose.This
programusesdihedralanglestoclassifytheconformationalstateof
amino acids as being in either a disordered or ordered state.
Throughoutthesimulation,thecategoriesofaminoacidsareupdated,
andtransitionsoftheaminoacidsbetweenthesestatesaretracked,
allowingthesoftwaretodeterminewhichgroupsofaminoacids
exhibit the kind of correlated motions we expect from allosteric
networks34.Usingthissoftwaretoanalyzethe5×2μssimulationsof
TRPV1runat310KrevealedthattheTRPboxservesasapathwayfor
the communication of allosteric information between the cytosolic
andtransmembranedomains(Fig.6B).Thiscentralimportanceofthe
TRPboxtoallostericcommunicationisconsistentwithfindingsfor
other members of the TRPV family35,36. In addition to the allosteric
communication between the cytosolic domains and the transmem-
brane domains, CARDS also revealed allosteric coupling between
adjacentcytosolicdomains.Thisisprimarilymediatedbyaminoacids
at the interface between the C-terminus and the ankyrin repeat
domains in the N-terminus. Close inspection of the ankyrin repeat
domain suggests that residue K240 and its neighbors are chiefly
responsible forthisinter-subunit allostericcommunication (Supple-
mentaryFig.S6B).
TocomparetheconformationalspaceofourMDsimulationsto
theresultsofthecryo-EMheterogeneityanalysismethods,weapplied
PCAtoourMDsimulations.Thisrevealeddeformationstotheprotein
similartothoseseenwith3DVAand3DFlex(SupplementaryMovies9,
10). In the first two principal components of the MD simulations,
opposing cytosolic domains exhibit correlated motion—when one
cytosolicdomainmovesinward,sodoestheoppositeone.Thisisthe
samebehaviorimpliedbytheresultsofthecryo-EMmethods.InPC1,
we see an asymmetric motion in the cytosolic domains, like those
foundinthefirstmodeof3DVAandthesecondin3DFlex.Whilein
PC2,evenwithDkTxbound,weseetheTMDstwistandthecytosolic
domains lengthen, as with the second mode of 3DVA and the first
mode of 3DFlex. This is consistent with an oscillation between the
toxin-boundopenstateandaconfigurationsimilartotheapoclosed
state(Fig.6D,SupplementaryMovie9).However,duetotheatomistic
nature of MD simulations, we have access to much more detailed
motionscomparedtoexistingcryo-EManalysismethods.Forexample,
calculatingthedistancebetweenthecenterofmassoftheproteinand
the cytosolicdomains reveals thatopposite subunits in factexhibit
weakly anti-correlated motion (Supplementary Fig. S6B, C). That is,
when one domain moves inward, the opposite one tends to move in
outward.AsPCAisadimensionalityreductiontechnique,theanalysis
ofindividualprincipalcomponentslosessomeoftheinformationfrom
the raw trajectories. The deficiencies of PCA in particular are well-
known in the MD simulation community37–39. Hence, we inspected
higher order PCs to test whether we could find any anti-correlated
eigenvectors.Indeed,PC3and4oftheMDsimulationsexhibitanti-
correlated motions between opposing cytosolic domains (Supple-
mentaryFig.S7A,B,SupplementaryMovie5).Accordingtothesizeof
theireigenvalues,theselattertwocomponentseachexplain7%ofthe
varianceinthetrajectory,whichissignificantcomparedtothefirstand
secondcomponents,whichaccountfor20%and9%ofthevariance,
respectively (Supplementary Fig. S7C). These higher-order
components resolve the apparent discrepancy between analysis of
therawMDtrajectoriesandtheinspectionofthefirsttwoPCs,high-
lighting how the use of dimensionality reduction may obscure
importantdetailsinthestudyofmolecularmotions.
AnothernotableaspectofourMDsimulationsisthatweobserved
thecollapseofthehydrophobiclowergatenearthebeginningofeach
simulation (Supplementary Fig. S5). However, this is a well-known
phenomenoninMDsimulationsofchannelswiththesixtransmem-
branehelix(6TM)architecture40–43,andtheothercharacteristicsofthe
ligand-bound open state were preserved during the simulation. We
summarizethisandtheresultsofthecryo-EMheterogeneitymethods
(Table1).Thistabledemonstratesthecorrespondencebetweenour
physics-basedMDsimulationsandtheimage-basedcryo-EMmethods.

Ensemblereweightingshowsadistributionofstructuresinthe
cryo-EMdata
Sofar,wehavebeencomparingtheresultsofheterogeneityanalysis
methods andMDsimulations byvisualinspectionoftheconforma-
tionalchanges ofTRPV1. However,with the MDsimulations nowin
hand,wehaveanavenueforamorequantitativeapproachtocom-
paringtheMDandcryo-EMimages—Bayesianensemblereweighting.
Inthismethod,thelikelihoodthatagivenstructurecamefromagiven
image is calculated12,44. The prior estimation for the populations of
differentstructures,whichcomesfromMDsimulations,isreweighted
toreflectthepopulationsofthisstructuralensembleinagivensetof
cryo-EMimages.Inthiscase,weusetheprincipalcomponentsfrom
MD,asdescribedabove,asourstatespace.Indeed,reweightingthe
MDdistributiontothatofthecryo-EMdataindicatesthatthemajority
of particles occupy a conformation close to the ligand-bound open
state,butwithasmallpopulationofparticlesclosertotheapoclosed
state(Fig.7A,B).
Oneadvantageofensemblereweightingisthatitrepresentsthe
distributiondirectlyintheinterpretableconformationspacedefined
bytheMDsimulations.Thisisincontrasttothelatentspacesofother
heterogeneity methods, which resemble simple Gaussian distribu-
tions. Another remarkable feature of the ensemble reweighting
methodisthatthedistributionsfromensemblereweightingseemto
convergewitharelativelysmallnumberofparticles:only5%ofthe
totalnumberofparticlesisnecessarytofindwhichstatesarepresent,
andwithonly23%ofthetotalnumberofparticles,theweightsofthe
differentstructureshaveconverged(SupplementaryFig.S2A).
However,thisinterpretationseemsinconsistentwiththeEMdata
analysis results, which indicate large conformational heterogeneity.
However,anuancedaspectoftheensemblereweightingmethodin
cryo-EMisthechoiceofwhichspecificstructurestousetorepresent
thefullMDensemble.Ideally,onerequiresrepresentativestructuresof
the entire conformational space, which are then reweighted by the
data.However,duetocomputationallimitations,itisoftendifficultto
determinewhetherMDsimulationshaveadequatelysampledallrele-
vant conformational states. To enhance the likelihood of sampling
structures that correspond to conformations observed in cryo-EM
data,weperformedsimulationsofTRPV1atthreedifferenttempera-
tures:150K,200K,and310K.Wenotethatimprovingconformational
sampling remains an important direction for future applications of
ensemblereweighting,withapotentialstrategyinvolvingthedirect
biasing of MD simulations using information derived from cryo-EM
imagesashasbeendonewithotherformsofexperimentaldata45,46.
Whenstructuresfromsimulationsatallthreetemperaturesareusedin
reweighting,the structures derived fromlower temperaturesimula-
tionsreceivedhigherweights(SupplementaryFig.S2B,C).According
toSupplementaryFig.S3,thesestructuresdeviatedmuchlessfromthe
startingconfiguration,whereTRPV1isclosertoanopenstatecom-
paredtothesimulationsat310K.Thisindicatesthatthemajorityof
particles are likely to be in a configuration closer to the cryo-EM
derived consensus structure, while a portion are in configurations
closertotheapo-closedstate.Theseresultsareconsistentwiththe
presenceoftheconformationalstatesexploredby3DFlexand3DVA,
butnowprojectedontotheinterpretableMDlatentspace.

ComparingheterogeneitymethodswithSingularValue
DecompositionandspectralanalysiswithAnaVox
Ensemble reweighting allows us to compare the MD simulations to
cryo-EMimagesdirectly.However,itisalsodesirabletocomparethe
resultsproducedbydifferentheterogeneitymethodswitheachother.
Wehavedesignedasimplepipelinetostudymotionscapturedfrom
differentmethodsdirectlyfromaseriesofvoxelmaps,whichwecall
AnaVox.Thistoolcanbeusedtocompareanytwovolumeseriesto
determineiftheycontainthesamemolecularheterogeneity.However,
herewewilluseittocomparevolumeseriesderivedfromtheeigen-
vectorsofthelatentspacesofdifferentmethods.Wefirstselectoneor
twovolumeseriespermethod,e.g.,atrajectoryofvolumesalongPC1
for 3DVA. Each volume is then masked, normalized, and the mean
volumeissubtractedfromeachmemberoftheseries.SingularValue
Decomposition(SVD) isthen performed onthe series of processed
volumes,resultinginabasisforeachindividualseries.Then,wetake
these eigenvectors and perform a second round of SVD to find a
commonbasisforalltheseries.BydoingtworoundsofSVD,wecan
compare the different series both quantitatively and qualitatively
(Fig.8A).Thefirstroundletsuscompareseriesquantitatively.When
weprojecttheeigenvectorsfromeachseriesontooneanother,wecan
calculateaquantityknownasProportionofCapturedVariance(PCV)13.
Thismeasureshowsimilartwogivenbasissetsare,withaPCVof1
indicatingthatthebasissetsareidentical,whileaPCVof0indicates
thatallthebasisvectorsinthetwosetsareorthogonal.Thesecond
round ofSVDallows us tocompare the seriesqualitatively.Bypro-
jectingvolumesintothecommonbasis,wecanseewhethertheymap
tosimilarregionsinthesharedspace.
Usingthispipeline,wecanclearlyseethatvolumessampledalong
thesecondlatentcoordinateof3DFlexarecloselyrelatedtovolumes
generated from PC1 of our MD simulations. Visually, these volume
seriesareverysimilar(SupplementaryMovie10).Whenprojectedinto
a common basis, they trace similar paths along the first twoeigen-
vectors,andtheirlog(PCV)is-0.75.Thisistheclosestmatchwehave
found between two series from different methods (Fig. 8B). This is
consistentwithourobservationsfromtheMDsimulations,confirming
thatboththesevolumeseriescaptureanasymmetricsqueezingofthe
cytosolicdomainsofTRPV1.
Apairwise comparisonofall methods using PCVdemonstrates
thatserieswhichcapturethesamequalitativechangestothemolecule
havemoresimilarSVDbases(Fig.8C).Ingeneral,serieswhichcapture
theasymmetricsqueezingofthecytosolicdomainssharehigherPCV
scoreswitheachother(MDmode1,3DVAmode1,and3DFlexmode2)
than with the series which capture changes to the height of the
molecule(MDmode2,3DVAmode2,and3DFlexmode1).Similarly,
we see that series that capture changes to the occupancy of the
double-knot toxin (cryoDRGN, ManifoldEM, and 3DVA) have higher
PCV scores with each other than with series that do not capture
changesinthetoxin(MDsimulationsand3DFlex).Thesetrendsare
alsoapparentfrominspectingthepathofallseriesinthecommon
basis,whereitisclearthattheseriesthatcapturethefadingofthe
toxingenerallyspanalargerregionofthecommonbasisthanthose
that describe only conformational heterogeneity (Supplementary
Fig. S10). Here, we can also see that linear eigenvectors for atomic
coordinatesmaybecomenon-linearwhenprojectedintovoxelspace
(SupplementaryFig.S4A,B).
InterpretingsummarystatisticslikePCVisnontrivial.Forvoxel
mapsinhigh-dimensionalspaces,itismorelikelyforrandomlychosen
unit vectors to be orthogonal rather than parallel and have a dot
productclosetozero.Hence,tohelpinterprettheabsolutemagnitude
ofthePCVmetric,wehavedesignedsimpleexperimentstoleexperimentstolearnhow
PCVchangesindifferentsituations.Byaddingwhitenoisetoaseriesof
volumes,wecanusePCVtocomparehowtheeigenvectorschangeasa
functionofthenoiselevel(SupplementaryFig.S15A).Wefindthata
log(PCV) of -2 would correspond to a signal-to-noise ratio (SNR) of
between10−1and10−2.ThisPCVscoreisthesameaswhatweobtained
when comparing modes 1 and 2 of 3DVA, a method that gives us
orthogonal components in volume space (Fig. 8C, Supplementary
Fig.S10),indicatingthatthereisnosignificantcorrelationbetweentwo
volumeserieswithlog(PCV)<−2.Forexample,inthecaseof3DVA
mode1(SupplementaryFig.S15B)ataSNRof10−1theoutlineofthe
original eigenvector is still visible; however, by 10−2 it is no longer
visible. It is also of interest what the contributions of different fre-
quenciesaretothePCV.Wefoundthatlow-passfilteringvolumesdoes
notgreatlychangetheirrelativePCVs,uptoalow-passfilteringof14Å.
Thisisconsistentwiththefindingthatlinearbasesofvoxelmapsare
most sensitive to low-resolution features (Supplementary
Fig.S15C,D)24.
WiththeSVDpipelinedemonstratedhere,itcanbedifficultto
distinguishbetweenconformationalandcompositionalheterogeneity.
To study compositional heterogeneity, one could filter the refined
map to detect contrast loss, but it is not simple to expand this
approachtoaseriesofvolumes47.Forourpurposes,wecaninspectthe
radiallyaveragedpowerspectraofavolumeseries.Sincetranslational
degreesoffreedomwillbedegenerateintheFourierbasis,aseries
whichonlycapturesconformationalheterogeneitywillexhibitsimilar
powerspectraforeachvolume(Fig.8D).Indeed,wecanseethatthe
methods which capture compositional heterogeneity from DkTx
(ManifoldEM,cryoDRGNand3DVA)havesignificantvariationintheir
powerspectracomparedtomethodswhichonlycaptureconforma-
tionalheterogeneity(MDsimulationsand3DFlex).
Togainanunderstandingofhowrobustheterogeneityanalysisis
toparticlecount,weapplied3DVAandcryoDRGNtoindependenthalf-
setsoftheparticleimages.Eachhalf-setcontained109,393particles.
Atfirst,wefoundthattheresultsbetweenthehalf-setsdifferedsig-
nificantly. By down-sampling the images and applying symmetry
expansion,theresultsbetweenthehalf-setsbecamemoresimilarto
eachotherandtheresultsobtainedfromusingalloftheparticlesinthe
dataset.Withthispreprocessing,weusedAnaVoxtodemonstratethat
eachmethoddelineatesthesamemolecularmotionsineachhalf-set.
The volume series for the same method in different half-sets are
visuallyverysimilar,andtheyfollowsimilarpathswhenprojectedinto
acommonbasis(SupplementaryFig.S14A).Thisconfirmsthatthey
capturethesamemolecularmotionsineachhalf-set.Calculatingthe
PCVforthevolumeseriesfromdifferenthalf-sets,wecanseethatthe
samemethodfromdifferenthalf-setssharesahighPCV,confirming
thattheycapturethesamemolecularheterogeneity(Supplementary
Fig.S14B).However,thePCVscoresindicatethatthebasesfordiffer-
ent half-sets are not the same when comparing across different
methods,indicatingsomedifferencesbetweenthem.Forexample,the
first3DVAmodeinthefirsthalf-setismoresimilartothecryoDRGN
resultsthanthefirst3DVAmodeinthesecondhalf-set.Visualanalysis
ofthevolumeseriesandinspectionoftheirpowerseriesindicatethat
these differences are most likely due to low-frequency noise in the
region surrounding the nanodisc (Supplementary Fig. S14C). This
highlightstheneedformorerobustmethodstofilterstructurednoise
fromadatasetduringheterogeneityanalysis.

ComplementaryanalysisofaRibosomedataset
While covering the diversity of heterogeneity in cryoEM datasets is
outsidethescopeofthiswork,wehaveincludedananalysisofasec-
onddatasettogetasenseofthegeneralizabilityofourconclusions:
theP80Sribosome.Wechosetheribosomeduetoitshistoricalsig-
nificance in cryo-EM methods development48. We again used
cryoDRGN,3DVA,3DFlex,andManifoldEMtoanalyzeasingleparticle
datasetofthematurePlasmodium80Sribosome(EMPIAR-10028)49.
With each method, we observed the greatest changes in the small
subunit.Thereappeartobetwoconcurrentheterogeneousmotions.
Thefirstisa“rocking"ofthesmallsubunit,consistentwithprevious
resultsusingManifoldEM6.Thesecondiscompositionalheterogeneity
ofpartsofthesmallsubunitcenteredonthees10andes12modulesin
the “beak" of the ribosome50. As in TRPV1, each method returned
slightly different results (Supplementary Fig. S11A, Supplementary
Movie11).However,inthiscasecryoDRGNand3DVAreturnedresults
thatwereremarkablysimilar.Visually,thevolumeseriesarehardto
distinguish.WhenembeddedinacommonbasisusingAnaVox,thereis
a clear correspondence between the paths drawn by the first3DVA
mode,cryoDRGN,andthetwo3DFlexmodes(SupplementaryFig.11B).
Thiscorrespondencebetweenthevolumeseriesisalsonicelycaptured
bythePCV(SupplementaryFig.S11C).Asbefore,thevarianceinthe
volumeseriesthatcapturecompositionalheterogeneitydisplaylarger
varianceintheirpowerspectra(3DVA,CryoDRGN,andManifoldEM)
comparedtothosethatmostlycaptureconformationalheterogeneity
(3DFlex)(SupplementaryFig.S11D).Thisconfirmsthatmonitoringthe
powerspectraisareliablemethodtoassessthepresenceofcompo-
sitionalheterogeneityinagivenvolumeseries.

Methods
Preprocessingthesingle-particlecryo-EMimagestackforhet-
erogeneityanalysis
Signalsubtractionofthenanodisc.Forthisstudy,weanalyzedthe
publiclydepositedrawelectronmicroscopydataatEMPIAR-1005921.
Thisdatasetconsistsofmicrographsandparticlestacksofananodisc-
embedded TRPV1 protein bound to two agonists: a peptide toxin
known as double knot tarantula toxin (DkTx) and a vanilloid com-
poundknownasResiniferatoxin (RTX).Since TPRV1isa membrane
protein, it was reconstituted into lipid nanodiscs with
MSP2N2scaffoldspriortoimaging.Thisformsnanodiscswithadia-
meterofroughly150Å,significantlylargerthanTRPV1itself.Nanodiscs
ingeneralcanformatdifferentsizesinsolution,leadingtocomposi-
tionalheterogeneitybetweenparticleimagesofthenanodiscitself73.
Additionally,itisthoughtthatvariationsinthelocationoftheprotein
within the nanodisc alsoresultin apparent nanodisc heterogeneity.
Consistentwiththeseaspectsofmembraneproteinreconstitutioninto
nanodiscs, in several cases our initial heterogeneity analyses were
dominatedbythecompositionalheterogeneityofthenanodisc(Sup-
plementaryMovies16–18)74.Bycreatingasignificantlydilatedmaskto
accountfortheparticleswithlargernanodiscs,wewereabletouse
signal subtraction to decrease the artifacts from this nanodisc het-
erogeneity. Of the algorithms we tested, cryoDRGN and 3DVA
appeared to be the most sensitive to nanodisc heterogeneity (Sup-
plementaryMovies17,18).Similartothistrend,wenoticedinMani-
foldEMthatsignalsubtractionandmaskingofimagesincreasedthe
amount of motion detected by the algorithm (Supplementary
Movie19).Forcompleteness,weconducted3DVAandcryoDRGNon
both the signal-subtracted and non-signal-subtracted data. When
viewing the volumes at high thresholds and carefully exploring the
latentspacesofeachmethod,wefoundthesameheterogeneitywithin
thetransmembraneregion oftheprotein,indicatingthatthesignal
subtractionprocedureislikelynotbiasingourresults(Supplementary
Movie20).Theseresultsreflecttheimportanceofcarefulpreproces-
singfortheanalysisofheterogeneityinsingleparticlecryo-EM.

HomogeneousreconstructionandparticlepickingofTRPV1.Forall
mainfigures,thesourceparticlestackwastheparticlestackof218,805
particlesdepositedinEMPIAR-1005921.Afteraroundof2Dclassification
toremovejunkparticles,wewereleftwith218,717,whichweusedfor
allourheterogeneityanalysis.Usinghomogeneousandlocalrefine-
ment with C4 symmetry, we were able to reconstruct this signal-
subtracted particle stack to a resolution of 3.0 Å. However, experi-
mentswerealsoperformedwithanenlargedparticlestack.Toexpand
thecountofparticles,wereferredbacktotheoriginalmicrographsin
theEMPIARdepositionandutilizedablobpickerwitharadiusbetween
200and300Å.Thisyielded598,273particles.With2Dclassification,
wediscardedparticlesthatwereobviouslyjunk,butkeptparticlesin
classesthatwouldusuallybediscardedbecauseoftheirlowquality.
Thisyielded525,910particlesintotal.Thisenlargedstackwasrecon-
structedtoaresolutionof2.8Åusinghomogeneousrefinementfol-
lowedbylocalrefinementunderC4symmetryconstraints.
Wealsoperformedexperimentsusingdrasticallyreducedparticle
countsusingcryoSieve53.Ourstartingparticlestackforthisalgorithm
wastheexpandedparticlestackmentionedabove,andtheconsensus
modelfromthisparticlestackwasusedasareferencemodelforthe
cryoSieve algorithm. Frequency marching began from 40 Å and 8
iterations of the cryoSieve algorithm were applied with a retention
ratioof80%betweeniterations.Thisresultedinaparticlestackwith
27,572particles,reconstructionresolutionof3.6ÅusingC4symmetry
constraintswithhomogeneousandlocalrefinementincryoSPARC.

Processing and heterogeneity analysis of the Plasmodium 80S
Ribosome.ThepubliclyavailableparticlestackforthePlasmodium
ribosome from EMPIAR-10028 was reconstructed to 3.1 Å with
homogeneousrefinementincryoSPARC49.Allheterogeneitymethods
wereappliedusingthesamehyperparametersasthoseforTRPV1,with
theexceptionsof3DFlexandManifoldEM.
For3DFlex,weused5latentcoordinates,128hiddenunits,anda
latentcenteringstrengthof30.Thiscombinationwasfoundtogivethe
bestlocalresolution.ForManifoldEM,weusedanaperturesizeof4,an
objectdiameterof200Å,andaminimumthresholdonthePrDsof200
images.NoPrDswerediscardedforhavingtoomanyimages.

Heterogeneitymethods
3Dflexiblerefinement.The signal-subtracted particles weredown-
sampledto128pixels,resultinginapixelsizeof1.83Å.Atetrahedral
meshof754pointsand2517cellswasgeneratedtocoverthedown-
sampled volume. The 3DFlex model itself consisted of 2 latent
dimensionswith64hiddenunitsinthemulti-perceptronarchitecture.
Ameshrigidity(λ)of1andlatentcenteringstrengthof5wereapplied
tothemodelduringtraining.Wefoundthatreducingtherigidityofthe
meshfrom2to1servedtoincreasetheresolutionofthereconstructed
maps.However,addinghiddenunitsandextralatentdimensionsdid
littletoimprovetheresults.Noconsiderationofsymmetrywasgiven
to the 3DFlex model during training unless otherwise specified. In
ordertoproduceaseriesofvolumes,wesampled50volumesalong
eachlatentcoordinatedefinedbythe3DFlexmodelinFig.2.

3Dvariabilityanalysis.Theparticleswereanalyzedusing3DVaria-
bilityanalysisaftertheywerepreparedusingsignalsubtractionand
localrefinement.Afilterresolutionof5Åwasusedduringall3DVA
calculations; changing this value was not observed to have a
meaningful effecton the results of the algorithm. Unless otherwise
stated,noconsiderationofsymmetrywasgivenfor3DVAcalculations.
Inordertoproduceaseriesofvolumes,wesampled50locationsalong
the first and second latent coordinates of the 3DVA model defined
inFig.3.

CryoDRGN.CryoDRGNwastrainedwiththedefaultparameters.For
theunderlyingVAE,alatentspaceof8dimensionswasused.Boththe
encoderanddecodernetworksconsistedof3hiddenlayerswith1024
nodeseach.TheactivationfunctionfortheVAEwastheReLUfunction.
100epochsoftrainingwereperformedforthemodelsanalyzedinthis
study. No consideration of symmetry was given to the cryoDRGN
modelduringtrainingunlessotherwisespecified.Inordertoproducea
series of volumes, we sampled 50 volumes along the first principal
componentdefinedbytheVAElatentspaceofthemodelfromFig.4.
Training the VAE beyond 50 epochs was onlyobserved to produce
marginalimprovementsinthedetectionofheterogeneityinthesys-
tem.IncreasingthenumberoflatentdimensionsintheVAEto10only
servedtoproducemorenoisyartifactsinthevolumeseriesproduced
bytheVAE.

Manifoldembedding.ManifoldEMwasrunwiththefollowingproto-
col.Anaperturesizeof4andanobjectdiameterof160Åresultedin
sufficiently populated projection directions with no noticeable arti-
facts from global angularrotations. To avoid the inclusion of noisy
projectiondirectionswithtoofewparticles,projectiondirectionswith
fewerthan400particleswerediscardedfromanalysis.Noprojection
directionswerediscardedforhavingtoomanyparticles.Thisresulted
in89PrDsbeingincludedintheanalysis.The’sense’ofthemotions
was determined by eye. No consideration of symmetry was made
duringtheuseofManifoldEMunlessotherwisespecified.Inorderto
produceaseriesofvolumes,wesampled50volumesalongthecon-
formational coordinate discovered by the ManifoldEM calculations
inFig.5.
Manyattemptsweremadetoimprovethequalityoftherecon-
structed volumes from ManifoldEM, including low-passfiltering the
inputparticles,symmetryexpandingtheparticlestacktoincreasethe
number of samples, and tweaking the aperture to include more or
fewerparticlesineachprojectiondirection.Intheend,wefoundthat
themostusefulstrategyistouseahigherminimumthresholdtoavoid
theinclusionofoverlynoisyPrDs,andtoprefertheuseofPrDswhere
themotionofinterestisclearlyvisible.Thisensuresthattheusercan
correctlyidentifythe‘sense’ordirectionofthemotionofinterest.
AnactivelymaintainedversionofManifoldEMcanbeaccessedat
thefollowinglinkhttps://github.com/flatironinstitute/ManifoldEM23.

Moleculardynamicsandensemblereweighting
SystembuildingandMDsimulationprotocols.Theproteinstructure
forthisworkwasanatomicmodel(PDB:5IRX)associatedwiththeraw
singleparticledatasetwehaveanalyzed(EMPIAR-10059)21.Amodelfor
the cytosolic Ankyrin Repeat Domains (ARDs) (PDB ID: 2PNN) was
dockedintothelow-resolutionelectrondensityofthecryo-EMdensity
map (EMDB: 8117). The missing N and C-termini were built using
ColabFoldmodelsandsplicedintothefullTRPV1proteinaswell75,76.
TheconformationoftheC-terminuswasmodeledontheconforma-
tion found in PDB ID 7LP9, as this structure resolved contacts the
C-terminalandneighboringcytosolicdomains.ThestructureofTRPV1
in this dataset was solved in a nanodisc and four 1-Palmitoyl-2-
oleoylphosphatidylcholine(POPC)lipidsareboundtotheinteriorof
thetransmembranedomains.However,theyhavemissingtails.Forour
simulations,thesetailswerebuiltmanuallyinVMDtoensurethatthese
lipids and their polar contacts with the protein were preserved77.
Parameters for the bound Resiniferatoxin co-factors were derived
from CGENFF78. The default protonation state for all residues was
appliedforthisstudy.
Thisprotein-ligandsystemwasinsertedintoaPOPClipidbilayer
andsolvatedusingtheCharmm-GUIwebserverwithTIP3Pwater79.The
system was neutralized with sodium and chloride ions at a con-
centration of 0.15 mol/L. The resulting simulation box was roughly
cubic,withasidelengthof16nm.Thecompletemoleculardynamics
systemcompositionissummarizedinTable2.
Minimizationviathesteepestdescentalgorithmwasperformed
untiltheforceonallatomswasbelow24kcal/mol/Å.Thiswasfollowed
by relaxation, wherea restraint of10 kcal/mol/Å2 was placed on all
heavyatoms,andthenthemagnitudeofthisrestraintwashalvedevery
200ps,in15iterations.Relaxationandproductionwererunwith1and
2 fs time steps, respectively. Relaxation was followed by 100 ns of
equilibration.Duringequilibration,watermoleculesandsodiumions
naturallymigratedintotheporeofthechannel.Duringrelaxation,a
Berendsenthermostatandbarostatwereapplied,andforproduction,
aNosé-Hooverthermostatwithacouplingtimeof1psandaParrinello-
Rahman barostat with a coupling time of 5 ps were applied during
productionrunsofoursimulations.Simulationswerecarriedoutat
atmospheric pressure. Each simulation replicate was minimized,
relaxedandequilibratedindependently.Productionrunswereexten-
ded up to 2 μs at 310 K. Inspection of the C-alpha RMSD in these
production runs is consistent with an equilibrated simulation (Sup-
plementaryFig.S17).GROMACS2022withtheCHARMM36mforce-
field with a cut-off distance of 1.2 nm was used for all molecular
dynamics simulations in this publication80,81. The trajectories them-
selvesarepubliclyavailableaswellathttps://doi.org/10.5281/zenodo.
15132985.

AnalysisofallostericnetworksusingCARDS.Networksofallostery
werestudiedinthesimulationsusingCorrelationofAllRotamericand
DynamicalStates(CARDS)withabufferwidthof15degrees20.Sym-
metry expansion was performed on these trajectories before this
analysiswasperformed.Thenatureofsymmetryexpansionisoutlined
inthenextsection.
CARDSmeasuresthecorrelationbetweenconformationaltransi-
tionsofpairsofdihedrals.Thiscanbeperformedforeverydihedralin
aproteintogainanunderstandingofwhichdihedralsarecontributing
to the communication of allostericinformation. To assign a mutual
informationquantitytoeachresidueasinFig.6B,weaddedthecon-
tributionsforeachdihedralinagivenresidueandtooktheaverage
mutualinformationacrossalldihedralsinthatresidue.Torankamino
acidsbyhowmuchallostericinformationtheypassthroughthesys-
tem,weperformedarowsumonthematrixinFig.S6A.Thisgivesa
spatialunderstandingofwhichaminoacidsarecontributingthemost
tothecommunicationofallostericinformation.

PCAonmoleculardynamicstrajectories.Typically,PCAofMDtra-
jectories is dominated by the motion of large disordered loops,
becausetheyhavealargevariance39.However,sincethemeanstruc-
ture of these disordered domains is not a good estimate for their
positioninanygivenframe,theseportionsofaproteincannotusually
beresolvedfromelectrondensities.Hence,tohelpignorethesedis-
orderedfluctuations,beforeperformingprincipalcomponentanaly-
sis, we normalized the data by subtracting the mean position from
eachatomanddividingtheirpositionineachframebytheirstandard
deviationoverthecourseofthetrajectory.Wefoundthatthisnor-
malization procedure emphasizes collective motions, meaning that
isolatedmotions,suchasthosefoundindisorderedloops,willnotbias
theanalysis82,83.Thissameprocedurehasbeenusedinanotherrecent
comparisonstudybetweencryo-EMandMDsimulations84.Inourcase,
a comparison between eigenvectors derived from normalized and
unnormalized coordinates reveals that motions in segments of the
proteinwithnosecondarystructure,liketheC-terminusandresidues
F720-W740 are deemphasized by this normalization procedure.
However,theoveralldifferenceisminor,andthesameconformational
changes are found for corresponding eigenvectors (Supplementary
Movie21).
Inthecontextofcryo-EM,theconventionalschemeforapplying
PCA to molecular dynamics has some deficiencies. Primarily, it
assumesthatthereisnosymmetryintheparticleofinterest.Hence,we
developedasimpleschemetoapplythenotionofsymmetryexpan-
siontothetrajectory.ThetrajectoryisclonedNtimes,whereNisthe
orderofthesymmetrygroupofinterest(inthiscase,4timesforC4
symmetry).Ineachoftheseclones,thechainlabelsareswappedto
simulatethesituationthatonechainratherthananotherunderwenta
largefluctuation.
Dimensionality reduction is a critical part of the analysis of
moleculardynamicstrajectories.PCAisapopularchoicebuthaswell-
knownlimitations37–39.OnepracticaldifficultywhenusingPCAsisthat
largedisorderedloopstypicallydominatetheanalysis,astheyproduce
the most variance. This is a particular problem when comparing
eigenvectorstocryo-EM,becausesuchloopsarenotvisible.Inorderto
limit the influence of disordered motion in our case, the following
procedureofnormalizationwasemployed:
TheCartesiancoordinatesoftheatomsinthesystemareplacedin
a3Nvectorattimeτ
xτ=ðxτ
1
,yτ
1
,zτ
1
,:::,xτ N,yτ N,zτ
N\Þ: ð1Þ
Thisletsusdefineavectorofmeansandvariancesforeachcoordi-
nate.
μ = ð (cid:3) hx q1 i ffitffiffi , ffiffi h ffiffi y ffiffi1ffiffi i ffitffiffi , ffiffi h ffiffi z ffiffi1ffiffi i ffiffitffiffi , ffiffi h ffiffi x ffiffiNffiffi i ffiffitffi , ffiffi h ffiffi y ffiffiffiNffiffi i ffiffitffiffi , ffiffi h ffiffi z ffiffiffiNffiffiffiffi i (cid:4)t Þ
σ = hx2 1 i t (cid:2)hx 1 i2 t ,:::, hz2 N i t (cid:2)hz N i2 t ∈[1,2.6]μm.
Eachxτ
i
isthustransformedmτ
i
=ðxτ
i
(cid:2)μ
i
Þ=σ i.Thesemτarestacked
intoamatrixMtocalculatethecorrespondingsingularvaluedecom-
position(SVD).
2 3
mτ, mτ, mτ, ::: mτ, mτ, mτ
1 1 1 N N N
M= 6 6 6 6 4 m2 1 . . . τ, m2 1 τ, m2 1 τ, . : . :: . m2 N τ, y2 N τ, m . . . 2 N τ7 7 7 7 5 ð2Þ
mTτ, mTτ, mTτ, ::: mTτ, mTτ, mTτ
1 1 1 3N 3N 3N
andthusweperformtheSVD
M=USVT ð3Þ
Following this procedure, the columns of VT will contain the
principal components of M. To visualize these eigenvectors, we
selectedthefirsttwocolumnsofthematrixV,whichcontaintheunit
eigenvectorsforprincipalcomponents1and2,andmultipliedthemby
theirrespectivesingularvalues.Thisallowsustovisualizetherelative
scaleofvariancecapturedbyagiveneigenvectoraswellasitsdirec-
tionofmotion.
To produce an interpretable sample of these eigenvectors in
Cartesian space, we created several atomic models by multiplying
thesescaledeigenvectorsbyascalarparameterbetween-1and1at50
equally spaced intervals. This allowed us to view motion along the
eigenvectorsasacontinuousdeformationintheinterpretableCarte-
sianbasis(SupplementaryMovie9).The50samplesofthesescaled
eigenvectorscanalsobeprojectedbackintotheprincipalcomponent
basisandplotted(SupplementaryFig.S4A).Codefortheconstruction
of interpretable eigenvectors in atomic space be found at https://
github.com/flatironinstitute/astore_natcomms_2025_scripts/blob/
main/pca.py.ForanalysisbyAnaVox,theseatomicmodelswerecon-
vertedintovolumesusingthemolmapcommandinChimeraX85witha
resolutionof8Å.

Calculatingmotionofthecytosolicdomains.Asaninterperatable
measurementoftheconformationalchangestoTRPV1,wecalculated
theverticalandhorizontaldisplacementofthecytosolicdomainsin
boththemoleculardynamicstrajectoriesandduringprojectionofthe
eigenvectors of the PCA (Fig. 6C, D). To track the horizontal dis-
placementofthesedomains,theradialdisplacementinthex-yplane
betweenthecentersofmassofresiduesinthecytosolicdomainswas
comparedtoareferencepointdefinedasthecenterofmassforthe
TMDs.The relevantresidues in the cytosolicdomains were114-183,
whiletheTMDsweredefinedasintervalsonresidues418-454,469-498,
511-532,536-554,563-598,630-642,and 656-712in order to exclude
disorderedloops.Similarly,totracktheverticaldisplacementofthe
cytosolic domains,wecalculated the difference in the z coordinate
betweenthecentersofmassinthecytosolicdomainsandtheTMDs,
witheachdomaindefinedwiththesameresiduesasbefore.Scriptsfor
performingthesecalculationsareavailableonGitHub:https://github.
com/flatironinstitute/astore_natcomms_2025_scripts.

Ensemble reweighting. Ensemble reweighting leverages experi-
mentaldatatoreweighapriorconformationalensemble,extracting
theconformationalprobabilitydensitysuggestedbythedata.Inthis
study, we employed Bayesian analysis and the cryo-EM particles, in
accordancewith the mathematicalframework detailed in ref. 12, to
reweightthe310KensemblefromtheMDsimulations.Wecalculated
theimage-to-structurelikelihoodusingBioEM86,87.Weconductedan
initialroundofBioEM,asdescribedinref.88,todeterminethebest
orientationsforeachimageusingthereferencePDBstructure,36000
quaternions uniformly distributed in SO(3) and variable defocus
∈[1,2.6]μm.Subsequently,foreachimage,wegeneratedalistof8000
quaternionswithafinerangularresolution(approximately2degrees)
around the top 300 orientations from the first round. For the 120
clustercentersobtainedfromallMDsimulations,weemployedthis
refinedlistofquaternionswithBioEM,incorporatingthebestdefocus
valuesfromthefirstroundandallowingforvariablecentertranslations
of 60 pixels along both axes. To prevent misalignment, a low-
passfilterof10Åwasappliedtoallstructures.WeusedtheBayesianframe-
work (specifically Eq. ref. 12) to infer the weight of each cluster
center.ThisinferencewasperformedusingMarkovchainMonteCarlo
(MCMC)with8replicasand20,000MCMCsteps,utilizingthecode
available at https://github.com/flatironinstitute/Ensemble-
reweighting-using-Cryo-EM-particles.Errorbarsarederivedfromthe
standarddeviationoftheMCMCsamples.Theinferredweightswere
then distributed to the structures from the 310K MD ensemble by
nearestneighborclusteringinPC1andPC2space.

ClusteringofMDtrajectories.Duringthereweightingprocedure,a
subset ofstructures is chosen fromthe full ensemble to limit com-
putationalexpense.Statesaretypicallychosentoberepresentativeof
the full ensemble; details can be found in ref. 5. For this study,
structureswerechosenusingK-medioidsclusteringonthePC1and
PC2variablesdefinedinthesection.100clustercenterswereselected
fromthesimulationsrunat310K,and150K,and200K,whileeachcluster
centerwasrepresentedbythecentroidofitscluster.

ComparingseriesofvolumeswithAnaVox
Embedding volumes in a common basis with singular value
decomposition.Singularvaluedecomposition(SVD)isalinearalgebra
technique.Inshort,SVDgivesusanequivalentdecompositionforboth
the rows and the columns of a given matrix. If we use SVD to
comparethesimilarityofvolumeseries,thenwecanrecoverPCA
fromdifferentheterogeneitymethodsbyembeddingallvolumesina
commonlow-dimensionalspace.
Torunourpipeline,westartbyflatteningandnormalizingall
volumes in the series and their mean-removed versions, where the
number of volumes in each series is the number of mean-removed
volume,withthevolumesinthereseriesareonthesamescale.Thenwe
remove the mean volume from each series, that is, for every method
we define the mean-removed matrix ½D i (cid:3) j =X ij (cid:2)N(cid:2) v 1 N j= v 1 X ij, with D i 2RN v×1923 .
Lastly, we compute an individual economy-SVD for each of these
matrices, that is, we obtain D i=U i S i VT i, where U i,S i 2RN v×N v,
andV i 2R1923×N v.

Comparing the SVD of different volume series. To compare the
different bases obtained by SVD, we adapt a metric introduced in
ref.13,namedtheproportionofcapturedvariance(PCV).Thismetric
quantifiestheproportionofcapturedvariancebetweenthesubspace
spanned by the left singular vectors (V i) of two volume series. Let
(U i,S i,V i)and(U j,S j,V j)betheSVDfromtheithandjthvolumeseries,
respectively.ThePCVisdefinedas
PCVði,i Þ= i 1 i 2 i 2 F = l,k i 2Pk i 1 l i 2 k ,
1 2 kVT i 2 V i 2 S i 2 k2 F k ½S i 2 (cid:3)2 k
wherekAk2= A2 denotestheFrobeniusnorm.Thismetrichasa
maximum value of 1, and quantifies how well one subspace covers
anotherone,takingintoaccounthowimportanttheeigenvectorsare.
Notice that this metric is asymmetric, meaning that instances exist
wherePCV(A,B)<PCV(B,A),meaningthatAcoversB,betterthanB
coversA.AnillustrativeexampleoftheasymmetricPCVmetricisgiven
inSupplementaryFig.S13.

Finding a common embedding for a collection of volume series
withAnaVox.WhenwecomputetheSVDforasinglevolumeseries,
we obtain the principal components ordered based on their
varianceacrosstheseries.Nowweareinterestedincomparingthe
features across different series. To do this, we perform a second
roundofSVDusingtheeigenvectors(V i)foundpreviously.Fromthe
first round, we have that we can represent any mean-removed
volumeas
½D i (cid:3) j = ½U i (cid:3) jk ½S i (cid:3) kk ½V i (cid:3) k ,
k=1
where ½V i (cid:3) k 2R1923 is the k-th left singular vector. Now, we can
perform a second round of SVD by weighting all the V i ’s by their
respective singular values ½S i (cid:3) kk , stacking them in a matrix Y, and
computing its SVD
Y=UeeSVe
. Then, we can represent any weighted
eigenvector(i.e.,arowofY)as
Ve
l=1
whereY q=½½S i (cid:3) kk V i (cid:3) k forsomei∈[N series ],andk∈[N V].Lastly,wecan
combine both series to represent the original volumes in the basis
definedbyVe
,whichgivesus
½D i (cid:3) j = XN v ½U i (cid:3) jk N v X×N series Ue kl eS ll Ve l
k=1 l=1
= series v ½U i (cid:3) jk Ue kl eS ll Ve l :
l=1 k=1
Here, the vectors
Ve
l form a common basis that we can use to
projectPany volume from our different series. The coefficients
bythesecondroundofSVD.

Detecting compositional heterogeneity by inspecting power
spectraofvolumeseries.DuetotheinsensitivityofaFourierbasis
totranslationaldegreesoffreedom,aserieswhichonlycaptures
conformationalheterogeneitywillexhibitsimilarpowerspectrafor
eachvolume.Indeed,wecanseethatthemethods which capture
compositional heterogeneity from DkTx (ManifoldEM,cryoDRGN
and3DVA)havesignificantvariationintheirpowerspectracompared
tomethodswhichonlycaptureconformationalheterogeneity(MD
simulationsand3DFlex).Tocomparethesechangesbetweenvolume
series,wefirstcalculatetheradiallyaveragedpowerspectrumPSi(k)
foreachvolumeiinaseriesi=1,…,N vwithktheradialreciprocal
component and N v the total number of volumes in the series.
Then, we normalize those power spectra according to the mean
of the values at k = 0, that is PSðkÞ i=PSðkÞ i =PSð0Þ i, and calculate
the vPariance of the series for each radial frequency
σ2ðkÞ=
i
ðPSi ðkÞ(cid:2)PSðkÞÞ2 where PSðkÞ is the mean of the power
spectraoftheseriesatradialfrequencyk.InFig.8D,wepresentthe
totalpowerspectruPmvariancethatisthesumofthevariancesover
frequencies,i.e., N 1 v k σ2ðkÞ,foreachvolumeseries,normalizedby
thenumberofvolumesintheseries.Byaddinguptheamountof
variancefortheseriesacrossitspowerspectrum,wecancompare
the relative amounts of compositional heterogeneity in different
volumeseries.
TheAnaVoxpipelinecanbeaccessedatthefollowinglinkhttps://
github.com/flatironinstitute/AnaVox.Volumeseriesfromeachanaly-
sismethodandtheMDtrajectoriescanbefoundonZenodo:https://
doi.org/10.5281/zenodo.13900835.

Synthetichomogeneousdatasetgenerationandevaluation
Tostudytheeffectsofabiasedposedistributiononthedetectionof
heterogeneitybyavailabletools,wecreatedasyntheticdataset.This
datasetcontainedimagesofasinglestructure(thestartingstructure
usedinourMDsimulations),butimageparametersweretakenfrom
the poses and defocus values derived from the consensus recon-
struction of the experimental TRPV1 dataset. We created 218,786
images with a signal-to-noise ratio of 0.01 using the cryoJAX
simulator89.
Leveraging the ground truth poses and defocus values used to
createtheimages,thissyntheticdatasetreconstructedto2.5Åwith
simplebackprojection(SupplementaryFig.S16A).Wethenran3DVA
withvariouslevelsoflowpassfilteronthisdataset.Wefoundthatalow
passfilterof8Åormoreproducedsignificanthallucinatedcomposi-
tional heterogeneity (Supplementary Fig. S16B, Supplementary
Movie15).

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In this manuscript, the authors applied four different CryoEM heterogeneity analysis methods, as well as MD simulation, on a TRPV1 dataset and compared their results. Conceptually, this provides a guide for structural biologists to choose the right method to analyze the structural heterogeneity within their systems. However, as a paper about the comparison between software, it seems necessary to include more than only one example, especially since the structural variability in TRPV1 is dominated by the lipid nanodisc fluctuation, which can cause difficulty in some methods, even after subtraction. The MD simulation and the SVD-based comparison sounds reasonable, but they lack details and validation and it is hard to judge how useful they are. Overall, I think the paper still needs some extra work to be a useful reference for the researchers in the field.\nHere are my main comments.\nFor 3DFlex, the author claimed to solve some cytosolic helices to below 3 A. However, there are no clear sidechain features visible on those helices in Figure 2. Some zoomed-in figures are needed to show those are real high resolution features instead of errors in local resolution measurement.\nMap quality in figure 4 is not as good as that in figures 2 and 3, and the map in figure 5 looks even worse. It is unclear whether it is just a difference in filtering, or an actual issue in software. It is recommended to filter the maps to the same (or at least similar) level so they can be compared visually.\nFor the CryoDRGN results, it is quite strange to observe toxin fading gradually, instead of simply being present/absent. Is the fading density caused by subtle movement of toxin, or the imperfection of the software? As a compositional heterogeneity, it can also be analyzed by traditional focused classification. It would be a good way to quantitatively analyze the latent space distribution.\nThe results from ManifoldEM are not very convincing. Does it also produce other conformational changes that are more visible, like the movement of cytosolic helices shown in figure 2? It would be more convincing if the authors showed some large-scale movements that are clearly visible first to build credit for the tool, before showing the less visible small motion.\nThe MD simulation part sounds surprisingly simple for a membrane protein. The text refers to figure S2, which is not very informative. Are the lipids from the nanodisc considered in the simulation? How are the water molecules included in the simulation, particularly the water inside the channel?\nThe MD simulation should provide much deeper insight into the dynamics of the system, instead of just being visually similar to the CryoEM analysis. For example, if the movement of the cytosolic helices from different subunits are truly correlated, what is the driving force that leads to the coordinated movement of helices from the opposite side, hundreds of Å apart? Since those helices are not interacting directly with each other, presumably there are residues at the transmembrane domain responsible for the coordination? Is it possible to identify them, maybe even mutate them in silico and show the correlation of cytosolic helical movement disappearing as a result?\nAdditionally, in table 1 authors claim the MD simulation shows the pore closing movement, which was not shown in figure 6. What is the scale of the movement? It has previously been reported the pore opening is related to simulation temperature. Is it also the case here? What is the evidence that the pore closing movement from MD simulation corresponds to the movement observed from ManifoldEM? From my understanding, the gating of most channels are controlled by sidechains and the change is only a few Å in radius, but the maps in figure 5 lack high resolution features. Or did the authors discover a pore closing movement at a much larger scale through their MD simulation?\nFor the SVD based comparison, there is a lack of a baseline to indicate the significance of the correlation. For example, what would be the correlation values if one generates many maps by grouping particles randomly, or just adding noise to the maps? Otherwise, it is impossible to know whether a correlation score of 0.06 means a significant agreement between two methods.\nFor SVD comparison, what has the author done to avoid the different filtering in the maps, as well as the SNR caused by particle count? Would the correlation values be different if we simply lowpass filter every map to 5 or 10Å? Would it matter if only a subset of particles were used in the reconstruction?\nAlso for the SVD based comparison, it seems that only movement trajectories along eigen-vectors are considered? What if one method describes a movement as a straight path and the other as a curve in the latent space?\nTo further validate the comparison method, maybe break the datasets to two halves and analyze the heterogeneity in each set separately. This will show whether the correlation is consistent and what is the noise level in the comparison. It still won’t overcome any potential systematic bias but should be a step forward.\nThe nanodisc subtraction strategy might be problematic. Typically for focused alignment, we remove the density of rigid parts by subtracting the projections of averaged structures from the particle images. In this case, we already know the nanodisc is highly variable, so the average structure does not reflect the actual nanodisc density in each particle, resulting in leftover density of different sizes and shapes. Even though it does not have a major impact on homogeneous reconstructions, they may be interpreted as structural variability, and might be the main cause of the heterogeneity at the transmembrane domain reported by CryoDRGN and ManifoldEM.", "reviewer_2": "This study provides a survey of tools designed to extract information on conformational heterogeneity from cryo-EM data sets, specifically evaluating five existing software packages using the data observing the TRPV1 protein complex. The authors describe similarities and differences in the output of each tool.\nIt is difficult to discern the main conclusions of the manuscript. While the authors discuss which aspects of heterogeneity each technique may capture more effectively, it remains unclear how these findings can guide the broader research community in analyzing other systems because this investigation focuses solely on TRPV1. In this regard, the cited paper 31 [CryoBench], co-authored by some of the present authors, appears to address a similar goal but on a wider scale. On the other hand, if the primary goal of the study is to provide a comprehensive overview of the software, the descriptions of the methods and procedures seem too brief to be of practical use to the reader.\nOther details:\n- There are mistakes in the caption of Figure 5.\n- L613: N_series and N_v not defined.\n- Eq (4): The meaning of [V_i]_k is unclear. Is it the k-th column vector?\n- Eq (5): What’s beta?\n- Section 5.3: In the main text, the term “the common basis” appears. Is it referring to Eq. (5)?\n- L639: What is the meaning of PS_i(k)(k) = PS_i?\n- GROMACS needs citations."}, "author_response": "We thank the reviewer for their detailed comments towards making this a truly useful reference for researchers in the field. We believe these revisions significantly improve the quality of the manuscript. Below we have point-by-point responses to their comments, including substantial additions to the MD simulation and SVD-based comparison portions of the manuscript.\nWe agree with the reviewer that cryo-EM datasets are diverse and present challenges beyond those offered by the analysis of TRPV1 alone. While covering the diversity of heterogeneity in cryo-EM datasets is outside the scope of this work, we have included an analysis of a second dataset that does not present some of the challenges of TRPV1: the P80S ribosome. These results reveal similar themes regarding the strengths and weaknesses of the various heterogeneity analysis methods that we saw for the TRPV1 analysis. A new subsection (Section 2.8) has been added to the end of the results section that includes a new supplemental figure (Supplementary Figure S10) and supplemental movie (Supplementary Movie 18).\nHere are my main comments.\nFor 3DFlex, the author claimed to solve some cytosolic helices to below 3 A. However, there are no clear sidechain features visible on those helices in Figure 2. Some zoomed-in figures are needed to show those are real high resolution features instead of errors in local resolution measurement.\nResponse: We thank the reviewer for this nice suggestion that adds clarity to the 3DFlex results. Additional panels have been added to Figure 2A showing the clear sidechain features in the cytosolic domains for the sharpened volumes produced by 3DFlex. The following sentence has been added to the caption of Figure 2A:\nZoomed panels of the B-factor sharpened volumes show that clear atomistic details have been recovered by 3D Flexible Refinement.\nMap quality in figure 4 is not as good as that in figures 2 and 3, and the map in figure 5 looks even worse. It is unclear whether it is just a difference in filtering, or an actual issue in software. It is recommended to filter the maps to the same (or at least similar) level so they can be compared visually.\nResponse: We thank the review for pointing this out, and we agree that the original versions of these figures were unclear. The maps in Figure 4C were reconstructions from a subset of particles. Hence there are higher levels of noise and dust. The maps from the cryoDRGN VAE in Figure 4A are of a quality comparable to those from 3DVA and 3DFlex. The maps from Figure 4A are the ones subsequently analyzed using AnaVox. We have now updated Figure 4 as follows: 1) volumes have been enlarged in Figure 4A so the detail and the similarity to the volumes in 4C are more obvious, and labels and captions have been updated to make this distinction clearer; 2) the caption in Figure 4A now reads “Volumes produced from the cryoDRGN VAE latent space display the most variation in the occupancy of the DkTx toxin, bound to the extracellular face of TRPV1, with relatively little conformational change in the protein itself.”; and 3) the caption for Figure 4C reads “Selection of particles surrounding the parts of the latent space around these clusters yields consensus reconstructions consistent with the visualizations of the cryoDRGN decoder network.”\nThe low map quality in Figure 5 is likely due to the filtering of particles during the ManifoldEM pipeline and is addressed in more detail in a later comment.\nFor the CryoDRGN results, it is quite strange to observe toxin fading gradually, instead of simply being present/absent. Is the fading density caused by subtle movement of toxin, or the imperfection of the software? As a compositional heterogeneity, it can also be analyzed by traditional focused classification. It would be a good way to quantitatively analyze the latent space distribution.\nResponse: We agree with the reviewer that heterogeneity methods currently have an imperfect treatment of compositional heterogeneity, showing volumes that fade gradually rather than simply being present or absent. Our understanding is that this is a current imperfection in heterogeneity methods – the tendency to produce maps that appear to be interpolations between two states, even if those intermediate states are not present in the data, as is the case for discrete heterogeneity. In general, it is a difficult problem to solve because the noise level in cryo-EM data makes it difficult to classify particles 100% accurately. However, we have added an additional figure to illustrate that the toxin binding is likely discrete. Supplementary Figure S11 has been added depicting the results of traditional 3D classification and text at line 535 onward in the Discussion section has been added to the manuscript:\nAn example of artifacts from these continuous latent spaces can be seen in the gradual fading of the toxin density in the volume series produced by explorations of the CryoDRGN and 3DVA latent spaces. This is despite the fact that binding of the toxin is likely discrete in nature as supported by the results of traditional 3D classification (Supplementary Figure S11).\nWe believe this entanglement of compositional and conformational heterogeneity is a critical point that users of these methods should be aware of, and we have elaborated on this more concretely in the Discussion section (mostly for the benefit of Reviewer 2, but it is also relevant here).\nThe results from ManifoldEM are not very convincing. Does it also produce other conformational changes that are more visible, like the movement of cytosolic helices shown in figure 2? It would be more convincing if the authors showed some large-scale movements that are clearly visible first to build credit for the tool, before showing the less visible small motion.\nResponse: We acknowledge that the volumes currently produced by ManifoldEM are not as convincing as those produced by the other tools in this paper. However, the 2D NLSA movies (Supplementary Movie 16) produced from the ManifoldEM pipeline gives insights into the heterogeneity of the sample. The lower quality of the volumes from ManifoldEM is likely due to the lower number of particles used during analysis because of the method’s intrinsic filtering of particles with low occupancies in parts of the orientation sphere. This is the likely source of the poor map quality alluded to in the reviewers previous comment. To emphasize the tools ability to detect large-scale motions, we have swapped the order of the side and top views in Figure 5. This emphasizes the fact that the same changes to the volumes can be seen in 3DVA and 3DFlex. We have also modified the text to make these points clearer and the paragraphs starting at line 229 now read:\nViewing the 2D NLSA eigenvectors derived by ManifoldEM reveals heterogeneous motions consistent with those found by other tools. When viewing PrDs derived from side views of TRPV1, changes to the conformations of the cytosolic domains are visible (Figure 5A). Whereas the PrDs from top and bottom views of TRPV1 allow one to see twists in the transmembrane domains (Figure 5B). These results are consistent with the conformational changes observed in 3DVA and 3DFlex. Additionally, the focus on projection directions appears to allow ManifoldEM to detect motions not visible to other methods. In top and bottom views of TRPV1, we could also detect changes to the size of the pore of the ion channel. This change to the size of the pore was not visible in the other methods.\nWhile the motions depicted in the 2D NLSA movies appear consistent with other methods, the volumes produced in the final stages of Manifold Embedding are not high enough resolution to capture them. It is likely that the reason for the lower quality volumes produced by this tool is because during preprocessing for this method, any PrDs with fewer than 400 images are discarded. This results in the loss of roughly half of the particle stack from the analysis. By increasing the minimum threshold to 400 as opposed to 200 it was found that more noisy PrDs were excluded leading to an increase in the quality of the reconstructed volumes. However, this thresholding procedure results in the loss of roughly half of the particle stack from the analysis.\nThe MD simulation part sounds surprisingly simple for a membrane protein. The text refers to figure S2, which is not very informative. Are the lipids from the nanodisc considered in the simulation? How are the water molecules included in the simulation, particularly the water inside the channel?\nResponse: We thank the reviewer for pointing this out, indeed the details of our MD simulation were too simplistic in the previous version of the manuscript. We have now added the relevant details to the methods section, including the construction of the lipid bilayer and inclusion of the solvent. Great care was taken to consider the lipids resolved in the original cryoEM density from the nanodisc. Lipids bound to the transmembrane subunits of TRPV1 were built by hand. A new panel A has been added to Figure 6, replacing Supplementary Figure 2, to give a clearer understanding of the molecular dynamics system we studied in this manuscript. The details added to the Methods section on line 685 are as follows:\nThe structure of TRPV1 in this dataset was solved in a nanodisc and four 1-Palmitoyl-2-oleoylphosphatidylcholine, palmitoyloleoylphosphatidylcholine (POPC) lipids are bound to the interior of the transmembrane domains. However, they have missing tails. For our simulations, these tails were built manually in VMD to ensure that these lipids and their polar contacts with the protein were preserved (Humphrey et al. 1996).\nThis protein-ligand system was inserted into a POPC lipid bilayer and solvated using the Charmm-GUI webserver with TIP3P water (Jo et al. 2008). The system was neutralized with sodium and chloride ions at a concentration of 0.15 mol/L. This molecular dynamics system included 78,962 water molecules, 239 chloride ions, 215 sodium ions and 424 POPC lipids. Before running molecular dynamics, the simulation box was roughly cubic, with a side length of 16 nm. The final system contained 382,773 atoms.\n…\nDuring equilibration, water molecules and sodium ions naturally migrated into the pore of the channel.\nThe MD simulation should provide much deeper insight into the dynamics of the system, instead of just being visually similar to the CryoEM analysis. For example, if the movement of the cytosolic helices from different subunits are truly correlated, what is the driving force that leads to the coordinated movement of helices from the opposite side, hundreds of Å apart? Since those helices are not interacting directly with each other, presumably there are residues at the transmembrane domain responsible for the coordination? Is it possible to identify them, maybe even mutate them in silico and show the correlation of cytosolic helical movement disappearing as a result?\nResponse: We thank the reviewer for these comments, and warmly welcomed the opportunity to look into molecular dynamics simulations more deeply. Indeed the analysis we performed showed several interesting insights into the allostery of TRPV1. We have included an analysis of the allosteric network in TRPV1 using the “Correlation of All Rotameric and Dynamical States” (CARDS) method (Singh, Bowman 2017). This method uses mutual information on dihedral angles to determine when conformational transitions are correlated. The results of this analysis are now shown in Figure 6B and Supplementary Figure S6A and described from line 261 onward:\nAllostery is known to be crucial for TRP channels (Hilton et al 2019, Astore et al 2024). Hence, Correlation of All Rotameric and Dynamical States (CARDS) was used to detect which amino acids were responsible for transmitting allosteric information during the 5x2µs simulations run at 310 K (Sing et al 2017). This analysis revealed that the TRP box serves as a pathway for the communication of allosteric information between the cytosolic and transmembrane domains (Figure 6B). This is consistent with findings for other members of the TRPV family (Nadezhdin et al 2021, Burns et al 2023). The CARDS analysis revealed that the region surrounding amino acid K240, at the interface of the cytosolic domain and the C-terminal of the adjacent subunit, is particularly important, suggesting that any correlated motion between neighboring subunits may be mediated by this cluster of residues (Supplementary Figure S6A).\nWe have also calculated the Pearson correlation coefficients for the fluctuations of all four subunits in TRPV1. These results can be seen in a new panel of Supplementary Figure S6 B-C.\nWe agree with the reviewer that mutating residues identified as important to the allostery of TRPV1 would be of considerable interest as there are several mutations known to affect the allosteric properties of TRPV1. However, such studies would be beyond the scope of the present manuscript as we are primarily concerned here with cryo-EM methods. We expect that our simulations and the insights from the CARDS analysis in particular will guide future mutagenesis studies. In keeping with this suggestion, we have chosen to make our simulation systems and trajectories public so that it is easier to perform mutational analysis in future. The following lines have been added to the manuscript on line 276:\nWe note that when the principal components of the MD simulations are visualized, the cytosolic domains appear to move in unison, but more careful study of the raw trajectories reveals that in fact the motions of these subunits are either independent or weakly anti-correlated (Supplementary Figure S6B-C).\nAdditionally, in table 1 authors claim the MD simulation shows the pore closing movement, which was not shown in figure 6. What is the scale of the movement? It has previously been reported the pore opening is related to simulation temperature. Is it also the case here? What is the evidence that the pore closing movement from MD simulation corresponds to the movement observed from ManifoldEM? From my understanding, the gating of most channels are controlled by sidechains and the change is only a few Å in radius, but the maps in figure 5 lack high resolution features. Or did the authors discover a pore closing movement at a much larger scale through their MD simulation?\nResponse: We agree with the reviewer, this section of the manuscript was indeed confusing. We hope the revisions below add clarity to our results. The reviewer is correct that changes to the pore size are difficult to see in the volumes produced by ManifoldEM, but these motions are clearly visible in Supplementary Movie 16 of the NLSA eigenvector. Notes have been added to Table 1 to clarify that some motions were only visible in the 2D NLSA movies and not the volumes.\nRegarding the MD simulations, we observe the TRPV1 pore collapse during the early stages of the trajectory, which is not uncommon for channels with hydrophobic gates of this architecture. The following has been added starting at line 286:\nOne notable aspect of our simulations is that we observed the collapse of the hydrophobic lower gate near the beginning of each simulation (Supplementary Figure S5). However, this is a well-known phenomenon in MD simulations of channels with the six transmembrane helix (6TM) architecture (Mo et al. 2012, Shafrir et al. 2008, Jensen et al. 2010, Anishkin et al. 2004) and the other characteristics of the ligand-bound open state were preserved during the simulation.\nFor the SVD based comparison, there is a lack of a baseline to indicate the significance of the correlation. For example, what would be the correlation values if one generates many maps by grouping particles randomly, or just adding noise to the maps? Otherwise, it is impossible to know whether a correlation score of 0.06 means a significant agreement between two methods.\nResponse: We would like to thank the reviewer for this suggestion in particular, as it made us think more carefully about how to interpret and improve the PCV metric (see details below). We agree that it is important to have a baseline to compare the values of the captured variance. As suggested by the reviewers, we generated a series of volumes with increasing noise levels and calculated the captured variance with the series with no added noise (Sup. Fig. 14). We find that a value of log(PCV)=-2 is equivalent to the PCV found when comparing a volume series to itself with a signal to noise ratio between 10-1 and 10-2. Since this is the same value obtained when comparing series coming from two orthogonal bases in volume space (3DVA Mode 1 vs. 3DVA Mode 2), we define this value as an empirical baseline for determining if the correlation between two methods is better than that obtained by orthogonal methods. We acknowledge that this is an initial baseline test for the PCV metric, and that more comprehensive and robust tests, such as p-value analysis, could be suggested. We will defer these to future work.\nWe would like to note that while addressing this comment we found an error in the code that calculates PCV – we were not squaring Equation 3 in our code. This has now been fixed and all related figures have been updated with the correct results. Note also that these results are now reported as base 10 logarithms. This error did not change the qualitative trends of our results. Additionally, we have chosen to rebrand cryoSVD as “AnaVox” (short for Analysis of voxel maps), to communicate the broader scope of this Python package, and improved the PCV metric definition as described below.\nThe following has been added to the manuscript at line 366:\nInterpreting summary statistics like PCV is nontrivial. For voxel maps in high dimensional spaces, it is more likely for randomly chosen unit vectors to be orthogonal rather than parallel and have a dot product close to zero. Hence, to help interpret the absolute magnitude of the PCV metric, we have designed simple experiments to learn how PCV changes in different situations. By adding white noise to a series of volumes, we can use PCV to compare how the eigenvectors change as a function of the noise level (Supplementary Figure S14A). We find that a log(PCV) of -2 would correspond to a signal to noise ratio (SNR) of between 10−1 and 10−2. This PCV score is the same as what we obtained when comparing modes 1 and 2 of 3DVA, a method that gives us orthogonal components in volume space (Figure 8C, Supplementary Figure S9), indicating that there is no significant correlation between two volume series with log(PCV)< −2. For example, in the case of 3DVA mode 1 (Supplementary Figure S14B) at a SNR of 10−1 the outline of the original eigenvector is still visible, however, by 10−2 it is no longer visible.\nIn addition, to make our PCV metric more informative, we have decided to compute this quantity in an asymmetric way. That is, in the previous revision we were computing PCV(A, B) = ½ d(U_A, V_B, S_B) + ½ d(V_B, U_A, S_A), and now we are have changed to simply computing PCV(A, B) = d(U_A, V_B, S_B). The reason being that this asymmetric form of PCV communicates how well the basis of one volume series spans another. Equation 3 now reflects this revision. An example of this asymmetry is given in Supplementary Figure S12 and line 822 now reads:\nNotice that this metric is asymmetric, meaning that instances exist where PCV (A, B) < PCV (B, A) or vice versa. PCV (A, B) is a measure of how well A covers B, and there can be cases where A covers B well but B does not cover A well. An illustrative example of the asymmetric PCV metric is given in Supplementary Figure S12.\nFor SVD comparison, what has the author done to avoid the different filtering in the maps, as well as the SNR caused by particle count? Would the correlation values be different if we simply lowpass filter every map to 5 or 10Å? Would it matter if only a subset of particles were used in the reconstruction?\nResponse: The reviewer poses valid concerns. Filtering of maps and changes in SNR will cause changes to SVD comparisons. Supplementary Figure S14A shows the effects of adding noise to volume series, and Supplementary Figure S14C shows the effects of low pass filtering on the volumes. Although the PCV matrix values exhibit some changes, we observe that some low-pass filtering or adding some levels of noise does not produce significant effects. To clarify this we have added the text at line 379:\nIt is also of interest what the contributions of different frequencies are to the PCV. We found that low pass filtering of the volumes does not greatly change their relative PCVs up to a low-pass filtering of 14 A. This is consistent with the finding that linear bases of voxel maps are most sensitive to low resolution features (Supplementary Figure S14C-D) (Sorzano, Carazo et al. 2021).\nThe consideration of different numbers of particles is beyond the scope of cryoSVD (now AnaVox) itself as it relates to the SNRs of the volumes. We believe this will be a powerful tool for comparing series of volumes, it is agnostic to how those volumes were produced. A case study of the effects of changing particle count can be found in supplementary Figure S13 where we analysed the volume series presented from independent half-sets of the particle images. This will be described in full in a later comment.\nAlso for the SVD based comparison, it seems that only movement trajectories along eigen-vectors are considered? What if one method describes a movement as a straight path and the other as a curve in the latent space?\nResponse: We share the reviewer's interest in exploring trajectories in high dimensional latent spaces. In fact, this is precisely the reason we developed cryoSVD (now AnaVox) — to compare volume series. If one method describes a straight path through a latent space, and another trajectory is curved, AnaVox is still capable of comparing both these volume series. It is agnostic to how those voxel series were produced. Our PCV metric compares the basis spanned by the eigenvectors of these series in voxel space, not the latent space of the method that produced them. Since this paper focuses on comparing the outputs of different heterogeneity methods, it is logical to compare the volumes produced by sampling in a straight line along the eigenvectors of the latent space of the source method, as this represents their typical outputs. However, one can use this pipeline to compare any two volume series to see if they contain the same heterogeneous motion. The following has been added on line 329:\nThis tool can be used to compare any two volume series to determine if they contain the same molecular heterogeneity. However, here we will use it to compare volume series derived from the eigenvectors from the latent spaces of different methods.\nTo further validate the comparison method, maybe break the datasets to two halves and analyze the heterogeneity in each set separately. This will show whether the correlation is consistent and what is the noise level in the comparison. It still won’t overcome any potential systematic bias but should be a step forward.\nResponse: This was a very useful suggestion by the reviewer and we learned a lot attempting to address this point. To further validate our comparison methods, we split the dataset into two halves and ran heterogeneity analysis on each set separately. We found that without downsampling the particle stacks our results showed poor correlation, indicating a low signal to noise ratio. We reran our analysis on independent half-sets of the data with cryoDRGN and 3DVA. These results show that the volume series from the different half-sets correlate well with each other. However, we did find some notable differences in their power spectrums which appear to be from different low frequency noise between the two half-sets. Supplementary Figure S13 has been added with these results and the following paragraph has been added to line 395 of the manuscript:\nTo gain an understanding of how robust heterogeneity analysis is to particle count, we applied 3DVA and cryoDRGN to independent half-sets of the particle images. Each half-set contained 109,393 particles. At first we found that the results between the half-sets differed significantly. By down sampling the images and applying symmetry expansion, the results between the half-sets became more similar. With this preprocessing, we used AnaVox to demonstrate that each method delineates the same molecular motions in each half-set. The volume series for the same method in different half-sets are visually very similar and they follow similar paths when projected into a common basis (Supplementary Figure S13A). This confirms that they capture the same molecular motions in each half-set. Calculating the PCV for the volume series from different half-sets, we can see that the same method from different half-sets shares a high PCV, confirming that they capture the same molecular heterogeneity (Supplementary Figure S13B). However, the PCV scores indicate that the bases for different half-sets are not the same when comparing across different methods, indicating some differences between them. For example, the first 3DVA mode in the first half-set is more similar to the cryoDRGN results than the first 3DVA mode in the second half-set. Visual analysis of the volume series and inspection of their power series indicate that these differences are most likely due to low-frequency noise in the region surrounding the nanodisc (Supplementary Figure S13C).\nThe nanodisc subtraction strategy might be problematic. Typically for focused alignment, we remove the density of rigid parts by subtracting the projections of averaged structures from the particle images. In this case, we already know the nanodisc is highly variable, so the average structure does not reflect the actual nanodisc density in each particle, resulting in leftover density of different sizes and shapes. Even though it does not have a major impact on homogeneous reconstructions, they may be interpreted as structural variability, and might be the main cause of the heterogeneity at the transmembrane domain reported by CryoDRGN and ManifoldEM.\nResponse: We share the reviewer’s concerns about proper treatment of nanodisc heterogeneity. Signal subtraction, indeed, is not an optimal solution, as it uses a static mask to treat the heterogeneous signal of the nanodisc. Thus, in the long run we believe that more rigorous forms of signal subtraction are necessary for heterogeneous ensembles. This influenced our addition to the AnaVox results section on line 413:\nThis highlights the need for more robust methods to filter structured noise from a dataset during heterogeneity analysis.\nTo ensure that, as the reviewer is concerned, we have not biased our analysis of protein conformational heterogeneity by the process of nanodisc signal subtraction, we have now conducted both the 3DVA and cryoDRGN analyses with and without signal subtraction. We found the same heterogeneity in the transmembrane domains in the volume series produced both with and without signal subtraction. The following has been added to the first Methods Section, elaborating on the nanodisc subtraction procedure on line 596:\nFor completeness, we conducted 3DVA and cryoDRGN on both the signal subtracted and non-signal subtracted data. When viewing the volumes at high thresholds and carefully exploring the latent spaces of each method, we found the same heterogeneity within the transmembrane region of the protein, indicating that the signal-subtraction procedure is likely not biasing our results (Supplementary Movie 17).\nFinally, the newly included heterogeneity analysis of the ribosome shows how the same qualitative trends between methods are reproduced even without signal subtraction.\nReviewer #2 (Remarks to the Author):\nThis study provides a survey of tools designed to extract information on conformational heterogeneity from cryo-EM data sets, specifically evaluating five existing software packages using the data observing the TRPV1 protein complex. The authors describe similarities and differences in the output of each tool.\nIt is difficult to discern the main conclusions of the manuscript. While the authors discuss which aspects of heterogeneity each technique may capture more effectively, it remains unclear how these findings can guide the broader research community in analyzing other systems because this investigation focuses solely on TRPV1.\nResponse: We thank the reviewer for their valuable feedback. We believe this work is useful as a guide to the broader research community in analyzing other systems, but also hope our findings highlight the current differences in the results of existing methods for analyzing conformational heterogeneity in cryo-EM. The reviewer’s suggestions have motivated us to 1) add a supplemental analysis of a ribosome dataset and 2) clarify the goals and findings of our manuscript, emphasizing what is novel and important about our work.\n1) In the additional analysis of the P80S ribosome we see similar trends to our TRPV1 analysis -- see new subsection 2.8, including a new supplemental figure (Figure S10) and supplemental movie.\n2) This work represents the most comprehensive study comparing the heterogeneity of a single experimental dataset with multiple methods to date. With this one dataset we are able to address how different algorithms confound compositional and conformational heterogeneity and the complicating issue of nanodisc heterogeneity for membrane proteins. While this paper starts with qualitative comparisons of different approaches to analyze heterogeneity, we conclude by presenting various metrics for assessing heterogeneity in cryo-EM data, namely comparisons with molecular dynamics simulations through Bayesian ensemble reweighting or by using the AnaVox tool newly presented here. Overall, while we are encouraged by all the new approaches to heterogeneity in cryo-EM, we want to emphasize to readers that heterogeneity in cryo-EM remains a hard problem, and methods still need to be pushed forward by comparing and critiquing methods as possible. We have added clarifying text to this effect in the Abstract, Introduction, and Discussion & Conclusion sections:\n(Introduction line 49): The increasing number of heterogeneity tools in the literature motivated us to perform the most rigorous study of heterogeneity in a biologically relevant sample to date and determine if these tools provide complementary or orthogonal information. … In this vein, we also developed a software framework called AnaVox to quantitatively determine whether two heterogeneity methods agree with each other.\n(Introduction line 100): Hence, we developed the software package AnaVox to quantitatively determine whether two volume series capture the same heterogeneity signal and also characterise the presence of compositional heterogeneity.\n(Discussion lines 453 onward): Overall we have found that for both the TRPV1 and the ribosome datasets we analyzed, different heterogeneity analysis methods more preferentially analyze compositional vs. conformational heterogeneity when both are present in a sample.\nThis study was motivated by the promise of computational tools to access heterogeneous molecular states in single particle cryo-EM data. However, our results demonstrate that while existing software tools are able to reveal the presence of this heterogeneity, they are not yet capable of enabling quantitative biophysics. This is apparent from the differing sensitivity and stochasticity of each tool. For example, the nature of the VAE architecture in cryoDRGN appears to make it especially sensitive to compositional heterogeneity. On the other hand, the volume preserving flow field of 3DFlex limits it to the detection of conformational heterogeneity. Outside these extremes, we found these tools tend to convolve the two kinds of heterogeneity. For this dataset, all of the methods we used produce clear evidence of a spectrum of states. However, the cryoDRGN result clearly indicates a population of channels in the apo closed state, but in combination with the results of the other methods, we also believe there is significant heterogeneity within the ligand bound states. Characterising the independent heterogeneous motion of the two states appears to be beyond the capability of the current generation of tools. This work establishes the strengths and weaknesses of different approaches to the study of heterogeneity. However, continued software development in this space is still needed before these tools can be used for quantitative biophysical investigations of macromolecules by cryo-EM practitioners.\n(Discussion lines 556 onward): Overall, this study of molecular heterogeneity in a real dataset serves as a proof of concept for more ambitious experiments that can exploit the heterogeneity of vitrified samples in cryo-EM. With the analysis of these datasets we have demonstrated the AnaVox pipeline as well as discussed issues with the entanglement of compositional and conformational heterogeneity, which we believe is a key issue for methods developers to tackle in the future. We have even suggested an initial solution in this direction, by presenting the mean power spectrum variance for each method.\nIn this regard, the cited paper 31 [CryoBench], co-authored by some of the present authors, appears to address a similar goal but on a wider scale.\nResponse: Although the motivations of CryoBench and this work are similar, there is a key difference in that the datasets analyzed in this work are experimentally collected datasets. CryoBench presents synthetic datasets as benchmarks for heterogeneity analysis in cryo-EM, while here we address the performance of a variety of heterogeneity methods on real datasets, with all the complications, nuances, and uncontrolled variables that implies, where no ground truth is available. In this work we consider practical considerations relevant to these analyses that cannot be captured from the analysis of synthetic data. Both the approach of CryoBench and the approach presented here are necessary and meaningful. To clarify this point, the following has been added to the Introduction on line 52:\nWhile a similar task has been performed using ground truth volumes with synthetic data (Jeon et al 2024), the study presented here provides comparison metrics and validation tools when using real data for which no ground truth is available.\nOn the other hand, if the primary goal of the study is to provide a comprehensive overview of the software, the descriptions of the methods and procedures seem too brief to be of practical use to the reader.\nResponse: Indeed the primary goal of this study is not to provide a comprehensive overview of the software, but to rather be a guide to the reader of our experiences of using these tools to tackle one experimental cryo-EM dataset with different methods and how we have gone about comparing the results. We hope that the additions to the reviewer’s first comment help to clarify this. We refer the reader to the original papers and documentation for the practical application of each method. We have also expanded the Methods sections where we have gone above and beyond the standard use of each method to help inform the readers on how we have found the specific results from the tools presented in this study.\nOther details:\n- There are mistakes in the caption of Figure 5.\n- L613: N_series and N_v not defined.\n- Eq (4): The meaning of [V_i]_k is unclear. Is it the k-th column vector?\n- Eq (5): What’s beta?\n- Section 5.3: In the main text, the term “the common basis” appears. Is it referring to Eq. (5)?\n- L639: What is the meaning of PS_i(k)(k) = PS_i?\n- GROMACS needs citations.\nResponse: These minor errors have been fixed and we thank the reviewer for catching them:\n- The caption for Figure 5C was missing and has been added.\n- N and N have now been defined, and they are now at line 812.\nseries v\n- [V] of equation 4 is now defined in line 833.\ni k\n- Equations 5 and 6 have now been altered for clarity. Indeed the definition of beta was originally unclear, and we have reformulated the equations to make the mathematics clearer. Originally, beta was simply the coefficients of the expansion in the common basis. We now write out these coefficients explicitly, showing how the original basis vectors from the individual volume series are projected into the new common basis.\n- The term ‘common basis’ is introduced in Section 2.7 and we now reference the correct methods section (now Section 4.4) for clarity, and have generally improved the description of what the common basis is within Section 4.4.\n- Indeed, the definition of the normalization of the power spectra was not well-defined, and should have been \\hat{\\mathrm{PS}(k)_i} = \\mathrm{PS}(k)_i/\\mathrm{PS}(0)_i this is now fixed at line 855.\n- A GROMACS citation (Van Der Spoel et al) has now been added to section 4.3."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The manuscript is greatly improved after revision. The CryoEM maps are now better presented, statistical analysis is more meaningful, and the detailed explanation of the MD simulation part helps bridge CryoEM heterogeneity analysis to the simulation. I am relatively happy with the technological aspects of the paper at its current version. My remaining concern is that in some cases, the paper might be over-claiming biological findings without sufficient evidence. Since this is mainly a technological paper, I think this can be solved by just down-weighting the biological part, which won’t really damage the impact of the paper.\n1. The claimed compositional heterogeneity, i.e. present vs faded toxin density, is still not convincing. In figure S8, the authors try to show the uniform angle distribution of the toxin faded class. However, the angle distribution of the entire dataset is not included. If the full dataset has preferred orientation and one class has uniform orientation distribution, it might indicate some interesting issue. Either binding to toxin changes molecule behavior at air water interface and induces preferred orientation, or the class without toxin contains more bad particles that get more random orientation assignment. In Figure S11, the authors try to reproduce the two classes by traditional classification methods. However, with similar particle numbers, protein features in the toxin faded classes are much worse. So this is hardly a supporting evidence of the existence of compositional heterogeneity. It also might be good to correlate the classes in S11 with the latent space in Figure 4B. It would be a good comparison between deep learning and traditional methods.\nWith all the evidence present in the paper, it is hard to argue there is a population without toxin rather than the toxin just binds somewhat flexibly. Since none of the methods reveals a concrete class without toxin, it might be better to not claim with certainty that this is a compositional heterogeneity.\n2. I would not interpret the change in movie 16 to be the “changes to the size of the pore of the ion channel”. It is a very strong biological conclusion without much evidence. Without well defined 3D volumes, it is hard to argue this is not just random fluctuation of leftover lipid from subtraction. It is better to not draw the conclusion without sufficient resolution structures.\n3. The anti-correlation between opposite cytosolic domains from MD simulation is interesting. Since now there is movement information from both observed data and from simulation, is it possible to do some statistical test to see whether they have the same distribution? It would be an exciting finding if we are observing statistically significantly more or less correlated motion in CryoEM than that is predicted by the simulation.\n4. Figure S9 shows an interesting comparison between the different software but it is very hard to distinguish the marks from the sample plot. It is better to make separate subplots.\n5. In Figure S13, can you also include the comparison of half sets with initial full particles that shows the differences? I think uncertainty of latent space construction by the machine learning methods is a useful thing to point out. Also I noticed that while the scale of motion by different methods vary a lot in Figure S9, the difference is much smaller in Figure S13. What is the cause of this? Some plotting issue or has to do with the downsampling?\n(Remarks on code availability)\nIt would be good if you can add more documentation in the notebook.", "reviewer_2": "With the modifications to the text, the motivation for this study is now clearer. However, the general conclusion remains rather simplistic: some methods are suitable for detecting compositional heterogeneity, while others are better for conformational heterogeneity.\nRegarding the related study, [15] CryoBench, I am curious whether the conclusions presented here align with those of CryoBench, which analysed synthetic data.\nBased on the revised manuscript, I have several comments concerning the data from the MD simulations.\nThe discussion in Section 2.5 regarding allostery is brief and difficult to follow. A brief explanation of what the CARDS analysis does would be helpful. Simply stating that K240 is “important” does not provide enough insight.\nAdditionally, volumes are generated using structures along PCA vectors. How was the maximum amplitude of displacement along these PCA vectors determined? The motions shown in the accompanying movie appear to be quite large. Did the structure actually undergo such significant movements during the relatively simple MD simulation used?\nIn relation to this, the description indicates that only the PCA vectors were used to compare motions with heterogeneity analysis from EM data. While the “direction” of motion from PCA may correspond to some of the EM data, what about the amplitude? Do the structural fluctuations observed in MD match the extent of heterogeneity indicated by the EM analysis?\nFurthermore, for the reweighting analysis, how can it be confirmed that the conformations captured in the cryo-EM images were indeed sampled during the MD simulations?\nL276: “We note that when the principal components of the MD simulations are visualized, the cytosolic domains appear to move in unison, but more careful study of the raw trajectories reveals that in fact the motions of these subunits are either independent or weakly anti-correlated (Supplementary Figure S6B-C).” How is this apparent discrepancy resolved? I suspect that the procedure used to “symmetrize” the MD trajectories may be responsible for this effect.\nL320: “Indicating that the sample’s ensemble is cooler than the temperature at which it was incubated before vitrification (4°C).” This statement is unclear. I assume this implies that the conformational ensemble is narrow, hence being described as “cooler.” However, this interpretation seems inconsistent with the EM data analysis results, which indicate large conformational heterogeneity.\nL746: What is meant by “a standardization transform”?"}, "author_response": "We would like to thank the reviewer for their positive feedback, rigor, and further suggestions to improve the manuscript. We agree that the main point of this paper is the technical insights on cryo-EM analysis rather than the biological insights on TRPV1 and have down-weighted our biological claims accordingly:\nAll references to compositional heterogeneity of DkTx have been changed to specifically refer to changes in the volume series themselves as opposed to implying that the sample itself exhibits compositional heterogeneity.\n● Figure 1 has been revised to no longer suggest that DkTx detaches from TRPV1 but merely document that we observe the toxin fading using some heterogeneity methods.\n● Section 2.2 is now titled “3DVA detects TMD twisting, cytosolic domain lengthening and variable density in DkTx.”\n● Section 2.3 is now titled “The CryoDRGN latent space encodes changes in DkTx density.”\n● Figure 4 C is now labeled “toxin present” and “toxin faded” rather than “toxin bound” and “toxin detached”\n● Table 1 has had the column labeled “toxin occupancy” changed to “Change in toxin density”.\nAll references to the pore opening and closing have been removed.\n● Figure 1 has had the depiction of changes to the pore removed.\n● Figure 3C now refers to “tall” and “short” conformations rather than “open” and “closed” ones.\n● The panels in figure 5 are now labeled “tall” and “short” as opposed to “closed” and “open”.\n● “Pore opening” in Table 1 has been changed to “Change in pore region”.\n1. The claimed compositional heterogeneity, i.e. present vs faded toxin density, is still not convincing. In figure S8, the authors try to show the uniform angle distribution of the toxin faded class. However, the angle distribution of the entire dataset is not included. If the full dataset has preferred orientation and one class has uniform orientation distribution, it might indicate some interesting issue. Either binding to toxin changes molecule behavior at air water interface and induces preferred orientation, or the class without toxin contains more bad particles that get more random orientation assignment. The reviewer’s inquiries into the compositional heterogeneity of DkTx in this dataset have led to several interesting insights and improvements to this manuscript. Indeed multiple phenomena could explain the class-dependent preferred orientation observed in Figure S9 (previously S8). With regards to this specific comment, we have added the full dataset to the first column of this figure as requested. These results show that the full dataset has preferred orientation similar but not identical to the toxin bound state, with many more side views in the full dataset. We believe it is unlikely that the toxin unbound state contains more bad particles than the other classes as its resolution with C1 symmetry is 3.5 Å which is comparable to that of the toxin bound state, 3.6 Å, while both classes have similar numbers of particles. These classes have only slightly worse resolution from the full dataset, down from 3.3 Å. These resolutions have also been added to Figure S9. Even though we find it unlikely that the toxin unbound class contains a large number of poor quality particles, we have remarked on another possibility (see below), which may be giving rise to the appearance of compositional heterogeneity in this dataset.\nThe following has been added to line 571:\nOn close inspection of the cryoDRGN-based reconstructions of the toxin-bound and unbound states, we found that the toxin-bound particles have a non-uniform distribution of viewing angles as compared to the unbound particles, which have a near uniform distribution (Supplementary Figure S9). This indicates that the air-water interface may be influencing the compositional state of the particles. However, an alternative explanation for the fading density of the toxin to true compositional heterogeneity could be the entanglement of compositional state with the pose distribution of the particles in the cryoDRGN latent space. This explanation is consistent with our results as particles that contain side views, in which the toxin is more visible, are over-represented (Supplementary Figure S9B). This suggests that what appears as compositional heterogeneity may in fact be an artifact of entangling pose with the cryoDRGN latent space. This entanglement has been observed with this method when applied to synthetic data (Jeon et al. 2024).\nIn Figure S11, the authors try to reproduce the two classes by traditional classification methods. However, with similar particle numbers, protein features in the toxin faded classes are much worse. So this is hardly a supporting evidence of the existence of compositional heterogeneity. It also might be good to correlate the classes in S11 with the latent space in Figure 4B. It would be a good comparison between deep learning and traditional methods.\nWe acknowledge that the protein features in the toxin faded classes are much worse than those of the toxin-bound classes despite similar particle numbers between both classes. By contrast, the class of toxin-faded structures found by cryoDRGN are remarkably high resolution for only 45 000 particles. Figure S12 (previously S11) now has the location of the classes found by traditional classification methods projected into the latent space from Figure 4B. Adding this was a nice suggestion as it provides an apples to apples comparison between deep learning and traditional methods. Line 605 has been changed such that it does not suggest evidence of compositional heterogeneity:\nThis is despite the fact that binding of the toxin is likely to be discrete in nature if this fading reflects true compositional heterogeneity. We attempted to use traditional 3D classification to reconstruct any apo states, but we found that classes without visible toxin density were significantly worse quality than those with clear density associated with DkTx (Supplementary Figure S12).\nWith all the evidence present in the paper, it is hard to argue there is a population without toxin rather than the toxin just binds somewhat flexibly. Since none of the methods reveals a concrete class without toxin, it might be better to not claim with certainty that this is a compositional heterogeneity.\nWe greatly appreciate the reviewer’s attention to detail regarding the apparent compositional heterogeneity of the toxin. This point in particular has led us to think very carefully about how to study compositional heterogeneity in experimental datasets. This line of inquiry has prompted us to perform what we have found to be an insightful computational experiment:\nWe produced a synthetic stack of particles with noise levels thought to approach experimental levels (a signal to noise ratio of 10^-2). The particle stack only contained a single structure, but its orientations came from the poses estimated from the reference experimental dataset. When heterogeneity analysis was applied to this synthetic particle stack, we found significant levels of hallucinated compositional heterogeneity in the toxin and other parts of TRPV1. This raises the possibility that the compositional heterogeneity in this dataset may in fact be hallucinated. The results of this experiment are now presented in the manuscript in Supplementary Movie 20 and Supplementary Figure S16 and all mentions of compositional heterogeneity in the manuscript have been appropriately edited to reflect this new possibility.\nIn the 3DVA results (Section 2.2) we have added the following on line 186\nHowever, we cannot be sure that this is true compositional heterogeneity, as experiments with synthetic data demonstrate that under some conditions, a non-uniform pose distribution in a dataset may cause 3DVA to detect compositional heterogeneity where none exists (Supplementary Figure S16, Supplementary Movie 20).\nIn the cryoDRGN results (Section 2.3), we have added the the following on line 219\n“However, whether these results represent true compositional heterogeneity or a hallucination due to the dataset’s non-uniform pose distribution is unclear (Supplementary Figure S9, S16).”\nFollowing on from the previous additions the following has been added to the discussion on line 586:\nTo explicitly test this possibility, we performed an experiment where we produced a synthetic dataset with only a single structure, but with the same distribution of poses as that found in the full dataset (Supplementary Figure S16 A). When we applied 3D Variability Analysis to this particle stack, we found significant amounts of compositional heterogeneity in the first mode (Supplementary Movie 20, Supplementary Figure S16 B). This raises the possibility that the presence of compositional heterogeneity we detected in this dataset may in fact be an artefact of its non-uniform pose distribution rather than a truly variable occupancy of DkTx.\n2. I would not interpret the change in movie 16 to be the “changes to the size of the pore of the ion channel”. It is a very strong biological conclusion without much evidence. Without well defined 3D volumes, it is hard to argue this is not just random fluctuation of leftover lipid from subtraction. It is better to not draw the conclusion without sufficient resolution structures.\nIn keeping with focusing the manuscript on technological rather than biological claims, mentions of the pore closing have been removed. However, the appearance of intricate changes to the ion channel in well populated PrDs of ManifoldEM is a possible unique advantage of using individual 2D projection directions to study heterogeneity. Hence, we have kept this observation in the manuscript but made the language more focused on the method rather than the biological findings. Figure 5 has had its text and caption updated to remove mention of the pore and line 249 and onward reads:\nFor example, the NLSA eigenvectors from the top and bottom views of TRPV1 show changes in intensity within the pore. However, while the larger motions depicted in the 2D NLSA movies appear consistent with other methods, the volumes produced in the final stages of ManifoldEM are not high enough resolution to capture changes within the pore, meaning we cannot discount the possibility that these changes may be an artefact leftover from signal subtraction.\n3. The anti-correlation between opposite cytosolic domains from MD simulation is interesting. Since now there is movement information from both observed data and from simulation, is it possible to do some statistical test to see whether they have the same distribution? It would be an exciting finding if we are observing statistically significantly more or less correlated motion in CryoEM than that is predicted by the simulation.\nWe share the reviewers excitement surrounding probing the agreement between the MD simulation and observed data from cryo-EM. A statistical test that measures whether the same distribution is extracted by both techniques would be extremely useful. The most direct method to compare the motions from the cryoEM datasets and the MD simulations would be to build atomic models from voxel maps sampled from the heterogeneity method’s latent space and then calculate the Pearson’s correlation between the subunits of these new models, and compare these statistics with the MD trajectories. However, building atomic models can be very time consuming and would be difficult in our case due to the low resolution nature of many of the volumes from some heterogeneity methods. We have added language to this effect in line 625:\nAs an example of a problem [AnaVox] is not yet able to help tackle, we noticed the anti-correlation between opposite cytosolic domains in our MD simulations, facilitated by the ability to accurately measure distances in atomic space. We would ideally be able to compare directly to the results from the cryo-EM analyses, however, this is currently difficult because of the time consuming task of constructing atomic models and the uncertainty associated with building models into lower resolution voxel maps. Recent work toward easily converting between atomic models and voxel maps may help make this process routine (Jamali et al. 2024, Chen. 2025).\n4. Figure S9 shows an interesting comparison between the different software but it is very hard to distinguish the marks from the sample plot. It is better to make separate subplots.\nSupplementary Figure 10 (previously S9) has been split into nine panels to highlight the similarities and differences between the paths traced by each volume series in the common linear basis. Splitting the plot in this way makes it easier to distinguish the datapoints from different software. We greatly appreciate this suggestion to improve the aesthetic and legibility of this figure.\n5. In Figure S13, can you also include the comparison of half sets with initial full particles that shows the differences? I think uncertainty of latent space construction by the machine learning methods is a useful thing to point out. Also I noticed that while the scale of motion by different methods vary a lot in Figure S9, the difference is much smaller in Figure S13. What is the cause of this? Some plotting issue or has to do with the downsampling?\nThis suggestion by the reviewer makes it easier to understand the uncertainty in the latent space construction by different methods. As requested, we have included the comparison of half sets with the volume series derived from the initial full particle stack which were depicted in Supplementary Figure 10 (previously S9). This makes this analysis easier to relate to the results from Supplementary Figure 10 and Supplementary Figure 14 (previously S13). Note that although we have downsampled the volumes from Figure S10 in order to compare them to the downsampled volumes from the half-sets in Figure S14, the range of the volumes from the full dataset is the same as the non-downsampled projections in Figure S10. This confirms that there are no plotting issues. Rather, the differences in the scale of motion are due to some methods having inherently less variance than others, and in Supplementary Figure 14 we have chosen to do our half-set analysis with methods which have similar ranges of motion to each other.\nReviewer #1 (Remarks on code availability):\nIt would be good if you can add more documentation in the notebook.\nDocumentation and technical explanations have been added to the notebook.\nReviewer #2 (Remarks to the Author):\nWith the modifications to the text, the motivation for this study is now clearer. However, the general conclusion remains rather simplistic: some methods are suitable for detecting compositional heterogeneity, while others are better for conformational heterogeneity.\nRegarding the related study, [15] CryoBench, I am curious whether the conclusions presented here align with those of CryoBench, which analysed synthetic data. Based on the revised manuscript, I have several comments concerning the data from the MD simulations.\nWe thank the reviewer for continuing to provide detailed feedback and assisting us with adding clarity to our manuscript with regards to how it compares to the current state of the literature. The presented study is an in-depth technical overview of the current generation of heterogeneity methods applied to real data. Beyond observing that some methods are better for detecting compositional or conformational heterogeneity, the presented manuscript gives a practitioner of cryoEM an insight into how they might use and compare the tools available in the literature and be aware of their potential issues.\nThe abstract has now been reworded to better reflect the limitations of the current generation of tools, while also relating our enthusiasm and optimism that biophysical studies will soon become routine with vitrified cryoEM samples.\nHowever, each method we tested produced different results, with some methods detecting either only compositional or conformational heterogeneity in the sample. In order to compare the diverse results of these methods, we developed a software framework, AnaVox, to quantitatively determine whether two heterogeneity methods agree with each other.\nAdditionally, we believe that our work has several notable aspects that make it a unique contribution which complements the conclusions of cryoBench. This is because we have different audiences for our studies. This study is presented for those who wish to study real data with existing methods-while cryoBench was developed for those who wish to create new methods and want to prototype those methods on simple datasets where ground truth is available.\nCryoBench only considers idealized synthetic datasets with uniform samples of poses, and datasets have been designed such that they do not mix compositional and conformational heterogeneity. Such clean data sets can rarely be captured experimentally. Most datasets will exhibit a non-uniform pose distribution due to interactions with the air-water interface, and membrane proteins in particular will exhibit both compositional and conformational heterogeneity due to the presence of a dynamic lipid membrane. This is now remarked on in the introduction in line 53:\nWhile recent work motivated by similar goals compared results from heterogeneity methods when applied to synthetic data (Jeon et al. 2024), the study presented here provides comparison metrics and validation tools when studying real data. Data collected in experiments typically exhibits significantly greater complications than idealised synthetic data created for benchmarks, such as a mixture of both compositional and conformational heterogeneity, a non-uniform distribution of poses and unknown ground truth structures.\nConsidering real datasets has encouraged us to create more practical metrics for studying heterogeneity of voxel maps than what is found in cryoBench. The metrics developed in cryoBench were intended to be applied to synthetic datasets, where ground truth embeddings in a latent space are known ahead of time. However, through the AnaVox pipeline we have developed a technique which allows one to compare a series of voxel maps without needing to know the ground truth they were drawn from. This will serve as a practical tool for users in studying real data, as opposed to only developers trying to create new methods. We have remarked on this in the conclusion on line 641:\nNotably, we have developed new metrics for studying heterogeneity methods without reference to ground truth data or latent embeddings, which was a significant limitation on previous iterations of metrics designed to assess the output of heterogeneity workflows (Jeon et al. 2024).\nThe discussion in Section 2.5 regarding allostery is brief and difficult to follow. A brief explanation of what the CARDS analysis does would be helpful. Simply stating that K240 is “important” does not provide enough insight.\nWe believe that the allosteric analysis was a valuable addition to the manuscript in the last round of revisions at the request of the reviewers and so we are pleased to add additional details. We have added further insights into our analysis and a short description of the CARDS method. Line 276 in Section 2.5 now reads as follows:\nHere we employed Correlation of All Rotameric and Dynamical States (CARDS) for this purpose. This program uses dihedral angles to classify the conformational state of amino acids as being in either a disordered or ordered state. Throughout the simulation, the categories of amino acids are updated, and transitions of the amino acids between these states are tracked, allowing the software to determine which groups of amino acids exhibit the kind of correlated motions we expect from allosteric networks (Singh, Bowman 2017).\n…\nIn addition to the allosteric communication between the cytosolic domains and the transmembrane domains, CARDS also revealed allosteric coupling between adjacent cytosolic domains. This is primarily mediated by amino acids at the interface between the C-terminus and the ankyrin repeat domains in the N-terminus. Close inspection of the ankyrin repeat domain suggests that residue K240 plays a central role in this inter-subunit allosteric communication (Supplementary Figure 6B).\nAdditionally, volumes are generated using structures along PCA vectors. How was the maximum amplitude of displacement along these PCA vectors determined? The motions shown in the accompanying movie appear to be quite large. Did the structure actually undergo such significant movements during the relatively simple MD simulation used? In relation to this, the description indicates that only the PCA vectors were used to compare motions with heterogeneity analysis from EM data. While the “direction” of motion from PCA may correspond to some of the EM data, what about the amplitude? Do the structural fluctuations observed in MD match the extent of heterogeneity indicated by the EM analysis?\nThese details are crucial to understanding the agreement between the molecular dynamics simulations and the cryoEM results, hence we are happy to elaborate.\nThe maximum amplitude of displacement along the PCA vectors was determined by multiplying unit eigenvectors from the Principal Component Analysis and scaling them by their respective singular values, to help visualize the magnitude of the different motions. To help visualize the extent of this motion, we have projected these scaled eigenvectors into the linear atomic basis in Supplementary Figure 4A and indicated their location on the figure. This confirms that our simulations undergo significant movements in our unbiased MD simulations.\nA new panel B has also been added to Supplementary Figure 4. Here we produced voxel maps from frames in the MD simulations and projected them into the common embedding derived in Section 2.7. This demonstrates that the structural fluctuations observed in MD match the extent of the heterogeneity indicated by some cryoEM analysis methods. However, we also note that methods that capture the toxin fading (ManifoldEM, cryoDRGN and 3DVA) tend to span a larger region of the common embedding than methods that only capture conformational heterogeneity (molecular dynamics and 3DFlex). This observation has been added on line 404:\nThese trends are also apparent from inspecting the path of all series in the common basis, where it is clear the series that capture the fading of the toxin generally span a larger region of the common basis than those that describe only conformational heterogeneity (Supplementary Figure S10).\nThe following details regarding the relative scales of the visualisation of the eigenvectors have also been added to the methods section on line 847:\nTo visualize these eigenvectors, we selected the first two columns of the matrix V^T, which contain the unit eigenvectors for principal components 1 and 2, and multiplied them by their respective singular values. This allows us to visualize the relative scale of variance captured by a given eigenvector as well as its direction of motion. To produce an interpretable sample of these eigenvectors in Cartesian space, we created several atomic models by multiplying the scaled eigenvectors by a scalar parameter between -1 and 1 at 50 equally spaced intervals. This allowed us to view motion along the eigenvectors as a continuous deformation in the interpretable Cartesian basis (Supplementary Movie 9). The 50 samples of these scaled eigenvectors can also be projected back into the principal component basis and plotted (Supplementary Figure S4A).\nFurthermore, for the reweighting analysis, how can it be confirmed that the conformations captured in the cryo-EM images were indeed sampled during the MD simulations?\nWe share the reviewer's concerns related to confirming that the conformations captured in the cryo-EM images are sampled by the MD simulations with respect to our ensemble reweighting analysis. In fact this is an active area of development for the ensemble reweighting method and hence addressing this point completely would be out of place for this manuscript. Here we merely aim to show that certain structures are upweighted, indicating that there may be some particles in the full dataset corresponding to structures more similar to the apo state. To reflect this, the paragraph on line 344 now reads:\nHowever, a nuanced aspect of the ensemble reweighting method in cryo-EM is the choice of which specific structures to use to represent the full MD ensemble. Ideally, one requires representative structures of the entire conformational space, which are then reweighted by the data. However, due to computational limitations, it is often difficult to determine whether MD simulations have adequately sampled all relevant conformational states. To enhance the likelihood of sampling structures that correspond to conformations observed in cryo-EM data, we performed simulations of TRPV1 at three different temperatures: 150 K, 200 K, and 310 K. We note that improving conformational sampling remains an important direction for future applications of ensemble reweighting, with a potential strategy involving the direct biasing of MD simulations using information derived from cryo-EM images as has been done with other forms of experimental data (Marinelli, Fraldo-Gómez. 2015, White et al. 2015).\nL276: “We note that when the principal components of the MD simulations are visualized, the cytosolic domains appear to move in unison, but more careful study of the raw trajectories reveals that in fact the motions of these subunits are either independent or weakly anti-correlated (Supplementary Figure S6B-C).” How is this apparent discrepancy resolved? I suspect that the procedure used to “symmetrize” the MD trajectories may be responsible for this effect.\nWe thank the reviewer for picking up on this subtle point. We believe the anti-correlated motion between opposing subunits arises from higher order principal components, not to analyzed in detail. These asymmetric modes account for a significant portion of the variance in the trajectory, with PCs 3 and 4 accounting for a combined 14% of the variance in the dataset. A new Supplementary Movie 5 has been uploaded to visualize these eigenvectors and they are also depicted in a new Supplementary Figure, Figure S7. We have also added a panel to this figure to depict the spectrum of eigenvalues and the amount of variance they each capture. For clarity, we have added a new panel D to Supplementary Figure 6 to explicitly outline the symmetry expansion procedure and present the results of the non-symmetry expanded trajectories. These results have been commented on in line 302 of the manuscript:\nThis anti-correlation may arise from the third and fourth principal components of the trajectory (Supplementary Figure S7A-B, Supplementary Movie 5), which according to the size of their eigenvalues, each explain 7% of the variance in the trajectory which is still significant compared to the first and second components which account for 20% and 9% of the variance, respectively (Supplementary Figure S7C).\nL320: “Indicating that the sample’s ensemble is cooler than the temperature at which it was incubated before vitrification (4°C).” This statement is unclear. I assume this implies that the conformational ensemble is narrow, hence being described as “cooler.” However, this interpretation seems inconsistent with the EM data analysis results, which indicate large conformational heterogeneity.\nWe agree with the reviewer that this statement was speculative. We have reworded it to be more general and better reflect our results. Line 355 now reads:\nWhen structures from simulations at all three temperatures are used in reweighting, the structures derived from lower temperature simulations received higher weights (Supplementary Figure S2B-C). According to Supplementary Figure S3, these structures deviated much less from the starting configuration where TRPV1 is closer to an open state compared to the simulations at 300K. This indicates that the majority of particles are likely to be in a configuration closer to the cryo-EM derived consensus structure, while a portion are in configurations closer to the apo-closed state. These results are consistent with the presence of the conformational states explored by 3DFlex and 3DVA, but now projected onto the interpretable MD latent space.\nL746: What is meant by “a standardization transform”?\n“Standardization” refers to the procedure of subtracting the mean from a series and dividing each value in the series by its standard deviation. This can also be called “normalization”. We have chosen to change to the latter term so that its meaning is clearer. Line 826 now reads:\nHence, to help ignore these disordered fluctuations, before performing principal component analysis, we normalized the data by subtracting the mean position from each atom and divided their position in each frame by the standard deviation over the course of the trajectory. This normalization procedure emphasizes large collective motions."}, {"round": 3, "reviewer_comments": {"reviewer_1": "I am happy with the revised article and response to reviewers, and I recommend the publication of this paper.\n(Remarks on code availability)", "reviewer_2": "Certain aspects of the PCA analysis conducted on the MD trajectories are still unclear.\n-Using “a version of principal component analysis” (L299), the authors observed “the cytosolic domains appear to move in unison” (L303).\n-“but more careful study of the raw trajectories reveals that in fact the motions of these subunits are either independent or weakly anti-correlated” (L304)\n-PC’1’ shown in Figure 6C looks “anti-correlated”.\n-“This anti-correlation may arise from the ‘third and fourth’ principal components of the trajectory.” (L305)\nThese points appear to contradict. 1) What do the authors mean by “move in unison” and “anti-correlated”? 2) Is the PC1 in Figure 6C derived from standard PCA or symmetrized PCA? 3) The symmetrized PCA does not describe the motions observed in the MD trajectory itself. The authors state that this approach was used to reveal the deformations “similar to those seen with 3DVA and 3DFlex” (L300), but is this approach justifiable?\nConventional PCA is typically used to identify large-scale motions. I understand that this approach may not be optimal for the authors’ analysis due to the presence of disordered domains. However, the described “normalization” procedure is not entirely clear to me. If all motions are normalized in the way described, it seems that the distinction between large-scale and small-scale motions would be lost. The manuscript states that “This normalization procedure emphasizes large collective motions” (L832), but since all atomic motions are normalized, it is unclear what types of motions are actually being captured.\nCould the authors clarify this point? Additionally, has this normalization strategy been employed in previous studies, and if so, could references be provided?\n(Remarks on code availability)"}, "author_response": "We appreciate the reviewer’s rigor and their assistance with making the results of our MD simulations clearer. We have endeavored to clarify any points of confusion by making changes to the manuscript where appropriate but we stress that the technical choice of normalizing the Cartesian coordinates prior to performing PCA is a minor one and we do not believe it changes anything significant about our results.\n-Using “a version of principal component analysis” (L299), the authors observed “the cytosolic domains appear to move in unison” (L303). See below, as we have rewritten the paragraph at L298 to clarify this and other points.\n-“but more careful study of the raw trajectories reveals that in fact the motions of these subunits are either independent or weakly anti-correlated” (L304) See below, as we have rewritten the paragraph at L298 to clarify this and other points.\n-PC’1’ shown in Figure 6C looks “anti-correlated”.\nIn Figure 6C we are measuring the distance from the protein’s center of mass for each subunit. To help emphasize this point, we have added legends to panels C and D. With this added legend, it is clearer that chains A and B move together, and chains C and D move together. When we use the term correlated or anti-correlated, we are specifically referring to the motions of opposite subunits. To add clarity to this point, line 301 now reads:\n“In the first two principal components of the MD simulations, opposing cytosolic domains exhibit correlated motion—when one cytosolic domain moves inward so does the opposite one”\nWhile line 311 now reads:\n“For example, calculating the distance between the center of mass of the protein and the cytosolic domains reveal that opposite subunits in fact exhibit weakly anti-correlated motion (Supplementary Figure 6B-C). That is, when one domain moves inward, the opposite one tends to move outward.”\n-“This anti-correlation may arise from the ‘third and fourth’ principal components of the trajectory.” (L305)\nThese points appear to contradict. 1) What do the authors mean by “move in unison” and “anti-correlated”?\nWe apologize for the confusion in terminology and appreciate the reviewer’s help to make the results of our analysis clearer. We have rewritten the paragraph on line 298 in order to more clearly emphasize the results from the MD simulations. We note that although counter-intuitive, there is no contradiction between variables “moving in unison” and being “anti-correlated”. By referring to the motions as “in unison” and “anti-correlated”, we meant to indicate that when one subunit moves inward, the other one moves outward. If they did not move in unison, the appropriate term for this kind of motion would be uncorrelated.\nTo alleviate any confusion, we have created a GitHub repository with all the code used to create Figure 6 and removed the phrase “move in unison” from the relevant paragraph on line 298. This paragraph now reads:\n“To compare the conformational space of our MD simulations to the results of the cryo-EM heterogeneity analysis methods, we applied PCA to our MD simulations. This revealed deformations to the protein similar to those seen with 3DVA and 3DFlex (Supplementary Movies 9-10). In the first two principal components of the MD simulations, opposing cytosolic domains exhibit correlated motion—when one cytosolic domain moves inward so does the opposite one. This is the same behavior implied by the results of the cryo-EM methods. In PC1, we see an asymmetric motion in the cytosolic domains like those found in the first mode of 3DVA and the second in 3DFlex. While in PC2, even with DkTx bound, we see the TMDs twist and the cytosolic domains lengthen, as with the second mode of 3DVA and the first mode of 3DFlex. This is consistent with an oscillation between the toxin-bound open state and a configuration similar to the apo closed state (Figure 6D, Supplementary Movie 9). However, due to the atomistic nature of MD simulations, we have access to much more detailed motions compared to existing cryo-EM analysis methods. For example, calculating the distance between the center of mass of the protein and the cytosolic domains reveal that opposite subunits in fact exhibit weakly anti-correlated motion (Supplementary Figure S6B-C). That is, when one domain moves inward, the opposite one tends to move outward. As PCA is a dimensionality reduction technique, analysis of individual principal components loses some of the information from the raw trajectories. The deficiencies of PCA in particular are well-known in the MD simulation community [36–38]. Hence, we inspected higher order PCs to test whether we could find any anti-correlated eigenvectors. Indeed, PC 3 and 4 of the MD simulations exhibit anti-correlated motions between opposing cytosolic domains (Supplementary Figure S7A-B, Supplementary Movie 5). According to the size of their eigenvalues, these latter two components each explain 7% of the variance in the trajectory which is significant compared to the first and second components which account for 20% and 9% of the variance respectively (Supplementary Figure S7C). These higher order components resolve the apparent discrepancy between analysis of the raw MD trajectories and the inspection of the first two PCs, highlighting how the use of dimensionality reduction may obscure important details in the study of molecular motions.”\n2) Is the PC1 in Figure 6C derived from standard PCA or symmetrized PCA?\nThe eigenvectors described in Figure 6C and D are derived from a symmetrized PCA. This is now indicated in the figure legend. We have also added a new paragraph to the methods section detailing how we have performed these calculations depicted in this figure and created a GitHub repository where one can download and run the scripts we used to analyze the motion of the cytosolic domains: https://github.com/flatironinstitute/astore_natcomms_2025_scripts\nLine 882 now reads:\nAs an interpretable measurement of the conformational changes to TRPV1, we calculated the vertical and horizontal displacement of the cytosolic domains in both the molecular dynamics trajectories and during projection of the eigenvectors of the PCA (Figure 6C-D). To track the horizontal displacement of these domains, the radial displacement in the x-y plane between two the centers of mass of residues in the cytosolic domains were compared to a reference point defined as the center of mass for the TMDs. The relevant residues in the cytosolic domains were 114-183 while the TMDs were defined as intervals on residues 418–454, 469–498, 511–532, 536–554, 563–598, 630–642 and 656–712 in order to exclude disordered loops. Similarly, to track the vertical displacement of the cytosolic domains, we calculated the difference in the z coordinate between the centers of mass in the cytosolic domains and the TMDs, with each domain defined with the same residue as before. Scripts for performing these calculations are available on GitHub:\nhttps://github.com/flatironinstitute/astore_natcomms_2025_scripts.\n3) The symmetrized PCA does not describe the motions observed in the MD trajectory itself. The authors state that this approach was used to reveal the deformations “similar to those seen with 3DVA and 3DFlex” (L300), but is this approach justifiable?\nAs mentioned in the new paragraph on line 298, while the MD trajectories contain more information than the PCA analyses, analyzing the eigenvectors of the trajectory allows us to make consistent comparisons between the voxel-based methods and our molecular simulations. Conventional PCA is typically used to identify large-scale motions. I understand that this approach may not be optimal for the authors’ analysis due to the presence of disordered domains. However, the described “normalization” procedure is not entirely clear to me. If all motions are normalized in the way described, it seems that the distinction between large-scale and small-scale motions would be lost. The manuscript states that “This normalization procedure emphasizes large collective motions” (L832), but since all atomic motions are normalized, it is unclear what types of motions are actually being captured. Could the authors clarify this point?\nThe normalization procedure employed in this manuscript is common practice in some areas of data science where it is also referred to as “standardization”. Importantly, we stress that this normalization procedure does not produce significantly different results compared to the usual PCA procedure. Our results would not qualitatively change if we had not employed it. For example, the latent embeddings obtained with and without normalization have the same geometry (see plots below).\nSimilarly, visualizations of the eigenvectors derived from both procedures capture the same conformational changes. There are only certain loops of the protein that move marginally more in the non-normalized eigenvectors, while the same large-scale conformational changes are captured by the first four eigenvectors.\nThese results have been rendered in Supplementary Movie 21 and they have been added to the manuscript on line 845:\nIn our case, a comparison between eigenvectors derived from normalized and unnormalized coordinates reveals that motions in segments of the protein with no secondary structure like the C-terminus and residues F720-W740 are deemphasized by this normalization procedure. However, the overall difference is minor and the same conformational changes are found for corresponding eigenvectors (Supplementary Movie 21).\nAdditionally, has this normalization strategy been employed in previous studies, and if so, could references be provided?\nThis normalization strategy was recently employed in Mancl et al. eLife 2025. Though here the authors did not describe its effects in detail (see the code in Supplementary File 5 of their manuscript).\nWe have also added two references to the manuscript specific to the normalization procedure used in this manuscript (Jolliffe and Cadima 2016, and Brereton 2025) on line 845. The former reference in particular details the effects of this preprocessing on the results of PCA. In summary, it ensures that isolated variables that change significantly in scale do not dominate the analysis. Instead, this preprocessing tends to emphasize collective changes between several variables.This is desirable in our case as we are more interested in collective domain motions rather than large-scale changes of poorly-ordered loops.\nThis is outlined in the text on line 844:\nWe found that this normalization procedure emphasizes collective motions, meaning that isolated motions, such as those found in disordered loops, will not bias the analysis (Jollife and Cadima 2016, Brereton 2025). This same procedure has been used in another recent comparison study between cryo-EM and MD simulations (Mancl et al. 2025)."}]

Decision letter