cryoSPHERE: Single-Particle HEterogeneous REconstruction from cryo EM

## Abstract

Abstract The three-dimensional structure of a protein plays a key role in determining its function. Methods like AlphaFold have revolutionized protein structure prediction based only on the amino-acid sequence. However, proteins often appear in multiple different conformations, and it is highly relevant to resolve the full conformational distribution.
Single-particle cryo-electron microscopy (cryo EM) is a powerful tool for capturing a large number of images of a given protein, frequently in different conformations (referred to as particles ). The images are, however, very noisy projections of the protein, and traditional methods for cryo EM reconstruction are limited to recovering a single, or a few, conformations.
In this paper, we introduce cryoSPHERE, a deep learning method that takes as input a nominal protein structure, e.g. from AlphaFold, learns how to divide it into segments, and how to move these as approximately rigid bodies to fit the different conformations present in the cryo EM dataset. This formulation is shown to provide enough constraints to recover meaningful reconstructions of single protein structures. This is illustrated in three examples where we show consistent improvements over the current state-of-the-art for heterogeneous reconstruction.

## 1 Introduction

Single-particle cryo-electron microscopy (cryo EM) is a powerful technique for determining the three-dimensional structure of biological macromolecules, including proteins. In a cryo EM experiment, millions of copies of the same protein are first frozen in a thin layer of vitreous ice and then imaged using an electron microscope. This yields a micrograph: a noisy image containing 2D projections of individual proteins. The protein projections are then located on this micrograph and cut out so that an experiment typically yields $10^{4}$ to $10^{7}$ images of size $N_{\textrm{pix}}\times N_{\textrm{pix}}$ of individual proteins, referred to as particles. Our goal is to reconstruct the possible structures of the proteins given these images.

Frequently, proteins are conformationally heterogeneous and each copy represents a different structure. Conventionally, this information has been discarded, and all of the sampled structures were assumed to be in only one or a few conformations (*homogeneous* reconstruction). Here, we would like to recover all of the structures in a *heterogeneous* reconstruction.

Structure reconstruction from cryo EM presents a number of challenges. Firstly, each image shows a particle in a different, unknown orientation. Secondly, because of the way the electrons interact with the protein, the spectrum of the images is flipped and reduced. Mathematically, this corresponds to a convolution of each individual image with the Point Spread Function (PSF). Thirdly, the images typically have a very low signal-to-noise ratio (SNR), usually around $0.1$, see e.g. Baxter et al. ().

For all these reasons, it is very challenging to perform de novo cryo EM reconstruction. Standard methods, produce electron densities averaged over many, if not all conformations (Scheres, ; Punjani et al., ), performing discrete heterogeneous reconstruction. More recent methods attempt to extract continuous conformational heterogeneity, e.g. by imposing constraints on the problem through an underlying structure deformed to fit the different conformations present in the dataset, see e.g Rosenbaum et al. (); Zhong et al. (); Li et al. (). AlphaFold (Jumper et al., ) and RosettaFold (Baek et al., ) can provide such a structure based on the primary sequence of the protein only. In spite of this strong prior, it is still difficult to recover meaningful conformations. The amount of noise and the fact that we observe only 2D projections creates local minima that are difficult to escape (Zhong et al., ; Rosenbaum et al., ), leading to unrealistic conformations.

To remedy this, we root our method in the observation that different conformations can often be explained by large scale movements of domains of the protein (Mardt et al., ). Specifically, we develop a variational auto-encoder (VAE) (Kingma & Welling, ) that, from a nominal structure and a set of cryo EM images:

- •
Learns how to divide the amino-acid chain into segments, given a user defined maximum number of segments; see Figure 2.
The nominal structure can for instance be obtained by AlphaFold (Jumper et al., ).
- •
For each image, learns approximately rigid transformations of the identified segments of the nominal structure, which effectively allows us to recover different conformations on an image-by-image (single particle) basis.

These two steps happen concurrently, and the model is end-to-end differentiable.
The model is illustrated in Figure 1.

Note that what we call a segment is conceptually different from a domain in the structural biology sense. The domains of a protein play a pivotal role in diverse functions, engaging in interactions with other proteins, DNA/RNA, or ligand, while also serving as catalytic sites that contribute significantly to the overall functionality of the protein, see e.g. Schulz & Schirmer (); Nelson et al. (). By comparison, the segments we learn do not necessarily have a biological function. However, while not strictly necessary for the function of the method, experiments in Section 4 show that our VAE often recovered the actual domains corresponding to different conformations.

