Molecular Dynamics Simulations Accelerated by GPU for Biological Macromolecules with a Non-Ewald Scheme for Electrostatic Interactions

Molecular dynamics simulation analysis of structural dynamic cross correlation induced by odorant hydrogen-bonding in mouse eugenol ol- factory receptor

Structural fluctuations and dynamic cross-correlations in the mouse eugenol olfactory receptor (Olfr73) were studied by molecular dynamics (MD) simulation to characterize the dynamic response of the protein upon ligand binding. The initial structure was generated by the artificial intelligence tool AlphaFold2 due to the current lack of experimental data. We focused on the hydrogen (H) bond of the odorant eugenol to Ser113, Asn207, and Tyr260 of the receptor protein, the importance of which has been suggested by previous experimental studies. The H-bond was not observed in docking simulations, but in subsequent MD simulations the H-bond to Ser113 was formed in 2-4 ns. The lifetime of the H-bond was in the range of 1-20 ns. On the trajectory with the most stable (20 ns) H-bond, the structural fluctuation of the α-carbon atoms of the receptor main chain was studied by calculating the root mean square fluctuations, the dynamic cross-correlation map, and the time-dependent dynamic cross-correlation. The analysis suggested a correlation transfer pathway Ser113 → Phe182 → (Leu259 or Tyr260) → Tyr291 induced by the ligand binding with a time scale of 4-6 ns.

A cutoff-based method with charge-distribution-data driven pair potentials for efficiently estimating electrostatic interactions in molecular systems

We introduce a simple cutoff-based method for precise electrostatic energy calculations in the molecular dynamics (MD) simulations of point-particle systems. Our method employs a theoretically derived smooth pair potential function to define electrostatic energy, offering stability and computational efficiency in MD simulations. Instead of imposing specific physical conditions, such as dielectric environments or charge neutrality, we focus on the relationship represented by a single summation formula of charge-weighted pair potentials. This approach allows an accurate energy approximation for each particle, enabling a straightforward error analysis. The resulting particle-dependent pair potential captures the charge distribution information, making it suitable for heterogeneous systems and ensuring an enhanced accuracy through distant information inclusion. Numerical investigations of the Madelung constants of crystalline systems validate the method's accuracy.

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.

Difference of binding modes among three ligands to a receptor mSin3B corresponding to their inhibitory activities

A preceding experiment suggested that a compound, which inhibits binding of the REST/NRSF segment to the cleft of a receptor protein mSin3B, can be a potential drug candidate to ameliorate many neuropathies. We have recently developed an enhanced conformational sampling method, genetic-algorithm-guided multi-dimensional virtual-system-coupled canonical molecular dynamics, and in the present study, applied it to three systems consisting of mSin3B and one of three compounds, sertraline, YN3, and acitretin. Other preceding experiments showed that only sertraline inhibits the binding of REST/NRSF to mSin3B. The current simulation study produced the spatial distribution of the compounds around mSin3B, and showed that sertraline and YN3 bound to the cleft of mSin3B with a high propensity, although acitretin did not. Further analyses of the simulation data indicated that only the sertraline-mSin3B complex produced a hydrophobic core similar to that observed in the molecular interface of the REST/NRSF-mSin3B complex: An aromatic ring of sertraline sunk deeply in the mSin3B's cleft forming a hydrophobic core contacting to hydrophobic amino-acid residues located at the bottom of the cleft. The present study proposes a step to design a compound that inhibits competitively the binding of a ligand to its receptor.

Enhanced Conformational Sampling Method Based on Anomaly Detection Parallel Cascade Selection Molecular Dynamics: ad-PaCS-MD

