Deflate to Understand Complex Molecular Kinetics

Brooke E. Husic & Frank Noé. arXiv: 1907.04101

Contributed by Jesper Madsen

Dimensionality reduction is at the core of understanding and making intuitive sense of complex dynamic phenomena in chemistry.  It is usually assumed that the slowest mode is the one of primary interest; however, it is critical to realize that this is not always so! A conceptual example hereof is a protein folding simulation (Lindorff-Larsen et al. Science 334, 517-520, 2011) where the slowest dynamical mode is not the folding itself (see Figure). What is the influence, then, of “non-slowest” modes in this process and how can it most appropriately be elucidated?

FIG: Figure 2 from the preprint: "(A) Sampled villin structures from the MD trajectory analyzed. Helical secondary structure is colored and coils are white. Each image represents five structures sampled from similar locations in TIC space as determined by a 250-center k-means model built upon the first three original TICs. The purple structure represents the folded state, and the blue structure represents the denatured state. The green structure is a rare helical misfolded state that we assert is an artifact. (B) Two-dimensional histograms for TICA transformations constructed from villin contact distances. Dashed lines indicate the regions corresponding to the sampled structures of the same color. The first TIC tracks the conversion to and from the rare artifact only. The second TIC tracks the majority of the folding process and correlates well with RMSD to the folded structure."

This work by Husic and Noé show how deflation can provide an answer to these questions. Technically speaking deflation is a collection of methods for how to modify a matrix after the largest eigenvalue is known in order to find the rest. In their provided example of the folding simulation, the dominant Time-lagged Independent Component (TIC) encapsulates the "artifact" variation that we are not really interested in. Thus, a constructed kinetic (Markov-state) model will be contaminated in several undesirable ways as discussed by the authors in great detail.  

In principle, this should be a very common problem since chemical systems have complex Hamiltonians. Perhaps the reason why we don’t see it discussed more is that ultra-rare events – real or artifact – may not usually be sampled during conventional simulations. So, with the increasing computational power available to us, and simulations approaching ever-longer timescales, this is likely something that we need to be able to handle. This preprint describes well how one can think about attacking these potential difficulties.   

Discovering chemistry with an ab initio nanoreactor.

Lee-Ping Wang, Alexey Titov, Robert McGibbon, Fang Liu, Vijay S. Pande and Todd J. Martínez

Nature Chemistry

6, 1044–1048 (2014)

Contributed by Alán Aspuru-Guzik

When is quantum chemistry going to take the primary discovery role, rather than an explanatory one? This one of the questions asked in the recent Nature Chemistry paper by the groups of Vijay Pande and Todd Martínez at Stanford.

With the advent of fast GPGPU quantum chemistry in codes such as Terachem, it is practical to carry out medium-scale quantum dynamics simulations for a considerable amount of time. In this paper, the authors employ a spherical "piston" to compress and decompress high-temperature reaction mixtures, which in turn result in an accelerated process for triggering chemical reactions. Automated algorithms are employed to recognize the nature of the chemical transformations involved. The authors use this approach to identify novel mechanisms for the Urey-Miller origins of life experiment.

In this modern world of Facebook, Twitter, Youtube, Facebook and blogs, a video is worth more than a thousand words. The video below summarizes the paper very well.

This paper is indeed, one of my top 5 favorite papers in the field of computational chemistry of 2014. 

Todd Martínez explaining the future of the nanoreactor in the QUITEL 2014 conference held in the Galápagos Islands in Ecuador. He suggested, in a joking fashion, that the nanoreactor could eventually evolve a flask with molecules into life ab initio. The virtual end product of the evolution seems to be the contributor of this post.

Will molecular dynamics simulations of proteins ever reach equilibrium

Samuel Genheden, Ulf Ryde. Phys. Chem. Chem.Phys., 2012, 14, 8662-8677

Contributed by Gerald Monard

When one runs molecular dynamics simulations, one always wonders whether the trajectories are “long enough”. The ergodic hypothesis states that with an infinite trajectory (in time), all accessible states should be reached. However, one cannot wait an infinite amount of time. In practice, trajectories are finite. Thus the question: did my molecular dynamics simulation have reach equilibrium?

S. Genheden and U. Ryde have analyzed conformational entropies of five protein and protein-ligand complexes with different methods for simulations up to 500 ns. They show that, for all cases, conformational entropies have not converged and that, therefore, all available states have not been reached.

Conformational entropies have been evaluated from the distribution of the dihedral angles, using different models: dihedral-distribution histogramming (DDH), von Mises approach (VMA), quasi-harmonic analysis (QHA), and normal-mode analysis (NMA). The latter was used in the estimation of ligand-protein binding free energies with the MM/GBSA method.