The implementation of the model is available on github (^1^11[h](h)ttps://github.com/Gabriel-Ducrocq/VAECryoEM).

Figure: Figure 1: Flow chart of our network. The learnable parts of the model are the encoder, the decoder and the Gaussian mixture. Note that even though the transformations predicted by the decoder are on a per image basis, that is not the case of the Gaussian mixture, which is shared across all particles.
Refer to caption: x1.png

## 2 Background

### 2.1 Notations and problem formulation

In what follows, $S=\{r_{i}\}_{i=1,\dots,R_{\mathrm{res}}}$ denotes a protein made of a number $R_{\mathrm{res}}\in\mathbb{N}^{\star}$ of residues $r_{i}$. Each residue is itself a set of atoms: $r_{i}=\{a_{k}\}_{k=1,\dots,A_{i}}$ where $A_{i}\in\mathbb{N}^{\star}$ is the number of atoms present in residue $i$. We can equivalently define the structure $S=\{a_{k}\}_{k=1,\dots,A_{\mathrm{total}}}$ in terms of its atoms, where $A_{\mathrm{total}}\in\mathbb{N}^{\star}$ is the total number of atoms in the structure. We denote by $a_{k}=(x_{k},y_{k},z_{k},\eta_{k})^{T}\in\mathbb{R}^{3}\times\{C_{\alpha},C_{
\beta},N\}$ the coordinates and atom type of $a_{k}$ for $k\in\{1,\dots,A_{\mathrm{total}}\}$.
This means that in this paper, we limit ourselves to the backbone of the protein: we consider only the $C_{\alpha},C_{\beta}$ and $N$ atoms.

The electron density map of a structure $S$, also called a volume, is a function $V_{S}:\mathbb{R}^{3}\rightarrow\mathbb{R}$, where $V_{S}(x)$ represents the probability density function of an electron of $S$ being present, evaluated at $x\in\mathbb{R}^{3}$. That is, if $B\subseteq\mathbb{R}^{3}$, the probability that the structure $S$ has an electron in $B$ is given by $\int_{B}V_{S}(x)\mathrm{d}x$.

Assume we have a set of 2D images $\{I_{i}\}_{i=1}^{N}$ of size $N_{\textrm{pix}}\times N_{\textrm{pix}}$, representing 2D projections of different copies of the same protein in different conformations. Traditionally, the goal of cryo EM heterogeneous reconstruction has been to recover, for each image $i$, the electron density map $V_{i}$ corresponding to the underlying conformation present in image $i$. See Subsection 2.2 for a review these methods.
However, following recent works, e.g., Rosenbaum et al. (); Zhong et al. (), we aim at recovering, for each image $i$, the underlying structure $S_{i}$ explaining the image. That is, we try to recover the precise position in $\mathbb{R}^{3}$ of each residue.

Since the different conformations of large proteins can often be explained by large scale movements of its domains (Mardt et al., ), our method learns to decompose the amino-acid chain of the protein into segments and, for each image $I_{i}$, to rigidly move the learnt segments of a base structure $S_{0}$ to match the conformation present in that image, as explained in Section 1.

AlphaFold (Jumper et al., ) and RosettaFold
(Baek et al., ) can provide such a structure $S_{0}$, based on the amino-acid sequence of the protein at hand. In Section 4, we further fit the AlphaFold predicted structure into a volume recovered by a custom backprojection algorithm provided by Zhong et al. ().

Figure: Figure 2: Example of segments recovered with a Gaussian mixture of 6 components.
Refer to caption: extracted/5629471/figures/experiments/phy_10k/mask_10k.png

### 2.2 Related work

Two of the most popular methods for cryo EM reconstruction, which are not based on deep learning, are
RELION (Scheres, ) and cryoSPARC (Punjani et al., ). Both methods perform volume reconstruction, hypothesize that $k$ conformations are present in the dataset and perform maximum a posteriori estimation over the $k$ density maps, thus performing discrete heterogeneous reconstruction. Both of these algorithms operate in Fourier space using an expectation-maximization algorithm (Dempster et al., ) and are non-amortized: the poses are learnt for each image, and adding new images to the dataset would require a new run of the algorithm. Other approaches perform continuous heterogeneous reconstruction. For example, 3DVA (Punjani & Fleet, ) uses a probabilistic principal component analysis model to learn a latent space while E2GMM (Chen & Ludtke, ) parameterizes the volume as a mixture of $N$ Gaussian p.d.f and learns their mean conditionally on an input image.

Another class of methods involve deep learning and typically performs continuous heterogeneous reconstruction using a VAE architecture. Some of them attempt to reconstruct a density map, while others target the atomic structures – or a backbone – directly. For example, cryoDRGN (Zhong et al., , ) and CryoAI (Levy et al., ) use a VAE acting on Fourier space to learn a latent space and a mapping that associates a 3D density map with each latent variable. They perform non-amortized and amortized inference over the poses, respectively. Other methods are defined in the image space, e.g. 3DFlex (Punjani & Fleet, ) and cryoPoseNet (Nashed et al., ). They both perform non-amortized inference over the poses. These methods either learn, for a given image $I_{i}$, $\{V_{i}(x_{k})\}$ the values at a set of $N_{\textrm{pix}}^{3}$ fixed 3D coordinates $\{x_{k}\}$, representing the volume on a grid (*explicit* parameterization), or an actual function $\hat{V}_{i}:\mathbb{R}^{3}\rightarrow\mathbb{R}$ in the form of a neural network that can be queried at chosen coordinates (*implicit* parameterization). The implicit parameterization requires smaller networks but needs an expensive $N_{\textrm{pix}}^{2}$ and $N_{\textrm{pix}}^{3}$ evaluation for image and volume generation, respectively.

Figure: Figure 3: Toy dataset. Leftmost: Histograms of the predicted and true angles of rotation in radians. The true angles are in green. The recovered distances are in blue. Left: Predicted against true angles in Ångström. The black line represent $x=y$. Middle: Base structure with the two domains, also used to generate the images. The domain in blue is rotated according to the axis $(0,1,0)$ with angles of rotations sampled from the green distribution on the leftmost figure. Note that the segments predicted by cryoSPHERE exactly match the two domains. The fourth segment is in blue and corresponds to the last 157 residues of chain B, matching exactly the ground truth domain. Right: First mode structure. Rightmost: Second mode structure.
Refer to caption: extracted/5629471/figures/experiments/toy_data_set/distribution_distance.png

Other deep learning methods attempt to directly reconstruct structures instead of volumes and share a common process: starting from a plausible base structure, obtained with e.g. AlphaFold (Jumper et al., ), for each image, they move each residue of the base structure to fit the conformation present in that specific image. These methods differ on how they parameterize the structure and in the prior they impose on the deformed structure or the motion of the residues. For example AtomVAE (Rosenbaum et al., ) considers only residues and penalizes the distances between two subsequent residues that deviate too much from an expected value. CryoFold (Zhong et al., ) considers the residue centers and their side-chain and also imposes a loss on the distances between subsequent residues and the distances between the residue centers and their side-chain. Unfortunately, due to the high level of noise and the fact that we observe only projections of the structures, these ”per-residue transformation” methods tend to be stuck in local minima, yielding unrealistic conformations unless the base structure is taken from the distribution of conformations present in the images (Zhong et al., ), limiting their applicability to real datasets. Even though AtomVAE (Rosenbaum et al., ) could roughly approximate the distribution of states of the protein, it was not able to recover the conformation given a specific image. In addition, to our knowledge, the code for these two methods is not available and non-trivial to reimplement. This prevents a comparison of our results with AtomVAE and CryoFold. To reduce the bias that the base structure brings, (Schwab et al., ) recently proposed to fit pseudo-atoms in a consensus map with a neural network directly. Similar to our work, several other methods constrain the atomic model to rigid body motions. For example (Chen et al., ), starting from a base structure $S_{\mathrm{big}}$ fitted in the consensus map (called the big GMM), learn a representation structure $S_{\mathrm{small}}$ of this consensus map using a smaller number using the method described in (Schwab et al., ). They subsequently learn a translation on each residue $S_{\mathrm{small}}$, that contribute to translate the residues of $S_{\mathrm{big}}$ weighted by how close they are from the residues of $S_{\mathrm{small}}$. This can be thought of as our segmentation GMM, except that it takes places in the 3D space and only the centers of the Gausian modes can be learnt, while ours is fitted on the residue indexes and all of the GMM parameters are learnable. In addition, they do not use this segmentation to perform rigid body motions. Instead, in a subsequent step, the user is in charge of identifying the rigidly moving segments of the protein by applying one or several masks to the protein. This is in contrast to our method where the segmentation used to perform rigid body motion is learnt by the network trying to move the domains to fit the images. Finally, this method involves a sequence of optimization problem, with potentially different architecture, while ours is end-to-end differentiable. Using the method in (Schwab et al., ), Chen et al. () develop a focused refinement on patches of the GMM representation of the protein. These patches are learnt using k-means on the location of residues and do not depend on the different conformations of the data set. This in contrast to cryoSPHERE which learns the segments of the protein are tightly linked to the change of conformation. Concurrently to our work, Li et al. () developed CryoSTAR which acts on each residue and enforces local rigidity of the motion of the protein by imposing a similarity loss between the base structure and the deformed structure.
The interested reader can see Donnat et al. () for an in-depth review of deep learning methods for cryo EM reconstruction.

The reconstruction methods relying on an atomic model, such as e.g cryoSTAR, DynaMight or cryoSPHERE offer the possibility to the user to provide prior information via this atomic model. They also offer the possibility of deforming the protein according to chemical force fields. This is not the case of the methods performing volume reconstruction without such an atomic model.

## 3 Method – cryoSPHERE

In this section, we present our method for single-particle heterogeneous reconstruction, denoted cryoSPHERE. The method focuses on structure instead of volume reconstruction. It differs from the previous works along this line in the way the movements of the residues are constrained: instead of deforming the base structure on a residue level and then imposing a loss on the reconstructed structure, it moves parts of the protein approximately rigidly by construction.

We use a typical VAE architecture, see Figure 1 for a depiction of how the information flows in our network. We map each image to a latent variable by a stochastic encoder. The latent variable is then passed to a decoder which outputs a transformation (rotation and translation) per segment. Based on these transformations and segment decomposition, the underlying structure $S_{0}$ is deformed and then turned into a volume. Finally, the volume is projected into a 2D image, according to its pose that can be estimated using existing software, e.g. (Scheres, ; Punjani et al., ). Finally, this image is convolved with the PSF
to give the reconstructed image that is compared to the input image. This completes the forward pass.
In the backward pass, the parameters of the encoder, decoder and Gaussian mixture are updated with a gradient descent based on the
evidence lower bound (see Section 3.2).
In the following, we describe details of our network: the encoder and decoder, the segment decomposition, how to turn a structure into a 2D projection of its electron density, and an efficient way to compare the true and predicted images.

### 3.1 Image formation model

To compute the 2D projection of the protein structure $S$, we first estimate its 3D electron density map $V$:

$$ $V_{S}(r):=\sum_{a\in S}\dfrac{1}{\sqrt{2\pi}\sigma}\exp\left({-\dfrac{||r-a||^ {2}}{2\sigma^{2}}}\right)$ (1) $$

where $r\in\mathbb{R}^{3}$ and $\sigma$ is a positive real number. In this definition the protein’s electron density is approximated as the sum of a Gaussian p.d.f centered on each atom. We use a fixed value of $\sigma=1$, following Rosenbaum et al. ().

From these density maps, we then compute an image projection $I\in\mathbb{R}^{N_{\textrm{pix}}\times N_{\textrm{pix}}}$ as:

$$ $I(R,t,S)(r_{x},r_{y})=g*\int_{\mathbb{R}}V_{RS+t}(r)dr_{z},$ (2) $$

where $(r_{x},r_{y})\in\mathbb{R}^{2}$ are the coordinates of a pixel, $r_{z}\in\mathbb{R}$ is the coordinate along the $z$ axis, $R\in SO(3)$ is a rotation matrix and $t\in\mathbb{R}^{3}$ is a translation vector. The abuse of notation $RS+t$ means that every atom of $S$ is rotated according to $R$ and then translated according to $t$. The image is finally convolved with the point spread function (PSF) $g$, which in Fourier space is the contrast transfer function (CTF), see Vulović et al. ().

In practice, we avoid approximating the integral in Equation (2). The density map, being a sum of Gaussian p.d.f, allows for an exact calculation by integrating over the last marginal. This significantly reduces the computing time.

Figure: Figure 4: Spike dataset. Left: Predicted angles. Middle left: FSC curves for the open and close conformations. Middle: Close conformation. Middle right: Open conformation. Right: Definition of the angles.
Refer to caption: extracted/5629471/figures/experiments/spike/spike_distrib.png

### 3.2 Maximum likelihood with variational inference

To learn a distribution of the different conformations, we hypothesize that the conformation seen in image $I_{i}$ depends on a latent variable $z_{i}\in\mathbb{R}^{L}$, with prior $p(z_{i})$.
Let $f_{\theta}(S_{0},z)$ be a function which, for a given base structure $S_{0}$ and latent variable $z$, outputs a new transformed structure $S$. This function depends on a set of learnable parameters $\theta$.
Then, the conditional likelihood of an image $I^{\star}\in\mathbb{R}^{N_{\textrm{pix}}\times N_{\textrm{pix}}}$ with a pose given by a rotation matrix $R$ and a translation vector $t$ is modeled as $p_{\theta}(I^{\star}|R,t,S_{0},z)=\mathcal{N}(I^{\star}|I(R,t,f_{\theta}(S_{0}
,z)),\sigma_{\mathrm{noise}}^{2})$, where
$\sigma_{\mathrm{noise}}^{2}$ is the variance of the observation noise.
The marginal likelihood is thus given by

$$ $p_{\theta}(I^{\star}|R,t,S_{0})=\int p_{\theta}(I^{\star}|R,t,S_{0},z)p(z)dz.$ (3) $$

In practice, the pose $(R,t)$ of a given image is unknown. However, following similar works (Zhong et al., ; Li et al., ), we suppose that we can estimate $R$ and $t$ to sufficient accuracy using off-the-shelf methods (Scheres, ; Punjani et al., ).

Directly maximizing the likelihood (3) is infeasible because one needs to marginalize over the latent variable. For this reason, we adopt the VAE framework, conducting variational inference on $p_{\theta}(z|I^{\star})\propto p_{\theta}(I^{\star}|z)p(z)$, and simultaneously performing maximum likelihood estimation on the parameters $\theta$.

Let $q_{\psi}(z|I^{\star})$ denote an approximate posterior distribution over the latent variables. We can then maximize the evidence lower-bound (ELBO):

$$ $\displaystyle\mathcal{L}(\theta,\psi)=\mathbb{E}_{q_{\psi}}[\log p_{\theta}(I^ {\star}|z)]-\mathrm{D}_{\mathrm{KL}}(q_{\psi}(z|I^{\star})||p(z))$ (4) $$

which lower bounds the log-likelihood $\log p_{\theta}(I^{\star})$. Here $\mathrm{D}_{\mathrm{KL}}$ denotes the Kullback-Leibler (KL) divergence. In this framework $f_{\theta}$ is called the decoder and $q_{\psi}(z|I^{\star})$ the encoder.

### 3.3 Segment decomposition

We fix a maximum number of segments $N_{\textrm{segm}}\in\{1,\dots,R_{\mathrm{res}}\}$ and we represent the decomposition of the protein by a stochastic matrix $G\in\mathbb{R}^{R_{\mathrm{res}}\times N_{\textrm{segm}}}$. The rows of $G$ represent ”how much of each residue belongs to each segment”, and our objective is to ensure that each residue mostly belongs to one segment, that is:

$$ $\begin{split}\forall i\in\{1,\dots,R_{\mathrm{res}}\},&\;\exists m^{\star}\in \{1,\dots,N_{\textrm{segm}}\}\\ \mathrm{s.t}&\sum_{m\neq m^{\star}}G_{im}\ll 1\end{split}$ (5) $$

We also aim for the segments to respect the sequential structure of the amino acid chain, and the model to be end-to-end differentiable. Without end-to-end differentiability, we could not apply the reparameterization trick and we would have to resort on Monte-Carlo estimation of the gradient of the segments, which has a higher variance, see e.g (Mohamed et al., ).

To meet these criteria, we fit a Gaussian mixture model (GMM) with $N_{\textrm{segm}}$ components on the real line supporting the residue indices, that is, on $\mathbb{R}$. Each component $m$ has a mean $\mu_{m}$, standard deviation $\sigma_{m}$ and a softmax weight $\alpha_{m}$. The $\{\alpha_{m}\}$ are passed into a softmax to obtain the weights $\{\pi_{m}\}$ of the GMM, ensuring they are positive and summing to one.
We define the probability that a residue $i$ belongs to segment $m$ as the probability that the index $i$ belongs to mode $m$, annealed by a temperature $\tau\in\mathbb{R}^{+\star}$:

$$ $G_{im}:=\dfrac{\{\phi(i|\mu_{m},\sigma^{2}_{m})\pi_{m})\}^{\tau}}{\sum_{k=1}^{ N_{\textrm{segm}}}\{\phi(i|\mu_{k},\sigma^{2}_{k})\pi_{k}\}^{\tau}}$ (6) $$