In the present study, we propose a rare-event sampling method called anomaly detection parallel cascade selection molecular dynamics (ad-PaCS-MD). The original PaCS-MD was designed to generate conformational transition pathways from a given reactant to a product when the latter is known a priori. As an extension of the original method, ad-PaCS-MD has been designed to efficiently search transition pathways from a given reactant without referring to a given product. In ad-PaCS-MD, rarely occurring but essential states (configurations) of proteins for the transitions are identified based on the degrees of an anomaly. In more detail, ad-PaCS-MD adopts an algorithm called an anomaly detection generative adversarial network (anoGAN) as a measure for detecting rarely occurring states to be resampled. Here, the essential configurations with higher degrees of the anomaly are selected with anoGAN and intensively resampled by restarting short-time MD simulations from the selected configurations. By repeating the detections and resampling of configurations with the higher degrees of the anomaly, ad-PaCS-MD automatically and efficiently promotes the rare events and gives a wide range of the free energy landscape by combining with the Markov state model construction. As demonstrations, open-closed transitions of two globular proteins (T4 lysozyme and maltose-binding protein) were promoted with ad-PaCS-MD by referring only to the given starting configurations. In each demonstration, ad-PaCS-MD promoted the large-amplitude open-closed transitions with nanosecond-order simulation times. In conclusion, our demonstrations showed a higher conformational sampling efficiency for ad-PaCS-MD than conventional MD (CMD) because CMD required computational costs of more than microsecond-order simulation times to promote the rare events.

GPU-Accelerated Discovery of Pathogen-Derived Molecular Mimics of a T-Cell Insulin Epitope

The strong links between (Human Leukocyte Antigen) HLA, infection and autoimmunity combine to implicate T-cells as primary triggers of autoimmune disease (AD). T-cell crossreactivity between microbially-derived peptides and self-peptides has been shown to break tolerance and trigger AD in experimental animal models. Detailed examination of the potential for T-cell crossreactivity to trigger human AD will require means of predicting which peptides might be recognised by autoimmune T-cell receptors (TCRs). Recent developments in high throughput sequencing and bioinformatics mean that it is now possible to link individual TCRs to specific pathologies for the first time. Deconvolution of TCR function requires knowledge of TCR specificity. Positional Scanning Combinatorial Peptide Libraries (PS-CPLs) can be used to predict HLA-restriction and define antigenic peptides derived from self and pathogen proteins. In silico search of the known terrestrial proteome with a prediction algorithm that ranks potential antigens in order of recognition likelihood requires complex, large-scale computations over several days that are infeasible on a personal computer. We decreased the time required for peptide searching to under 30 min using multiple blocks on graphics processing units (GPUs). This time-efficient, cost-effective hardware accelerator was used to screen bacterial and fungal human pathogens for peptide sequences predicted to activate a T-cell clone, InsB4, that was isolated from a patient with type 1 diabetes and recognised the insulin B-derived epitope HLVEALYLV in the context of disease-risk allele HLA A*0201. InsB4 was shown to kill HLA A*0201⁺ human insulin producing β-cells demonstrating that T-cells with this specificity might contribute to disease. The GPU-accelerated algorithm and multispecies pathogen proteomic databases were validated to discover pathogen-derived peptide sequences that acted as super-agonists for the InsB4 T-cell clone. Peptide-MHC tetramer binding and surface plasmon resonance were used to confirm that the InsB4 TCR bound to the highest-ranked peptide agonists derived from infectious bacteria and fungi. Adoption of GPU-accelerated prediction of T-cell agonists has the capacity to revolutionise our understanding of AD by identifying potential targets for autoimmune T-cells. This approach has further potential for dissecting T-cell responses to infectious disease and cancer.

Free-energy landscape of molecular interactions between endothelin 1 and human endothelin type B receptor: fly-casting mechanism

The free-energy landscape of interaction between a medium-sized peptide, endothelin 1 (ET1), and its receptor, human endothelin type B receptor (hETB), was computed using multidimensional virtual-system coupled molecular dynamics, which controls the system's motions by introducing multiple reaction coordinates. The hETB embedded in lipid bilayer was immersed in explicit solvent. All molecules were expressed as all-atom models. The resultant free-energy landscape had five ranges with decreasing ET1-hETB distance: completely dissociative, outside-gate, gate, binding pocket, and genuine-bound ranges. In the completely dissociative range, no ET1-hETB interaction appeared. In the outside-gate range, an ET1-hETB attractive interaction was the fly-casting mechanism. In the gate range, the ET1 orientational variety decreased rapidly. In the binding pocket range, ET1 was in a narrow pathway with a steep free-energy slope. In the genuine-bound range, ET1 was in a stable free-energy basin. A G-protein-coupled receptor (GPCR) might capture its ligand from a distant place.

Theoretical equations of state for a charged fluid

