New parallel computing algorithm of molecular dynamics for extremely huge scale biological systems

Supercomputing in the biological sciences: Toward Zettascale and Yottascale simulations

Molecular simulations of biological systems tend to be significantly more compute-intensive than those in materials science and astrophysics, due to important contributions of long-range electrostatic forces and large numbers of time steps (>1E9) required. Simulations of biomolecular complexes of microseconds to milliseconds are considered state-of-the-art today. However, these time scales are miniscule in comparison to physiological time scales relevant to molecular machine activity, drug action, and elongation cycles for protein synthesis, RNA synthesis, and DNA synthesis (seconds to days). While an exascale supercomputer has simulated an entire virus for nanoseconds, this supercomputer would need to be 10 billion times faster to simulate that virus for 3 hours of physiological time, demonstrating the insatiable need for computing power. With growing interest in computational drug design from the pharmaceutical sector, the biological sciences are positioned to be an industry driver in computing.

Biomolecular simulations at the exascale: From drug design to organelles and beyond

The rapid advancement in computational power available for research offers to bring not only quantitative improvements, but also qualitative changes in the field of biomolecular simulation. Here, we review the state of biomolecular dynamics simulations at the threshold to exascale resources becoming available. Both developments in parallel and distributed computing will be discussed, providing a perspective on the state of the art of both. A main focus will be on obtaining binding and conformational free energies, with an outlook to macromolecular complexes and (sub)cellular assemblies.

A methodology of quantifying membrane permeability based on returning probability theory and molecular dynamics simulation

We propose a theoretical approach to estimate the permeability coefficients of substrates (permeants) for crossing membranes from donor (D) phase to acceptor (A) phase by means of molecular dynamics (MD) simulation. A fundamental aspect of our approach involves reformulating the returning probability (RP) theory, a rigorous bimolecular reaction theory, to describe permeation phenomena. This reformulation relies on the parallelism between permeation and bimolecular reaction processes. In the present method, the permeability coefficient is represented in terms of the thermodynamic and kinetic quantities for the reactive (R) phase that exists within the inner region of a membrane. One can evaluate these quantities using multiple MD trajectories starting from phase R. We apply the RP theory to the permeation of ethanol and methylamine at different concentrations (infinitely dilute and 1 mol % conditions of permeants). Under the 1 mol% condition, the present method yields a larger permeability coefficient for ethanol (0.12 ± 0.01 cm s-1) than for methylamine (0.069 ± 0.006 cm s-1), while the values of the permeability coefficient are satisfactorily close to those obtained from the brute-force MD simulations (0.18 ± 0.03 and 0.052 ± 0.005 cm s-1 for ethanol and methylamine, respectively). Moreover, upon analyzing the thermodynamic and kinetic contributions to the permeability, we clarify that a higher concentration dependency of permeability for ethanol, as compared to methylamine, arises from the sensitive nature of ethanol's free-energy barrier within the inner region of the membrane against ethanol concentration.

GENESIS 2.1: High-Performance Molecular Dynamics Software for Enhanced Sampling and Free-Energy Calculations for Atomistic, Coarse-Grained, and Quantum Mechanics/Molecular Mechanics Models

GENeralized-Ensemble SImulation System (GENESIS) is a molecular dynamics (MD) software developed to simulate the conformational dynamics of a single biomolecule, as well as molecular interactions in large biomolecular assemblies and between multiple biomolecules in cellular environments. To achieve the latter purpose, the earlier versions of GENESIS emphasized high performance in atomistic MD simulations on massively parallel supercomputers, with or without graphics processing units (GPUs). Here, we implemented multiscale MD simulations that include atomistic, coarse-grained, and hybrid quantum mechanics/molecular mechanics (QM/MM) calculations. They demonstrate high performance and are integrated with enhanced conformational sampling algorithms and free-energy calculations without using external programs except for the QM programs. In this article, we review new functions, molecular models, and other essential features in GENESIS version 2.1 and discuss ongoing developments for future releases.

Conformational changes of ribosomes during translation elongation resolved by molecular dynamics simulations