where $\phi(x|\mu,\sigma^{2})$ is the unidimensional Gaussian probability density function with mean $\mu$ and variance $\sigma^{2}$ and $\tau$ is a fixed hyperparameter. If $\tau$ is sufficiently large, we can expect condition (5) to be verified. See Figure 2 for an example of a segment decomposition using a Gaussian mixture.

In this ”soft” decomposition of the protein, each residue can belong to more than one segment, allowing for smooth deformations. In addition, the differentiable architecture is amenable to gradient descent methods, and a well chosen $\tau$ can approximate a ”hard” decomposition of the protein. We set $\tau=20$ in the experiment section. In our experience, this segmentation procedure is very robust to different initialization and converge in only a few epochs.

### 3.4 Decoder architecture

The decoder describes the distribution of the images given the latent variables. The latent variables include:

- 1.
One latent variable $z_{i}\in\mathbb{R}^{L}$ per image, parameterizing the conformation.
- 2.
The global parameters $\{\mu_{m},\sigma_{m},\alpha_{m}\}_{m=1,\dots,N_{\textrm{segm}}}$ of the GMM describing the segment decomposition.

Given these latent variables and a base structure $S_{0}$, we parameterize the decoder $f_{\theta}$ in three steps.

First, a neural network with parameters $\theta$ maps $z_{i}\in\mathbb{R}^{L}$ to one 3D translation per segment $\{\vec{t}_{m}\}_{m=1,\dots,N_{\textrm{segm}}}$ and one rotation per segment. The rotation is represented by a unit quaternion $\{\vec{q}_{m}\}_{m=1,\dots,N_{\textrm{segm}}}$ (Vicci, ), which can be equivalently described as a set of axes of rotation $\{\vec{\phi}_{m}\}_{m=1,\dots,N_{\textrm{segm}}}$ and corresponding angles of rotations $\{\delta_{m}\}_{m=1,\dots,N_{\textrm{segm}}}$.

Second, given the parameters of the GMM, we can compute the matrix $G$.
Finally, for each residue $i$ of $S_{0}$, we update the coordinates of all its atoms $\{a_{ik}\}_{k=1,\dots,A_{i}}$:

- 1.
First, $a_{ik}$ is successively rotated around the axis $\vec{\phi}_{m}$ with an angle $G_{im}\delta_{m}$ for $m\in\{1,\dots,N_{\textrm{segm}}\}$ to obtain updated coordinates $a_{ik}^{\prime}$.
- 2.
Second, it is translated according to:
$a_{ik}^{\prime\prime}=a_{ik}^{\prime}+\sum_{j=m}^{N}G_{im}\vec{t}_{m}.$
(7)

This way, the rotation angles and translations for a residue incorporate contributions from all segments, proportionally on how much they belong to the segments.
If condition (5) is met, a roughly rigid motion for each segment can be expected.

### 3.5 Encoder and priors

Figure: Figure 5: MD dataset. Left: Histograms of the distances of the two upper domains. The true distances are in green. The recovered distances are in blue. Right: Predicted against true distances in Ångström. The black line represent $x=y$.
Refer to caption: extracted/5629471/figures/experiments/phy_flexible_20/distribution_distance.png

We follow the classical variational-autoencoder framework. The distribution $q_{\psi}(y|I^{\star})$ is given by a normal distribution $\mathcal{N}(\mu(I^{\star}),\mathrm{diag}(\sigma^{2}(I^{\star})))$ where $\mu\in\mathbb{R}^{L}$ and $\sigma\in\mathbb{R}_{+}^{L}$ are generated by a neural network with parameters $\psi$, taking an image $I^{\star}$ as input.
Additionaly, the approximate posterior distribution on the parameters of the GMM is chosen to be Gaussian and independent of the input image:

$$ $\displaystyle\mu_{m}$ $\displaystyle\sim\mathcal{N}(\nu_{\mu_{m}},\beta^{2}_{\mu_{m}})$ $\displaystyle\sigma_{m}$ $\displaystyle\sim\mathcal{N}(\nu_{\sigma_{m}},\beta^{2}_{\sigma_{m}})$ $\displaystyle\alpha_{m}$ $\displaystyle\sim\mathcal{N}(\nu_{\alpha_{m}},\beta^{2}_{\alpha_{m}})$ $$

where $\{\nu_{\mu_{m}},\beta_{\mu_{m}},\nu_{\sigma_{m}},\beta_{\sigma_{m}},\nu_{
\alpha_{m}},\beta_{\alpha_{m}}\}_{m=1,\dots,N_{\textrm{segm}}}$ are parameters that are directly optimized and in practice we feed $\sigma_{m}$ to an ELU+1 layer in order to avoid a negative or null standard deviation.

Finally, we assign standard Gaussian priors to both the local latent variable $z_{i}$
$z_{i}\sim\mathcal{N}(0,I_{L}),$
and the global GMM parameters $\{\mu_{m},\sigma_{m},\alpha_{m},\}_{m=1,\dots,N_{\textrm{segm}}}$.

Note that applying the reparameterization trick (Kingma & Welling, ) is straightforward for a Gaussian distribution. Calculating the KL-divergence between two Gaussian distributions as in Equation (4), is also straightforward.

## 4 Experiments

In this section, we test cryoSPHERE on different datasets and compare the results to cryoDRGN (Zhong et al., ). CryoDRGN is a state-of-the-art method for continuous heterogeneous reconstruction, in which the refinement occurs at the level of electron densities. To synthesize data for all experiments, we applied CTF corruption to the images and add Gaussian noise to achieve $\mathrm{SNR}\approx 0.1$. See Appendix A for a detailed description of the CTF. Furthermore, we make the assumption that we know the exact poses. However, since we use a structure $S_{0}$, which is different from the structures used to generate the datasets, these poses can only be an approximation. This will not be the case of cryoDRGN, for which these poses will indeed be exact, as this method does not use a base structure. Consequently, in this context, the comparison may introduce a bias in favour of cryoDRGN.

### 4.1 Toy dataset: a continuous motion.

For an initial assessment of cryoSPHERE, we examine a simulated dataset based on
a bacterial phytochrome with UniProt (The UniProt Consortium, ) entry Q9RZA4. This protein is a dimer, with 755 residues per chain. To simulate a continuous heterogeneous set of images, we predicted the structure with AlphaFold multimer (Evans et al., ) and divided the protein into two domains. The first one consists of the full chain A and 598 residues of chain B, 1353 residues in total. The second domain comprises the last 157 residues of the chain B, see Figure 3. We then sample $10^{4}$ rotation angles from a two-components Gaussian mixture. Each conformation is generated by rotating the second domain according to the axis $(0,1,0)$ and one of the sampled rotation angles. This yields $10^{4}$ different conformations. Finally, we uniformly sample $15$ different poses for each structure, yielding a total of $150$k images. Appendix A.1 provides more details, and Figure 3 shows examples of structures.

Testing the segment decomposition, we then run cryoSPHERE by requesting division into $N_{\textrm{segm}}=4$. The program learnt a first and third segment with 0 0 residues, a second segment with $1353$ residues and a fourth segment with $157$ residues (Figure 3). Thus, cryoSPHERE learnt segments according to the ground truth.

Moreover, Figure 3 shows that most of the predicted angles of rotation for the fourth segment are in excellent agreement with the ground truth structural changes. However, approximately $12$ % of the dataset is predicted as belonging to the wrong mode. Figure 8 in the appendix shows that the norms of the predicted translations for both segments are negligible, also indicating that cryoSPHERE has recovered the ground truth. The same figure shows that we recover the right axis of rotation. In Appendix A.1, we plot the predicted angles against the mean of the latent variable distributions output by the encoder. We observe that cryoSPHERE has recovered the circular motion: the greater the latent mean, the lower the predicted angle, with a quick transition from one mode to the other. This indicates that cryoSPHERE is able to recover the structure of single-particles from cryo EM data.

### 4.2 A discrete conformation dataset: the SARS-CoV-2 spike protein

We now test cryoSPHERE on a dataset comprising two discrete conformations of the SARS-CoV-2 spike protein: a close conformation with PDB entry 6ZP5 and an open conformation with PDB entry 6ZP7, see Melero et al. (). We generate $150$k images of size $N_{\textrm{pix}}=300$ with uniformly sampled poses and neglected translations, with seventy percent corresponding to the close conformation and thirty percent to the open one.

Prior to running cryoSPHERE, we obtain a base structure by first running the multimer version of AlphaFold (Evans et al., ) on the amino acid sequence of the spike protein, selecting the best (according to AlphaFold) structures. Second, we utilize a custom backprojection algorithm provided by (Zhong et al., ) to get a homogeneous reconstruction. This yields a volume with low resolution in the mobile parts but decent resolution in the less mobile ones. To reduce a potential bias introduced by using an AlphaFold structure, we subsequently fit this structure into the volume, initially performing rigid-body fitting in ChimeraX (Meng et al., ) and then using ISOLDE (Croll, ).
Subsequently, we run cryoSPHERE with this fitted structure as the base structure for 48 hours on a single GPU, totalling 174 epochs.

We define the two conformations in terms of the angle between two domains, see Figure 4. The same figure demonstrates that we recover the angles within roughly $1$ degree. After clustering the structures predicted by cryoSPHERE based on their angle defined in Figure 4, only $4$ of them get predicted to be in the wrong conformation. Additionally, by turning the structures predicted by cryoSPHERE into volumes using our operator (1), we observe excellent agreement with the ground truth volumes, as shown by the Fourier shell correlation (FSC) curves (Klaholz, ). This illustrates that cryoSPHERE can effectively recover discrete conformational disorder.

### 4.3 Molecular dynamics dataset: bacterial phytochrome.

We now test cryoSPHERE with a more difficult dataset. This time, we simulate a continuous motion of a bacterial phytochrome, with PDB entry 4Q0J (Burgie et al., ). This dimer has 503 residues per chain. We define two upper domains corresponding to the residues 321-502 from both chain A and B. Then, we run a Molecular Dynamics (MD) simulation to simulate a dissociation of these two domains, as depicted in Figure 9 in Appendix A.2. We take $10^{4}$ structures along the trajectory and sample $15$ poses uniformly for each structure. Refer to Appendix A.2 for more details.

