Scaling of Multimillion-Atom Biological Molecular Dynamics Simulation on a Petascale Supercomputer

Large-Scale Molecular Dynamics Simulation Based on Heterogeneous Many-Core Architecture

As the application of molecular dynamics (MD) simulations continues to evolve, the demand for accelerating large-scale simulation systems and handling of enormous simulation tasks is steadily increasing. We propose a parallel acceleration method for large-scale MD simulations based on Sunway heterogeneous many-core processors. This method integrates task scheduling, simulation calculations, and data storage, effectively tackling issues related to large-scale simulations and numerous simulation tasks. The task scheduling strategy flexibly handles tasks on various scales and enables parallel execution of multiple tasks. During the simulation calculations, we ported GROMACS to the Sunway architecture and accelerated the calculation of short-range forces through a heterogeneous processor. Our method achieves approximately 10-fold acceleration and 90% scalability when executing a single simulation task. When handling numerous simulation tasks, our method achieves parallel execution of all of the tasks with 90% scalability. By employing our method, we carried out 50 ns simulations on over 3000 distinct conotoxin structures individually within just 5 h. Additionally, we evaluated more than 200 protein-ligand complexes, and the simulation efficiency significantly exceeded that of midsized to small GPU clusters.

Design of potential anti-cancer agents as COX-2 inhibitors, using 3D-QSAR modeling, molecular docking, oral bioavailability proprieties, and molecular dynamics simulation

Modeling the structural properties of novel morpholine-bearing 1, 5-diaryl-diazole derivatives as potent COX-2 inhibitor, two proposed models based on CoMFA and CoMSIA were evaluated by external and internal validation methods. Partial least squares analysis produced statistically significant models with Q 2 values of 0.668 and 0.652 for CoMFA and CoMSIA, respectively, and also a significant non-validated correlation coefficient R² with values of 0.882 and 0.878 for CoMFA and CoMSIA, respectively. Both models met the requirements of Golbraikh and Tropsha, which means that both models are consistent with all validation techniques. Analysis of the CoMFA and CoMSIA contribution maps and molecular docking revealed that the R1 substituent has a very significant effect on their biological activity. The most active molecules were evaluated for their thermodynamic stability by performing MD simulations for 100 ns; it was revealed that the designed macromolecular ligand complex with 3LN1 protein exhibits a high degree of structural and conformational stability. Based on these results, we predicted newly designed compounds, which have acceptable oral bioavailability properties and would have high synthetic accessibility.

Emerging Trends of Computational Chemistry and Molecular Modeling in Froth Flotation: A Review

Froth flotation is the most versatile process in mineral beneficiation, extensively used to concentrate a wide range of minerals. This process comprises mixtures of more or less liberated minerals, water, air, and various chemical reagents, involving a series of intermingled multiphase physical and chemical phenomena in the aqueous environment. Today's main challenge facing the froth flotation process is to gain atomic-level insights into the properties of its inherent phenomena governing the process performance. While it is often challenging to determine these phenomena via trial-and-error experimentations, molecular modeling approaches not only elicit a deeper understanding of froth flotation but can also assist experimental studies in saving time and budget. Thanks to the rapid development of computer science and advances in high-performance computing (HPC) infrastructures, theoretical/computational chemistry has now matured enough to successfully and gainfully apply to tackle the challenges of complex systems. In mineral processing, however, advanced applications of computational chemistry are increasingly gaining ground and demonstrating merit in addressing these challenges. Accordingly, this contribution aims to encourage mineral scientists, especially those interested in rational reagent design, to become familiarized with the necessary concepts of molecular modeling and to apply similar strategies when studying and tailoring properties at the molecular level. This review also strives to deliver the state-of-the-art integration and application of molecular modeling in froth flotation studies to assist either active researchers in this field to disclose new directions for future research or newcomers to the field to initiate innovative works.

Local Molecular Field Theory for Coulomb Interactions in Aqueous Solutions

Coulomb interactions play a crucial role in a wide array of processes in aqueous solutions but present conceptual and computational challenges to both theory and simulations. We review recent developments in an approach addressing these challenges─local molecular field (LMF) theory. LMF theory exploits an exact and physically suggestive separation of intermolecular Coulomb interactions into strong short-range and uniformly slowly varying long-range components. This allows us to accurately determine the averaged effects of the long-range components on the short-range structure using effective single particle fields and analytical corrections, greatly reducing the need for complex lattice summation techniques used in most standard approaches. The simplest use of these ideas in aqueous solutions leads to the short solvent (SS) model, where both solvent-solvent and solute-solvent Coulomb interactions have only short-range components. Here we use the SS model to give a simple description of pairing of nucleobases and biologically relevant ions in water.

Design of novel anti-cancer drugs targeting TRKs inhibitors based 3D QSAR, molecular docking and molecular dynamics simulation