Molecular dynamics simulations have emerged as a powerful set of tools to unravel the intricate dynamics of ribosomes during protein synthesis. Recent advancements in this field have enabled simulations to delve deep into the conformational rearrangements of ribosomes and associated factors, providing invaluable insights into the intricacies of translation. Emphasis on simulations has recently been on translation elongation, such as tRNA selection, translocation, and ribosomal head-swivel motions. These studies have offered crucial structural interpretations of how genetic information is faithfully translated into proteins. This review outlines recent discoveries concerning ribosome conformational changes occurring during translation elongation, as elucidated through molecular dynamics simulations.

Explainable chemical artificial intelligence from accurate machine learning of real-space chemical descriptors

Machine-learned computational chemistry has led to a paradoxical situation in which molecular properties can be accurately predicted, but they are difficult to interpret. Explainable AI (XAI) tools can be used to analyze complex models, but they are highly dependent on the AI technique and the origin of the reference data. Alternatively, interpretable real-space tools can be employed directly, but they are often expensive to compute. To address this dilemma between explainability and accuracy, we developed SchNet4AIM, a SchNet-based architecture capable of dealing with local one-body (atomic) and two-body (interatomic) descriptors. The performance of SchNet4AIM is tested by predicting a wide collection of real-space quantities ranging from atomic charges and delocalization indices to pairwise interaction energies. The accuracy and speed of SchNet4AIM breaks the bottleneck that has prevented the use of real-space chemical descriptors in complex systems. We show that the group delocalization indices, arising from our physically rigorous atomistic predictions, provide reliable indicators of supramolecular binding events, thus contributing to the development of Explainable Chemical Artificial Intelligence (XCAI) models.

Probing the geometric and electronic structures of the transition metal oxides RhOn-1/0 (n = 1-4) clusters

The photoelectron spectroscopies of RhOn- (n = 1-2) were obtained via using the photoelectron velocity-map imaging (PE-VMI) approach. The experimental values of the adiabatic detachment energy (ADE) and vertical detachment energy (VDE) for RhO- were reported to be 1.58 ± 0.02 eV. The experimental AED and VDE values of RhO2- were reported to be 2.70 ± 0.02 eV and 2.79 ± 0.02 eV, respectively. The vibrational frequencies of RhO- and RhO2- measured from photoelectron spectra (PES) were 817(76) cm-1 and 932(55) cm-1, respectively. Based on the density functional theory (DFT), the RhOn-1/0 (n = 1-4) clusters were investigated. The optimized configurations of corresponding ground states and low-lying clusters were discovered. Meanwhile, the simulated photoelectron spectroscopy (PES) of RhOn- (n = 1-4) and the theoretical ADE and VDE values of RhOn- (n = 1-4) clusters were unveiled to assist future experimental studies of Rhodium oxide clusters. Moreover, the associated molecular orbitals (MOs), natural population analysis (NPA) and bond order analysis have been utilized to investigate the chemical bonding in these groups.

Solvation dynamics on the diffusion timescale elucidated using energy-represented dynamics theory

Photoexcitation of a solute alters the solute-solvent interaction, resulting in the nonequilibrium relaxation of the solvation structure, often called a dynamic Stokes shift or solvation dynamics. Thanks to the local nature of the solute-solvent interaction, the characteristics of the local solvent environment dissolving the solute can be captured by the observation of this process. Recently, we derived the energy-represented Smoluchowski-Vlasov (ERSV) equation, a diffusion equation for molecular liquids, which can be used to analyze the solvation dynamics on the diffusion timescale. This equation expresses the time development for the solvent distribution on the solute-solvent pair interaction energy (energy coordinate). Since the energy coordinate can effectively treat the solvent flexibility in addition to the position and orientation, the ERSV equation can be utilized in various solvent systems. Here, we apply the ERSV equation to the solvation dynamics of 6-propionyl-2-dimethylamino naphthalene (Prodan) in water and different alcohol solvents (methanol, ethanol, and 1-propanol) for clarifying the differences of the relaxation processes among these solvents. Prodan is a solvent-sensitive fluorescent probe and is thus widely utilized for investigating heterogeneous environments. On the long timescale, the ERSV equation satisfactorily reproduces the relaxation time correlation functions obtained from the molecular dynamics (MD) simulations for these solvents. We reveal that the relaxation time coefficient on the diffusion timescale linearly correlates with the inverse of the translational diffusion coefficients for the alcohol solvents because of the Prodan-solvent energy distributions among the alcohols. In the case of water, the time coefficient deviates from the linear relationship for the alcohols due to the difference in the extent of importance of the collective motion between the water and alcohol solvents.