Figure: Figure 6: MD dataset. Left: Recovered segments on the fitted AlphaFold structure. The colors denotes different contiguous domains: if two segments share the same color, it does not mean their residues belong to each others. Middle left to middle right: cryoSPHERE predicted structures in green, corresponding ground truth in blue. Right: Fourier shell correlation curves. The mean across structures for cryoSPHERE is in plain blue. Dotted blue represents the mean $\pm$ one standard deviation. The red, plain line represents the mean for cryoDRGN across volumes, with the dotted red line representing the mean $\pm$ one standard deviation.
Refer to caption: extracted/5629471/figures/experiments/phy_flexible_20/mask_20_flexible.png

Analogously to Section 4.2, we obtain an AlphaFold structure that we fit into a volume reconstructed by backprojection. Further details can be found in Appendix A.2. Finally, we run our method with the fitted AlphaFold structure as our base structure and $N_{\textrm{segm}}=20$.

Since cryoSPHERE outputs structures, we can compute the predicted distances between the above-mentioned domains. We plot the true and recovered histograms of the distances between the two domains in Figure 5. Additionally, we plot the predicted distances against the true distances in the same figure. These two plots show that the predicted structures are in excellent agreement with the ground truth structures. This is even true for the most dissociated structures, even though the majority of the structures are in a closed conformation. Figure 6 shows the recovered segments as well as a selection of predicted structures and the corresponding ground truths. Our model has learnt to separate the two mobile top domains from the fix bottom one.

For comparison, we train cryoDRGN for the same duration using the same single GPU. We use the default settings of cryoDRGN, as detailed in Appendix A.2.

We plot the mean of the FSC curves between the predicted volumes and the corresponding ground truth volumes in Figure 6, for both cryoSPHERE and cryoDRGN. This plot shows that we perform better than cryoDRGN at both the $0.5$ and $0.143$ cutoffs. There are two reasons behind this.

First, we fit our base structure into a volume reconstructed from a homogeneous reconstruction method. This aligns the non-moving parts of the protein with the right places, as provided by the dataset. Without this step, it is likely that some medium-scale part of the protein would have been misplaced, potentially reducing the performance of our algorithm at the $0.5$ cutoff; see Appendix A.2 for evidence of this. It is important to note, that this step significantly improves our results, and it comes at almost no extra cost since it relies on the backprojection algorithm.

Secondly, deforming a good AlphaFold structure offers finer details than cryoDRGN. This accounts for the better performance of cryoSPHERE over cryoDRGN at the $0.143$ cutoff. Refer to Appendix A.2 for examples of reconstructed structures and volumes.

Moreover, our method demonstrates robustness to a limited number of images. Even with only one pose per structure (totalling $10^{4}$ images), cryoSPHERE maintains similar performance, while cryoDRGN suffers a significant drop in the FSC curves. Further details can be found in Appendix A.2.

## 5 Discussion

CryoSPHERE presents several advantages compared to other methods for volume and structure reconstruction:

##### Efficiency in Deformation:

Deforming a base structure into a density map avoids the computationally expensive $N_{\textrm{pix}}^{2}$ evaluation required by a decoder neural network in methods implicitly parameterising the grid, such as Zhong et al. (); Levy et al. (). The neural network size is vastly reduced compared to state of the art methods, which explicitly parameterise the density grid of size $N_{\textrm{pix}}^{3}$, see e.g. Punjani & Fleet (). Furthermore, direct deformation of a structure directly avoids the need for subsequent fitting into the recovered density map.

##### Reduced Dimensionality and Noise Resilience:

Learning one rigid transformation per segment, where the number of segments is much smaller than the number of residues reduces the dimensionality of the problem. This results in a smaller neural network size compared to approaches acting on each residues, such as Rosenbaum et al. (). Rigidly moving large portions of the protein corresponds to low-frequency movements, less prone to noise pollution than the high-frequency movements associated with moving each residue independently.
In addition, since our goal is to learn one rotation and one translation per segment, a latent variable of dimension $6\times N_{\textrm{segm}}$ is, in principle, a sufficiently flexible choice to model any transformation of the base structure. Choosing the latent dimension is more difficult for volume reconstruction methods such as (Zhong et al., ).

##### Interpretability:

CryoSPHERE outputs segments along with one rotation and one translation per segment, providing valuable and interpretable information. Practitioners can easily interpret how different parts are moving based on the transformations the network outputs. This interpretability is often challenging for deep learning models such as Zhong et al. (); Rosenbaum et al. ().

Section 4 demonstrates cryoSPHERE’s capability to recover continuous and discrete heterogeneity while performing structure reconstruction. It also demonstrate that even though the choice of $N_{\textrm{segm}}$ impact the FSC to the ground truth, cryoSPHERE is able to recover the correct motion for the entire range of $N_{\textrm{segm}}$ values. It is also able to keep the minimum necessary number of domains when the user sets an unnecessarily high $N_{\textrm{segm}}$, as seen in Section 4.1.

While using a base structure $S_{0}$ provides detailed information missed by cryoDRGN, especially evidenced by our better $0.143$ cutoffs, it can introduce bias on medium scale elements, especially if $N_{\textrm{segm}}$ is too small. We alleviate this problem by fitting $S_{0}$ into a volume reconstructed using a simple backprojection algorithm. This, combined with cryoSPHERE, achieves better $0.5$ cutoffs than cryoDRGN, showcasing its effectiveness in debiasing the results. If such a volume is unavailable, simply increasing $N_{\textrm{segm}}$ can reduce the bias.

Our study opens up for significant advancements in predicting protein ensembles and dynamics, critically important for unraveling the complexity of biological systems. By predicting all-atom structures from cryo EM datasets through more realistic deformations, our work lays the foundation for extracting direct insights into thermodynamic and kinetic properties. In the future, we anticipate the ability to predict rare and high-energy intermediate states, along with their kinetics, a feat beyond the reach of conventional methods such as molecular dynamics simulations.

It would be interesting to assess how much our segmentation correlates with bottom-up segmentation into domains conducted on the “omics” scale, see e.g (Lau et al., ). To achieve this quantitatively, we would need many examples of moving segments from cryo EM investigations to match the millions of segments from the “omics” studies. Therefore, we leave this investigation to later work.

Finally, this work is an important milestone in showing that one can learn a segmentation of the protein that is intimately linked to the change of conformation of the underlying protein, in an end-to-end fashion. This is in contrast to other works that perform the segmentation based on residue distances of the underlying atomic model only, usually not in an end-to-end way. The application of cryoSPHERE to real dataset requires a number of careful choices that we leave to future work.

## Acknowledgement

G.D is financially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) and the Data-Driven Life Science (DDLS) program through the WASP-DDLS collaboration.

This project would not have been possible without the computing resources provided by the Knut and Alice Wallenberg Foundation at the National Supercomputer Centre (Berzelius resource).

We also thank Claudio Mirabello at the National Bioinformatics Infrastructure Sweden at SciLifeLab for providing us access to his AlphaFold installation on Berzelius.

## Appendix A Experiments

In this appendix, we provide more details on the experiments of Section 4.

In all the experiments of Section 4, each posed structure is converted into a volume, which is then projected into 2D images according to our image formation model in (2), with $\sigma=1$. The same CTF corruption is applied with parameters specified in Table 1 and finally noise is added to achieve a SNR equal to 0.1. Here, SNR is defined as the ratio of the variance of the images to the variance of the noise.

**Table 1: Table of the parameters used to CTF corrupt the generated images. The same values were used for all images of all datasets.**
| Parameter | Value |
| --- | --- |
| dfU | 15301.1   Å |
| dfV | 14916.4   Å |
| dfang | 5.28   degrees |
| spherical aberration | 2.7  mm |
| accelerating voltage | 300  keV |
| amplitude contrast ratio | 0.07 |

### A.1 Toy dataset

For this experiment, we predict the phytochrome structure using AlphaFold multimer(Evans et al., ) on its
amino acid sequence with the UniProt (The UniProt Consortium, ) entry Q9RZA4. This protein forms a dimer with 755 residues on each chain. We define two domains for simulation purposes. The first domain comprises the first chain and the first 598 residues of the second chain. The second domain consists of the remaining 157 residues of the second chain. In Figure 3, we present the base structure, the decomposition into the two domains as well as the structures corresponding to the deformed base. These deformations represent the mean rotation of each mode. We rotate the second domain around the $(0,1,0)$ axis, sampling $10^{4}$ rotation angles from the Gaussian mixture:

$$ $0.5\times\mathcal{N}(-\pi/3,0.04)+0.5\times\mathcal{N}(-2\pi/3,0.04)$ (8) $$

This gives $10^{4}$ structures. For each image, we uniformly sample 15 rotation poses and 15 translation poses on $[-10,10]^{2}$. The structure undergoes rotation, translation, and is then turned into an image according to image formation model 2. Subsequently, the images undergo CTF corruption, and noise is added to achieve $\mathrm{SNR}\approx 0.1$. This process generates a total of $150$k images, each of size $N_{\textrm{pix}}=220$.

We run cryoSPHERE with $N_{\textrm{segm}}=4$ for 48 hours on a single NVIDIA A100 GPU, equivalent to 779 epochs.

Due to computational constraints, the plots in Subsection 4.1 are based on only 10000 images, one per conformation.

Figure 7 illustrates the predicted angles against the latent means, demonstrating that the model effectively learns rotational motion.

Additionally, in Figure 8, we plot the distribution of predicted translation norms for the two segments, the predicted angles of rotation for segment 2 and the distribution of the dot product between the predicted axis of rotation for segment 4 and the true axis of rotation $(0,1,0)$. Our model predicts no translation for each segment and no rotation for the second segment. The axis of rotation is almost exactly recovered.

Figure: Figure 7: Toy dataset. Predicted angle against latent mean for cryoSPHERE. Note that for clarity $0.3$ percent of the points were removed.
Refer to caption: extracted/5629471/figures/experiments/toy_data_set/predAngle_vs_latentMean.png

Figure: Figure 8: Toy dataset. Leftmost: Boxplot of the norm of the translations predicted for the fourth segment. Middle left: Boxplot of the norm of the translations predicted for the second segment. Middle right: Boxplot of the predicted angle of rotation for the second segment, in radians. Rightmost: Boxplot of the dot product between the predicted axis of rotation for the fourth segment and the true angle of rotation. CryoSPHERE recovers the right axis of rotation almost perfectly.
Refer to caption: extracted/5629471/figures/experiments/toy_data_set/moving_trans.png

### A.2 Molecular dynamics dataset

We take the structure of a phytochrome with PDB ID 4Q0J (Burgie et al., ) and define two domains: residues 321 to 502 of the first chain and residues 321 to 502 of the second chain. To simulate the dissociation process of the two upper domains, we perform Metadynamics(Barducci et al., ) simulations in GROMACS
(Abraham et al., ; Pronk et al., ) with the PLUMED 2 implementation (Tribello et al., ). The collective variable chosen is the distance
between the self-defined centers of mass (COMs) of the upper domains (residues 321-502 of chain A and
B). A 100 ns simulation is conducted using the NpT ensemble, maintaining pressure control through the
Parrinello-Rahman barostat. Gaussian deposition occurs every 5000 steps, featuring a height of 0.1
kJ/mol and a width of 0.05 nm. Afterwards, we extract $10^{4}$ structures along the dry trajectory. See Figure 9 for examples of structures. From this, two datasets are generated.

Figure: Figure 9: From left to right. 1/ The starting close conformation of the MD simulation. 2/ A medium-open structure arising from the MD simulation. 3/ An open conformation towards the end of the MD simulation. 4/ The AlphaFold structure used as a base structure.
Refer to caption: extracted/5629471/figures/experiments/phy_150k/close.png

The first dataset consists $15$ rotation poses uniformly sampled per structure. Translation poses are neglected in this data set, resulting in a total of $150$k images.

The same process is followed to generate the second dataset, with the exception that only one rotation pose per structure is assumed. This results in $10$k images.

For both datasets, we use $N_{\textrm{pix}}=190$, and each pixel is of size $1$Å. The encoder is a 4-hidden-layer neural network with fully connected hidden layers of dimension of $2048,1024,512,512$. The decoder is a 2-hidden-layer neural network with fully connected hidden layers of dimension $350,350$. We set the batch size to $100$ and use the ADAM optimizer (Kingma & Ba, ) with a learning rate of $0.0003$, unless otherwise stated. The latent dimension is set to $40$ with $N_{\textrm{segm}}=6$ for all runs, unless stated otherwise. For the 150k dataset, we train for 24 hours on a single NVIDIA A100 GPU, corresponding to 642 epochs. The 10k dataset is trained for 12 hours on the same single GPU, corresponding to 4114 epochs.

To compare cryoSPHERE to cryoDRGN, we train cryoDRGN using default settings - a batch size of 8 and a latent dimension equal to 8 - for the same duration, on the same single GPU used for our algorithm. To maintain consistency when comparing volumes generated from cryoDRGN, our methods, and the ground truth, we convert the structures (both ground truth and predicted with cryoSPHERE) into volumes using the same image formation model employed to generate the dataset.

For computational efficiency, the FSC plots and distances plots are not based on all structures. In the case of the 150k dataset, only one image per structure is used to compute the distances. All images in the 10k dataset are utilized. Consequently, the distance plots are based on 10k images in both cases. For the computation of the FSC curves, we select only 1000 images evenly distributed among these 10000 structures.

#### A.2.1 150k dataset

Now, we test the robustness of our method to the choice of a base structure.

First, we perform two runs with $N_{\textrm{segm}}=6$. For the first run, we use the same AlphaFold structure as in section 4.3, before it was fitted into the volume obtained by backprojection. In the second run, we use the same AlphaFold structure but fitted as described in Section 4.3.

Figure 10 displays the results for the unfitted run. Despite cryoSPHERE showing excellent agreement with the ground truth in terms of distances compared to cryoDRGN, it misses some medium scale elements, while maintaining a better level of small scale details. We observe underperformance at the $0.5$ cutoff but outperformance at the $0.143$ cutoff. For example, considering the ground truth volumes of the first MD structure and the corresponding reconstructed volumes from cryoSPHERE and cryoDRGN in Figure 13, we note that when only using an AlphaFold structure without postprocessing, our method can be slightly off on some medium-scale elements, especially at the bottom of the protein. As hypothesizing $N_{\textrm{segm}}=6$ is not sufficient to break these elements in specific segments, cryoSPHERE does not move them, retaining the features of the AlphaFold structure, potentially misplaced. In contrast, cryoDRGN seems to be closer to the ground truth on these elements. However, upon closer inspection of the reconstructed volumes, we observe that cryoSPHERE gives much finer details throughout the volume, a charecteristic not replicated by cryoDRGN.

Figure 11 showcases the results for the fitted run. CryoSPHERE now outperforms cryoDRGN at both the $0.5$ and $0.143$ cutoffs while maintaining an excellent agreement with the ground truth. The fitting of the AlphaFold structure has a significantly positive impact on the performances of our method.

For both runs, we compare the reconstructed volumes of the same image from our method to the ground truth in Figure 14. We observe that throughout the volume, the parts that were off compared to the ground truth with the unfitted AlphaFold structure are much better placed with the fitted AlphaFold structure.

Figure: Figure 10: MD 150k dataset unfitted AlphaFold structure. Left: Histograms of the distances between the two upper domains. The true distances are in green. The recovered distances are in blue. Middle: Predicted against true distances in Ångström. The black line corresponds to $x=y$. Right: Plot of the Fourier shell correlation curves. The mean across structures for cryoSPHERE is in plain blue. Dotted blue represents the mean $\pm$ one standard deviation. The red, plain line represents the mean of cryoDRGN across volumes, with the dotted red line representing the mean $\pm$ one standard deviation.
Refer to caption: extracted/5629471/figures/experiments/phy_150k/distribution_distance.png

Figure: Figure 11: MD 150k dataset fitted AlphaFold structure. Left: Histograms of the distances between the two upper domains. The true distances are in green. The recovered distances are in blue. Middle: predicted against true distances in Ångström. The black line corresponds to $x=y$. Right: Plot of the Fourier shell correlation curves. The mean across structures is in plain blue. Dotted blue represents the mean $\pm$ one standard deviation. The red, plain line represents the mean of cryoDRGN across volumes, with the dotted red line representing the mean $\pm$ one standard deviation.
Refer to caption: extracted/5629471/figures/experiments/phy_150k_fitted/distribution_distance.png