Tropomyosin receptor kinase (TRK) enzymes are responsible for different types of tumors caused by neurotrophic tyrosine receptor kinase gene fusion and have been identified as an effective target for anticancer therapy. The study of the mechanism between polo-like kinase (PLKs) and pyrazol inhibitors was performed using 3D-QSAR modeling, molecular docking, and MD simulations in order to design high-activity inhibitors. The HQSAR (Q2 = 0.793, R2 = 0.917, R2ext = 0.961), CoMFA (Q2 = 0.582, R2 = 0.722, R2ext = 0.951), CoMSIA/SE (Q2 = 0.603, R2 = 0.801, R2ext = 0.849), and Topomer CoMFA (Q2 = 0.726, R2 = 0.992, R2ext = 0.717) showed good reliability and predictability. All models have been successfully tested by external validation, so all five established models are reliable. The analysis of the different contour maps of different models gives structural information to improve the inhibitory function. Molecular docking results show that the amino acids Met 592, GLU 590, LEU 657, VAL 524, and PHE 589 are the active sites of the tropomyosin receptor TRKs. The results obtained by MD showed that compound 19i could form a more stable complex protein (PDB id: 5KVT). Based on these results, we developed new compounds and their expected inhibitory activities. The results of physicochemical and ADME-Tox properties showed that the four proposed molecules are orally bioavailable, and they are not toxic in the Ames test. Thus, these results would provide modeling information that could help experimental researchers find TRK type I inhibitors more efficiently.Communicated by Ramaswamy H. Sarma.

Non-Ewald methods for evaluating the electrostatic interactions of charge systems: similarity and difference

In molecular simulations, it is essential to properly calculate the electrostatic interactions of particles in the physical system of interest. Here we consider a method called the non-Ewald method, which does not rely on the standard Ewald method with periodic boundary conditions, but instead relies on the cutoff-based techniques. We focus on the physicochemical and mathematical conceptual aspects of the method in order to gain a deeper understanding of the simulation methodology. In particular, we take into account the reaction field (RF) method, the isotropic periodic sum (IPS) method, and the zero-multipole summation method (ZMM). These cutoff-based methods are based on different physical ideas and are completely distinguishable in their underlying concepts. The RF and IPS methods are "additive" methods that incorporate information outside the cutoff region, via dielectric medium and isotropic boundary condition, respectively. In contrast, the ZMM is a "subtraction" method that tries to remove the artificial effects, generated near the boundary, from the cutoff sphere. Nonetheless, we find physical and/or mathematical similarities between these methods. In particular, the modified RF method can be derived by the principle of neutralization utilized in the ZMM, and we also found a direct relationship between IPS and ZMM.

Quantum machine learning for chemistry and physics

Machine learning (ML) has emerged as a formidable force for identifying hidden but pertinent patterns within a given data set with the objective of subsequent generation of automated predictive behavior. In recent years, it is safe to conclude that ML and its close cousin, deep learning (DL), have ushered in unprecedented developments in all areas of physical sciences, especially chemistry. Not only classical variants of ML, even those trainable on near-term quantum hardwares have been developed with promising outcomes. Such algorithms have revolutionized materials design and performance of photovoltaics, electronic structure calculations of ground and excited states of correlated matter, computation of force-fields and potential energy surfaces informing chemical reaction dynamics, reactivity inspired rational strategies of drug designing and even classification of phases of matter with accurate identification of emergent criticality. In this review we shall explicate a subset of such topics and delineate the contributions made by both classical and quantum computing enhanced machine learning algorithms over the past few years. We shall not only present a brief overview of the well-known techniques but also highlight their learning strategies using statistical physical insight. The objective of the review is not only to foster exposition of the aforesaid techniques but also to empower and promote cross-pollination among future research in all areas of chemistry which can benefit from ML and in turn can potentially accelerate the growth of such algorithms.

Influence of Long-Range Forces on the Transition States and Dynamics of NaCl Ion-Pair Dissociation in Water

We study NaCl ion-pair dissociation in a dilute aqueous solution using computer simulations both for the full system with long-range Coulomb interactions and for a well-chosen reference system with short-range intermolecular interactions. Analyzing results using concepts from Local Molecular Field (LMF) theory and the recently proposed AI-based analysis tool "State predictive information bottleneck" (SPIB), we show that the system with short-range interactions can accurately reproduce the transition rate for the dissociation process, the dynamics for moving between the underlying metastable states, and the transition state ensemble. Contributions from long-range interactions can be largely neglected for these processes because long-range forces from the direct interionic Coulomb interactions are almost completely canceled (>90%) by those from solvent interactions over the length scale where the transition takes place. Thus, for this important monovalent ion-pair system, short-range forces alone are able to capture detailed consequences of the collective solvent motion, allowing the use of physically suggestive and computationally efficient short-range models for the dissociation event. We believe that the framework here should be applicable to disentangling mechanisms for more complex processes such as multivalent ion disassociation, where previous work has suggested that long-range contributions may be more important.

Combining Computational Chemistry and Crystallography for a Better Understanding of the Structure of Cellulose

The approaches in this article seek to enhance understanding of cellulose at the molecular level, independent of the source and the particular crystalline form of cellulose. Four main areas of structure research are reviewed. Initially, the molecular shape is inferred from the crystal structures of many small molecules that have β-(1→4) linkages. Then, conformational analyses with potential energy calculations of cellobiose are covered, followed by the use of Atoms-In-Molecules theory to learn about interactions in experimental and theoretical structures. The last section covers models of cellulose nanoparticles. Controversies addressed include the stability of twofold screw-axis conformations, the influence of different computational methods, the predictability of crystalline conformations by studies of isolated molecules, and the twisting of model cellulose crystals.

