Feedback to your DynaMight preprint

[biochem-fan]

I enjoyed reading your preprint "DynaMight: estimating molecular motions with improved reconstruction from cryo-EM images". Although I could have sent this by private email, now that the preprint is public, I thought it would be valuable to make my comments public as well. (And GitHub Markdown gives nice formatting.) Feel free to reply by email if you prefer to do so.

Several factual errors

You say "Flex3D" in page 13 and 14 but the right name is "3DFlex". Their paper title is "3DFlex: determining structure and motion of flexible proteins from cryo-EM".

Page 2

the Zernike3D approach expresses the cryo-EM map in a basis of 3D Zernike polynomials

This is not correct. Zernike3D expresses the deformation field of a cryo-EM map using 3D Zernike polynomials.

Page 14

It is possible that iterative real-space methods like stochastic gradient descent or algebraic reconstruction techniques, as implemented in Flex3D [12] or Zernike 3D [13], may yield better results.

3DFlex does not use stochastic gradient descent for 3D reconstruction. They use full-batch L-BFGS.

They say in their paper:

We also found that minibatch stochastic gradient descent methods for directly optimizing V did not yield high quality results, potentially because noise in minibatch gradient estimates is problematic for this task. Instead, we were able to solve the problem using full-batch L-BFGS (ref. 16) in real space

Something unclear in the manuscript

Page 6

Additional penalties that prevent Gaussians coming too close to each other, or moving too far away from other Gaussians, also exist

Can you add the formula for this?

We train a neural network as a regression function to estimate a deformation field that coincides on the given sampling points Γ(ci ), but can be evaluated on arbitrary positions

Why do you need another neural network for this? Cannot you average shift vectors of nearest N Gaussians, perhaps weighted by the distance from the point to be interpolated? This doesn't add any additional parameters to model and probably is safer and more robust.

Page 7

I don't understand the equation 7 and discussion following it...

How is the equation 7 derived? What is the explicit form of the matrix D?

We approximate equation (7) by using the filter that would correspond to the homogeneous case.

What is the explicit form of "the filter"?

Although even in the optimal scenario of having complete data of clean projection images, this method does not yield a minimum of functional (6), it still allows to correct for the deformation to some degree.

This is again vague. Is the approximation derived theoretically and valid, for example, up to the first order? Or is it an arbitrary heuristics?

Suggestions (for future research; not necessarily this paper)

Why is the resolution worse than MultiBody refinement? Is the problem in the estimation of the deformation model or the new reconstruction algorithm?

• You can test this by running your new reconstruction algorithm using deformations from the MultiBody refinement job.
• You can also compare output deformations of your decoder (or your regression neutral network) with motion from the good MultiBody refinement job.
• What happens if you replace the regression neural network with a simpler interpolation scheme I proposed above? The regression network might be behaving wild.

You can also perform above tests using a synthetic dataset with known deformations.

[dkimanius]

Thanks Takanori!

As for the factual errors: "Flex3D" is a typo and "stochastic gradient descent" should be "gradient based".

D in eq. (7) is diagonal (in Fourier space) because of the Fourier slice theorem (without the defomrations). The diagonal is one over CTF^2. See eq. (9) in our old paper (Kimanius et al. 2021). In real-space (with deformations) it becomes a non-linear projection. The approximation is that we don't deform it.

We'll clearify these.

[biochem-fan]

Thanks for clarification.

D in eq. (7) is diagonal (in Fourier space) because of the Fourier slice theorem (without the defomrations)

This is not correct in the strict sense, because Fourier components mix due to trilinear interpolation in the reciprocal space. Gridding correction compensates for this but we abandoned it. Fortunately, the fact that this didn't make maps any worse means that this effect is negligible.

In real-space (with deformations) it becomes a non-linear projection. The approximation is that we don't deform it.

I guess the problem is that we don't know how to deform CTF in the real space. If there were no CTF, we can accumulate interpolation weights along curved trajectories in the real space. Is my understanding correct?

We'll clearify these

One way to make this more approachable is to explain the algorithm in a concrete way. e.g. The Fourier transform of a particle is multiplied with CTF in the reciprocal space and inverse Fourier transformed. This image is backprojected in real space along a curved trajectory.

Also note that most readers probably have never heard of adjoint operators. $P^{*}$ is explained but $C^*$ not. For real valued CTF (anti-symmetric aberrations pre-corrected), this is the same as $C$, right?

[schwabjohannes]

Hi Takanori,

It is true that the matrix after discretization is not diagonal (because of interpolation). But in a continuous (function space) setting the underlying operator is diagonal and is multiplication in Fourier space.

I guess the problem is that we don't know how to deform CTF in the real space. If there were no CTF, we can accumulate interpolation weights along curved trajectories in the real space. Is my understanding correct?

Even if there were no CTF, it would not be that easy. If we define (linear) operators $\Gamma_i$ (discretized this would be a matrix in $\mathbf{R}^{N^3\times N^3}$, where $N$ is the box size) that perform the deformations of a volume, then the operator $D$ from equation (7) is $\sum_i \Gamma_i^\ast P_i^\ast P_i \Gamma_i$. I don't know if there is a simple way to compute these opertors, but it is not a diagonal operator (can not be computed by multiplication in Fourier space) and would correspond a full matrix of size $N^3\times N^3$.

One way to make this more approachable is to explain the algorithm in a concrete way. e.g. The Fourier transform of a particle is multiplied with CTF in the reciprocal space and inverse Fourier transformed. This image is backprojected in real space along a curved trajectory.

I agree, that we could have explained that better and in a concrete way.

Also note that most readers probably have never heard of adjoint operators. is explained but not. For real valued CTF (anti-symmetric aberrations pre-corrected), this is the same as , right?

Yes thats right.

[biochem-fan]

It is true that the matrix after discretization is not diagonal (because of interpolation). But in a continuous (function space) setting the underlying operator is diagonal and is multiplication in Fourier space.

I agree. Probably the algorithm should be described in two steps: first explain mathematics and then describe how it is actually implemented in a concrete way (i.e. discretization, interpolation, Fourier transform etc).

Even if there were no CTF, it would not be that easy.

Reading this paragraph, I can follow your logic. But at the same time I cannot see why my naive algorithm (i.e. accumulate values and interpolation weights along curved trajectories in the real space and divide) does not work. Probably I should study real space Filtered Back Projection theories again.

Thanks for explanations!

[biochem-fan]

In the Zernike3D paper, their ART update formula (Eq. 9) is

$$ V(r)^{k+1} = V(r)^k + \lambda P_H^* \left(P_H V(r + g_L(r)) - I_k(s) \right).$$

$P_H V(r + g_L(r))$ is a projection of a deformed reference volume and $I_k(s)$ is an observed particle image. Because they apply $P_H^*$, i.e. the non-deformed backprojector, to the residual, they seem to use a similar approximation as yours.