In this article, we present a molecular thermodynamic study of a system of N particles contained within a volume V and interacting via a hard-core pair potential with an attractive interaction according to the Wolf model for charged systems. This variable-range potential is characterized by three parameters: the repulsive hard-core diameter σ, the energy-well depth ϵ, and the inverse range α; a fourth parameter of the model is a cut-off distance x_c that depends on α according to the relation x_c = 2/α. Two equations of state (EOSs) are presented and derived from thermodynamic perturbation theory and Monte Carlo (MC) simulation data. The first EOS is given by the standard Zwanzig's high-temperature expansion of the Helmholtz free energy, where the first three perturbation terms a₁, a₂, and a₃ were obtained from MC simulations in the canonical ensemble (NVT) and parameterized as functions of α and the reduced density of particles ρ^* = Nσ³/V. The second EOS was obtained from the discrete perturbation theory applied to a discrete representation of the Wolf potential. Results for pressures, internal energies, and isochoric heat capacities are compared to the MC computer simulation data of the Wolf system, including vapor-liquid coexistence curves, for different values of α. Overall, both EOSs give a very good representation of the thermodynamic properties of the Wolf fluid when 0.3 ≤ α ≤ 1.0 and 0.05 ≤ ρ^* ≤ 0.8. Since the Yukawa fluid can reproduce information of screened ionic interactions, we discuss the equivalence between the Wolf and Yukawa fluids in the context of equivalent systems in liquid theory.

Thermal stability of single-domain antibodies estimated by molecular dynamics simulations

Single‐domain antibodies (sdAbs) function like regular antibodies, however, consist of only one domain. Because of their low molecular weight, sdAbs have advantages with respect to production and delivery to their targets and for applications such as antibody drugs and biosensors. Thus, sdAbs with high thermal stability are required. In this work, we chose seven sdAbs, which have a wide range of melting temperature (T_m) values and known structures. We applied molecular dynamics (MD) simulations to estimate their relative stability and compared them with the experimental data. High‐temperature MD simulations at 400 K and 500 K were executed with simulations at 300 K as a control. The fraction of native atomic contacts, Q, measured for the 400 K simulations showed a fairly good correlation with the T_m values. Interestingly, when the residues were classified by their hydrophobicity and size, the Q values of hydrophilic residues exhibited an even better correlation, suggesting that stabilization is correlated with favorable interactions of hydrophilic residues. Measuring the Q value on a per‐residue level enabled us to identify residues that contribute significantly to the instability and thus demonstrating how our analysis can be used in a mutant case study.

Structural Dynamics of the PET-Degrading Cutinase-like Enzyme from Saccharomonospora viridis AHK190 in Substrate-Bound States Elucidates the Ca2+-Driven Catalytic Cycle

A cutinase-type polyesterase from Saccharomonospora viridis AHK190 (Cut190) has been shown to degrade the inner block of polyethylene terephthalate. A unique feature of Cut190 is that its function and stability are regulated by Ca²⁺ binding. Our previous crystal structure analysis of Cut190S226P showed that one Ca²⁺ binds to the enzyme, which induces large conformational changes in several loop regions to stabilize an open conformation [Miyakawa, T., et al. (2015) Appl. Microbiol. Biotechnol. 99, 4297]. In this study, to analyze the substrate recognition mechanism of Cut190, we determined the crystal structure of the inactive form of a Cut190 mutant, Cut190*S176A, in complex with calcium ions and/or substrates. We found that three calcium ions bind to Cut190*S176A, which is supported by analysis using native mass spectrometry experiments and 3D Reference Interaction Site Model calculations. The complex structures with the two substrates, monoethyl succinate and monoethyl adipate (engaged and open forms), presumably correspond to the pre- and post-reaction states, as the ester bond is close to the active site and pointing outward from the active site, respectively, for the two complexes. Ca²⁺ binding induces the pocket to open, enabling the substrate to access the pocket more easily. Molecular dynamics simulations suggest that a post-reaction state in the engaged form presumably exists between the experimentally observed forms, indicating that the substrate would be cleaved in the engaged form and then requires the enzyme to change to the open form to release the product, a process that Ca²⁺ can greatly accelerate.