Machine Learning Force Fields

In recent years, the use of machine learning (ML) in computational chemistry has enabled numerous advances previously out of reach due to the computational complexity of traditional electronic-structure methods. One of the most promising applications is the construction of ML-based force fields (FFs), with the aim to narrow the gap between the accuracy of ab initio methods and the efficiency of classical FFs. The key idea is to learn the statistical relation between chemical structure and potential energy without relying on a preconceived notion of fixed chemical bonds or knowledge about the relevant interactions. Such universal ML approximations are in principle only limited by the quality and quantity of the reference data used to train them. This review gives an overview of applications of ML-FFs and the chemical insights that can be obtained from them. The core concepts underlying ML-FFs are described in detail, and a step-by-step guide for constructing and testing them from scratch is given. The text concludes with a discussion of the challenges that remain to be overcome by the next generation of ML-FFs.

Molecular Insight into the Cosolvent Effect on Lignin-Cellulose Adhesion

Understanding and controlling the physical adsorption of lignin compounds on cellulose pulp are key parameters in the successful optimization of organosolv processes. The effect of binary organic-aqueous solvents on the coordination of lignin to cellulose was studied with molecular dynamics simulations, considering ethanol and acetonitrile to be organic cosolvents in aqueous solutions in comparison to their monocomponent counterparts. The structures of the solvation shells around cellulose and lignin and the energetics of lignin-cellulose adhesion indicate a more effective disruption of lignin-cellulose binding by binary solvents. The synergic effect between solvent components is explained by their preferential interactions with lignin-cellulose complexes. In the presence of pure water, long-lasting H-bonds in the lignin-cellulose complex are observed, promoted by the nonfavorable interactions of lignin with water. Ethanol and acetonitrile compete with water and lignin for cellulose oxygen binding sites, causing a nonlinear decrease in the lignin-cellulose interactions with the amount of the organic component. This effect is modulated by the water exclusion from the cellulose solvation shell by the organic solvent component. The amount and rate of water exclusion depend on the type of organic cosolvent and its concentration.

Heterogeneous parallelization and acceleration of molecular dynamics simulations in GROMACS

The introduction of accelerator devices such as graphics processing units (GPUs) has had profound impact on molecular dynamics simulations and has enabled order-of-magnitude performance advances using commodity hardware. To fully reap these benefits, it has been necessary to reformulate some of the most fundamental algorithms, including the Verlet list, pair searching, and cutoffs. Here, we present the heterogeneous parallelization and acceleration design of molecular dynamics implemented in the GROMACS codebase over the last decade. The setup involves a general cluster-based approach to pair lists and non-bonded pair interactions that utilizes both GPU and central processing unit (CPU) single instruction, multiple data acceleration efficiently, including the ability to load-balance tasks between CPUs and GPUs. The algorithm work efficiency is tuned for each type of hardware, and to use accelerators more efficiently, we introduce dual pair lists with rolling pruning updates. Combined with new direct GPU-GPU communication and GPU integration, this enables excellent performance from single GPU simulations through strong scaling across multiple GPUs and efficient multi-node parallelization.

Local Molecular Field Theory for Nonequilibrium Systems

We provide a framework for extending equilibrium local molecular field (LMF) theory to a statistical ensemble evolving under a time-dependent applied field. In this context, the self-consistency of the original LMF equation is achieved dynamically, which provides an efficient method for computing the equilibrium LMF potential, in addition to providing the nonequilibrium generalization. As a concrete example, we investigate water confined between hydrophobic or charged walls, systems that are very sensitive to the treatment of long-ranged electrostatics. We then analyze confined water in the presence of a time-dependent applied electric field, generated by a sinusoidal or abrupt variation of the magnitudes of uniform charge densities on each wall. Very accurate results are found from the time-dependent LMF formalism even for strong static fields and for time-dependent systems that are driven far from equilibrium where linear response methods fail.

Complex Behavior of Ordered and Icelike Water in Carbon Nanotubes near Its Bulk Boiling Point

We report the results of extensive molecular dynamics (MD) simulation of water in a carbon nanotube (CNT) with a specific diameter over a wide range of temperatures from 343 to 423 K. In order to characterize the nature of water, we have computed the Kirkwood g-factor, the ten Wolde parameter, the radial distribution, the cage correlation, the intermediate scattering functions, the mean-square displacements of the water molecules, and the connectivity of the oxygen atoms. The computed properties provide evidence for complex behavior. Some of the properties indicate an icelike structure, while others point to ordered (but not necessarily frozen) water. The connectivity is close to 9. The ordered water exists both below and above its bulk boiling point. The order is identified based on the ten Wolde parameter and may explain, along with the dynamic slow down, the recent discovery of "ice" in CNTs near the bulk boiling point in a certain range of CNT diameters, not seen in tubes of other sizes.

Alchemical free energy calculations and umbrella sampling with local molecular field theory