Perspective on Integrative Simulations of Bioenergetic Domains

Bioenergetic processes in cells, such as photosynthesis or respiration, integrate many time and length scales, which makes the simulation of energy conversion with a mere single level of theory impossible. Just like the myriad of experimental techniques required to examine each level of organization, an array of overlapping computational techniques is necessary to model energy conversion. Here, a perspective is presented on recent efforts for modeling bioenergetic phenomena with a focus on molecular dynamics simulations and its variants as a primary method. An overview of the various classical, quantum mechanical, enhanced sampling, coarse-grained, Brownian dynamics, and Monte Carlo methods is presented. Example applications discussed include multiscale simulations of membrane-wide electron transport, rate kinetics of ATP turnover from electrochemical gradients, and finally, integrative modeling of the chromatophore, a photosynthetic pseudo-organelle.

Time dependent dihedral angle oscillations of the spike protein of SARS-CoV-2 reveal favored frequencies of dihedral angle rotations

The spike protein of SARS-CoV-2 is critical to viral infection of human host cells which ultimately results in COVID-19. In this study we analyzed the behavior of dihedral angles (phi and psi) of the wild-type spike protein over time from molecular dynamics and identified that their oscillations are dominated by a few discrete, relatively low frequencies in the 23-63 MHz range with 42.969 MHz being the most prevalent frequency sampled by the oscillations. We thus observed the spike protein to favor certain frequencies more than others. Gaps in the tally of all observed frequencies for low-abundance amino acids also suggests that the frequency components of dihedral angle oscillations may be a function of position in the primary structure since relatively more abundant amino acids lacked gaps. Lastly, certain residues identified in the literature as constituting the inside of a druggable pocket, as well as others identified as allosteric sites, are observed in our data to have distinctive time domain profiles. This motivated us to propose additional residues with similar time domain profiles, which may be of potential interest to the vaccine and drug design communities for further investigation. Thus these findings indicate that there is a particular frequency domain profile for the spike protein hidden within the time domain data and this information, perhaps with the suggested residues, might provide additional insight into therapeutic development strategies for COVID-19 and beyond.

Molecular Dynamics and Other HPC Simulations for Drug Discovery

High performance computing (HPC) is taking an increasingly important place in drug discovery. It makes possible the simulation of complex biochemical systems with high precision in a short time, thanks to the use of sophisticated algorithms. It promotes the advancement of knowledge in fields that are inaccessible or difficult to access through experimentation and it contributes to accelerating the discovery of drugs for unmet medical needs while reducing costs. Herein, we report how computational performance has evolved over the past years, and then we detail three domains where HPC is essential. Molecular dynamics (MD) is commonly used to explore the flexibility of proteins, thus generating a better understanding of different possible approaches to modulate their activity. Modeling and simulation of biopolymer complexes enables the study of protein-protein interactions (PPI) in healthy and disease states, thus helping the identification of targets of pharmacological interest. Virtual screening (VS) also benefits from HPC to predict in a short time, among millions or billions of virtual chemical compounds, the best potential ligands that will be tested in relevant assays to start a rational drug design process.

Molecular simulation for food protein-ligand interactions: A comprehensive review on principles, current applications, and emerging trends

In recent years, investigations on molecular interaction mechanisms between food proteins and ligands have attracted much interest. The interaction mechanisms can supply much useful information for many fields in the food industry, including nutrient delivery, food processing, auxiliary detection, and others. Molecular simulation has offered extraordinary insights into the interaction mechanisms. It can reflect binding conformation, interaction forces, binding affinity, key residues, and other information that physicochemical experiments cannot reveal in a fast and detailed manner. The simulation results have proven to be consistent with the results of physicochemical experiments. Molecular simulation holds great potential for future applications in the field of food protein-ligand interactions. This review elaborates on the principles of molecular docking and molecular dynamics simulation. Besides, their applications in food protein-ligand interactions are summarized. Furthermore, challenges, perspectives, and trends in molecular simulation of food protein-ligand interactions are proposed. Based on the results of molecular simulation, the mechanisms of interfacial behavior, enzyme-substrate binding, and structural changes during food processing can be reflected, and strategies for hazardous substance detection and food flavor adjustment can be generated. Moreover, molecular simulation can accelerate food development and reduce animal experiments. However, there are still several challenges to applying molecular simulation to food protein-ligand interaction research. The future trends will be a combination of international cooperation and data sharing, quantum mechanics/molecular mechanics, advanced computational techniques, and machine learning, which contribute to promoting food protein-ligand interaction simulation. Overall, the use of molecular simulation to study food protein-ligand interactions has a promising prospect.