Phosphorylation of an intrinsically disordered region of Ets1 shifts a multi-modal interaction ensemble to an auto-inhibitory state

Multi-modal interactions are frequently observed in intrinsically disordered regions (IDRs) of proteins upon binding to their partners. In many cases, post-translational modifications in IDRs are accompanied by coupled folding and binding. From both molecular simulations and biochemical experiments with mutational studies, we show that the IDR including a Ser rich region (SRR) of the transcription factor Ets1, just before the DNA-binding core domain, undergoes multi-modal interactions when the SRR is not phosphorylated. In the phosphorylated state, the SRR forms a few specific complex structures with the Ets1 core, covering the recognition helix in the core and drastically reducing the DNA binding affinities as the auto-inhibitory state. The binding kinetics of mutated Ets1 indicates that aromatic residues in the SRR can be substituted with other hydrophobic residues for the interactions with the Ets1 core.

Molecular dynamics coupled with a virtual system for effective conformational sampling

An enhanced conformational sampling method is proposed: virtual-system coupled canonical molecular dynamics (VcMD). Although VcMD enhances sampling along a reaction coordinate, this method is free from estimation of a canonical distribution function along the reaction coordinate. This method introduces a virtual system that does not necessarily obey a physical law. To enhance sampling the virtual system couples with a molecular system to be studied. Resultant snapshots produce a canonical ensemble. This method was applied to a system consisting of two short peptides in an explicit solvent. Conventional molecular dynamics simulation, which is ten times longer than VcMD, was performed along with adaptive umbrella sampling. Free-energy landscapes computed from the three simulations mutually converged well. The VcMD provided quicker association/dissociation motions of peptides than the conventional molecular dynamics did. The VcMD method is applicable to various complicated systems because of its methodological simplicity. © 2018 Wiley Periodicals, Inc.

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.

Molecular mechanisms of cooperative binding of transcription factors Runx1-CBFβ-Ets1 on the TCRα gene enhancer

Ets1 is an essential transcription factor (TF) for several important physiological processes, including cell proliferation and differentiation. Its recognition of the enhancer region of the TCRα gene is enhanced by the cooperative binding of the Runx1–CBFβ heterodimer, with the cancelation of phosphorylation-dependent autoinhibition. The detailed mechanism of this interesting cooperativity between Ets1 and the Runx1–CBFβ heterodimer is still largely unclear. Here, we investigated the molecular mechanisms of this cooperativity, by using molecular dynamics simulations. Consequently, we detected high flexibility of the loop region between the HI2 and H1 helices of Ets1. Upon Runx1–CBFβ heterodimer binding, this loop transiently adopts various sub-stable conformations in its interactions with the DNA. In addition, a network analysis suggested an allosteric pathway in the molecular assembly and identified some key residues that coincide with previous experimental studies. Our simulations suggest that the cooperative binding of Ets1 and the Runx1–CBFβ heterodimer alters the DNA conformation and induces sub-stable conformations of the HI2–H1 loop of Ets1. This phenomenon increases the flexibility of the regulatory module, including the HI2 helix, and destabilizes the inhibitory form of this module. Thus, we hypothesize that this effect facilitates Ets1–DNA binding and prevents the phosphorylation-dependent DNA binding autoinhibition.

Efficient Conformational Search Based on Structural Dissimilarity Sampling: Applications for Reproducing Structural Transitions of Proteins

Structural Dissimilarity Sampling (SDS) is proposed as an efficient conformational search method to promote structural transitions essential for the biological functions of proteins. In SDS, initial structures are selected based on structural dissimilarity, and conformational resampling is repeated. Conformational resampling is performed as follows: (I) arrangement of initial structures for a diverse distribution at the edge of a conformational subspace and (II) promotion of the structural transitions with multiple short-time molecular dynamics (MD) simulations restarting from the diversely distributed initial structures. Cycles of (I) and (II) are repeated to intensively promote structural transitions because conformational resampling from the initial structures would quickly expand conformational distributions toward unvisited conformational subspaces. As a demonstration, SDS was first applied to maltodextrin binding protein (MBP) in explicit water to reproduce structural transitions between the open and closed states of MBP. Structural transitions of MBP were successfully reproduced with SDS in nanosecond-order simulation times. Starting from both the open and closed forms, SDS successfully reproduced the structural transitions within 25 cycles (a total of 250 ns of simulation time). For reference, a conventional long-time (500 ns) MD simulation under NPT (300 K and 1 bar) starting from the open form failed to reproduce the structural transition. In addition to the open-closed motions of MBP, SDS was applied to folding processes of the fast-folding proteins (chignolin, Trp-cage, and villin) and successfully sampled their native states. To confirm how the selections of initial structures affected conformational sampling efficiency, numbers of base sets for characterizing structural dissimilarity of initial structures were addressed in distinct trials of SDS. The parameter searches showed that the conformational sampling efficiency was relatively insensitive with respect to the numbers of base sets, indicating the robustness of SDS for actual applications.