Understanding the thermodynamic driving forces underlying any chemical process requires a description of the underlying free energy surface. However, computation of free energies is difficult, often requiring advanced sampling techniques. Moreover, these computations can be further complicated by the evaluation of any long-ranged interactions in the system of interest, such as Coulomb interactions in charged and polar media. Local molecular field theory is a promising approach to avoid many of the conceptual and computational difficulties associated with long-ranged interactions. We present frameworks for performing alchemical free energy calculations and non-Boltzmann sampling with local molecular field theory. We demonstrate that local molecular field theory can be used to perform these free energy calculations with accuracy comparable to traditional methodologies while eliminating the need for explicit treatment of long-ranged interactions in simulations.

Light-induced dilation in nanosheets of charge-transfer complexes

We report the observation of a sizable photostrictive effect of 5.7% with fast, submillisecond response times, arising from a light-induced lattice dilation of a molecular nanosheet, composed of the molecular charge-transfer compound dibenzotetrathiafulvalene (DBTTF) and C60 An interfacial self-assembly approach is introduced for the thickness-controlled growth of the thin films. From photoabsorption measurements, molecular simulations, and electronic structure calculations, we suggest that photostriction within these films arises from a transformation in the molecular structure of constituent molecules upon photoinduced charge transfer, as well as the accommodation of free charge carriers within the material. Additionally, we find that the photostrictive properties of the nanosheets are thickness-dependent, a phenomenon that we suggest arises from surface-induced conformational disorder in the molecular components of the film. Moreover, because of the molecular structure in the films, which results largely from interactions between the constituent π-systems and the sulfur atoms of DBTTF, the optoelectronic properties are found to be anisotropic. This work enables the fabrication of 2D molecular charge-transfer nanosheets with tunable thicknesses and properties, suitable for a wide range of applications in flexible electronic technologies.

Dynamic Neutron Scattering by Biological Systems

Dynamic neutron scattering directly probes motions in biological systems on femtosecond to microsecond timescales. When combined with molecular dynamics simulation and normal mode analysis, detailed descriptions of the forms and frequencies of motions can be derived. We examine vibrations in proteins, the temperature dependence of protein motions, and concepts describing the rich variety of motions detectable using neutrons in biological systems at physiological temperatures. New techniques for deriving information on collective motions using coherent scattering are also reviewed.

Graphics processing units in bioinformatics, computational biology and systems biology

Several studies in Bioinformatics, Computational Biology and Systems Biology rely on the definition of physico-chemical or mathematical models of biological systems at different scales and levels of complexity, ranging from the interaction of atoms in single molecules up to genome-wide interaction networks. Traditional computational methods and software tools developed in these research fields share a common trait: they can be computationally demanding on Central Processing Units (CPUs), therefore limiting their applicability in many circumstances. To overcome this issue, general-purpose Graphics Processing Units (GPUs) are gaining an increasing attention by the scientific community, as they can considerably reduce the running time required by standard CPU-based software, and allow more intensive investigations of biological systems. In this review, we present a collection of GPU tools recently developed to perform computational analyses in life science disciplines, emphasizing the advantages and the drawbacks in the use of these parallel architectures. The complete list of GPU-powered tools here reviewed is available at http://bit.ly/gputools.

Nanoscale Engineering of Designer Cellulosomes

Biocatalysts showcase the upper limit obtainable for high-speed molecular processing and transformation. Efforts to engineer functionality in synthetic nanostructured materials are guided by the increasing knowledge of evolving architectures, which enable controlled molecular motion and precise molecular recognition. The cellulosome is a biological nanomachine, which, as a fundamental component of the plant-digestion machinery from bacterial cells, has a key potential role in the successful development of environmentally-friendly processes to produce biofuels and fine chemicals from the breakdown of biomass waste. Here, the progress toward so-called "designer cellulosomes", which provide an elegant alternative to enzyme cocktails for lignocellulose breakdown, is reviewed. Particular attention is paid to rational design via computational modeling coupled with nanoscale characterization and engineering tools. Remaining challenges and potential routes to industrial application are put forward.

Long-ranged contributions to solvation free energies from theory and short-ranged models

Long-standing problems associated with long-ranged electrostatic interactions have plagued theory and simulation alike. Traditional lattice sum (Ewald-like) treatments of Coulomb interactions add significant overhead to computer simulations and can produce artifacts from spurious interactions between simulation cell images. These subtle issues become particularly apparent when estimating thermodynamic quantities, such as free energies of solvation in charged and polar systems, to which long-ranged Coulomb interactions typically make a large contribution. In this paper, we develop a framework for determining very accurate solvation free energies of systems with long-ranged interactions from models that interact with purely short-ranged potentials. Our approach is generally applicable and can be combined with existing computational and theoretical techniques for estimating solvation thermodynamics. We demonstrate the utility of our approach by examining the hydration thermodynamics of hydrophobic and ionic solutes and the solvation of a large, highly charged colloid that exhibits overcharging, a complex nonlinear electrostatic phenomenon whereby counterions from the solvent effectively overscreen and locally invert the integrated charge of the solvated object.

Efficient Multistate Reactive Molecular Dynamics Approach Based on Short-Range Effective Potentials