Accelerated preprocessing of large numbers of brain images by parallel computing on supercomputers

"Preprocessing" is the first step required in brain image analysis that improves the overall quality and reliability of the results. However, it is computationally demanding and time-consuming, particularly to handle and parcellate complicatedly folded cortical ribbons of the human brain. In this study, we aimed to shorten the analysis time for data preprocessing of 1410 brain images simultaneously on one of the world's highest-performing supercomputers, "Fugaku." The FreeSurfer was used as a benchmark preprocessing software for cortical surface reconstruction. All the brain images were processed simultaneously and successfully analyzed in a calculation time of 17.33 h. This result indicates that using a supercomputer for brain image preprocessing allows big data analysis to be completed shortly and flexibly, thus suggesting the possibility of supercomputers being used for expanding large data analysis and parameter optimization of preprocessing in the future.

Elucidating protein-ligand binding kinetics based on returning probability theory

The returning probability (RP) theory, a rigorous diffusion-influenced reaction theory, enables us to analyze the binding process systematically in terms of thermodynamics and kinetics using molecular dynamics (MD) simulations. Recently, the theory was extended to atomistically describe binding processes by adopting the host-guest interaction energy as the reaction coordinate. The binding rate constants can be estimated by computing the thermodynamic and kinetic properties of the reactive state existing in the binding processes. Here, we propose a methodology based on the RP theory in conjunction with the energy representation theory of solution, applicable to complex binding phenomena, such as protein-ligand binding. The derived scheme of calculating the equilibrium constant between the reactive and dissociate states, required in the RP theory, can be used for arbitrary types of reactive states. We apply the present method to the bindings of small fragment molecules [4-hydroxy-2-butanone (BUT) and methyl methylthiomethyl sulphoxide (DSS)] to FK506 binding protein (FKBP) in an aqueous solution. Estimated binding rate constants are consistent with those obtained from long-timescale MD simulations. Furthermore, by decomposing the rate constants to the thermodynamic and kinetic contributions, we clarify that the higher thermodynamic stability of the reactive state for DSS causes the faster binding kinetics compared with BUT.

Toward a structural identification of metastable molecular conformations

Interpreting high-dimensional data from molecular dynamics simulations is a persistent challenge. In this paper, we show that for a small peptide, deca-alanine, metastable states can be identified through a neural net based on structural information alone. While processing molecular dynamics data, dimensionality reduction is a necessary step that projects high-dimensional data onto a low-dimensional representation that, ideally, captures the conformational changes in the underlying data. Conventional methods make use of the temporal information contained in trajectories generated through integrating the equations of motion, which forgoes more efficient sampling schemes. We demonstrate that EncoderMap, an autoencoder architecture with an additional distance metric, can find a suitable low-dimensional representation to identify long-lived molecular conformations using exclusively structural information. For deca-alanine, which exhibits several helix-forming pathways, we show that this approach allows us to combine simulations with different biasing forces and yields representations comparable in quality to other established methods. Our results contribute to computational strategies for the rapid automatic exploration of the configuration space of peptides and proteins.

Connected Component Analysis of Dynamical Perturbation Contact Networks

Describing protein dynamical networks through amino acid contacts is a powerful way to analyze complex biomolecular systems. However, due to the size of the systems, identifying the relevant features of protein-weighted graphs can be a difficult task. To address this issue, we present the connected component analysis (CCA) approach that allows for fast, robust, and unbiased analysis of dynamical perturbation contact networks (DPCNs). We first illustrate the CCA method as applied to a prototypical allosteric enzyme, the imidazoleglycerol phosphate synthase (IGPS) enzyme from Thermotoga maritima bacteria. This approach was shown to outperform the clustering methods applied to DPCNs, which could not capture the propagation of the allosteric signal within the protein graph. On the other hand, CCA reduced the DPCN size, providing connected components that nicely describe the allosteric propagation of the signal from the effector to the active sites of the protein. By applying the CCA to the IGPS enzyme in different conditions, i.e., at high temperature and from another organism (yeast IGPS), and to a different enzyme, i.e., a protein kinase, we demonstrated how CCA of DPCNs is an effective and transferable tool that facilitates the analysis of protein-weighted networks.

