Energy conservation in molecular dynamics simulations of classical systems

Non-unique Hamiltonians for discrete symplectic dynamics

An outstanding property of any Hamiltonian system is the symplecticity of its flow, namely, the continuous trajectory preserves volume in phase space. Given a symplectic but discrete trajectory generated by a transition matrix applied at a fixed time-increment (τ > 0), it was generally believed that there exists a unique Hamiltonian producing a continuous trajectory that coincides at all discrete times (t = nτ with n integers) as long as τ is small enough. However, it is now exactly demonstrated that, for any given discrete symplectic dynamics of a harmonic oscillator, there exist an infinite number of real-valued Hamiltonians for any small value of τ and an infinite number of complex-valued Hamiltonians for any large value of τ. In addition, when the transition matrix is similar to a Jordan normal form with the supradiagonal element of 1 and the two identical diagonal elements of either 1 or -1, only one solution to the Hamiltonian is found for the case with the diagonal elements of 1, but no solution can be found for the other case.

Energy stable scheme for random batch molecular dynamics

The computational bottleneck of molecular dynamics is pairwise additive long-range interactions between particles. The random batch Ewald (RBE) method provides a highly efficient and superscalable solver for long-range interactions, but the stochastic nature of this algorithm leads to unphysical self-heating effect during the simulation. We propose an energy stable scheme (ESS) for particle systems by employing a Berendsen-type energy bath. The scheme removes the notorious energy drift, which exists due to the force error even when a symplectic integrator is employed. Combining the RBE with the ESS, the new method provides a perfect solution to the computational bottleneck of molecular dynamics at the microcanonical ensemble. Numerical results for a primitive electrolyte and all-atom pure water systems demonstrate the attractive performance of the algorithm, including its dramatically high accuracy, linear complexity, and overcoming the energy drift for long-time simulations.

Energy, temperature, and heat capacity in discrete classical dynamics

Simulations of objects with classical dynamics are in fact a particular version of discrete dynamics, since almost all the classical dynamics simulations in natural science are performed with the use of the simple "leapfrog" or "Verlet" algorithm. It was, however, Newton who in Principia, Proposition I in 1687 first formulated the discrete algorithm, which much later in 1967 was rederived by L. Verlet. Verlet also formulated a first-order approximation for the velocity v(t) at time t, which has been used in simulations since then. The approximated expressions for v(t) and the kinetic energy lead to severe errors in the thermodynamics at high densities, temperatures, strong repulsive forces, or for large discrete time increments used in discrete "molecular dynamics" (MD) simulations. Here we derive the exact expressions for the discrete dynamics, and show by simulations of a Lennard-Jones system that these expressions now result in equality between temperatures determined from the kinetic energies and the corresponding configurational temperatures determined from the expresssion of Landau and Lifshitz, derived from the forces.

An in-silico scaffold- hopping approach to design novel inhibitors against gp130: A potential therapeutic application in cancer and Covid-19

An upregulation of the gp130-signalling cascade has been reported in multiple cancers, making gp130 an attractive target for the development of anticancer drugs. An inverted-funnel-like approach was utilised along with various structure-based drug designing strategies to discover and optimise novel potential inhibitors of gp130. The study resulted in the discovery of 2 ligands- 435 and 510, both of which exhibit a very high-binding affinity towards the gp130 D1 domain which controls cytokine recognition and interaction thus being involved in complexation. The two resulting complexes remained stable over time with the ligands maintaining a steady interaction with the target. This inference is drawn from their RMSD, Rg, SASA and RMSF analysis. We also tested the protein folding patterns based on their principal component analysis, energy of surface and landscape. The leads also displayed a more favourable ADMET profile than their parent compounds. The two lead candidates show a better therapeutic profile in comparison to the two existing drugs- bazedoxifene and raloxifene. Both these potential leads can be addressed for their activity in-vitro and can be used as a potential anti-cancer treatment as well as to combat Covid-19 related cytokine storm.

Energy-conserving molecular dynamics is not energy conserving