Nonbonded interactions between molecules usually include the van der Waals force and computationally expensive long-range electrostatic interactions. This article develops a more efficient approach: the effective-interaction multistate empirical-valence-bond (EI-MS-EVB) model. The EI-MS-EVB method relies on a mapping of all interactions onto a short-range and thus, computationally efficient effective potential. The effective potential is tabulated by matching its force to known trajectories obtained from the full-potential empirical multistate empirical-valence-bond (MS-EVB) model. The effective pairwise interaction depends on and is uniquely determined by the atomic configuration of the system, varying only with respect to the hydrogen-bonding topology. By comparing the EI-MS-EVB and full MS-EVB calculations of several equilibrium and dynamic properties important to hydrated excess proton solvation and transport, we show that the EI-MS-EVB model produces very accurate results for the specific system in which the tabulated potentials were generated. The EI-MS-EVB potential also transfers reasonably well to similar systems with different temperatures and box sizes. The EI-MS-EVB method also reduces the computational cost of the nonbonded interactions by about 1 order of magnitude in comparison with the full algorithm.

Mechanism of lignin inhibition of enzymatic biomass deconstruction

BACKGROUND

The conversion of plant biomass to ethanol via enzymatic cellulose hydrolysis offers a potentially sustainable route to biofuel production. However, the inhibition of enzymatic activity in pretreated biomass by lignin severely limits the efficiency of this process.

RESULTS

By performing atomic-detail molecular dynamics simulation of a biomass model containing cellulose, lignin, and cellulases (TrCel7A), we elucidate detailed lignin inhibition mechanisms. We find that lignin binds preferentially both to the elements of cellulose to which the cellulases also preferentially bind (the hydrophobic faces) and also to the specific residues on the cellulose-binding module of the cellulase that are critical for cellulose binding of TrCel7A (Y466, Y492, and Y493).

CONCLUSIONS

Lignin thus binds exactly where for industrial purposes it is least desired, providing a simple explanation of why hydrolysis yields increase with lignin removal.

The zero-multipole summation method for estimating electrostatic interactions in molecular dynamics: analysis of the accuracy and application to liquid systems

In the preceding paper [I. Fukuda, J. Chem. Phys. 139, 174107 (2013)], the zero-multipole (ZM) summation method was proposed for efficiently evaluating the electrostatic Coulombic interactions of a classical point charge system. The summation takes a simple pairwise form, but prevents the electrically non-neutral multipole states that may artificially be generated by a simple cutoff truncation, which often causes large energetic noises and significant artifacts. The purpose of this paper is to judge the ability of the ZM method by investigating the accuracy, parameter dependencies, and stability in applications to liquid systems. To conduct this, first, the energy-functional error was divided into three terms and each term was analyzed by a theoretical error-bound estimation. This estimation gave us a clear basis of the discussions on the numerical investigations. It also gave a new viewpoint between the excess energy error and the damping effect by the damping parameter. Second, with the aid of these analyses, the ZM method was evaluated based on molecular dynamics (MD) simulations of two fundamental liquid systems, a molten sodium-chlorine ion system and a pure water molecule system. In the ion system, the energy accuracy, compared with the Ewald summation, was better for a larger value of multipole moment l currently induced until l ≲ 3 on average. This accuracy improvement with increasing l is due to the enhancement of the excess-energy accuracy. However, this improvement is wholly effective in the total accuracy if the theoretical moment l is smaller than or equal to a system intrinsic moment L. The simulation results thus indicate L ∼ 3 in this system, and we observed less accuracy in l = 4. We demonstrated the origins of parameter dependencies appearing in the crossing behavior and the oscillations of the energy error curves. With raising the moment l we observed, smaller values of the damping parameter provided more accurate results and smoother behaviors with respect to cutoff length were obtained. These features can be explained, on the basis of the theoretical error analyses, such that the excess energy accuracy is improved with increasing l and that the total accuracy improvement within l ⩽ L is facilitated by a small damping parameter. Although the accuracy was fundamentally similar to the ion system, the bulk water system exhibited distinguishable quantitative behaviors. A smaller damping parameter was effective in all the practical cutoff distance, and this fact can be interpreted by the reduction of the excess subset. A lower moment was advantageous in the energy accuracy, where l = 1 was slightly superior to l = 2 in this system. However, the method with l = 2 (viz., the zero-quadrupole sum) gave accurate results for the radial distribution function. We confirmed the stability in the numerical integration for MD simulations employing the ZM scheme. This result is supported by the sufficient smoothness of the energy function. Along with the smoothness, the pairwise feature and the allowance of the atom-based cutoff mode on the energy formula lead to the exact zero total-force, ensuring the total-momentum conservations for typical MD equations of motion.

Enhanced sampling techniques in molecular dynamics simulations of biological systems

BACKGROUND

Molecular dynamics has emerged as an important research methodology covering systems to the level of millions of atoms. However, insufficient sampling often limits its application. The limitation is due to rough energy landscapes, with many local minima separated by high-energy barriers, which govern the biomolecular motion.

SCOPE OF REVIEW

In the past few decades methods have been developed that address the sampling problem, such as replica-exchange molecular dynamics, metadynamics and simulated annealing. Here we present an overview over theses sampling methods in an attempt to shed light on which should be selected depending on the type of system property studied.

MAJOR CONCLUSIONS