Next, we conduct two additional runs with the 150k dataset: one using the first structure from the MD simulation and the other using the last structure of the MD simulation. Figure 16 displays the results for the first MD structure as base structure, showcasing an excellent agreement with the ground truth. In addition, the FSC curves experience a drop compared to the fitted AlphaFold structure, but they are similar to the unfitted AlphaFold structure. This is likely because using a base structure belonging to the dataset introduces less bias than using an AlphaFold-generated structure with no fitting. Figure 17 presents the same plot, starting from the last MD structures. The conclusions are the same, although we tend to overestimate the distances to some extent. However, the FSC curves are significantly improved and now match the fitted AlphaFold structure run.

Figure 18 shows the segment decompositions recovered by cryoSPHERE for all these runs, revealing their overall similarity. Notably, two small segments are identified in the helices linking the upper part to the lower one. This is not problematic since our decomposition is soft, allowing for smooth transitions between no motion and the motion of the upper domains.

As expected, the quality of the results depends on the quality of the base structure used. Fortunately, AlphaFold can provide reasonable structures that we can further fit into density maps obtained with traditional methods. This approach combines the strength of both world: the more static parts of a protein can exhibit good resolution with traditional tools, and we can fit in the corresponding part of an AlphaFold structure, thereby debiasing the method to some extent. On the more mobile parts where traditional methods may have lower resolution, AlphaFold can provide improved resolution.

In a final experiment, we attempt to increase the number of segments to $N_{\textrm{segm}}=16$ using the unfitted AlphaFold structure as the base structure. We run our algorithm for 24 hours on the same single GPU used in previous experiments. Figure 12 demonstrates an excellent agreement with the ground truth in terms of distances. The FSC curves indicate that increasing the number of segments actually improves the results, approaching the performance of the run with the fitted AlphaFold structure, despite using the unfitted AlphaFold structure. Figure 15 visually confirms this improvement. As expected, more segments allow the model to adjust better to the data, reducing the bias introduced by the AlphaFold structure. However, this improvement comes at a cost: in 24 hours, cryoSPHERE with the unfitted AlphaFold structure and 6 segments makes 642 epochs, whereas the run with 16 segments completes 421 epochs.

Figure: Figure 12: 150k dataset with the unfitted AlphaFold structure and 16 segments. Left: Histograms of the distances of the two upper domains. The true distances are in green. The recovered distances are in blue. Middle: predicted against true distances in Ångström. The black line corresponds to $x=y$. Right: Fourier shell correlation curves for cryoSPHERE. The mean of of the unfitted AlphaFold with 6 segments is in blue. Dotted blue represents the mean $\pm$ one standard deviation. The red, plain line represents the mean of the unfitted AlphaFold run with 16 segments, with the dotted red line representing the mean $\pm$ one standard deviation. The green, plain line represents the mean of the run with the fitted AlphaFold structure and 6 segments.
Refer to caption: extracted/5629471/figures/experiments/phy_16/distribution_distance.png

Figure: Figure 13: Left: comparison of a ground truth volume and the corresponding reconstructed volume from cryoSPHERE. Gray is cryoSPHERE, purple is the ground truth. Middle left: comparison of a ground truth volume and the corresponding reconstructed volume from cryoDRGN. Pink is cryoDRGN, purple is the ground truth. Middle right: cryoSPHERE. Right: cryoDRGN.
Refer to caption: extracted/5629471/figures/comparison_volumes/true_unfitted.png

Figure: Figure 14: left: comparison of a ground truth volume and the corresponding reconstructed volume from cryoSPHERE with the unfitted AlphaFold structure. Gray is cryoSPHERE, purple is the ground truth. Middle left: comparison of a ground truth volume and the corresponding reconstructed volume from cryoSPHERE with the fitted Alphafols structure. Blue is fitted, purple is the ground truth. Middle right: cryoSPHERE without fitting. Right: cryoSPHERE with fitting.
Refer to caption: extracted/5629471/figures/comparison_volumes/true_unfitted.png

Figure: Figure 15: Comparison of volumes for an image of the MD 150k dataset. Gray is cryoSPHERE with unfitted AlphaFold structure and 6 segments. Tan is the same with 16 segments. Blue is ground truth. Having 16 segments leads to a volume matching the ground truth much better than 6 segments. This is true for medium scale elements throughout the volume, even at the bottom of the protein, which tends to stay still. Have a look at the coil on the middle right of the figure. Having 16 segments allows the method to break this coil into a segment and to adjust its position. Which is not the case with only 6 segments.
Refer to caption: extracted/5629471/figures/experiments/phy_16/comparison_16_fitted.png

Figure: Figure 16: Plots for the MD 150k dataset starting from the first MD structure. Left: histograms of the true and predicted distances. Middle:predicted against true distances. Right: mean FSC curves for the AlphaFold base structure and the first MD structure as a base structure.
Refer to caption: extracted/5629471/figures/experiments/phy_150k/close/distribution_distance.png

Figure: Figure 17: Plots for the MD 150k dataset with the last MD structure. Left: histograms of the true and predicted distances. Middle: Predicted against true distances. Right: Mean FSC curves for the AlphaFold base structure and the last MD structure as a base structure.
Refer to caption: extracted/5629471/figures/experiments/phy_150k/open/distribution_distance.png

Figure: Figure 18: Examples of recovered segments. From left to right: 10k dataset, 6 segments. 150k dataset unfitted AlphaFold structure, 6 segments. 150k dataset with the first MD structure as the base structure, 6 segments. 150k dataset with the last MD structure as the base structure, 6 segments. 150k dataset unfitted AlphaFold structure, 16 segments.
Refer to caption: extracted/5629471/figures/experiments/phy_10k/mask_10k.png

#### A.2.2 10k dataset

To assess the performance under data conditions, we run our method with only $10^{4}$ images, $N_{\textrm{segm}}=6$ and the same unfitted AlphaFold base structure and settings as in subsection A.2.1. We train our network and cryoDRGN for only 12 hours only.

Figure 19 demonstrates that cryoSPHERE maintains an excellent agreement with the ground truth and is not significantly impacted by the scarcity of data. On the contrary, the FSC curves indicate that the performance of cryoDRGN is affected. This is further confirmed by Figure 20, where cryoSPHERE has roughly the same FSC curves between the 10k and 150k datasets, while cryoDRGN experiences a significant drop in performance. This is likely because the AlphaFold structure provides strong prior knowledge.

Additionaly, we plot the recovered segments in Figure 18. Despite having fewer images, our method succesfully divides the protein into meaningful parts for the given task.

Figure: Figure 19: For the MD 10k dataset. Left: Histograms of the true and predicted distances. Middle: Predicted against true distances. Right: mean FSC curves for the cryoSPHERE with the same base AlphaFold structure as the 150k dataset and cryoDRGN ran on this 10k dataset. The blue lines are for cryoSPHERE, the red lines for cryoDRGN.
Refer to caption: extracted/5629471/figures/experiments/phy_10k/distribution_distance.png

Figure: Figure 20: For the MD 10k dataset. Left: comparison of the FSC curves for cryoSPHERE on the 150k and 10k datasets. Right: comparison of the FSC curves for cryoDRGN on the 150k and 10k datasets. Blue represents the 150k results and red the 10k results.
Refer to caption: extracted/5629471/figures/experiments/phy_10k/comparison_10k_150k.png

## Peer Reviews