Structural effects of spike protein D614G mutation in SARS-CoV-2

A single mutation from aspartate to glycine at position 614 has dominated all circulating variants of the severe acute respiratory syndrome coronavirus 2. D614G mutation induces structural changes in the spike (S) protein that strengthen the virus infectivity. Here, we use molecular dynamics simulations to dissect the effects of mutation and 630-loop rigidification on S-protein structure. The introduction of the mutation orders the 630-loop structure and thereby induces global structural changes toward the cryoelectron microscopy structure of the D614G S-protein. The ordered 630-loop weakens local interactions between the 614th residue and others in contrast to disordered structures in the wild-type protein. The mutation allosterically alters global interactions between receptor-binding domains, forming an asymmetric and mobile down conformation and facilitating transitions toward up conformation. The loss of salt bridge between D614 and K854 upon the mutation generally stabilizes S-protein protomer, including the fusion peptide proximal region that mediates membrane fusion. Understanding the molecular basis of D614G mutation is crucial as it dominates in all variants of concern, including Delta and Omicron.

Analyzing the Geometry and Dynamics of Viral Structures: A Review of Computational Approaches Based on Alpha Shape Theory, Normal Mode Analysis, and Poisson-Boltzmann Theories

The current SARS-CoV-2 pandemic highlights our fragility when we are exposed to emergent viruses either directly or through zoonotic diseases. Fortunately, our knowledge of the biology of those viruses is improving. In particular, we have more and more structural information on virions, i.e., the infective form of a virus that includes its genomic material and surrounding protective capsid, and on their gene products. It is important to have methods that enable the analyses of structural information on such large macromolecular systems. We review some of those methods in this paper. We focus on understanding the geometry of virions and viral structural proteins, their dynamics, and their energetics, with the ambition that this understanding can help design antiviral agents. We discuss those methods in light of the specificities of those structures, mainly that they are huge. We focus on three of our own methods based on the alpha shape theory for computing geometry, normal mode analyses to study dynamics, and modified Poisson-Boltzmann theories to study the organization of ions and co-solvent and solvent molecules around biomacromolecules. The corresponding software has computing times that are compatible with the use of regular desktop computers. We show examples of their applications on some outer shells and structural proteins of the West Nile Virus.

An exa-scale high-performance molecular dynamics simulation program: MODYLAS

A new version of the highly parallelized general-purpose molecular dynamics (MD) simulation program MODYLAS with high performance on the Fugaku computer was developed. A benchmark test using Fugaku indicated highly efficient communication, single instruction, multiple data (SIMD) processing, and on-cache arithmetic operations. The system's performance deteriorated only slightly, even under high parallelization. In particular, a newly developed minimum transferred data method, requiring a significantly lower amount of data transfer compared to conventional communications, showed significantly high performance. The coordinates and forces of 101 810 176 atoms and the multipole coefficients of the subcells could be distributed to the 32 768 nodes (1 572 864 cores) in 2.3 ms during one MD step calculation. The SIMD effective instruction rates for floating-point arithmetic operations in direct force and fast multipole method (FMM) calculations measured on Fugaku were 78.7% and 31.5%, respectively. The development of a data reuse algorithm enhanced the on-cache processing; the cache miss rate for direct force and FMM calculations was only 2.74% and 1.43%, respectively, on the L1 cache and 0.08% and 0.60%, respectively, on the L2 cache. The modified MODYLAS could complete one MD single time-step calculation within 8.5 ms for the aforementioned large system. Additionally, the program contains numerous functions for material research that enable free energy calculations, along with the generation of various ensembles and molecular constraints.

Acceleration of generalized replica exchange with solute tempering simulations of large biological systems on massively parallel supercomputer

Generalized replica exchange with solute tempering (gREST) is one of the enhanced sampling algorithms for proteins or other systems with rugged energy landscapes. Unlike the replica-exchange molecular dynamics (REMD) method, solvent temperatures are the same in all replicas, while solute temperatures are different and are exchanged frequently between replicas for exploring various solute structures. Here, we apply the gREST scheme to large biological systems containing over one million atoms using a large number of processors in a supercomputer. First, communication time on a multi-dimensional torus network is reduced by matching each replica to MPI processors optimally. This is applicable not only to gREST but also to other multi-copy algorithms. Second, energy evaluations, which are necessary for the multistate bennet acceptance ratio (MBAR) method for free energy estimations, are performed on-the-fly during the gREST simulations. Using these two advanced schemes, we observed 57.72 ns/day performance in 128-replica gREST calculations with 1.5 million atoms system using 16,384 nodes in Fugaku. These schemes implemented in the latest version of GENESIS software could open new possibilities to answer unresolved questions on large biomolecular complex systems with slow conformational dynamics.