Enhanced sampling methods have been employed for a broad range of biological systems and the choice of a suitable method is connected to biological and physical characteristics of the system, in particular system size. While metadynamics and replica-exchange molecular dynamics are the most adopted sampling methods to study biomolecular dynamics, simulated annealing is well suited to characterize very flexible systems. The use of annealing methods for a long time was restricted to simulation of small proteins; however, a variant of the method, generalized simulated annealing, can be employed at a relatively low computational cost to large macromolecular complexes.

GENERAL SIGNIFICANCE

Molecular dynamics trajectories frequently do not reach all relevant conformational substates, for example those connected with biological function, a problem that can be addressed by employing enhanced sampling algorithms. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

Recent advances in QM/MM free energy calculations using reference potentials

Background

Recent years have seen enormous progress in the development of methods for modeling (bio)molecular systems. This has allowed for the simulation of ever larger and more complex systems. However, as such complexity increases, the requirements needed for these models to be accurate and physically meaningful become more and more difficult to fulfill. The use of simplified models to describe complex biological systems has long been shown to be an effective way to overcome some of the limitations associated with this computational cost in a rational way.

Scope of review

Hybrid QM/MM approaches have rapidly become one of the most popular computational tools for studying chemical reactivity in biomolecular systems. However, the high cost involved in performing high-level QM calculations has limited the applicability of these approaches when calculating free energies of chemical processes. In this review, we present some of the advances in using reference potentials and mean field approximations to accelerate high-level QM/MM calculations. We present illustrative applications of these approaches and discuss challenges and future perspectives for the field.

Major conclusions

The use of physically-based simplifications has shown to effectively reduce the cost of high-level QM/MM calculations. In particular, lower-level reference potentials enable one to reduce the cost of expensive free energy calculations, thus expanding the scope of problems that can be addressed.

General significance

As was already demonstrated 40 years ago, the usage of simplified models still allows one to obtain cutting edge results with substantially reduced computational cost. This article is part of a Special Issue entitled Recent developments of molecular dynamics.

•

We present some of the advances to accelerate high-level QM/MM calculations.

•

Quantitative limitations of low-level methods can be overcome by these approaches.

•

Reference potentials make free energy simulations feasible for large systems.

•

Automated fitting reduces the need of expensive sampling of high-level approaches.

•

Application of reference potentials can be extended to a wide range of processes.

Epitope fluctuations in the human papillomavirus are under dynamic allosteric control: a computational evaluation of a new vaccine design strategy

The dynamic properties of the capsid of the human papillomavirus (HPV) type 16 were examined using classical molecular dynamics simulations. By systematically comparing the structural fluctuations of the capsid protein, a strong dynamic allosteric connection between the epitope containing loops and the h4 helix located more than 50 Å away is identified, which was not recognized thus far. Computer simulations show that restricting the structural fluctuations of the h4 helix is key to rigidifying the epitopes, which is thought to be required for eliciting a proper immune response. The allostery identified in the components of the HPV is nonclassical because the mean structure of the epitope carrying loops remains unchanged, but as a result of allosteric effect the structural fluctuations are altered significantly, which in turn changes the biochemical reactivity profile of the epitopes. Exploiting this novel insight, a new vaccine design strategy is proposed wherein a relatively small virus capsid fragment is deposited on a silica nanoparticle in such a way that the fluctuations of the h4 helix are suppressed. The structural and dynamic properties of the epitope carrying loops on this hybrid nanoparticle match the characteristics of epitopes found on the full virus-like particle precisely, suggesting that these nanoparticles may serve as potent, cost-effective, and safe alternatives to traditionally developed vaccines. The structural and dynamic properties of the hybrid nanoparticle are examined in detail to establish the general concepts of the proposed new design.

Relative free enthalpies for point mutations in two proteins with highly similar sequences but different folds

Enveloping distribution sampling was used to calculate free-enthalpy changes associated with single amino acid mutations for a pair of proteins, GA95 and GB95, that show 95% sequence identity yet fold into topologically different structures. Of the L → A, I → F, and L → Y mutations at positions 20, 30, and 45, respectively, of the 56-residue sequence, the first and the last contribute the most to the free-enthalpy difference between the native and non-native sequence-structure combinations, in agreement with the experimental findings for this protein pair. The individual free-enthalpy changes are almost sequence-independent in the four-strand/one-helix structure, the stable form of GB95, while in the three-helix bundle structure, the stable form of GA95, an interplay between residues 20 and 45 is observed.

Hierarchical Multiscale Modeling of Macromolecules and their Assemblies

Soft materials (e.g., enveloped viruses, liposomes, membranes and supercooled liquids) simultaneously deform or display collective behaviors, while undergoing atomic scale vibrations and collisions. While the multiple space-time character of such systems often makes traditional molecular dynamics simulation impractical, a multiscale approach has been presented that allows for long-time simulation with atomic detail based on the co-evolution of slowly-varying order parameters (OPs) with the quasi-equilibrium probability density of atomic configurations. However, this approach breaks down when the structural change is extreme, or when nearest-neighbor connectivity of atoms is not maintained. In the current study, a self-consistent approach is presented wherein OPs and a reference structure co-evolve slowly to yield long-time simulation for dynamical soft-matter phenomena such as structural transitions and self-assembly. The development begins with the Liouville equation for N classical atoms and an ansatz on the form of the associated N-atom probability density. Multiscale techniques are used to derive Langevin equations for the coupled OP-configurational dynamics. The net result is a set of equations for the coupled stochastic dynamics of the OPs and centers of mass of the subsystems that constitute a soft material body. The theory is based on an all-atom methodology and an interatomic force field, and therefore enables calibration-free simulations of soft matter, such as macromolecular assemblies.

GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit

MOTIVATION

Molecular simulation has historically been a low-throughput technique, but faster computers and increasing amounts of genomic and structural data are changing this by enabling large-scale automated simulation of, for instance, many conformers or mutants of biomolecules with or without a range of ligands. At the same time, advances in performance and scaling now make it possible to model complex biomolecular interaction and function in a manner directly testable by experiment. These applications share a need for fast and efficient software that can be deployed on massive scale in clusters, web servers, distributed computing or cloud resources.

RESULTS

Here, we present a range of new simulation algorithms and features developed during the past 4 years, leading up to the GROMACS 4.5 software package. The software now automatically handles wide classes of biomolecules, such as proteins, nucleic acids and lipids, and comes with all commonly used force fields for these molecules built-in. GROMACS supports several implicit solvent models, as well as new free-energy algorithms, and the software now uses multithreading for efficient parallelization even on low-end systems, including windows-based workstations. Together with hand-tuned assembly kernels and state-of-the-art parallelization, this provides extremely high performance and cost efficiency for high-throughput as well as massively parallel simulations.

AVAILABILITY

GROMACS is an open source and free software available from http://www.gromacs.org.

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Behavior of solvent-exposed hydrophobic groove in the anti-apoptotic Bcl-XL protein: clues for its ability to bind diverse BH3 ligands from MD simulations

Bcl-XL is a member of Bcl-2 family of proteins involved in the regulation of intrinsic pathway of apoptosis. Its overexpression in many human cancers makes it an important target for anti-cancer drugs. Bcl-XL interacts with the BH3 domain of several pro-apoptotic Bcl-2 partners. This helical bundle protein has a pronounced hydrophobic groove which acts as a binding region for the BH3 domains. Eight independent molecular dynamics simulations of the apo/holo forms of Bcl-XL were carried out to investigate the behavior of solvent-exposed hydrophobic groove. The simulations used either a twin-range cut-off or particle mesh Ewald (PME) scheme to treat long-range interactions. Destabilization of the BH3 domain-containing helix H2 was observed in all four twin-range cut-off simulations. Most of the other major helices remained stable. The unwinding of H2 can be related to the ability of Bcl-XL to bind diverse BH3 ligands. The loss of helical character can also be linked to the formation of homo- or hetero-dimers in Bcl-2 proteins. Several experimental studies have suggested that exposure of BH3 domain is a crucial event before they form dimers. Thus unwinding of H2 seems to be functionally very important. The four PME simulations, however, revealed a stable helix H2. It is possible that the H2 unfolding might occur in PME simulations at longer time scales. Hydrophobic residues in the hydrophobic groove are involved in stable interactions among themselves. The solvent accessible surface areas of bulky hydrophobic residues in the groove are significantly buried by the loop LB connecting the helix H2 and subsequent helix. These observations help to understand how the hydrophobic patch in Bcl-XL remains stable in the solvent-exposed state. We suggest that both the destabilization of helix H2 and the conformational heterogeneity of loop LB are important factors for binding of diverse ligands in the hydrophobic groove of Bcl-XL.

Compartmentalization and metabolic channeling for multienzymatic biosynthesis: practical strategies and modeling approaches

: The construction of efficient enzyme complexes for multienzymatic biosynthesis is of increasing interest in order to achieve maximum yield and to minimize the interference due to shortcomings that are typical for straightforward one-pot multienzyme catalysis. These include product or intermediate feedback inhibition, degeneration, and diffusive losses of reaction intermediates, consumption of co-factors, and others. The main mechanisms in nature to tackle these effects in transient or stable protein associations are the formation of metabolic channeling and microcompartments, processes that are desirable also for multienzymatic biosynthesis in vitro. This chapter provides an overview over two main aspects. First, numerous recent strategies for establishing compartmentalized multienzyme associations and constructed synthetic enzyme complexes are reviewed. Second, the computational methods at hand to investigate and optimize such associations systematically, especially with focus on large multienzyme complexes and metabolic channeling, are discussed. Perspectives on future studies of multienzymatic biosynthesis concerning compartmentalization and metabolic channeling are presented.

Force fields for classical molecular dynamics

In this chapter we review the basic features and the principles underlying molecular mechanics force fields commonly used in molecular modeling of biological macromolecules. We start by summarizing the historical background and then describe classical pairwise additive potential energy functions. We introduce the problem of the calculation of nonbonded interactions, of particular importance for charged macromolecules. Different parameterization philosophies are then presented, followed by a section on force field validation. We conclude with a brief overview on future perspectives for the development of classical force fields.

Simple and accurate scheme to compute electrostatic interaction: zero-dipole summation technique for molecular system and application to bulk water

