Symmetry-Driven Discovery of Dynamical Variables in Molecular Simulations

Thank you for your comments. WE have updated the figures in the paper and are working on fixing typos. Our detailed response:

1. Only validated on "toy" example Alanine dipeptide: While we agree Alanine dipeptide is toy baseline (the MNIST of comp bio), we want to emphasize that, to our knowledge, no prior work has discovered the torsion angles in the interpretable way our method does, and certainly not without requiring data. Therefore, we felt that it was absolutely necessary to validate our method on Alanine dipeptide. We do agree that larger systems would be the true test and we are working on longer peptide chains and small proteins. The challenge is that a nice phase space like the one for Ala dipeptide does not exist for these larger systems, so it would take a lot more computational biology discussions and simulations to have presentable results. This would not leave enough room to discuss our methodology and therefore we think larger systems should be a follow-up work.
2. Not clear the potential performance of the model on more general MD simulations, in which more degrees of freedom are standard: One advantage of our method is that it can discover many DoF easily. In fact, the second step (maximization of symmetry-breaking, eq 15) is there to sift out the many DoFs and just pick a few. Most steps in our method are very cheap and just require a simple eigendecomposition. For example, using Hessian backbone degenerate subspace, we easily found 6 degenerate vectors. Any pair of these vectors can be used to define an $\\mathbf{L}$ as a rotation in this 2D subspace, yielding 15 (6 choose 2) DoFs.
3. Not compare with other models…not clear what is the potential advantage over existing models: We have added a general comment above where we discuss related work. We think that our method provides an orthogonal contribution to the works mentioned above, because:
   1. Conceptual leap: Other works think of DoF (or collective variables) as vectors, while we map them to transformations
   2. Data-free: All other methods require many real or simulated conformations, but our method just needs a decent energy function (forcefield).
4. Code not available: We have added the code in the Supplementary Materials.
5. Interesting algorithm for analyzing MD, it does not fall into the general machine/deep learning category…more suitable for a computational chemistry/biology journal: We should have chosen the category more carefully. But, we believe because this is an ML algorithm, an ML venue is much better qualified to judge its validity than a bio venue. Additionally, comp bio, drug discovery and physical simulation are among the fastest growing and most important use cases of ML.
6. Typos and writing quality: We apologize for this. We are doing our best to polish the paper. We have updated the figures, and will send another note when the typos have been fixed.

We have now shared the code in the supplement.

10. What optimization algorithm? How costly, in terms of number of energy gradients? Since this system was small, we directly solved the quadratic problem by doing an eigendecomposition on the matrix. For larger systems one can use linear solvers or gradient-based methods. In the end finding Ls is like a linear regression problem with some regularizers. The only subtlety is whether the large gap in the magnitude of gradients in different directions (i.e. hierarchy of interaction strengths) can cause an issue. This is hierarchy certainly exists in our experiments, but it didn’t cause any issues. We suspect this is because we initialized from an already decent conformer, where bond and angles energies were already minimized. In practice any primitive initial layout for the molecule should work, e.g. using RDKit, which also gets bond lengths and angles right. For the Hessian-based method, we require only one Hessian evaluation, while the optimization-based approach needs us to evaluate O(n^2) force values (energy gradients) on points sampled near the initial conformation.

11. What are the min/max bounds on the 31x31 gridpoints? The transformations learned in our method correspond to Lie algera elements. We get the 31 grid points by repeatedly applying the learned transform on the starting configuration. x_t = x_0 . L^(2*pi*(t/31)) where L is the transformation matrix.

12. suggested papers: Thank you very much. We will review them shortly and write about them.

We really appreciate your helpful comments. Detailed response:

1. Qualitative conformation: We are adding visualizations of molecules transformed using our DoF overlaid onto initial conformation. We'll notify when ready.
2. DoF invariant in solvent and vacuum: We are rewording this. We meant to say that the DoF were transferable to other conditions: We observe the DOFs discovered in vacuum, and for one starting configurations, allowed us to reach various conformations of alanine dipeptide even in the presence of solvent. Note that the DOFs only to provide starting conformer from which we run OpenMM simulations. These starting points covered a wide enough range of the phi-psi angles to allow us to discover all the conformers with short simulations.
3. What does the beta-sheet conformation look like: It's as in top left of Fig 4 ($\beta$). It is nearly planar, with $\phi,\psi = -150,150$ falling at the upper left corner of the phase plot. We are updating the images and are using a star to show the location of the $\beta$ sheet conformer in the phase space. Other conformers are the hotspots of the density plots. For each of the 3D conformers in the old Fig 4 we will put a marker on the phi-psi plot.
4. In Figs 1 and 3 gridlines differ for each method, do not cover the entire phi-psi: That is correct. Each method has some pros and cons. Also, we noticed we had bug in plotting the gridlines, which is fixed now. Please see the general comment above for detailed findings regarding each method. TLDR; Only the Hessian subspace rotation method covers all of phi-psi, but all other method do a better job at finding more rare conformations.
5. Fig 4 do not resemble the leftmost column: Note that phi and psi are only the angles of the two bonds around the central “stem” pointing down. These angles need to match, while the rest of the molecule need not. For instance, the methyl groups on the left and right of the molecule can rotate up or down without costing much energy. That being said, we are improving this plot, and experimenting with overlaying the conformers for better comparison. We will update you on this soon.
6. stochastic - do multiple trials from the same init converge to similar effective DoF? Note that our Hessian approaches are not stochastic. But regarding our direct symm loss approach, you are right, they do not converge to the same DoF. Nevertheless, moving along them still allowed us to find all known conformers. We are adding a figure to the appendix showing multiple runs of the direct loss and their trajectories in phi-psi. We will reference it in the coming days. The only method that almost perfectly covers the entire phi-psi is the Hessian backbone degenerate subspace.
7. Method locally searches for flat directions, would not be able to discover new minima behind large energy barriers? Actually, it can. First, our flatness is approximate. More importantly, the second maximization step (eq 14 and 15) actively works against flatness. It finds directions where the gradient increases the most, thereby moving straight toward energy barriers or walls. The way the flat directions should be understood is by considering the hierarchy of forces here. For example, the quadratic bond length and angle potentials introduce directions in the Hessian which are constant and very large. Compared to these large components, directions corresponding to weak van der Waals forces are extremely flat, but eventually show curvature when moving far enough from the current point.
8. Method only requires calculating the Hessian once? Yes we only compute the Hessian once for the current experiments. However, the experiments can be done in an iterative manner where newly discovered conformers serve as a starting point for computing local DOFs. For our experiments only one round of Hessian computation suffices to find conformers of alanine dipeptide.
9. What, if any, challenges are expected when applied to neural force fields? In principle, our direct optimization method should directly work with any forcefield, including neural. In this work, we directly use the OpenMM forcefields. Our Hessian-based methods, on the other hand, require the force fields to provide accurate second order information. This might not hold true for ReLU based networks which have 0 Hessian almost everywhere. Our direct optimization approach should still work as it only requires an energy function and the curvature is estimated from perturbed samples instead of analytically.