Variation of free-energy landscape of the p53 C-terminal domain induced by acetylation: Enhanced conformational sampling

The C-terminal domain (CTD) of tumor suppressor protein p53 is an intrinsically disordered region that binds to various partner proteins, where lysine of CTD is acetylated/nonacetylated and histidine neutralized/non-neutralized. Because of the flexibility of the unbound CTD, a free-energy landscape (FEL) is a useful quantity for determining its statistical properties. We conducted enhanced conformational sampling of CTD in the unbound state via virtual system coupled multicanonical molecular dynamics, in which the lysine was acetylated or nonacetylated and histidine was charged or neutralized. The fragments were expressed by an all-atom model and were immersed in an explicit solvent. The acetylation and charge-neutralization varied FEL greatly, which might be convenient to exert a hub property. The acetylation slightly enhanced alpha-helix structures that are more compact than sheet/loop conformations. The charge-neutralization produced hairpins. Additionally, circular dichroism experiments confirmed the computational results. We propose possible binding mechanisms of CTD to partners by investigating FEL. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Revisiting antibody modeling assessment for CDR-H3 loop

The antigen-binding site of antibodies, also known as complementarity-determining region (CDR), has hypervariable sequence properties. In particular, the third CDR loop of the heavy chain, CDR-H3, has such variability in its sequence, length, and conformation that ordinary modeling techniques cannot build a high-quality structure. At Stage 2 of the Second Antibody Modeling Assessment (AMA-II) held in 2013, the model structures of the CDR-H3 loops were submitted by the seven modelers and were critically assessed. After our participation in AMA-II, we rebuilt one of the long CDR-H3 loops with 13 residues (A52 antibody) by a more precise method, using enhanced conformational sampling with the explicit water model, as compared to our previous method employed at AMA-II. The current stable models obtained from the free energy landscape at 300 K include structures similar to the X-ray crystal structures. Those models were not built in our previous work at AMA-II. The current free energy landscape suggested that the CDR-H3 loop structures in the crystal are not stable in solution, but they are stabilized by the crystal packing effect.

myPresto/omegagene: a GPU-accelerated molecular dynamics simulator tailored for enhanced conformational sampling methods with a non-Ewald electrostatic scheme

Molecular dynamics (MD) is a promising computational approach to investigate dynamical behavior of molecular systems at the atomic level. Here, we present a new MD simulation engine named "myPresto/omegagene" that is tailored for enhanced conformational sampling methods with a non-Ewald electrostatic potential scheme. Our enhanced conformational sampling methods, e.g., the virtual-system-coupled multi-canonical MD (V-McMD) method, replace a multi-process parallelized run with multiple independent runs to avoid inter-node communication overhead. In addition, adopting the non-Ewald-based zero-multipole summation method (ZMM) makes it possible to eliminate the Fourier space calculations altogether. The combination of these state-of-the-art techniques realizes efficient and accurate calculations of the conformational ensemble at an equilibrium state. By taking these advantages, myPresto/omegagene is specialized for the single process execution with Graphics Processing Unit (GPU). We performed benchmark simulations for the 20-mer peptide, Trp-cage, with explicit solvent. One of the most thermodynamically stable conformations generated by the V-McMD simulation is very similar to an experimentally solved native conformation. Furthermore, the computation speed is four-times faster than that of our previous simulation engine, myPresto/psygene-G. The new simulator, myPresto/omegagene, is freely available at the following URLs: http://www.protein.osaka-u.ac.jp/rcsfp/pi/omegagene/ and http://presto.protein.osaka-u.ac.jp/myPresto4/.