Computing the Volume, Surface Area, Mean, and Gaussian Curvatures of Molecules and Their Derivatives

Geometry is crucial in our efforts to comprehend the structures and dynamics of biomolecules. For example, volume, surface area, and integrated mean and Gaussian curvature of the union of balls representing a molecule are used to quantify its interactions with the water surrounding it in the morphometric implicit solvent models. The Alpha Shape theory provides an accurate and reliable method for computing these geometric measures. In this paper, we derive homogeneous formulas for the expressions of these measures and their derivatives with respect to the atomic coordinates, and we provide algorithms that implement them into a new software package, AlphaMol. The only variables in these formulas are the interatomic distances, making them insensitive to translations and rotations. AlphaMol includes a sequential algorithm and a parallel algorithm. In the parallel version, we partition the atoms of the molecule of interest into 3D rectangular blocks, using a kd-tree algorithm. We then apply the sequential algorithm of AlphaMol to each block, augmented by a buffer zone to account for atoms whose ball representations may partially cover the block. The current parallel version of AlphaMol leads to a 20-fold speed-up compared to an independent serial implementation when using 32 processors. For instance, it takes 31 s to compute the geometric measures and derivatives of each atom in a viral capsid with more than 26 million atoms on 32 Intel processors running at 2.7 GHz. The presence of the buffer zones, however, leads to redundant computations, which ultimately limit the impact of using multiple processors. AlphaMol is available as an OpenSource software.

Diffusion theory of molecular liquids in the energy representation and application to solvation dynamics

The generalized Langevin equation (GLE) formalism is a useful theoretical fundament for analyzing dynamical phenomena rigorously. Despite the systematic formulation of dynamics theories with practical approximations, however, the applicability of GLE-based methods is still limited to simple polyatomic liquids due to the approximate treatment of molecular orientations involved in the static molecular liquid theory. Here, we propose an exact framework of dynamics based on the GLE formalism incorporating the energy representation theory of solution, an alternative static molecular liquid theory. A fundamental idea is the projection of the relative positions and orientations of solvents around a solute onto the solute-solvent interaction, namely the energy coordinate, enabling us to describe the dynamics on a one-dimensional coordinate. Introducing systematic approximations, such as the overdamped limit, leads to the molecular diffusion equation in the energy representation that is described in terms of the distribution function of solvents on the energy coordinate and the diffusion coefficients. The present theory is applied to the solvation dynamics triggered by the photoexcitation of benzonitrile. The long-time behavior of the solvation time correlation function is in good agreement with that obtained by the molecular dynamics simulation.

Novel Calcium-Binding Ablating Mutations Induce Constitutive RET Activity and Drive Tumorigenesis

Distinguishing oncogenic mutations from variants of unknown significance (VUS) is critical for precision cancer medicine. Here, computational modeling of 71,756 RET variants for positive selection together with functional assays of 110 representative variants identified a three-dimensional cluster of VUSs carried by multiple human cancers that cause amino acid substitutions in the calmodulin-like motif (CaLM) of RET. Molecular dynamics simulations indicated that CaLM mutations decrease interactions between Ca2+ and its surrounding residues and induce conformational distortion of the RET cysteine-rich domain containing the CaLM. RET-CaLM mutations caused ligand-independent constitutive activation of RET kinase by homodimerization mediated by illegitimate disulfide bond formation. RET-CaLM mutants possessed oncogenic and tumorigenic activities that could be suppressed by tyrosine kinase inhibitors targeting RET. This study identifies calcium-binding ablating mutations as a novel type of oncogenic mutation of RET and indicates that in silico-driven annotation of VUSs of druggable oncogenes is a promising strategy to identify targetable driver mutations.

SIGNIFICANCE

Comprehensive proteogenomic and in silico analyses of a vast number of VUSs identify a novel set of oncogenic and druggable mutations in the well-characterized RET oncogene.

Extending machine learning beyond interatomic potentials for predicting molecular properties