More in next comment

We have added Figs 7, 8 to show how the DoF gridlines change when running direct optimization multiple times with different perturbed samples. We show three runs for different setting with or without solvents. Each column shows one setting. Overall, the number of conformers found in each run varies. We will check if increasing the number of samples can make this approach more robust.

Thank you for the very pertinent comments. Our response:

1. They don’t mention similar work: We are improving the related work. Please see general comments above.

2. Figs need improvement. Which is a? b? c? etc in Figure 1: We agree and apologize. We are updating the figures. Please see the updated paper. We have remade all figures:

   1. Added abc to Figs 1, 3
   2. Aligned layouts of Figs 2, 4 to make comparison easier
   3. Fixed a bug in drawing gridlines in Figs 1, 3
   4. Labeled location of known conformers on the density plots.

   Please also read our summary in the general comment above about what each of our four algorithms discovers.

3. Plot the distribution of energies: Thanks, a very good point. We have added Figs 5, 6 showing the energies of the discovered conformations. Please also see our general comment above titled Energies of conformations which shows tables of energies of prominent conformations with or without solvent using each of our methods. In short, our energies are very good. We are adding the tables to the appendix.

4. Code not available: We have added the code in the Supplementary Materials.

5. ** $maybe$ the Ramachandran plots look decent, but the energies are completely off and non-physical:** We should clarify that we take steps to avoid unphysical states. Our procedure is:

   1. Use discovered DoF to get a starting configuration for short simulations.
   2. Use openmm optimize to resolve clashes
   3. Run short openMM simulation (1000 steps)

   It is true that the starting point we get from our DoF often has very high energies because we only look for approximate symmetries and the maximization step (eq 15) explicitly chooses the DoFs that maximally change the energy (i.e. violate the symmetry condition the most). However, one call of openMM.minimize and a short simulation quickly resolves allosteric clashes and yields good physical states. In the energy distributions in Fig 5, 6 we remove the first 100 steps of simulations to avoid the steps where state may have been unphysical.

6. They mention that 'we see that almost all the Hessian-based methods recover all the major conformations .. with short simulations.': Our new Figs 5, 6 now show this more clearly. The distributions shown below each contour plot are the different steps of the 1000 steps (each 1 femtoseconds) openMM simulations, excluding the first 100 steps to remove initial non-physical states. The colors of points show their energies. We observe that all methods, *except* Hessian degenerate subspace, explore all conformations, including rare ones not explored by the long simulation (baseline).

7. They do not compare to other baselines in recent literature. As discussed in our general comment on related work, we feel our method provides an orthogonal contribution to the works mentioned above, because:

   1. Our main objective is to data-free exploration of the collective variables (CV) space and
   2. Unlike the given papers where conformation or trajectory between them is known, our task is to discover new conformers using locally-approximated CVs computed only using the atom positions and the local energy function.

   It may be possible to use the points in the long simulation to apply one of the Cv methods in the papers above, but we can say immediately that our method is faster and cheaper than that because it doesn’t require many simulation points or trajectories. Our method could also be used for the tasks discussed in those papers. Future directions can be understanding/sampling transition paths or providing better estimations of the free-energy surface using the afore-mentioned methods.

Energies are in units $kJ/mol$. The direct optimization method has two $\epsilon$ parameters (units nanometers), one for sampling around the starting point to estimate the Hessian (eqs 5 and 7), and one for the maximization step to choose best $\mathbf{L}$'s (eq 15).

Vacuum Results

The last three columns are for the three prominent conformations.

Solvent Results

These are for simulations in water. The last three columns are for the three prominent conformations. "Solvent Direct Optim" means using the symmetry loss eq 3 and 5 with the energy function of a forcefield that includes explicit solvent terms (using openMM).

Please see the supplement zip file for the code.