Molecular dynamics (MD) is a widely-used tool for simulating molecular and materials properties. It is common wisdom that molecular dynamics simulations should obey physical laws and, hence, lots of effort is put into ensuring that molecular dynamics simulations are energy conserving. The emergence of machine learning (ML) potentials for MD leads to a growing realization that monitoring conservation of energy during simulations is of low utility because the dynamics is often unphysically dissociative. Other ML methods for MD are not based on a potential and provide only forces or trajectories which are reasonable but not necessarily energy-conserving. Here we propose to clearly distinguish between the simulation-energy and true-energy conservation and highlight that the simulations should focus on decreasing the degree of true-energy non-conservation. We introduce very simple, new criteria for evaluating the quality of molecular dynamics by estimating the degree of true-energy non-conservation and we demonstrate their practical utility on an example of infrared spectra simulations. These criteria are more important and intuitive than simply evaluating the quality of the ML potential energies and forces as is commonly done and can be applied universally, e.g., even for trajectories with unknown or discontinuous potential energy. Such an approach introduces new standards for evaluating MD by focusing on the true-energy conservation and can help in developing more accurate methods for simulating molecular and materials properties.

Shadow energy functionals and potentials in Born-Oppenheimer molecular dynamics

In Born-Oppenheimer molecular dynamics (BOMD) simulations based on the density functional theory (DFT), the potential energy and the interatomic forces are calculated from an electronic ground state density that is determined by an iterative self-consistent field optimization procedure, which, in practice, never is fully converged. The calculated energies and forces are, therefore, only approximate, which may lead to an unphysical energy drift and instabilities. Here, we discuss an alternative shadow BOMD approach that is based on backward error analysis. Instead of calculating approximate solutions for an underlying exact regular Born-Oppenheimer potential, we do the opposite. Instead, we calculate the exact electron density, energies, and forces, but for an underlying approximate shadow Born-Oppenheimer potential energy surface. In this way, the calculated forces are conservative with respect to the approximate shadow potential and generate accurate molecular trajectories with long-term energy stabilities. We show how such shadow Born-Oppenheimer potentials can be constructed at different levels of accuracy as a function of the integration time step, δt, from the constrained minimization of a sequence of systematically improvable, but approximate, shadow energy density functionals. For each energy functional, there is a corresponding ground state Born-Oppenheimer potential. These pairs of shadow energy functionals and potentials are higher-level generalizations of the original "zeroth-level" shadow energy functionals and potentials used in extended Lagrangian BOMD [Niklasson, Eur. Phys. J. B 94, 164 (2021)]. The proposed shadow energy functionals and potentials are useful only within this extended dynamical framework, where also the electronic degrees of freedom are propagated as dynamical field variables together with the atomic positions and velocities. The theory is quite general and can be applied to MD simulations using approximate DFT, Hartree-Fock, or semi-empirical methods, as well as to coarse-grained flexible charge models.

Graph-based quantum response theory and shadow Born-Oppenheimer molecular dynamics

Graph-based linear scaling electronic structure theory for quantum-mechanical molecular dynamics simulations [A. M. N. Niklasson et al., J. Chem. Phys. 144, 234101 (2016)] is adapted to the most recent shadow potential formulations of extended Lagrangian Born-Oppenheimer molecular dynamics, including fractional molecular-orbital occupation numbers [A. M. N. Niklasson, J. Chem. Phys. 152, 104103 (2020) and A. M. N. Niklasson, Eur. Phys. J. B 94, 164 (2021)], which enables stable simulations of sensitive complex chemical systems with unsteady charge solutions. The proposed formulation includes a preconditioned Krylov subspace approximation for the integration of the extended electronic degrees of freedom, which requires quantum response calculations for electronic states with fractional occupation numbers. For the response calculations, we introduce a graph-based canonical quantum perturbation theory that can be performed with the same natural parallelism and linear scaling complexity as the graph-based electronic structure calculations for the unperturbed ground state. The proposed techniques are particularly well-suited for semi-empirical electronic structure theory, and the methods are demonstrated using self-consistent charge density-functional tight-binding theory both for the acceleration of self-consistent field calculations and for quantum-mechanical molecular dynamics simulations. Graph-based techniques combined with the semi-empirical theory enable stable simulations of large, complex chemical systems, including tens-of-thousands of atoms.

Field-Pulse-Induced Annealing of 2D Colloidal Polycrystals