Molmil: a molecular viewer for the PDB and beyond

We have developed a new platform-independent web-based molecular viewer using JavaScript and WebGL. The molecular viewer, Molmil, has been integrated into several services offered by Protein Data Bank Japan and can be easily extended with new functionality by third party developers. Furthermore, the viewer can be used to load files in various formats from the user’s local hard drive without uploading the data to a server. Molmil is available for all platforms supporting WebGL (e.g. Windows, Linux, iOS, Android) from http://gjbekker.github.io/molmil/. The source code is available at http://github.com/gjbekker/molmil under the LGPLv3 licence.

Elastic properties of dynein motor domain obtained from all-atom molecular dynamics simulations

Dyneins are large microtubule motor proteins that convert ATP energy to mechanical power. High-resolution crystal structures of ADP-bound cytoplasmic dynein have revealed the organization of the motor domain, comprising the AAA(+) ring, the linker, the stalk/strut and the C sequence. Recently, the ADP.vanadate-bound structure, which is similar to the ATP hydrolysis transition state, revealed how the structure of dynein changes upon ATP binding. Although both the ADP- and ATP-bound state structures have been resolved, the dynamic properties at the atomic level remain unclear. In this work, we built two models named 'the ADP model' and 'the ATP model', where ADP and ATP are bound to AAA1 in the AAA(+) ring, respectively, to observe the initial procedure of the structural change from the unprimed to the primed state. We performed 200-ns molecular dynamics simulations for both models and compared their structures and dynamics. The motions of the stalk, consisting of a long coiled coil with a microtubule-binding domain, significantly differed between the two models. The elastic properties of the stalk were analyzed and compared with the experimental results.

mDCC_tools: characterizing multi-modal atomic motions in molecular dynamics trajectories

We previously reported the multi-modal Dynamic Cross Correlation (mDCC) method for analyzing molecular dynamics trajectories. This method quantifies the correlation coefficients of atomic motions with complex multi-modal behaviors by using a Bayesian-based pattern recognition technique that can effectively capture transiently formed, unstable interactions. Here, we present an open source toolkit for performing the mDCC analysis, including pattern recognitions, complex network analyses and visualizations. We include a tutorial document that thoroughly explains how to apply this toolkit for an analysis, using the example trajectory of the 100 ns simulation of an engineered endothelin-1 peptide dimer.

AVAILABILITY AND IMPLEMENTATION

The source code is available for free at http://www.protein.osaka-u.ac.jp/rcsfp/pi/mdcctools/, implemented in C ++ and Python, and supported on Linux.

CONTACT

kota.kasahara@protein.osaka-u.ac.jp

SUPPLEMENTARY INFORMATION

Supplementary data are available at Bioinformatics online.

Enhanced sampling simulations to construct free-energy landscape of protein-partner substrate interaction

Molecular dynamics (MD) simulations using all-atom and explicit solvent models provide valuable information on the detailed behavior of protein-partner substrate binding at the atomic level. As the power of computational resources increase, MD simulations are being used more widely and easily. However, it is still difficult to investigate the thermodynamic properties of protein-partner substrate binding and protein folding with conventional MD simulations. Enhanced sampling methods have been developed to sample conformations that reflect equilibrium conditions in a more efficient manner than conventional MD simulations, thereby allowing the construction of accurate free-energy landscapes. In this review, we discuss these enhanced sampling methods using a series of case-by-case examples. In particular, we review enhanced sampling methods conforming to trivial trajectory parallelization, virtual-system coupled multicanonical MD, and adaptive lambda square dynamics. These methods have been recently developed based on the existing method of multicanonical MD simulation. Their applications are reviewed with an emphasis on describing their practical implementation. In our concluding remarks we explore extensions of the enhanced sampling methods that may allow for even more efficient sampling.

Simple, yet powerful methodologies for conformational sampling of proteins