The zero-dipole summation method was extended to general molecular systems, and then applied to molecular dynamics simulations of an isotropic water system. In our previous paper [I. Fukuda, Y. Yonezawa, and H. Nakamura, J. Chem. Phys. 134, 164107 (2011)], for evaluating the electrostatic energy of a classical particle system, we proposed the zero-dipole summation method, which conceptually prevents the nonzero-charge and nonzero-dipole states artificially generated by a simple cutoff truncation. Here, we consider the application of this scheme to molecular systems, as well as some fundamental aspects of general cutoff truncation protocols. Introducing an idea to harmonize the bonding interactions and the electrostatic interactions in the scheme, we develop a specific algorithm. As in the previous study, the resulting energy formula is represented by a simple pairwise function sum, enabling facile applications to high-performance computation. The accuracy of the electrostatic energies calculated by the zero-dipole summation method with the atom-based cutoff was numerically investigated, by comparison with those generated by the Ewald method. We obtained an electrostatic energy error of less than 0.01% at a cutoff length longer than 13 Å for a TIP3P isotropic water system, and the errors were quite small, as compared to those obtained by conventional truncation methods. The static property and the stability in an MD simulation were also satisfactory. In addition, the dielectric constants and the distance-dependent Kirkwood factors were measured, and their coincidences with those calculated by the particle mesh Ewald method were confirmed, although such coincidences are not easily attained by truncation methods. We found that the zero damping-factor gave the best results in a practical cutoff distance region. In fact, in contrast to the zero-charge scheme, the damping effect was insensitive in the zero-charge and zero-dipole scheme, in the molecular system we treated. We discussed the origin of this difference between the two schemes and the dependence of this fact on the physical system. The use of the zero damping-factor will enhance the efficiency of practical computations, since the complementary error function is not employed. In addition, utilizing the zero damping-factor provides freedom from the parameter choice, which is not trivial in the zero-charge scheme, and eliminates the error function term, which corresponds to the time-consuming Fourier part under the periodic boundary conditions.

Multiscale macromolecular simulation: role of evolving ensembles

Multiscale analysis provides an algorithm for the efficient simulation of macromolecular assemblies. This algorithm involves the coevolution of a quasiequilibrium probability density of atomic configurations and the Langevin dynamics of spatial coarse-grained variables denoted order parameters (OPs) characterizing nanoscale system features. In practice, implementation of the probability density involves the generation of constant OP ensembles of atomic configurations. Such ensembles are used to construct thermal forces and diffusion factors that mediate the stochastic OP dynamics. Generation of all-atom ensembles at every Langevin time step is computationally expensive. Here, multiscale computation for macromolecular systems is made more efficient by a method that self-consistently folds in ensembles of all-atom configurations constructed in an earlier step, history, of the Langevin evolution. This procedure accounts for the temporal evolution of these ensembles, accurately providing thermal forces and diffusions. It is shown that efficiency and accuracy of the OP-based simulations is increased via the integration of this historical information. Accuracy improves with the square root of the number of historical timesteps included in the calculation. As a result, CPU usage can be decreased by a factor of 3-8 without loss of accuracy. The algorithm is implemented into our existing force-field based multiscale simulation platform and demonstrated via the structural dynamics of viral capsomers.

Non-Ewald methods: theory and applications to molecular systems

Several non-Ewald methods for calculating electrostatic interactions have recently been developed, such as the Wolf method, the reaction field method, the pre-averaging method, and the zero-dipole summation method, for molecular dynamics simulations of various physical systems, including biomolecular systems. We review the theories of these approaches and their potential applications to molecular simulations, and discuss their relationships.

GPU/CPU Algorithm for Generalized Born/Solvent-Accessible Surface Area Implicit Solvent Calculations

Molecular dynamics methodologies comprise a vital research tool for structural biology. Molecular dynamics has benefited from technological advances in computing, such as multi-core CPUs and graphics processing units (GPUs), but harnessing the full power of hybrid GPU/CPU computers remains difficult. The generalized Born/solvent-accessible surface area implicit solvent model (GB/SA) stands to benefit from hybrid GPU/CPU computers, employing the GPU for the GB calculation and the CPU for the SA calculation. Here, we explore the computational challenges facing GB/SA calculations on hybrid GPU/CPU computers and demonstrate how NAMD, a parallel molecular dynamics program, is able to efficiently utilize GPUs and CPUs simultaneously for fast GB/SA simulations. The hybrid computation principles demonstrated here are generally applicable to parallel applications employing hybrid GPU/CPU calculations.

Carbohydrate force fields

Carbohydrates present a special set of challenges to the generation of force fields. First, the tertiary structures of monosaccharides are complex merely by virtue of their exceptionally high number of chiral centers. In addition, their electronic characteristics lead to molecular geometries and electrostatic landscapes that can be challenging to predict and model. The monosaccharide units can also interconnect in many ways, resulting in a large number of possible oligosaccharides and polysaccharides, both linear and branched. These larger structures contain a number of rotatable bonds, meaning they potentially sample an enormous conformational space. This article briefly reviews the history of carbohydrate force fields, examining and comparing their challenges, forms, philosophies, and development strategies. Then it presents a survey of recent uses of these force fields, noting trends, strengths, deficiencies, and possible directions for future expansion.