Two-dimensional colloidal crystals are of considerable fundamental and practical importance. However, their quality is often low due to the widespread presence of domain walls and defects. In this work, we explored the annealing process undergone by monolayers of superparamagnetic colloids adsorbed onto fluid interfaces in the presence of magnetic field pulses. These systems present the extraordinary peculiarity that both the extent and the character of interparticle interactions can be adjusted at will by simply varying the strength and orientation of the applied field so that the application of field pulses results in a sudden input of energy. Specifically, we have studied the effect of polycrystal size, pulse duration, slope and frequency on the efficiency of the annealing process and found that (i) this strategy is only effective when the polycrystal consists of less than approximately 10 domains; (ii) that the pulse duration should be of the order of magnitude of the time required for the outer particles to travel one diameter during the heating step; (iii) that the quality of larger polycrystals can be slightly improved by applying tilted pulses. The experimental results were corroborated by Brownian dynamics simulations.

Substrate spectrum of PPM1D in the cellular response to DNA double-strand breaks

PPM1D is a p53-regulated protein phosphatase that modulates the DNA damage response (DDR) and is frequently altered in cancer. Here, we employed chemical inhibition of PPM1D and quantitative mass spectrometry-based phosphoproteomics to identify the substrates of PPM1D upon induction of DNA double-strand breaks (DSBs) by etoposide. We identified 73 putative PPM1D substrates that are involved in DNA repair, regulation of transcription, and RNA processing. One-third of DSB-induced S/TQ phosphorylation sites are dephosphorylated by PPM1D, demonstrating that PPM1D only partially counteracts ATM/ATR/DNA-PK signaling. PPM1D-targeted phosphorylation sites are found in a specific amino acid sequence motif that is characterized by glutamic acid residues, high intrinsic disorder, and poor evolutionary conservation. We identified a functionally uncharacterized protein Kanadaptin as ATM and PPM1D substrate upon DSB induction. We propose that PPM1D plays a role during the response to DSBs by regulating the phosphorylation of DNA- and RNA-binding proteins in intrinsically disordered regions.

•

MS-based phosphoproteomic profiling of PPM1D substrates in U2OS and HCT116 cells

•

PPM1D counteracts ATM in the cellular response to DNA double-strand breaks

•

PPM1D target sites localize to glutamic acid-rich regions with high intrinsic disorder

•

Kanadaptin is a putative DNA damage response factor regulated by ATM and PPM1D

Optimization of the Ewald method for calculating Coulomb interactions in molecular simulations

Practical implementations of the Ewald method used to compute Coulomb interactions in molecular dynamics simulations are hampered by the requirement to truncate its reciprocal space series. It is shown that this can be mitigated by representing the contributions from the neglected reciprocal lattice vector terms as a simple modification of the real space expression in which the real and reciprocal space series have slightly different charge spreading parameters. This procedure, called the α′ method, enables significantly fewer reciprocal lattice vectors to be taken than is currently typical for Ewald, with negligible additional computational cost, which is validated on model systems representing different classes of charged system, a CsI crystal and melt, water, and a room temperature ionic liquid. A procedure for computing accurate energies and forces based on a periodic sampling of an additional number of reciprocal lattice vectors is also proposed and validated by the simulations. The convergence characteristics of expressions for the pressure based on the forces and the potential energy are compared, which is a useful assessment of the accuracy of the simulations in reproducing the Coulomb interaction. The techniques developed in this work can reduce significantly the total computer simulation times for medium sized charged systems, by factors of up to ∼5 for those in the classes studied here.

Anticancer potential of β-sitosterol and oleanolic acid as through inhibition of human estrogenic 17beta-hydroxysteroid dehydrogenase type-1 based on an in silico approach

The human estrogenic enzyme 17beta-hydroxysteroid dehydrogenase type-1 (HSD17B1) provides biosynthesis regulation of active estrogen in stimulating the development of breast cancer through cell proliferation. The β-sitosterol is classified as a steroid compound and is actually a type of triterpenoid compound that has a similar structure to a steroid. This similarity provides a great opportunity for the inhibitor candidate to bind to the HDS17B1 enzyme because of the template similarity on the active site. Several in silico approaches have been applied in this study to examine the potential of these two inhibitor candidates. Pharmacokinetic studies showed positive results by meeting several drug candidate criteria, such as drug-likeness, bioavailability, and ADMET properties. A combination of molecular docking and MD simulation showed good conformational interaction of the inhibitors and HSD17B1. Prediction of binding free energy (ΔG bind) using the Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) approach shows ΔG bind (kcal mol-1) of C1-HSD17B1: -49.31 ± 0.23 and C2-HSD17B1: -33.54 ± 0.34. Meanwhile, decomposition energy analysis (ΔG residue bind) suggested several key residues that were also responsible for the interaction with inhibitors, such as C1-HSD17B1 (six residues: Leu96, Leu149, Pro187, Met193, Val225, and Phe226) and C2-HSD17B1 (four residues: Ile14, Gly94, Pro187, and Val188). Hopefully, the obtained results from this research could be considered for the mechanistic inhibition of the HSDS17B1 enzyme at molecular and atomistic levels.