Several biological functions, such as molecular recognition, enzyme catalysis, signal transduction, allosteric regulation, and protein folding, are strongly related to conformational transitions of proteins. These conformational transitions are generally induced as slow dynamics upon collective motions, including biologically relevant large-amplitude fluctuations of proteins. Although molecular dynamics (MD) simulation has become a powerful tool for extracting conformational transitions of proteins, it might still be difficult to reach time scales of the biological functions because the accessible time scales of MD simulations are far from biological time scales, even if straightforward conventional MD (CMD) simulations using massively parallel computers are employed. Thus, it is desirable to develop efficient methods to achieve canonical ensembles with low computational costs. From this perspective, we review several enhanced conformational sampling techniques of biomolecules developed by us. In our methods, multiple independent short-time MD simulations are employed instead of single straightforward long-time CMD simulations. Our basic strategy is as follows: (i) selection of initial seeds (initial structures) for the conformational sampling in restarting MD simulations. Here, the seeds should be selected as candidates with high potential to transit. (ii) Resampling from the selected seeds by initializing velocities in restarting short-time MD simulations. A cycle of these simple protocols might drastically promote the conformational transitions of biomolecules. (iii) Once reactive trajectories extracted from the cycles of short-time MD simulations are obtained, a free energy profile is evaluated by means of umbrella sampling (US) techniques with the weighted histogram analysis method (WHAM) as a post-processing technique. For the selection of the initial seeds, we proposed four different choices: (1) Parallel CaScade molecular dynamics (PaCS-MD), (2) Fluctuation Flooding Method (FFM), (3) Outlier FLOODing (OFLOOD) method, and (4) TaBoo SeArch (TBSA) method. We demonstrate applications of our methods to several biological systems, such as domain motions of proteins with large-amplitude fluctuations, conformational transitions upon ligand binding, and protein folding/refolding to native structures of proteins. Finally, we show the conformational sampling efficiencies of our methods compared with those by CMD simulations and other previously developed enhanced conformational sampling methods.

Miscellaneous Topics in Computer-Aided Drug Design: Synthetic Accessibility and GPU Computing, and Other Topics

BACKGROUND

Computer-aided drug design is still a state-of-the-art process in medicinal chemistry, and the main topics in this field have been extensively studied and well reviewed. These topics include compound databases, ligand-binding pocket prediction, protein-compound docking, virtual screening, target/off-target prediction, physical property prediction, molecular simulation and pharmacokinetics/pharmacodynamics (PK/PD) prediction. Message and Conclusion: However, there are also a number of secondary or miscellaneous topics that have been less well covered. For example, methods for synthesizing and predicting the synthetic accessibility (SA) of designed compounds are important in practical drug development, and hardware/software resources for performing the computations in computer-aided drug design are crucial. Cloud computing and general purpose graphics processing unit (GPGPU) computing have been used in virtual screening and molecular dynamics simulations. Not surprisingly, there is a growing demand for computer systems which combine these resources. In the present review, we summarize and discuss these various topics of drug design.

Protein-protein docking on hardware accelerators: comparison of GPU and MIC architectures

Background

The hardware accelerators will provide solutions to computationally complex problems in bioinformatics fields. However, the effect of acceleration depends on the nature of the application, thus selection of an appropriate accelerator requires some consideration.

Results

In the present study, we compared the effects of acceleration using graphics processing unit (GPU) and many integrated core (MIC) on the speed of fast Fourier transform (FFT)-based protein-protein docking calculation. The GPU implementation performed the protein-protein docking calculations approximately five times faster than the MIC offload mode implementation. The MIC native mode implementation has the advantage in the implementation costs. However, the performance was worse with larger protein pairs because of memory limitations.

Conclusion

The results suggest that GPU is more suitable than MIC for accelerating FFT-based protein-protein docking applications.

A novel approach of dynamic cross correlation analysis on molecular dynamics simulations and its application to Ets1 dimer-DNA complex