confidence: 5
contribution: 2
initial_rating: 6
presentation: 2
questions: Can the authors elaborate how the density maps (volumes) are reconstructed in cryoSPHERE (like in Fig. 6)? I understand how backprojection in cryo-EM works, but still find it difficult to understand how this is performed exactly.
soundness: 2
strengths: Instead of finding the deformation for each residue, cryoSPHERE learns a segmentation of the amino acid or nucleotide chains and deforms the segment to fit the heterogenous cryo-EM data. This is indeed a valid assumption in many cases especially for large complexes  like the spliceosome (EMPIAR-10180).
summary: The authors introduced a new method, cryoSPHERE, that resolves the continuous heterogeneity problem in cryo-EM reconstruction. Similar to some other methods, cryoSPHERE uses an atomic model as the input, makes segmentations, and uses the results as the regularization in finding the flexible movement in the cryo-EM dataset.
weaknesses: Major:
1. CryoSPHERE uses an atomic model as the reference and formulate the conformational heterogeneity in cryo-EM dataset as the deformation of the segments of the reference model. In finding the optimal segments, the approach is similar to what used in e2gmm, where a GMM with N_{segm} components is fitted. This segmentation result is then used in a cryoSTAR/DynaMight-like setting, replacing the regularization losses and compute predicted projections to compare with the particle images. The author should also compare the result with other similar methods, including DynaMight, e2gmm, and maybe 3DFlex, and discuss the outcomes.

2. In cryo-EM, density maps are extremely important in evaluating the validity of the results. However, the authors only presented the density maps for one of the synthetic dataset. Especially in the case of under the hypothesis that the heterogeneity is modeled as the deformations of a canonical density, it is very important that the reconstructed canonical density is free from any obvious artifact.

3. The FSC comparison in Sec. 5.2 between cryoSPHERE and cryoDRGN and cryoSTAR may not be very fair, since cryoSPHERE's density is reconstructed under the assumption that there exist a canonical density, while cryoDRGN does not hold the same hypothesis, and cryoSTAR explicitly avoids this to reduce the unwanted bias from the reference model. A fair comparison would be made against DynaMight and 3DFlex.

4. Lack of ablation study. The authors should discuss how N_{segm}, which is a key hyperparameter that is related to the degree of freedom in the segmentation, affects the results in both synthetic and real data.

5. The authors compare cryoSPHERE to cryoDRGN and cryoSTAR, but from the experiments, it is unclear what the advantage that cryoSPHERE brings to the field. The only advantage seems to be the increase FSC resolution, but it is only shown in one synthetic dataset.

Minor:
1. The sequence of Fig.6 and Fig.7 in the paper is reversed.
2. In the related works, e2gmm should also be considered as a method using deep learning.
3. In the FSC comparison, if not comparing two half maps, the cutoff of 0.5 should be used and 0.143 is not meaningful. See [1] appendix for the reason (and also why 0.143 is the number widely used for half maps comprison). Therefore the discussion at the end of Sec. 5.2 about 0.5 and 0.143 cutoff improvement is not very correct.
4. The writing could be better polished.
5. I do not think the "Ethics Statement" content is proper as it lacks evidence.


[1] Rosenthal, Peter B., and Richard Henderson. "Optimal determination of particle orientation, absolute hand, and contrast loss in single-particle electron cryomicroscopy." Journal of molecular biology 333.4 (2003): 721-745.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 3
questions: Questions for clarification 
- In section 3.6 the authors say “Similar to Li et al. (2023), instead of using a mean squared error loss between the predicted and ground truth image, we use a correlation loss between the true and predicted image”. However, they do not motivate this choice. Would you please care to elaborate?
- In section 5.1 the authors write “Testing the segment decomposition, we then run cryoSPHERE by requesting division into Nsegm = 4. The program learnt a first and third segment with 0 residues, a second segment with 1353 residues and a fourth segment with 157 residues (Figure 5). Thus, cryoSPHERE learnt segments according to the ground truth.” Is there a reason why we should expect that the model prefers to learn two segments?

Additional feedback
- It reads as if equation 7 is the only loss that is being used. I assume that the VAE latent distribution is also trainable? It would benefit the readability if the authors clarify whether the loss is just the decoder loss or the whole loss.“
- In section 5.3 the authors write “The first two principal components explain more than 96 percent of the variance.” In general, I would be very reluctant to make any claims about PCA in the latent space as you could come up with a transformation that scrambles the latent space and still gives the same decoding, but covers much less variance. I can imagine that Figure 7 is a nice illustration, but I would recommend to leave it with that.
soundness: 3
strengths: The work is of high quality and overall very well-written. I’d like to highlight the following strong points of the paper:
- The authors seem very well aware of the difficulties atomistic models have in terms of the optimisation landscape and address the problem adequately by regularising while keeping sequential information intact
- The authors give useful biological remarks such as that their methods seems to recover domains, which is useful to know as a practitioner
- The authors show that their method works on real data, which is actually the first time I’ve read this for atomistic models
summary: - The authors propose a VAE and GMM-based approach to atomistic heterogeneous Cryo-EM reconstruction
- The GMM part learns how to divide the amino acid chan in segments
- The VAE learns how to encode images into latent variables and decode these variables along with the GMM output into segment-wise conformational changes
- The authors claim that their method overcomes current challenges that atomistic models face, while also achieving improved interpretability by learning these segments (that ideally correspond to domains in the biological sense)
weaknesses: There is one main weak point in the paper in my regard:
- The authors do not address limitations of their method very thoroughly:
    - whether or not their method is able to overcome poor initialisation of the initial structure S_0. This is one of the major risks of working with (pseudo)-atomistic models as shown in works such as DynaMight (Schwab et al. 2024) 
    - Whether or not local physical information such as bond lengths between atoms at the boundaries of the segments is preserved. This might be challenging

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 2
questions: Could you include additional comparisons with methods like CryoSTAR, or Dynamight (e.g., computational cost, novelty, etc)
soundness: 3
strengths: 1. The approach of learning to decompose the protein’s amino acid chain into segments to represent various conformations is novel.
2. The interpretability of moving part is highly beneficial for practitioners, offering a clearer understanding of conformational dynamics.
summary: This paper introduces cryoSPHERE, a deep learning model that uses nominal protein structures (e.g., from AlphaFold) to break proteins into segments and fit these segments to various conformations in cryo-EM data.
weaknesses: 1. No meaningful improvements were observed over CryoStar
    - Figures 3, 12, 13: It’s unclear what meaningful differences exist compared to CryoStar.
    - Figures 14, 21: There doesn’t appear to be a significant difference between CryoStar and CryoSphere (CryoStar may even appear better).
    - Figure 29: It would be great to compare CryoStar with CryoSphere rather than with CryoDRGN.
2. Should change Figure 6:
    - Instead of showing only CryoDRGN results and G.T, a comparison with CryoSTAR or Dynamight would be more appropriate (perhaps include Figure 32 and the CryoSphere result, as there are no CryoSphere results with SNR 0.001).

Meta-review: The paper proposes cryoSPHERE, a method for reconstructing proteins that can occur in different configurations. 

The paper received three reviews, and after the rebuttal all three reviewers are supportive of publication.

* Reviewer RJxo was initially not supportive of publication; the reviewer noted that that some comparisons are not fair, and pointed out that it is unclear what the advantage of the method is over methods like cryoSTAR and others. In the rebuttal, the authors highlight that cryoSPHERE can provide physically plausible reconstructions for examples where cryoSTAR cannnot. After the rebuttal, the reviewer still has concerns about the performance of the method, but finds the novelty acceptable to recommend weak accept.

* Reviewer NciV is supportive of publication and highlights that the method is evaluated on real data, contrarily to other works in the area that are based on atomistic models. 

* Reviewer JGA5 finds elements of the approach to be novel and the interpretability of the moving part beneficial for understanding conformational dynamics of proteins. However, the reviewer also notet that there are no meaningful improvements over CryoSTAR. In the rebuttal, the authors provided new figures that partially addressed the reviewers concerns, and while the reviewer still has concerns about the novelty relative to CryoSTAR, the reviewer raised the score to weak accept.

The paper provides a new method for reconstructing proteins in different configuration, and demonstrates that the method has, in some cases, a slight benefit over existing methods, in particular over cryoSPHERE. I agree with the reviewers that while the method uses many elements from prior methods (such as the way the VAE is used, the atomistic representations, etc) and while the benefit over existing methods is relativly marginal, the method and evaluations contribute to the important problem of reconstructing a single protein in different configurations, and I therefore recommend to accept the paper.

Decision: Accept (Poster)