Shadow Hamiltonian in classical NVE molecular dynamics simulations involving Coulomb interactions

Microcanonical ensemble (NVE) Molecular Dynamics (MD) computer simulations are performed with negligible energy drift for systems incorporating Coulomb interactions and complex constraint schemes. In principle, such systems can now be simulated in the NVE ensemble for millisecond time scales, with no requirement for system thermostatting. Numerical tools for assessing drift in MD simulations are outlined, and drift rates of 10^-6 K/μs are demonstrated for molten salts, polar liquids, and room temperature ionic liquids. Such drift rates are six orders of magnitude smaller than those typically quoted in the literature. To achieve this, the standard Ewald method is slightly modified so the first four derivatives of the real space terms go smoothly to zero at the truncation distance, r_c. New methods for determining standard Ewald errors and the new perturbation errors introduced by the smoothing procedure are developed and applied, these taking charge correlation effects explicitly into account. The shadow Hamiltonian, E_s, is shown to be the strictly conserved quantity in these systems, and standard errors in the mean of one part in 10¹⁰ are routinely calculated. Expressions for the shadow Hamiltonian are improved over previous work by accounting for O(h⁴) terms, where h is the MD time step. These improvements are demonstrated by means of extreme out-of-equilibrium simulations. Using the new methodology, the very low diffusion coefficients of room temperature 1-hexyl-3-methyl-imidazolium chloride are determined from long NVE trajectories in which the equations of motion are known to be integrated correctly, with negligible drift.

Neural networks for protein structure and function prediction and dynamic analysis

Hardware and software advancements along with the accumulation of large amounts of data in recent years have together spurred a remarkable growth in the application of neural networks to various scientific fields. Machine learning based on neural networks with multiple (hidden) layers is becoming an extremely powerful approach for analyzing data. With the accumulation of large amounts of protein data such as structural and functional assay data, the effects of such approaches within the field of protein informatics are increasing. Here, we introduce our recent studies based on applications of neural networks for protein structure and function prediction and dynamic analysis involving: (i) inter-residue contact prediction based on a multiple sequence alignment (MSA) of amino acid sequences, (ii) prediction of protein-compound interaction using assay data, and (iii) detection of protein allostery from trajectories of molecular dynamic (MD) simulation.

Shadow Hamiltonian in classical NVE molecular dynamics simulations: A path to long time stability

The shadow energy, E_s, is the conserved quantity in microcanonical ensemble (NVE) molecular dynamics simulations carried out with the position Verlet central-difference algorithm. A new methodology for calculating precise and accurate values of E_s is presented. It is shown for the first time that E_s rather than E is constant during structural changes occurring within a supercooled liquid. It is also explained how to prepare and conduct microsecond range bulk-phase NVE simulations with essentially zero energy drift without the need for thermostating. The drift is analyzed with block averaging and new drift functions of the shadow energy. With such minimal drift, extremely small and accurate standard errors in the mean for quantities like E_s, E, and temperature, T, can be obtained. Values of the standard error for E_s of ≈10^-10 in molecule-based reduced units can be routinely achieved for simulations of 10⁸ time steps. This corresponds to a simulation temperature drift of ≈10^-6 K/μs, six orders of magnitude smaller than generally considered to be acceptable for protein simulations. We also show for the first time how these treatments can be extended with no loss of accuracy to polyatomic systems with both flexible degrees of freedom and arbitrary geometric constraints imposed via the SHAKE algorithm. As a bonus, estimates of simulation-average kinetic and total energies from high order velocity expressions can be obtained to a good approximation from 2nd order velocities and the average mean square force (for polyatomics, this refers to per site, including any constraint forces).