The dynamic cross correlation (DCC) analysis is a popular method for analyzing the trajectories of molecular dynamics (MD) simulations. However, it is difficult to detect correlative motions that appear transiently in only a part of the trajectory, such as atomic contacts between the side-chains of amino acids, which may rapidly flip. In order to capture these multi-modal behaviors of atoms, which often play essential roles, particularly at the interfaces of macromolecules, we have developed the "multi-modal DCC (mDCC)" analysis. The mDCC is an extension of the DCC and it takes advantage of a Bayesian-based pattern recognition technique. We performed MD simulations for molecular systems modeled from the (Ets1)2-DNA complex and analyzed their results with the mDCC method. Ets1 is an essential transcription factor for a variety of physiological processes, such as immunity and cancer development. Although many structural and biochemical studies have so far been performed, its DNA binding properties are still not well characterized. In particular, it is not straightforward to understand the molecular mechanisms how the cooperative binding of two Ets1 molecules facilitates their recognition of Stromelysin-1 gene regulatory elements. A correlation network was constructed among the essential atomic contacts, and the two major pathways by which the two Ets1 molecules communicate were identified. One is a pathway via direct protein-protein interactions and the other is that via the bound DNA intervening two recognition helices. These two pathways intersected at the particular cytosine bases (C110/C11), interacting with the H1, H2, and H3 helices. Furthermore, the mDCC analysis showed that both pathways included the transient interactions at their intermolecular interfaces of Tyr396-C11 and Ala327-Asn380 in multi-modal motions of the amino acid side chains and the nucleotide backbone. Thus, the current mDCC approach is a powerful tool to reveal these complicated behaviors and scrutinize intermolecular communications in a molecular system.

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.

Monte Carlo computer simulation of sedimentation of charged hard spherocylinders

In this article we present a NVT Monte Carlo computer simulation study of sedimentation of an electroneutral mixture of oppositely charged hard spherocylinders (CHSC) with aspect ratio L/σ = 5, where L and σ are the length and diameter of the cylinder and hemispherical caps, respectively, for each particle. This system is an extension of the restricted primitive model for spherical particles, where L/σ = 0, and it is assumed that the ions are immersed in an structureless solvent, i.e., a continuum with dielectric constant D. The system consisted of N = 2000 particles and the Wolf method was implemented to handle the coulombic interactions of the inhomogeneous system. Results are presented for different values of the strength ratio between the gravitational and electrostatic interactions, Γ = (mgσ)/(e(2)/Dσ), where m is the mass per particle, e is the electron's charge and g is the gravitational acceleration value. A semi-infinite simulation cell was used with dimensions Lx ≈ Ly and Lz = 5Lx, where Lx, Ly, and Lz are the box dimensions in Cartesian coordinates, and the gravitational force acts along the z-direction. Sedimentation effects were studied by looking at every layer formed by the CHSC along the gravitational field. By increasing Γ, particles tend to get more packed at each layer and to arrange in local domains with an orientational ordering along two perpendicular axis, a feature not observed in the uncharged system with the same hard-body geometry. This type of arrangement, known as tetratic phase, has been observed in two-dimensional systems of hard-rectangles and rounded hard-squares. In this way, the coupling of gravitational and electric interactions in the CHSC system induces the arrangement of particles in layers, with the formation of quasi-two dimensional tetratic phases near the surface.

Structure of the entire stalk region of the Dynein motor domain

Dyneins are large microtubule-based motor complexes that power a range of cellular processes including the transport of organelles, as well as the beating of cilia and flagella. The motor domain is located within the dynein heavy chain and comprises an N-terminal mechanical linker element, a central ring of six AAA+ modules of which four bind or hydrolyze ATP, and a long stalk extending from the AAA+ring with a microtubule-binding domain (MTBD) at its tip. A crucial mechanism underlying the motile activity of cytoskeletal motor proteins is precise coupling between the ATPase and track-binding activities. In dynein, a stalk region consisting of a long (~15nm) antiparallel coiled coil separates these two activities, which must facilitate communication between them. This communication is mediated by a small degree of helix sliding in the coiled coil. However, no high-resolution structure is available of the entire stalk region including the MTBD. Here, we have reported the structure of the entire stalk region of mouse cytoplasmic dynein in a weak microtubule-binding state, which was determined using X-ray crystallography, and have compared it with the dynein motor domain from Dictyostelium discoideum in a strong microtubule-binding state and with a mouse MTBD with its distal portion of the coiled coil fused to seryl-tRNA synthetase from Thermus thermophilus. Our results strongly support the helix-sliding model based on the complete structure of the dynein stalk with a different form of coiled-coil packing. We also propose a plausible mechanism of helix sliding together with further analysis using molecular dynamics simulations. Our results present the importance of conserved proline residues for an elastic motion of stalk coiled coil and imply the manner of change between high-affinity state and low-affinity state of MTBD.