Machine learning (ML) is becoming a method of choice for modelling complex chemical processes and materials. ML provides a surrogate model trained on a reference dataset that can be used to establish a relationship between a molecular structure and its chemical properties. This Review highlights developments in the use of ML to evaluate chemical properties such as partial atomic charges, dipole moments, spin and electron densities, and chemical bonding, as well as to obtain a reduced quantum-mechanical description. We overview several modern neural network architectures, their predictive capabilities, generality and transferability, and illustrate their applicability to various chemical properties. We emphasize that learned molecular representations resemble quantum-mechanical analogues, demonstrating the ability of the models to capture the underlying physics. We also discuss how ML models can describe non-local quantum effects. Finally, we conclude by compiling a list of available ML toolboxes, summarizing the unresolved challenges and presenting an outlook for future development. The observed trends demonstrate that this field is evolving towards physics-based models augmented by ML, which is accompanied by the development of new methods and the rapid growth of user-friendly ML frameworks for chemistry.

Quantitative molecular simulations

All-atom simulations can provide molecular-level insights into the dynamics of gas-phase, condensed-phase and surface processes. One important requirement is a sufficiently realistic and detailed description of the underlying intermolecular interactions. The present perspective provides an overview of the present status of quantitative atomistic simulations from colleagues' and our own efforts for gas- and solution-phase processes and for the dynamics on surfaces. Particular attention is paid to direct comparison with experiment. An outlook discusses present challenges and future extensions to bring such dynamics simulations even closer to reality.

Pre-exascale HPC approaches for molecular dynamics simulations. Covid-19 research: A use case

Exascale computing has been a dream for ages and is close to becoming a reality that will impact how molecular simulations are being performed, as well as the quantity and quality of the information derived for them. We review how the biomolecular simulations field is anticipating these new architectures, making emphasis on recent work from groups in the BioExcel Center of Excellence for High Performance Computing. We exemplified the power of these simulation strategies with the work done by the HPC simulation community to fight Covid-19 pandemics. This article is categorized under:Data Science > Computer Algorithms and ProgrammingData Science > Databases and Expert SystemsMolecular and Statistical Mechanics > Molecular Dynamics and Monte-Carlo Methods.

The inherent flexibility of receptor binding domains in SARS-CoV-2 spike protein

Spike (S) protein is the primary antigenic target for neutralization and vaccine development for the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). It decorates the virus surface and undergoes large motions of its receptor binding domains (RBDs) to enter the host cell. Here, we observe Down, one-Up, one-Open, and two-Up-like structures in enhanced molecular dynamics simulations, and characterize the transition pathways via inter-domain interactions. Transient salt-bridges between RBDA and RBDC and the interaction with glycan at N343B support RBDA motions from Down to one-Up. Reduced interactions between RBDA and RBDB in one-Up induce RBDB motions toward two-Up. The simulations overall agree with cryo-electron microscopy structure distributions and FRET experiments and provide hidden functional structures, namely, intermediates along Down-to-one-Up transition with druggable cryptic pockets as well as one-Open with a maximum exposed RBD. The inherent flexibility of S-protein thus provides essential information for antiviral drug rational design or vaccine development.

Protein assembly and crowding simulations

Proteins encounter frequent molecular interactions in biological environments. Computer simulations have become an increasingly important tool in providing mechanistic insights into how such interactions in vivo relate to their biological function. The review here focuses on simulations describing protein assembly and molecular crowding effects as two important aspects that are distinguished mainly by how specific and long-lived protein contacts are. On the topic of crowding, recent simulations have provided new insights into how crowding affects protein folding and stability, modulates enzyme activity, and affects diffusive properties. Recent studies of assembly processes focus on assembly pathways, especially for virus capsids, amyloid aggregation pathways, and the role of multivalent interactions leading to phase separation. Also, discussed are technical challenges in achieving increasingly realistic simulations of interactions in cellular environments.

Atomistic description of molecular binding processes based on returning probability theory