Transforming the Accuracy and Numerical Stability of ReaxFF Reactive Force Fields

Molecular dynamics (MD) simulations provide an important link between theories and experiments. While ab initio methods can be prohibitively costly, the ReaxFF force field has facilitated in silico studies of chemical reactivity in complex, condensed-phase systems. However, the relatively poor energy conservation in ReaxFF MD has either limited the applicability to short time scales, in cases where energy propagation is important, or has required a continuous coupling of the system to a heat bath. In this study, we reveal the root cause of the unsatisfactory energy conservation, and offer a straightforward solution. The new scheme results in orders of magnitude improvement in energy conservation, numerical stability, and accuracy of ReaxFF force fields, compared to the previous state-of-the-art, at no additional cost. We anticipate that these improvements will open new avenues of research for more accurate reactive simulations in complex systems on long time scales.

Autoencoder-Based Detection of Dynamic Allostery Triggered by Ligand Binding Based on Molecular Dynamics

Dynamic allostery on proteins, triggered by regulator binding or chemical modifications, transmits information from the binding site to distant regions, dramatically altering protein function. It is accompanied by subtle changes in side-chain conformations of the protein, indicating that the changes in dynamics, and not rigid or large conformational changes, are essential to understand regulation of protein function. Although a lot of experimental and theoretical studies have been dedicated to investigate this issue, the regulation mechanism of protein function is still being debated. Here, we propose an autoencoder-based method that can detect dynamic allostery. The method is based on the comparison of time fluctuations of protein structures, in the form of distance matrices, obtained from molecular dynamics simulations in ligand-bound and -unbound forms. Our method detected that the changes in dynamics by ligand binding in the PDZ2 domain led to the reorganization of correlative fluctuation motions among residue pairs, which revealed a different view of the correlated motions from the PCA and DCCM. In addition, other correlative motions were also found as a result of the dynamic perturbation from the ligand binding, which may lead to dynamic allostery. This autoencoder-based method would be usefully applied to the signal transduction and mutagenesis systems involved in protein functions and severe diseases.

A Method for Molecular Dynamics on Curved Surfaces

Dynamics simulations of constrained particles can greatly aid in understanding the temporal and spatial evolution of biological processes such as lateral transport along membranes and self-assembly of viruses. Most theoretical efforts in the field of diffusive transport have focused on solving the diffusion equation on curved surfaces, for which it is not tractable to incorporate particle interactions even though these play a crucial role in crowded systems. We show here that it is possible to take such interactions into account by combining standard constraint algorithms with the classical velocity Verlet scheme to perform molecular dynamics simulations of particles constrained to an arbitrarily curved surface. Furthermore, unlike Brownian dynamics schemes in local coordinates, our method is based on Cartesian coordinates, allowing for the reuse of many other standard tools without modifications, including parallelization through domain decomposition. We show that by applying the schemes to the Langevin equation for various surfaces, we obtain confined Brownian motion, which has direct applications to many biological and physical problems. Finally we present two practical examples that highlight the applicability of the method: 1) the influence of crowding and shape on the lateral diffusion of proteins in curved membranes; and 2) the self-assembly of a coarse-grained virus capsid protein model.

Mathematical modeling of blood clot fragmentation during flow-mediated thrombolysis

A microscale mathematical model of blood clot dissolution based on coarse-grained molecular dynamics is presented. In the model, a blood clot is assumed to be an assembly of blood cells interconnected with elastic fibrin bonds, which are cleaved either biochemically (bond degradation) or mechanically (bond overstretching) during flow-mediated thrombolysis. The effect of a thrombolytic agent on biochemical bond degradation was modeled phenomenologically by assuming that the decay rate of an individual bond is a function of the remaining noncleaved bonds in the vicinity of that bond (spatial corrosion) and the relative stretching of the bond (deformational corrosion). The results of simulations indicate that the blood clot dissolution process progresses by a blood-flow-promoted removal of clot fragments, the sizes of which are flow-dependent. These findings are in good agreement with the results of our recent optical-microscopy experimental studies on a model of blood clot dissolution, as well as with clinical observations. The findings of this study may contribute to a better understanding of the clot fragmentation process and may therefore also help in designing new, safer thrombolytic approaches.