What is interesting in this article is that, whatever the entropy model used, no simulation converged and that entropies are still increasing at the end of the simulations. To explain this, the authors used a very simple model of protein containing only dihedral angles. In this case, the entropy of the model protein increases logarithmically with the simulation time, and the conformational entropy does not converge until the simulation time is larger than all the normal mode mean lifetimes (from 0.2 ps for 1 kJ/mol barrier to 12 h for 100 kJ/mol barrier). Using such an approach, the authors state that after 500 ns of simulation, most of their systems have sampled less than 70% of the conformational spaces and that the entropy will continue to increase with simulation time. Moreover, in the case of BPTI, after 1 ms of simulation, entropy has still not converged!

So, does this mean that we are all wrong when performing MD simulations? Will MD simulations ever reach equilibrium? In the strict sense, the answer is no: MD simulations do not converge, at least when regarding the entropy term. But, do we really need full convergence? As the authors say, reaching all available states for a protein would mean to fold and unfold the protein several times during the MD simulation. Do we (always) need this? In most cases, no.

Another interesting point is the MM/GBSA calculations. Using the NMA model to compute the entropy, the ligand-protein binding free energies converge in the 500 ns simulation time, showing a variation of less than a few kJ/mol. However, this is mainly due to error cancellation in the sum that defines the binding free energy in the MM/GBSA approach. In addition, the NMA approach produces much more stable estimates of the entropy than the other three methods (DDH, VMA, and QHA). This is due to the NMA approach that assumes a single conformation for the protein and estimates the entropy form the vibrational spectrum in this conformation. As long as the protein does not drastically change its conformation and its vibrational spectrum, the NMA entropy does not vary much (but then one could argue that this means that the NMA approach may not be very accurate).

In conclusion, I would say that this article shows that it is hopeless to see MD simulations of proteins to ever reach equilibrium (at least using a single trajectory). But that does not mean we should stop running MD simulations. We have to be (very) careful about the quantities that we are looking at. Different properties show varying sensitivity to this lack of strict equilibration.

Atomic-level simulations of current-voltage relationships in single-file ion channels

M.Ø. Jensen; V. Jogini; M. P. Eastwood and D.E. Shaw. J.Gen. Physiol, 141: 619-632 (2013)

Contributed by Ben Corry

There has been much discussion in these pages about papers assessing the accuracy of classical molecular dynamics in reproducing structural aspects of proteins, but less about whether such simulations can be used to understand physiological processes. Many such processes take place over time scales of ms to s and simulating them directly has long been beyond the scope of atomistic methods. But, the situation is rapidly changing with advances in computer hardware and software. One example of an important physiological process close to my heart is the conduction of ions through trans-membrane channels which underlies electrical signalling in nerve cells and takes place on a ms time scale. I am sure that many of us working in the field have thought that direct atomistic simulations have the potential to elucidate the physical mechanisms underlying this process, and have wondered how accurately such simulations could reproduce this phenomena. Yet again, D.E. Shaw Research provide us with an answer.

Making use of extensive computational resources and more than 1ms of simulation time, Jensen et al¹ measure the ionic current passing through a voltage gated potassium channel and the simpler gramicidin A channel under a range of voltages, allowing for a direct comparison to one of the most fundamental experimentally measureable properties. Sadly, the results are disappointing. Currents are about 40 times less and 300 times less than equivalent (highly accurate) experimental measurements in the potassium channel and gramicidin A respectively.

What are the reasons for this poor performance? Jensen et al point the finger at the most likely culprit, the accuracy of the non-polarisable biomolecular force field. Altering parameters such as the interaction strength between permeating ions and the protein has only a very small influence on the calculated current highlighting that simple modifications are unlikely to resolve the discrepancy between simulation and experiment. Deficiencies in the lipid model, such as the overestimation of the membrane dipolar potential, are also discussed; and while these are suggested to be a significant factor in the poor performance of the simulation they do not appear to be the major reason for the underestimation of ion currents. The final suggestion is that polarisable force fields may be required to accurately reproduce permeation rates under experimental conditions. Indeed, in the simple case of gA in which ions permeate one at a time, it would appear plausible that the inclusion of polarisability would improve the results by stabilising ions in the pore and reducing the barriers to permeation. It is not discussed is how much the structure of the proteins change during the simulations, especially under the influence of the electric field. If the structures of the proteins deviate from reality (again due to force field limitations) it is possible that this could also contribute to the poor performance of the simulations, in addition to problems with the ion-protein interactions and lack of polarisability. Given that it has only just become feasible to directly simulate ion currents with non polarisable force fields, it may be some time yet before we know if the use of polarisable ones will make the brute force simulation of physiological processes reliable.

References

1.       M.Ø. Jensen; V. Jogini; M. P. Eastwood and D.E. Shaw. J. Gen. Physiol, 141: 619-632