The efficiency of molecular binding such as host-guest binding is commonly evaluated in terms of kinetics, such as rate coefficients. In general, to compute the coefficient of the overall binding process, we need to consider both the diffusion of reactants and barrier crossing to reach the bound state. Here, we develop a methodology of quantifying the rate coefficient of binding based on molecular dynamics simulation and returning probability (RP) theory proposed by Kim and Lee [J. Chem. Phys. 131, 014503 (2009)]. RP theory provides a tractable formula of the rate coefficient in terms of the thermodynamic stability and kinetics of the intermediate state on a predefined reaction coordinate. In this study, the interaction energy between reactants is utilized as the reaction coordinate, enabling us to effectively describe the reactants' relative position and orientation on one-dimensional space. Application of this method to the host-guest binding systems, which consist of β-cyclodextrin and small guest molecules, yields the rate coefficients consistent with the experimental results.

CATBOSS: Cluster Analysis of Trajectories Based on Segment Splitting

Molecular dynamics (MD) simulations are an exceedingly and increasingly potent tool for molecular behavior prediction and analysis. However, the enormous wealth of data generated by these simulations can be difficult to process and render in a human-readable fashion. Cluster analysis is a commonly used way to partition data into structurally distinct states. We present a method that improves on the state of the art by taking advantage of the temporal information of MD trajectories to enable more accurate clustering at a lower memory cost. To date, cluster analysis of MD simulations has generally treated simulation snapshots as a mere collection of independent data points and attempted to separate them into different clusters based on structural similarity. This new method, cluster analysis of trajectories based on segment splitting (CATBOSS), applies density-peak-based clustering to classify trajectory segments learned by change detection. Applying the method to a synthetic toy model as well as four real-life data sets–trajectories of MD simulations of alanine dipeptide and valine dipeptide as well as two fast-folding proteins–we find CATBOSS to be robust and highly performant, yielding natural-looking cluster boundaries and greatly improving clustering resolution. As the classification of points into segments emphasizes density gaps in the data by grouping them close to the state means, CATBOSS applied to the valine dipeptide system is even able to account for a degree of freedom deliberately omitted from the input data set. We also demonstrate the potential utility of CATBOSS in distinguishing metastable states from transition segments as well as promising application to cases where there is little or no advance knowledge of intrinsic coordinates, making for a highly versatile analysis tool.

Exploring the Minimum-Energy Pathways and Free-Energy Profiles of Enzymatic Reactions with QM/MM Calculations

Understanding molecular mechanisms of enzymatic reactions is of vital importance in biochemistry and biophysics. Here, we introduce new functions of hybrid quantum mechanical/molecular mechanical (QM/MM) calculations in the GENESIS program to compute the minimum-energy pathways (MEPs) and free-energy profiles of enzymatic reactions. For this purpose, an interface in GENESIS is developed to utilize a highly parallel electronic structure program, QSimulate-QM (https://qsimulate.com), calling it as a shared library from GENESIS. Second, algorithms to search the MEP are implemented, combining the string method (E et al. J. Chem. Phys. 2007, 126, 164103) with the energy minimization of the buffer MM region. The method implemented in GENESIS is applied to an enzyme, triosephosphate isomerase, which converts dihyroxyacetone phosphate to glyceraldehyde 3-phosphate in four proton-transfer processes. QM/MM-molecular dynamics simulations show performances of greater than 1 ns/day with the density functional tight binding (DFTB), and 10-30 ps/day with the hybrid density functional theory, B3LYP-D3. These performances allow us to compute not only MEP but also the potential of mean force (PMF) of the enzymatic reactions using the QM/MM calculations. The barrier height obtained as 13 kcal mol-1 with B3LYP-D3 in the QM/MM calculation is in agreement with the experimental results. The impact of conformational sampling in PMF calculations and the level of electronic structure calculations (DFTB vs B3LYP-D3) suggests reliable computational protocols for enzymatic reactions without high computational costs.

A compression strategy for particle mesh Ewald theory

Particle Mesh Ewald (PME) has become a standard method for treating long-range electrostatics in molecular simulations. Although the method has inferior asymptotic computational complexity to its linear scaling competitors, it remains enormously popular due to its high efficiency, which stems from the use of fast Fourier transforms (FFTs). This use of FFTs provides great challenges for scaling the method up to massively parallel systems, in large part because of the need to transfer large amounts of data. In this work, we demonstrate that this data transfer volume can be greatly reduced as a natural consequence of the structure of the PME equations. We also suggest an alternative algorithm that supplants the FFT with a linear algebra approach, which further decreases communication costs at the expense of increased asymptotic computational complexity. This linear algebra based approach is demonstrated to have great potential for latency hiding by interleaving communication and computation steps of the short- and long-range electrostatic terms.