Publications2017-2021

Fundamental vibrations of the chromophore in the membrane protein bacteriorhodopsin (BR), a protonated Schiff base retinal, have been studied for decades, both by resonance Raman and by infrared (IR) difference spectroscopy. Such studies started comparing vibrational changes between the initial BR state (all-trans retinal) and the K intermediate (13-cis retinal), being later extended to the rest of intermediates. They contributed to our understanding of the proton-pumping mechanism of BR by exploiting the sensitivity of fundamental vibrational transitions of the retinal to its conformation. Here, we report on new bands in the 2,500 to 1,800 cm−1 region of the K-BR difference FT-IR spectrum. We show that the bands between 2,500 and 2,300 cm−1 originate from overtone and combination transitions from C-C stretches of the retinal. We assigned bands below 2,300 cm−1 to the combination of retinal C-C stretches with methyl rocks and with hydrogen-out-of-plane vibrations. Remarkably, experimental C-C overtone bands appeared at roughly twice the wavenumber of their fundamentals, with anharmonic mechanical constants ≤3.5 cm−1, and in some cases of ∼1 cm−1. Comparison of combination and fundamental bands indicates that most of the mechanical coupling constants are also very small. Despite the mechanical quasi-harmonicity of the C-C stretches, the area of their overtone bands was only ∼50 to ∼100 times smaller than of their fundamental bands. We concluded that electrical anharmonicity, the second mechanism giving intensity to overtone bands, must be particularly high for the retinal C-C stretches. We corroborated the assignments of negative bands in the K-BR difference FT-IR spectrum by ab initio anharmonic vibrational calculations of all-trans retinal in BR using a quantum-mechanics/molecular mechanics approach, reproducing reasonably well the small experimental anharmonic and coupling mechanical constants. Yet, and in spite accounting for both mechanical and electrical anharmonicities, the intensity of overtone C-C transitions was underestimated by a factor of 4–20, indicating room for improvement in state-of-the-art anharmonic vibrational calculations. The relatively intense overtone and combination bands of the retinal might open the possibility to detect retinal conformational changes too subtle to significantly affect fundamental transitions but leaving a footprint in overtone and combination transitions.

The spectroscopic features of protonated water species in dilute acid solutions have been long sought after for understanding the microscopic behavior of the proton in water with gas-phase water clusters H+(H2O)n extensively studied as bottom-up model systems. We present a new protocol for the calculation of the infrared (IR) spectra of complex systems, which combines the fragment-based Coupled Cluster method and anharmonic vibrational quasi-degenerate perturbation theory, and demonstrate its accuracy towards the complete and accurate assignment of the IR spectrum of the H+(H2O)21 cluster. The site-specific IR spectral signatures reveal two distinct structures for the internal and surface four-coordinated water molecules, which are ice-like and liquid-like, respectively. The effect of inter-molecular interaction between water molecules is addressed, and the vibrational resonance is found between the O-H stretching fundamental and the bending overtone of the nearest neighboring water molecule. The revelation of the spectral signature of the excess proton offers deeper insight into the nature of charge accommodation in the extended hydrogen-bonding network underpinning this aqueous cluster.

Sarcoplasmic reticulum (SR) Ca2+-ATPase transports two Ca2+ ions from the cytoplasm to the SR lumen against a large concentration gradient. X-ray crystallography has revealed the atomic structures of the protein before and after the dissociation of Ca2+, while biochemical studies have suggested the existence of intermediate states in the transition between E1P⋅ADP⋅2Ca2+ and E2P. Here, we explore the pathway and free energy profile of the transition using atomistic molecular dynamics simulations with the mean-force string method and umbrella sampling. The simulations suggest that a series of structural changes accompany the ordered dissociation of ADP, the A-domain rotation, and the rearrangement of the transmembrane (TM) helices. The luminal gate then opens to release Ca2+ ions toward the SR lumen. Intermediate structures on the pathway are stabilized by transient sidechain interactions between the A- and P-domains. Lipid molecules between TM helices play a key role in the stabilization. Free energy profiles of the transition assuming different protonation states suggest rapid exchanges between Ca2+ ions and protons when the Ca2+ ions are released toward the SR lumen.

An efficient anharmonic vibrational method is developed exploiting the locality of molecular vibration. Vibrational coordinates localized to a group of atoms are employed to divide the potential energy surface (PES) of a system into intra- and inter-group contributions. Then, the vibrational Schrödinger equation is solved based on a PES, in which the inter-group coupling is truncated at the harmonic level while accounting for the intra-group anharmonicity. The method is applied to a pentagonal hydrogen bond network (HBN) composed of internal water molecules and charged residues in a membrane protein, bacteriorhodopsin. The PES is calculated by the quantum mechanics/molecular mechanics (QM/MM) calculation at the level of B3LYP-D3/aug-cc-pVDZ. The infrared (IR) spectrum is computed using a set of coordinates localized to each water molecule and amino acid residue by second-order vibrational quasi-degenerate perturbation theory (VQDPT2). Benchmark calculations show that the proposed method yields the N–D/O–D stretching frequencies with an error of 7 cm–1 at the cost reduced by more than five times. In contrast, the harmonic approximation results in a severe error of 150 cm–1. Furthermore, the size of QM regions is carefully assessed to find that the QM regions should include not only the pentagonal HBN itself but also its HB partners. VQDPT2 calculations starting from transient structures obtained by molecular dynamics simulations have shown that the structural sampling has a significant impact on the calculated IR spectrum. The incorporation of anharmonicity, sufficiently large QM regions, and structural samplings are of essential importance to reproduce the experimental IR spectrum. The computational spectrum paves the way for decoding the IR signal of strong HBNs and helps elucidate their functional roles in biomolecules.

In recent years, molecular dynamics (MD) simulations with larger time steps have been performed with the hydrogen-mass-repartitioning (HMR) scheme, where the mass of each hydrogen atom is increased to reduce high-frequency motion while the mass of a non-hydrogen atom bonded to a hydrogen atom is decreased to keep the total molecular mass unchanged. Here, we optimize the scaling factors in HMR and combine them with previously developed accurate temperature/pressure evaluations. The heterogeneous HMR scaling factors are useful to avoid the structural instability of amino acid residues having a five- or six-membered ring in MD simulations with larger time steps. It also reproduces kinetic properties, namely translational and rotational diffusions, if the HMR scaling factors are applied to only solute molecules. The integration scheme is tested for biological systems that include soluble/membrane proteins and lipid bilayers for about 200 μs MD simulations in total and give consistent results in MD simulations with both a small time step of 2.0 fs and a large, multiple time step integration with time steps of 3.5 fs (for fast motions) and 7.0 fs (for slower motions). We also confirm that the multiple time step integration scheme used in this study provides more accurate energy conservations than the RESPA/C1 and is comparable to the RESPA/C2 in NAMD. In summary, the current integration scheme combining the optimized HMR with accurate temperature/pressure evaluations can provide stable and reliable MD trajectories with a larger time step, which are computationally more than 2-fold efficient compared to the conventional methods.

Structural modeling of proteins from cryo-electron microscopy (cryo-EM) density maps is one of the challenging issues in structural biology. De novo modeling combined with flexible fitting refinement (FFR) has been widely used to build a structure of new proteins. In de novo prediction, artificial conformations containing local structural errors such as chirality errors, cis peptide bonds, and ring penetrations are frequently generated and cannot be easily removed in the subsequent FFR. Moreover, refinement can be significantly suppressed due to the low mobility of atoms inside the protein. To overcome these problems, we propose an efficient scheme for FFR, in which the local structural errors are fixed first, followed by FFR using an iterative simulated annealing (SA) molecular dynamics protocol with the united atom (UA) model in an implicit solvent model; we call this scheme “SAUA-FFR”. The best model is selected from multiple flexible fitting runs with various biasing force constants to reduce overfitting. We apply our scheme to the decoys obtained from MAINMAST and demonstrate an improvement of the best model of eight selected proteins in terms of the root-mean-square deviation, MolProbity score, and RWplus score compared to the original scheme of MAINMAST. Fixing the local structural errors can enhance the formation of secondary structures, and the UA model enables progressive refinement compared to the all-atom model owing to its high mobility in the implicit solvent. The SAUA-FFR scheme realizes efficient and accurate protein structure modeling from medium-resolution maps with less overfitting.

The inside of a cell is highly crowded with proteins and other biomolecules. How proteins express their specific functions together with many off-target proteins in crowded cellular environments is largely unknown. Here, we investigate an inhibitor binding with c-Src kinase using atomistic molecular dynamics (MD) simulations in dilute as well as crowded protein solution. The populations of the inhibitor, 4-amino-5-(4-methylphenyl)−7-(t-butyl)pyrazolo[3,4-d]pyrimidine (PP1), in bulk solution and on the surface of c-Src kinase are reduced as the concentration of crowder bovine serum albumins (BSAs) increases. This observation is consistent with the reduced PP1 inhibitor efficacy in experimental c-Src kinase assays in addition with BSAs. The crowded environment changes the major binding pathway of PP1 toward c-Src kinase compared to that in dilute solution. This change is explained based on the population shift mechanism of local conformations near the inhibitor binding site in c-Src kinase.

Enzymes catalyzing [4+2] cycloaddition have attracted increasing attention because of their key roles in natural product biosynthesis. Here, we solved the X-ray crystal structures of a pair of decalin synthases, Fsa2 and Phm7, that catalyze intramolecular [4+2] cycloadditions to form enantiomeric decalin scaffolds during biosynthesis of the HIV-1 integrase inhibitor equisetin and its stereochemical opposite, phomasetin. Computational modeling, using molecular dynamics simulations as well as quantum chemical calculations, demonstrates that the reactions proceed through synergetic conformational constraints assuring transition state-like substrates folds and their stabilization by specific protein-substrate interactions. Site-directed mutagenesis experiments verified the binding models. Intriguingly, the flexibility of bound substrates is largely different in two enzymes, suggesting the distinctive mechanism of dynamics regulation behind these stereoselective reactions. The proposed reaction mechanism herein deepens the basic understanding how these enzymes work but also provides a guiding principle to create artificial enzymes.

Large-scale conformational transitions in multi-domain proteins are often essential for their functions. To investigate the transitions, it is necessary to explore multiple potential pathways, which involve different intermediate structures. Here, we present a multi-basin (MB) coarse-grained (CG) structure-based Go̅ model for describing transitions in proteins with more than two moving domains. This model is an extension of our dual-basin Go̅ model in which system-dependent parameters are determined systematically using the multistate Bennett acceptance ratio method. In the MB Go̅ model for multi-domain proteins, we assume that intermediate structures may have partial inter-domain native contacts. This approach allows us to search multiple transition pathways that involve distinct intermediate structures using the CG molecular dynamics (MD) simulations. We apply this scheme to an enzyme, adenylate kinase (AdK), which has three major domains and can move along two different pathways. Using the optimized mixing parameters for each pathway, AdK shows frequent transitions between the Open, Closed, and the intermediate basins and samples a wide variety of conformations within each basin. The explored multiple transition pathways could be compared with experimental data and examined in more detail by atomistic MD simulations.

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

Conformational changes of proteins upon ligand binding are usually explained in terms of several mechanisms including the induced fit, conformational selection, or their mixtures. Due to the slow time scales, conventional molecular dynamics (cMD) simulations based on the atomistic models cannot easily simulate the open-to-closed conformational transition in proteins. In our previous study, we have developed an enhanced sampling scheme (generalized replica exchange with solute tempering selected surface charged residues: gREST_SSCR) for multidomain proteins and applied it to ligand-mediated conformational changes in the G134R mutant of ribose-binding protein (RBPG134R) in solution. The free-energy landscape (FEL) of RBPG134R in the presence of a ribose at the binding site included the open and closed states and two intermediates, open-like and closed-like forms. Only the open and open-like forms existed in the FEL without a ribose. In the current study, the coupling between the conformational changes and ligand binding is further investigated using coarse-grained MD, multiple atomistic cMD, and free-energy calculations. The ribose is easily dissociated from the binding site of wild-type RBP and RBPG134R in the cMD simulations starting from the open and open-like forms. In contrast, it is stable at the binding site in the simulations from the closed and closed-like forms. The free-energy calculations provide the binding affinities of different structures, supporting the results of cMD simulations. Importantly, cMD simulations from the closed-like structures reveal transitions toward the closed one in the presence of a bound ribose. On the basis of the computational results, we propose a molecular mechanism in which conformational selection and induced fit happen in the first and second halves of the open-to-closed transition in RBP, respectively.

In this paper, we address high performance extreme‐scale molecular dynamics (MD) algorithm in the GENESIS software to perform cellular‐scale molecular dynamics (MD) simulations with more than 100,000 CPU cores. It includes (1) the new algorithm of real‐space nonbonded interactions maximizing the performance on ARM CPU architecture, (2) reciprocal‐space nonbonded interactions minimizing communicational cost, (3) accurate temperature/pressure evaluations that allows a large time step, and (4) effective parallel file inputs/outputs (I/O) for MD simulations of extremely huge systems. The largest system that contains 1.6 billion atoms was simulated using MD with a performance of 8.30 ns/day on Fugaku supercomputer. It extends the available size and time of MD simulations to answer unresolved questions of biomacromolecules in a living cell.

Ongoing COVID-19 pandemic caused by new coronavirus, SARS-CoV-2, calls for urgent developments of vaccines and antiviral drugs. The spike protein of SARS-CoV-2 (S-protein), which consists of trimeric polypeptide chains with glycosylated residues on the surface, triggers the virus entry into a host cell. Extensive structural and functional studies on this protein have rapidly advanced our understanding of the S-protein structure at atomic resolutions, while most of structural studies overlook the effect of glycans attached to S-protein on the conformational stability and functional motions between the inactive Down and the active Up forms. Here, we performed all-atom molecular dynamics (MD) simulations of both Down and Up forms of a fully glycosylated S-protein in solution as well as targeted MD (TMD) simulations between them to elucidate key inter-domain interactions for stabilizing each form and inducing the large-scale conformational transitions. The residue-level interaction analysis of the simulation trajectories detects distinct amino-acid residues and N-glycans as determinants on conformational stability of each form. During the conformational transitions between them, inter-domain interactions mediated by glycosylated residues are switched to play key roles on the stabilization of another form. Electrostatic interactions as well as hydrogen bonds between the three receptor binding domains work as driving forces to initiate the conformational transitions toward the active form. This study sheds light on the mechanisms underlying conformational stability and functional motions of S-protein, which are relevant for vaccines and antiviral drugs developments.

Structure determination using cryo-electron microscopy (cryo-EM) medium-resolution density maps is often facilitated by flexible fitting. Avoiding overfitting, adjusting force constants driving the structure to the density map, and emulating complex conformational transitions are major concerns in the fitting. To address them, we develop a new method based on a three-step multi-scale protocol. First, flexible fitting molecular dynamics (MD) simulations with coarse-grained structure-based force field and replica-exchange scheme between different force constants replicas are performed. Second, fitted Cα atom positions guide the all-atom structure in targeted MD. Finally, the all-atom flexible fitting refinement in implicit solvent adjusts the positions of the side chains in the density map. Final models obtained via the multi-scale protocol are significantly better resolved and more reliable in comparison with long all-atom flexible fitting simulations. The protocol is useful for multi-domain systems with intricate structural transitions as it preserves the secondary structure of single domains.

Short, structured fragments of non-coding mRNA may act as molecular switches upon binding specific ligands, regulating the translation of proteins encoded downstream this mRNA sequence. One switch, called riboswitch N1, is regulated by aminoglycosides such as neomycin. Nucleobase mutations in the apical loop, although distant from the binding pocket, significantly affect neomycin affinity and riboswitch regulatory efficiency. To explain this influence, we conducted molecular dynamics simulations using generalized replica exchange with solute tempering (gREST). Translation assay of a reporter protein in a yeast system shows that mutating A17 to G in the riboswitch apical loop reduces 6-fold the translation regulation efficiency of the mutant. Indeed, simulations of the unbound riboswitch show that G17 frequently stacks with base 7, while base 8 is stabilized towards the binding site in a way that it may interfere with the conformational selection mechanism and decrease riboswitch regulatory activity. In the riboswitch complexes, this single-point A to G mutation disrupts a strong hydrogen bond between nucleotides 5 and 17 and, instead, a new hydrogen bond between residue 17 and neomycin is created. This change forces neomycin to occupy a slightly shifted position in the binding pocket, which increases neomycin flexibility. Our simulations of the U14C mutation suggest that the riboswitch complex with neomycin is more stable if cytosine 14 is protonated. A hydrogen bond between the RNA phosphate and protonated cytosine appears as the stabilizing factor. Also, based on the cell-free translation assay and isothermal titration calorimetry experiments, mutations of nucleotides 14 and 15 affect only slightly the riboswitch ability to bind the ligand and its activity. Indeed, the simulation of the unbound U15A mutant suggests conformations preformed for ligand binding, which may explain slightly higher regulatory activity of this mutant. Overall, our results corroborate the in vivo and in vitro experiments on the N1 riboswitch-neomycin system, detail the relationship between nucleobase mutations and RNA dynamics, and reveal the conformations playing the major role in the conformational selection mechanism.

Conformational transitions in multidomain proteins are essential for biological functions. The Apo conformations are typically open and flexible, while the Holo states form more compact conformations stabilized by protein-ligand interactions. Unfortunately, the atomically detailed mechanisms for such open-closed conformational changes are difficult to be accessed experimentally as well as computationally. To simulate the transitions using atomistic molecular dynamics (MD) simulations, efficient conformational sampling algorithms are required. In this work, we propose a new approach based on generalized replica-exchange with solute tempering (gREST) for exploring the open-closed conformational changes in multidomain proteins. Wherein, selected surface charged residues in a target protein are defined as the solute region in gREST simulation and the solute temperatures are different in replicas and exchanged between them to enhance the domain motions. This approach is called gREST selected surface charged residues (gREST_SSCR) and is applied to the Apo and Holo states of ribose binding protein (RBP) in solution. The conformational spaces sampled with gREST_SSCR are much wider than those with the conventional MD, sampling open-closed conformational changes while maintaining RBP domains’ stability. The free-energy landscapes of RBP in the Apo and Holo states are drawn along with twist and hinge angles of the two moving domains. The inter-domain salt-bridges that are not observed in the experimental structures are also important in the intermediate states during the conformational changes.

Recently, we proposed novel temperature and pressure evaluations in molecular dynamics (MD) simulation to preserve the accuracy up to the third order of a time step, δt (Jung et al., J. Chem. Theory Comput. 15, 84-94 (2019) and Jung et al., J. Chem. Phys. 148, 164109 (2018)). These approaches allow us to extend δt of MD simulation in isothermal-isobaric condition up to 5 fs with velocity Verlet integrator. Here, we further improve the isothermal-isobaric MD integration by introducing the group-based evaluations of system temperature and pressure to our previous approach. The group-based scheme increases the accuracy even using inaccurate temperature and pressure evaluations by neglecting the high-frequency vibrational motions of hydrogen atoms. It also improves the overall performance by avoiding iterations in thermostat and barostat updates and by allowing a multiple time step (MTS) integration such as r-REPSA (reversible reference system propagation algorithm) with our proposed high-precision evaluations of temperature and pressure. Now, the improved integration scheme conserves physical properties of lipid-bilayer systems up to δt = 5 fs with velocity Verlet as well as δt = 3.5 fs for fast motions in r-RESPA, respectively.

Alchemical free energy simulations have long been utilized to predict free energy changes for binding affinity and solubility of small molecules. However, while the theoretical foundation of these methods is well established, seamlessly handling many of the practical aspects regarding the preparation of the different thermodynamic end states of complex molecular systems and the numerous processing scripts often remains a burden for successful applications. In this work, we present CHARMM-GUI Free Energy Calculator (http://www.charmm-gui.org/input/fec) that provides various alchemical free energy perturbation molecular dynamics (FEP/MD) systems with input and post-processing scripts for NAMD and GENESIS. Four submodules are available: Absolute Ligand Binder (for absolute ligand binding FEP/MD), Relative Ligand Binder (for relative ligand binding FEP/MD), Absolute Ligand Solvator (for absolute ligand solvation FEP/MD), and Relative Ligand Solvator (for relative ligand solvation FEP/MD). Each module is designed to build multiple systems of a set of selected ligands at once for high-throughput FEP/MD simulations. The capability of Free Energy Calculator is illustrated by absolute and relative solvation FEP/MD of a set of ligands and absolute and relative binding FEP/MD of a set of ligands for T4-lysozyme in solution and the adenosine A2A receptor in a membrane. The calculated free energy values are overall consistent with the experimental and published free energy results (within ∼1 kcal/mol). We hope that Free Energy Calculator is useful to carry out high-throughput FEP/MD simulations in the field of biomolecular sciences and drug discovery.

The molecular function of Atg9, the sole transmembrane protein in the autophagosome-forming machinery, remains unknown. Atg9 colocalizes with Atg2 at the expanding edge of the isolation membrane (IM), where Atg2 receives phospholipids from the endoplasmic reticulum (ER). Here we report that yeast and human Atg9 are lipid scramblases that translocate phospholipids between outer and inner leaflets of liposomes in vitro. Cryo-EM of fission yeast Atg9 reveals a homotrimer, with two connected pores forming a path between the two membrane leaflets: one pore, located at a protomer, opens laterally to the cytoplasmic leaflet; the other, at the trimer center, traverses the membrane vertically. Mutation of residues lining the pores impaired IM expansion and autophagy activity in yeast and abolished Atg9’s ability to transport phospholipids between liposome leaflets. These results suggest that phospholipids delivered by Atg2 are translocated from the cytoplasmic to the luminal leaflet by Atg9, thereby driving autophagosomal membrane expansion.

The accurate prediction of protein–ligand binding affinity is a central challenge in computational chemistry and in-silico drug discovery. The free energy perturbation (FEP) method based on molecular dynamics (MD) simulation provides reasonably accurate results only if a reliable structure is available via high-resolution X-ray crystallography. To overcome the limitation, we propose a sequential prediction protocol using generalized replica exchange with solute tempering (gREST) and FEP. At first, ligand binding poses are predicted using gREST, which weakens protein–ligand interactions at high temperatures to sample multiple binding poses. To avoid ligand dissociation at high temperatures, a flat-bottom restraint potential centered on the binding site is applied in the simulation. The binding affinity of the most reliable pose is then calculated using FEP. The protocol is applied to the bindings of ten ligands to FK506 binding proteins (FKBP), showing the excellent agreement between the calculated and experimental binding affinities. The present protocol, which is referred to as the gREST+FEP method, would help to predict the binding affinities without high-resolution structural information on the ligand-bound state.

Conventional torsion angle potentials used in molecular dynamics (MD) have a singularity problem when three bonded particles are collinearly aligned. This problem is often encountered in coarse-grained (CG) simulations. Here, we propose a new form of the torsion angle potential, which introduces an angle-dependent modulating function. By carefully tuning the parameters for this modulating function, our method can eliminate the problematic angle-dependent singularity while being combined with existing models. As an example, we optimized the modulating function of the torsion angle potential for popular CG models of biomolecules based on the statistics over experimental structures deposited in the Protein Data Bank. By applying our method to designed and natural biomolecules, we show that the new torsion angle potential is able to eliminate the singularity problem while maintaining the structural features in the original models. Furthermore, by comparing our design with previous methods, we found that our new potential has advantages in computational efficiency and numerical stability. We strongly recommend the usage of our new potential in the CG simulations of flexible molecules.

Three-dimensional structural diversity (conformational diversity) of glycans is essential for their biological functions. Molecular dynamics (MD) simulations provide insight into the three-dimensional (3D-) structures of biomolecules in atomic resolution which are difficult to obtain by experiment. However, their highly complex and flexible structures make glycans difficult to simulate. This minireview summarizes recent challenges in MD simulations of glycans using enhanced sampling techniques. Using MD simulations, various 3D-structures are revealed at atomic resolution, not only for isolated glycans but also for glycoconjugates and glycan-bound complexes. The role of MD simulations is not limited for interpretation of experimental results. The simulations are used to predict experiments, leading the integrated computational and experimental studies on the structure-function relationship of glycans. We hope that this article will help prompt various integrated studies of glycans in the future.

Soluble prion proteins contingently encounter foreign prion aggregates, leading to cross-species prion transmission. However, how its efficiency is regulated by structural fluctuation of the host soluble prion protein remains unsolved. In the present study, through the use of two distantly related yeast prion Sup35 proteins, we found that a specific conformation of a short disordered segment governs interspecies prion transmissibility. Using a multidisciplinary approach including high-resolution NMR and molecular dynamics simulation, we identified critical residues within this segment that allow interspecies prion transmission in vitro and in vivo, by locally altering dynamics and conformation of soluble prion proteins. Remarkably, subtle conformational differences caused by a methylene group between asparagine and glutamine sufficed to change the short segment structure and substantially modulate the cross-seeding activity. Thus, our findings uncover how conformational dynamics of the short segment in the host prion protein impacts cross-species prion transmission. More broadly, our study provides mechanistic insights into cross-seeding between heterologous proteins.

Atomically detailed description of conformational dynamics in biomolecules is often essential to understand biological functions. Combining experimental measurements with molecular simulations significantly improves the outcome. Ensemble refinements, where the simulations are utilized to refine ensemble averaged data in NMR, SAXS, or cryo-EM, are a popular approach in integrative structural biology. Single-molecule time-series data contain rich temporal information of biomolecular dynamics. However, direct usage of the time-series data together with molecular simulations is just beginning. Here, we review data-assimilation approaches linking molecular simulations with experimental time-series data and discuss current limitations and potential applications of this approach in integrative structural biology.

The heme ATP-binding cassette (ABC) transporter BhuUV-T of bacterial pathogen Burkholderia cenocepacia is required to transport heme across the inner cell membrane. The current hypothesis is that the binding of two ATPs to the nucleotide-binding domains of the transporter drives the initial steps of the transport cycle in which the empty transport sites are reoriented from the cytosol to the periplasm. Molecular details are missing because the structure of a key occluded intermediate remains hypothetical. Here we perform molecular simulations to analyze the free energy surface (FES) of the first step of the reorientation, namely the transition from an open inward-facing (IF) transport site to an occluded (Occ) conformation. We have modeled the latter structure in silico in a previous study. A simple annealing procedure removes residual bias originating from non-equilibrium targeted molecular dynamics. The calculated FES reveals the role of the ATPs in inducing the IF → Occ conformational change and validates the modeled Occ conformation.

Fibroblast growth factor receptor 3 (FGFR3) is a member of receptor tyrosine kinases, which is involved in skeletal cell growth, differentiation, and migration. FGFR3 transduces biochemical signals from the extracellular ligand-binding domain to the intracellular kinase domain through the conformational changes of the transmembrane (TM) helix dimer. Here, we apply generalized replica exchange with solute tempering method to wild type (WT) and G380R mutant (G380R) of FGFR3. The dimer interface in G380R is different from WT and the simulation results are in good agreement with the solid-state nuclear magnetic resonance (NMR) spectroscopy. TM helices in G380R are extended more than WT, and thereby, G375 in G380R contacts near the N‐termini of the TM helix dimer. Considering that both G380R and G375C show the constitutive activation, the formation of the N‐terminal contacts of the TM helices can be generally important for the activation mechanism.

Slow polypeptide conformational changes on time scales of >1 s are generally assumed to be highly cooperative two-state transitions, reflecting the high energy barrier. However, few experimental characterizations have tested the validity of this assumption. We performed residue-specific NMR thermodynamic analysis of the 27-residue lantibiotic peptide, nukacin ISK-1, to characterize the isomerization between two topological states on the second time scale. Unexpectedly, the thermal transition behaviors were distinct among peptide regions, indicating that the topological isomerization process is a mosaic of different degrees of cooperativity. The conformational change path between the two NMR structures was deduced by a targeted molecular dynamics simulation. The unique side-chain threading motions through the monosulfide rings are the structural basis of the high energy barrier, and the nonlocal interactions in the hydrophobic core are the structural basis of the cooperativity. Taken together, we provide an energetic description of the topological isomerization of nukacin ISK-1.

Dynamic interactions between proteins and nucleic acids are involved in many essential biological processes such as DNA replication, repair, and transcription. Recent advances in experimental and theoretical techniques have greatly increased our understanding of gene transcription at the molecular and cellular levels, and there is a need for integrative approaches to combine the data and knowledge obtained from multiple techniques. In this session, we discussed how hybrid approaches using single-molecule imaging, X-ray crystallography, cryo-electron microscopy (cryo-EM), and computer simulations could contribute to the comprehensive understanding of gene expression.

Reverse osmosis membranes based on aromatic polyamide (ar-PA) are widely used in desalination of seawater, yet the microscopic mechanism of water diffusion through a polyamide layer remains elusive. Here, we study the structure and dynamics of polymer chains and water molecules in ar-PA in comparison to nylon 6 (one of aliphatic polyamides) under various water contents (0.0–15.9 wt%). The infrared (IR) difference spectrum between dry and moist ar-PA shows little change in amide A bands, in contrast to that of nylon 6 which yields a prominent dip. Theoretical analyses using molecular dynamics simulations and quantum electronic and vibrational calculations reveal that the dip in nylon 6 is caused by breaking of hydrogen bonds (HBs) among amide groups. The incoming water molecules that break amide-amide HBs are bound to polyamide chains nearby and diffuse slowly. On the other hand, the amide-amide HBs of ar-PA are kept upon hydration. Such polymer structure facilitates growth of large water clusters with more than 100 water molecules and rapid diffusion of water molecules. The amide A band serves as a fingerprint to characterize the water permeability of polyamide materials.

Surfactant micelles are often utilized as membrane mimetics for structure determination and functional analysis of membrane proteins. The curved-surface effects of the micelle can perturb membrane protein structure. However, it is difficult to assess such effects and membrane mimetic artifacts by experimental and theoretical methods. Here, we propose an implicit micelle model (IMIC) to be used in molecular dynamics (MD) simulations of membrane proteins. IMIC is an extension of the IMM1 implicit membrane model and additionally introduces a superellipsoid approximation to represent the curved-surface effects. Most of the IMIC parameters are obtained from all-atom explicit solvent MD simulations of 12 membrane proteins in various micelles. The HIV envelope protein gp41, M13 major coat protein gp8, and amyloid precursor protein (APP) dimer are simulated via MD simulations with IMIC. These simulations clearly show how the micelle influences membrane protein structures compared to the bilayer environments. The MD simulations with IMIC provide reliable membrane protein structures in various micelle environments quickly with smaller computational cost than that for an explicit solvent/micelle model.

Human leukocyte Ig-like receptors (LILR) LILRB1 and LILRB2 are immune checkpoint receptors that regulate a wide range of physiological responses by binding to diverse ligands, including HLA-G. HLA-G is exclusively expressed in the placenta, some immunoregulatory cells, and tumors and has several unique isoforms. However, the recognition of HLA-G isoforms by LILRs is poorly understood. In this study, we characterized LILR binding to the β2-microglobulin (β2m)-free HLA-G1 isoform, which is synthesized by placental trophoblast cells and tends to dimerize and multimerize. The multimerized β2m-free HLA-G1 dimer lacked detectable affinity for LILRB1, but bound strongly to LILRB2. We also determined the crystal structure of the LILRB1 and HLA-G1 complex, which adopted the typical structure of a classical HLA class I complex. LILRB1 exhibits flexible binding modes with the α3 domain, but maintains tight contacts with β2m, thus accounting for β2m-dependent binding. Notably, both LILRB1 and B2 are oriented at suitable angles to permit efficient signaling upon complex formation with HLA-G1 dimers. These structural and functional features of ligand recognition by LILRs provide novel insights into their important roles in the biological regulations.

Atomistic molecular dynamics simulations of concentrated protein solutions in the presence of a phospholipid bilayer are presented to gain insights into the dynamics and interactions at the cytosol–membrane interface. The main finding is that proteins that are not known to specifically interact with membranes are preferentially excluded from the membrane, leaving a depletion zone near the membrane surface. As a consequence, effective protein concentrations increase, leading to increased protein contacts and clustering, whereas protein diffusion becomes faster near the membrane for proteins that do occasionally enter the depletion zone. Since protein–membrane contacts are infrequent and short-lived in this study, the structure of the lipid bilayer remains largely unaffected by the crowded protein solution, but when proteins do contact lipid head groups, small but statistically significant local membrane curvature is induced, on average.

The dual-basin Gō-model is a structural-based coarse-grained model for simulating a conformational transition between two known structures of a protein. Two parameters are required to produce a dual-basin potential mixed using two single-basin potentials, although the determination of mixing parameters is usually not straightforward. Here, we have developed an efficient scheme to determine the mixing parameters using the Multistate Bennett Acceptance Ratio (MBAR) method after short simulations with a set of parameters. In the scheme, MBAR allows us to predict observables at various unsimulated conditions, which are useful to improve the mixing parameters in the next round of iterative simulations. The number of iterations that are necessary for obtaining the converged mixing parameters are significantly reduced in the scheme. We applied the scheme to two proteins, the glutamine binding protein and the ribose binding protein, for showing the effectiveness in the parameter determination. After obtaining the converged parameters, both proteins show frequent conformational transitions between open and closed states, providing the theoretical basis to investigate structure-dynamics-function relationships of the proteins.

In this study, a replica-exchange method was developed to overcome conformational sampling difficulties in computer simulations of spin glass or other systems with rugged free-energy landscapes. This method was then applied to the protein-folding problem in combination with molecular dynamics (MD) simulation. Owing to its simplicity and sampling efficiency, the replica-exchange method has been applied to many other biological problems and has been continuously improved. The method has often been combined with other sampling techniques, such as umbrella sampling, free-energy perturbation, metadynamics, and Gaussian accelerated MD (GaMD). In this chapter, we first summarize the original replica-exchange molecular dynamics (REMD) method and discuss how new algorithms related to the original method are implemented to add new features. Heterogeneous and flexible structures of an N-glycan in a solution are simulated as an example of applications by REMD, replica exchange with solute tempering, and GaMD. The sampling efficiency of these methods on the N-glycan system and the convergence of the free-energy changes are compared. REMD simulation protocols and trajectory analysis using the GENESIS software are provided to facilitate the practical use of advanced simulation methods.

All-atom molecular dynamics (MD) simulation has become a useful tool to investigate protein dynamics in solution or in membranes. Such simulation can provide atomic detailed information of protein dynamics on microseconds or longer timescales these days. Due to advances of simulation methodologies and computer power, MD simulation is now applied to crowded protein systems to study the relationship between the structure, dynamics, and function of proteins in the environment. Multiple proteins in the system are simulated in a solvated box at a high concentration, where proteins interact with each other via weak and non-specific molecular interactions. Protein stability and translational/rotational diffusion are compared to NMR chemical shifts and relaxation parameters, respectively. In this chapter, we review the current status of all-atom MD simulations of crowded protein systems and discuss how simulation can contribute to understanding the effects of macromolecular crowding on proteins. All-atom models of the cytoplasm in Mycoplasma genitalium were built at the atomic resolution, which includes more than a hundred million atoms and consists of proteins, RNAs, ribosomes, metabolites, ions, and water molecules. All-atom MD simulations of the cytoplasmic model were carried out using highly parallelized MD software on supercomputers and have given us several predictions about cellular crowding effects, which should be examined by future experiments.

α-Amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA)-type glutamate receptors produce postsynaptic current by transmitting an agonist-induced structural change in the ligand-binding domain (LBD) to the transmembrane channel. Receptors carrying T686S/A substitutions in their LBDs produce weaker glutamate-evoked currents than wild-type (WT) receptors. However, the substitutions induce little differences in the crystal structures of their LBDs. To understand the structural mechanism underlying reduced activities of these AMPAR variants, we analyzed the structural dynamics of WT, T686S, and T686A variants of LBD using nuclear magnetic resonance. The HD exchange studies of the LBDs showed that the kinetic step where the ligand-binding cleft closes was changed by the substitutions, and the substitution-induced population shift from cleft-closed to cleft-open structures is responsible for the reduced activities of the variants. The chemical shift analyses revealed another structural equilibrium between cleft-locked and cleft-partially-open conformations. The substitution-induced population shift in this equilibrium may be related to slower desensitization observed for these variants.

We have developed an enhanced conformational sampling method combining replica-exchange umbrella sampling (REUS) with Gaussian accelerated molecular dynamics (GaMD). REUS enhances the sampling along predefined reaction coordinates, while GaMD accelerates the conformational dynamics by adding a boost potential to the system energy. The method, which we call GaREUS (Gaussian accelerated replica-exchange umbrella sampling), enhances the sampling more efficiently than REUS or GaMD, while the computational resource for GaREUS is the same as that required for REUS. The two-step reweighting procedure using the multistate Bennett acceptance ratio method and the cumulant expansion for the exponential average is applied to the simulation trajectories for obtaining the unbiased free-energy landscapes. We apply GaREUS to the calculations of free-energy landscapes for three different cases: conformational equilibria of N-glycan, folding of chignolin, and conformational change of adenyl kinase. We show that GaREUS speeds up the convergences of free-energy calculations using the same amount of computational resources as REUS. The free-energy landscapes reweighted from the trajectories of GaREUS agree with previously reported ones. GaREUS is applicable to free-energy calculations of various biomolecular dynamics and functions with reasonable computational costs.

Modern drug discovery increasingly focuses on the drug-target binding kinetics which depend on drug (un)binding pathways. The conventional molecular dynamics simulation can observe only a few binding events even using the fastest supercomputer. Here, we develop 2D gREST/REUS simulation with enhanced flexibility of the ligand and the protein binding site. Simulation (43 μs in total) applied to an inhibitor binding to c-Src kinase covers 100 binding and unbinding events. On the statistically converged free-energy landscapes, we succeed in predicting the X-ray binding structure, including water positions. Furthermore, we characterize hidden semibound poses and transient encounter complexes on the free-energy landscapes. Regulatory residues distant from the catalytic core are responsible for the initial inhibitor uptake and regulation of subsequent bindings, which was unresolved by experiments. Stabilizing/blocking of either the semibound poses or the encounter complexes can be an effective strategy to optimize drug-target residence time.

Molecular recognition underpins all specific protein-ligand interactions and is essential for biomolecular functions. Prediction of canonical binding poses and distinguishing binders from non-binders are much sought after goals. Here, we apply the generalized replica exchange with solute tempering method, gREST, combined with a flat-bottom potential to evaluate binder and non-binder interaction with T4 lysozyme L99A mutant. The buried hydrophobic cavity and possibility of coupled conformational changes in this protein make binding predictions difficult. The present gREST simulations, enabling enhanced flexibilities of the ligand and protein residues near the binding site, sample bindings in multiple poses and correctly portray X-ray structures. The free-energy profiles of binders (benzene, ethylbenzene, and n-hexylbenzene) are distinct from those of non-binders (phenol and benzaldehyde). Bindings of the two larger molecules seem to be associated with a structural change toward an excited conformation of the protein, which agrees with experimental findings. The protocol is generally applicable to various proteins having buried cavities with limited access for ligands with different shapes, sizes and chemical properties.

The heme importer from pathogenic bacteria is a member of the ATP-binding cassette (ABC) transporter family, which uses the energy of ATP-binding and hydrolysis for extensive conformational changes. Previous studies have indicated that conformational changes after heme translocation are triggered by ATP-binding to nucleotide binding domains (NBDs) and then, in turn, induce conformational transitions of the transmembrane domains (TMDs). In this study, we applied a template-based iterative all-atom molecular dynamics (MD) simulation to predict the ATP-bound outward-facing conformation of the Burkholderia cenocepacia heme importer BhuUV-T. The resulting model showed a stable conformation of the TMD with the cytoplasmic gate in the closed state and the periplasmic gate in the open state. Furthermore, targeted MD simulation predicted the intermediate structure of an occluded form (Occ) with bound ATP, in which both ends of the heme translocation channel are closed. The MD simulation of the predicted Occ revealed that Ser147 on the ABC signature motifs (LSGG[Q/E]) of NBDs occasionally flips and loses the active conformation required for ATP-hydrolysis. The flipping motion was found to be coupled to the inter-NBD distance. Our results highlight the functional significance of the signature motif of ABC transporters in regulation of ATPase and chemo-mechanical coupling mechanism.

A weight averaged method is reviewed, which combines molecular dynamics (MD) simulations and vibrational quasi‐degenerate perturbation theory (VQDPT) for computing the vibrational spectrum of complex systems. The main idea of the method is (i) the vibrational states are well localized to a relatively narrow region in space (e.g. a water molecule, chemical groups, etc.), (ii) such regions give rise to a set of local spectra, and (iii) the total spectrum is described by a weighted sum of these spectra. Applications to a pentapeptide (Ser‐Ile‐Val‐Ser‐Phe), a lipid (sphingomyelin) bilayer, and a hydration of nylon 6 demonstrate that the method reproduces the experimental spectrum accurately and gives an assignment of the observed bands. The developed method makes it feasible to establish a firm connection between spectral and structural information, which provides new insight into the molecular mechanism of complex systems.

Comprehensive data about the composition and structure of cellular components have enabled the construction of quantitative whole-cell models. While kinetic network–type models have been established, it is also becoming possible to build physical, molecular-level models of cellular environments. This review outlines challenges in constructing and simulating such models and discusses near- and long-term opportunities for developing physical whole-cell models that can connect molecular structure with biological function.

The growing interest in the complexity of biological interactions is continuously driving the need to increase system size in biophysical simulations, requiring not only powerful and advanced hardware but adaptable software that can accommodate a large number of atoms interacting through complex forcefields. To address this, we developed and implemented strategies in the GENESIS molecular dynamics package designed for large numbers of processors. Long‐range electrostatic interactions were parallelized by minimizing the number of processes involved in communication. A novel algorithm was implemented for nonbonded interactions to increase single instruction multiple data (SIMD) performance, reducing memory usage for ultra large systems. Memory usage for neighbor searches in real‐space nonbonded interactions was reduced by approximately 80%, leading to significant speedup. Using experimental data describing physical 3D chromatin interactions, we constructed the first atomistic model of an entire gene locus (GATA4). Taken together, these developments enabled the first billion‐atom simulation of an intact biomolecular complex, achieving scaling to 65,000 processes (130,000 processor cores) with 1 ns/day performance.

Quantum mechanics/molecular mechanics (QM/MM) calculations are applied for anharmonic vibrational analyses of biomolecules and solvated molecules. The QM/MM method is implemented into a molecular dynamics (MD) program, GENESIS, by interfacing with external electronic structure programs. Following the geometry optimization and the harmonic normal-mode analysis based on a partial Hessian, the anharmonic potential energy surface (PES) is generated from QM/MM energies and gradients calculated at grid points. The PES is used for vibrational self-consistent field (VSCF) and post-VSCF calculations to compute the vibrational spectrum. The method is first applied to a phosphate ion in solution. With both the ion and neighboring water molecules taken as a QM region, IR spectra of representative hydration structures are calculated by the second-order vibrational quasi-degenerate perturbation theory (VQDPT2) at the level of B3LYP/cc-pVTZ and TIP3P force field. A weight-average of IR spectra over the structures reproduces the experimental spectrum with a mean absolute deviation of 16 cm^–1. Then, the method is applied to an enzyme, P450 nitric oxide reductase (P450nor), with the NO molecule bound to a ferric (Fe^III) heme. Starting from snapshot structures obtained from MD simulations of P450nor in solution, QM/MM calculations have been carried out at the level of B3LYP-D3/def2-SVP(D). The spin state of Fe^III(NO) is likely a closed-shell singlet state based on a ratio of N–O and Fe–NO stretching frequencies (ν_N–O and ν_Fe–NO) calculated for closed- and open-shell singlet states. The calculated ν_N–O and ν_Fe–NO overestimate the experimental ones by 120 and 75 cm^–1, respectively. The electronic structure and solvation of Fe^III(NO) affect the structure around the heme of P450nor leading to an increase in ν_N–O and ν_Fe–NO.

The rotational diffusion of a protein in the presence of protein crowder molecules was analyzed via computer simulations. Cluster formation as a result of transient intermolecular contacts was identified as the dominant effect for reduced rotational diffusion upon crowding. The slow-down in diffusion was primarily correlated with direct protein–protein contacts rather than indirect interactions via shared hydration layers. But increased solvent viscosity due to crowding contributed to a lesser extent. Key protein–protein contacts correlated with a slow-down in diffusion involve largely interactions between charged and polar groups suggesting that the surface composition of a given protein and the resulting propensity for forming interactions with surrounding proteins in a crowded cellular environment may be the major determinant of its diffusive properties.

α-amino-3-hydroxy-5-methyl-4-isoaxazolepropionic acid (AMPA) ionotropic glutamate receptors mediate fast excitatory neurotransmission in the central nervous system, and their dysfunction is associated with neurological diseases. Glutamate binding to ligand-binding domains (LBDs) of AMPA receptors induces channel opening in the transmembrane domains of the receptors. The T686A mutation reduces glutamate efficacy so that the glutamate behaves as a partial agonist. The crystal structures of wild-type and mutant LBDs are very similar and cannot account for the observed behavior. To elucidate the molecular mechanism inducing partial agonism of the T686A mutant, we computed the free-energy landscapes governing GluA2 LBD closure using replica-exchange umbrella sampling simulations. A semiclosed state, not observed in crystal structures, appears in the mutant during simulation. In this state, the LBD cleft opens slightly because of breaking of interlobe hydrogen bonds, reducing the efficiency of channel opening. The energy difference between the LBD closed and semiclosed states is small, and transitions between the two states would occur by thermal fluctuations. Evidently, glutamate binding to the T686A mutant induces a population shift from a closed to a semiclosed state, explaining the partial agonism in the AMPA receptor.

In molecular dynamics (MD) simulations, an accurate evaluation of temperature is essential for controlling temperature as well as pressure in the isothermal–isobaric conditions. According to the Tolman’s equipartition theorem, all motions of all particles should share a single temperature. However, conventional temperature estimation from kinetic energy does not include Hessian terms properly, and thereby, the equipartition theorem is not satisfied with a large time step. In this paper, we show how to evaluate temperature the most accurately without increasing computational cost. We define two kinds of kinetic energies, evaluated at full- and half-time steps that underestimate or overestimate temperature, respectively. A combination of these two kinetic energies provides an optimal instantaneous temperature up to the third order of the time step. The method is tested for a one-dimensional harmonic oscillator, pure water molecules, a Bovine pancreatic trypsin inhibitor (BPTI) protein in water molecules, and a hydrated 1,2-dispalmitoyl-sn-phosphatidylcholine (DPPC) lipid bilayer in water molecules. In all tests, the optimal temperature estimator fulfills the equipartition theorem better than existing methods and reproduces well the usual physical properties for time steps up to and including 5 fs.

Flexible fitting is a powerful technique to build the 3D structures of biomolecules from cryoelectron microscopy (cryo-EM) density maps. One popular method is a cross-correlation coefficient-based approach, where the molecular dynamics (MD) simulation is carried out with the biasing potential that includes the cross-correlation coefficient between the experimental and simulated density maps. Here, we propose efficient parallelization schemes for the calculation of the cross-correlation coefficient to accelerate flexible fitting. Our schemes are tested for small, medium, and large biomolecules using CPU and hybrid CPU + GPU architectures. The scheme for the atomic decomposition MD is suitable for small proteins such as Ca²⁺-ATPase with the all-atom Go model, while that for the domain decomposition MD is better for larger systems such as ribosome with the all-atom Go or the all-atom explicit solvent models. Our methods allow flexible fitting for various biomolecules with reasonable computational cost. This approach also connects high-resolution structure refinements with investigation of protein structure-function relationship.

A synthetic riboswitch N1, inserted into the 5′-untranslated mRNA region of yeast, regulates gene expression upon binding ribostamycin and neomycin. Interestingly, a similar aminoglycoside, paromomycin, differing from neomycin by only one substituent (amino versus hydroxyl), also binds to the N1 riboswitch, but without affecting gene expression, despite NMR evidence that the N1 riboswitch binds all aminoglycosides in a similar way. Here, to explore the details of structural dynamics of the aminoglycoside-N1 riboswitch complexes, we applied all-atom molecular dynamics (MD) and temperature replica-exchange MD simulations in explicit solvent. Indeed, we found that ribostamycin and neomycin affect riboswitch dynamics similarly but paromomycin allows for more flexibility because its complex lacks the contact between the distinctive 6′ hydroxyl group and the G9 phosphate. Instead, a transient hydrogen bond of 6′-OH with A17 is formed, which partially diminishes interactions between the bulge and apical loop of the riboswitch, likely contributing to riboswitch inactivity. In many ways, the paromomycin complex mimics the conformations, interactions, and Na⁺ distribution of the free riboswitch. The MD-derived interaction network helps understand why riboswitch activity depends on aminoglycoside type, whereas for another aminoglycoside-binding site, aminoacyl-tRNA site in 16S rRNA, activity is not discriminatory.

Replica-exchange molecular dynamics (REMD) and their variants have been widely used in simulations of the biomolecular structure and dynamics. Replica exchange with solute tempering (REST) is one of the methods where temperature of a pre-defined solute molecule is exchanged between replicas, while solvent temperatures in all the replicas are kept constant. REST greatly reduces the number of replicas compared to the temperature REMD, while replicas at low temperatures are often trapped under their conditions, interfering with the conformational sampling. Here, we introduce a new scheme of REST, referred to as generalized REST (gREST), where the solute region is defined as a part of a molecule or a part of the potential energy terms, such as the dihedral-angle energy term or Lennard-Jones energy term. We applied this new method to folding simulations of a β-hairpin (16 residues) and a Trp-cage (20 residues) in explicit water. The protein dihedral-angle energy term is chosen as the solute region in the simulations. gREST reduces the number of replicas necessary for good random walks in the solute-temperature space and covers a wider conformational space compared to the conventional REST2. Considering the general applicability, gREST should become a promising tool for the simulations of protein folding, conformational dynamics, and an in silico drug design.

Single-molecule experiments and molecular dynamics (MD) simulations are indispensable tools for investigating protein conformational dynamics. The former provide timeseries data, such as donor-acceptor distances, whereas the latter give atomistic information, although this information is often biased by model parameters. Here, we devise a machine-learning method to combine the complementary information from the two approaches and construct a consistent model of conformational dynamics. It is applied to the folding dynamics of the formin-binding protein WW domain. MD simulations over 400 μs led to an initial Markov state model (MSM), which was then “refined” using single-molecule Förster resonance energy transfer (FRET) data through hidden Markov modeling. The refined or data-assimilated MSM reproduces the FRET data and features hairpin one in the transition-state ensemble, consistent with mutation experiments. The folding pathway in the data-assimilated MSM suggests interplay between hydrophobic contacts and turn formation. Our method provides a general framework for investigating conformational transitions in other proteins.

A data-driven modeling scheme is proposed for conformational dynamics of biomolecules based on molecular dynamics (MD) simulations and experimental measurements. In this scheme, an initial Markov State Model (MSM) is constructed from MD simulation trajectories, and then, the MSM parameters are refined using experimental measurements through machine learning techniques. The second step can reduce the bias of MD simulation results due to inaccurate force-field parameters. Either time-series trajectories or ensemble-averaged data are available as a training data set in the scheme. Using a coarse-grained model of a dye-labeled polyproline-20, we compare the performance of machine learning estimations from the two types of training data sets. Machine learning from time-series data could provide the equilibrium populations of conformational states as well as their transition probabilities. It estimates hidden conformational states in more robustways compared to that from ensemble-averaged data although there are limitations in estimating the transition probabilities between minor states. We discuss how to use the machine learning scheme for various experimental measurements including single-molecule time-series trajectories.

In molecular dynamics (MD) simulations, a proper definition of kinetic energy is essential for controlling pressure as well as temperature in the isothermal-isobaric condition. The virial theorem provides an equation that connects the average kinetic energy with the product of particle coordinate and force. In this paper, we show that the theorem is satisfied in MD simulations with a larger time step and holonomic constraints of bonds, only when a proper definition of kinetic energy is used. We provide a novel definition of kinetic energy, which is calculated from velocities at the half-time steps (t − Δt/2 and t + Δt/2) in the velocity Verlet integration method. MD simulations of a 1,2-dispalmitoyl-sn-phosphatidylcholine (DPPC) lipid bilayer and a water box using the kinetic energy definition could reproduce the physical properties in the isothermal-isobaric condition properly. We also develop a multiple time step (MTS) integration scheme with the kinetic energy definition. MD simulations with the MTS integration for the DPPC and water box systems provided the same quantities as the velocity Verlet integration method, even when the thermostat and barostat are updated less frequently.

The 99 amino acid C-terminal fragment of Amyloid Precursor Protein APP-C99 (C99) is cleaved by γ-secretase to form Aβ peptide, which plays a critical role in the etiology of Alzheimer’s Disease (AD). The structure of C99 consists of a single transmembrane domain flanked by intra and intercellular domains. While the structure of the transmembrane domain has been well characterized, little is known about the structure of the flanking domains and their role in C99 processing by γ-secretase. To gain insight into the structure of full-length C99, REMD simulations were performed for monomeric C99 in model membranes of varying thickness. We find equilibrium ensembles of C99 from simulation agree with experimentally-inferred residue insertion depths and protein backbone chemical shifts. In thin membranes, the transmembrane domain structure is correlated with extra-membrane structural states and the extra-membrane domain structural states become less correlated to each other. Mean and variance of the transmembrane and G₃₇G₃₈ hinge angles are found to increase with thinning membrane. The N-terminus of C99 forms β-strands that may seed aggregation of Aβ on the membrane surface, promoting amyloid formation. In thicker membranes the N-terminus forms α-helices that interact with the nicastrin domain of γ-secretase. The C-terminus of C99 becomes more α-helical as the membrane thickens, forming structures that may be suitable for binding by cytoplasmic proteins, while C-terminal residues essential to cytotoxic function become α-helical as the membrane thins. The heterogeneous but discrete extra-membrane domain states analyzed here open the path to new investigations of the role of C99 structure and membrane in amyloidogenesis.

Anharmonic vibrational state calculations were performed for PtCO and Ar–PtCO via the direct vibrational configuration interaction (VCI) method based on CCSD(T) energies and CCSD dipole moments at tens of thousands of grids, to get insights into the anomalous effect of a solid argon matrix on the vibrational spectra of PtCO. It was shown that, through the binding of Ar to PtCO via a strong van der Waals interaction, the Pt–C–O bending fundamental level drastically loses the infrared intensity although the corresponding overtone band shows a relatively large intensity. The origin of this phenomenon was analyzed based on the dipole moment surfaces and electron densities around the equilibrium structure. The present computations have solved the inconsistency between the gas-phase and the matrix-isolation experiments for PtCO.

Ion mobility mass spectrometry (IM-MS) is a technique capable of investigating structural changes of biomolecules based on their collision cross section (CCS). Recent advances in IM-MS allow us to separate carbohydrate isomers with subtle conformational differences, but the relationship between CCS and atomic structure remains elusive. Here, we characterize conformational ensembles of gasphase N-glycans under the electrospray ionization condition using molecular dynamics simulations with enhanced sampling. We show that the separation of CCSs between isomers reflects folding features of N-glycans, which are determined both by chemical compositions and protonation states. Providing a physicochemical basis of CCS for N-glycans helps not only to interpret IM-MS measurements but also to estimate CCSs of complex glycans.

For a long time, the effect of a crowded cellular environment on protein dynamics has been largely ignored. Recent experiments indicate that proteins diffuse more slowly in a living cell than in a diluted solution, and further studies suggest that the diffusion depends on the local surroundings. Here, detailed insight into how diffusion depends on protein–protein contacts is presented based on extensive all-atom molecular dynamics simulations of concentrated villin headpiece solutions. After force field adjustments in the form of increased protein–water interactions to reproduce experimental data, translational and rotational diffusion was analyzed in detail. Although internal protein dynamics remained largely unaltered, rotational diffusion was found to slow down more significantly than translational diffusion as the protein concentration increased. The decrease in diffusion is interpreted in terms of a transient formation of protein clusters. These clusters persist on sub-microsecond time scales and follow distributions that increasingly shift toward larger cluster size with increasing protein concentrations. Weighting diffusion coefficients estimated for different clusters extracted from the simulations with the distribution of clusters largely reproduces the overall observed diffusion rates, suggesting that transient cluster formation is a primary cause for a slow-down in diffusion upon crowding with other proteins.

Atomic structures of proteins, nucleic acids, and their complexes are determined using X-ray crystallography, NMR, or cryo-Electron Microscopy. These structures are essential to understand their structure-function relationships. However, the experimental conditions are totally different from the actual cellular environments and it is hard to understand how biomolecules behave in such cellular environments, just using the atomic structures. We have recently built protein crowding systems in computers and carried out all-atom molecular dynamics (MD) simulations of the systems to understand biomolecular dynamics in the crowded environments. The largest simulations we have ever performed were the all-atom MD simulations of a bacterial cytoplasm using K computer. By analyzing the simulation trajectories, we observed that non-specific protein-protein and protein-metabolite interactions play important roles in biomolecular dynamics and stability in a cell. The new insight from the simulations is useful not only for basic life science in molecular and cellular biology but also drug discovery in future for introducing the effect of non-specific protein-drug interactions.

膜タンパク質は細胞膜近傍で生じる生命現象の多くを司る重要な生体分子である．X線結晶構造解析を中心とする構造生物学の発展により，膜タンパク質の立体構造が原子解像度で決定され，その立体構造と機能の関係が明らかになってきた．さらに膜タンパク質の機能を詳しく理解するためには立体構造に加えて分子ダイナミクスの詳細を知る必要がある．構造解析だけでは理解できない膜タンパク質の分子ダイナミクスと機能の関係を調べる手法として，近年，分子動力学計算などの計算機シミュレーションが注目されている．とくにアメリカで開発された分子動力学専用計算機Antonや，統計力学の原理に基づく効率のよい構造探索手法を用いることで，従来は不可能であった膜タンパク質の機能にかかわる遅い分子運動のシミュレーションが可能となり，実験結果と直接比較することができるようになった．さらに，膜タンパク質の構造ダイナミクスを考慮した創薬応用や膜タンパク質がかかわる疾患の分子機構の解明につながることが期待される．

Nitric oxide (NO) plays diverse and significant roles in biological processes despite its cytotoxicity, raising the question of how biological systems control the action of NO to minimize its cytotoxicity in cells. As a great example of such a system, we found a possibility that NO-generating nitrite reductase (NiR) forms a complex with NO-decomposing membrane-integrated NO reductase (NOR) to efficiently capture NO immediately after its production by NiR in anaerobic nitrate respiration called denitrification. The 3.2-Å resolution structure of the complex of one NiR functional homodimer and two NOR molecules provides an idea of how these enzymes interact in cells, while the structure may not reflect the one in cells due to the membrane topology. Subsequent all-atom molecular dynamics (MD) simulations of the enzyme complex model in a membrane and structure-guided mutagenesis suggested that a few interenzyme salt bridges and coulombic interactions of NiR with the membrane could stabilize the complex of one NiR homodimer and one NOR molecule and contribute to rapid NO decomposition in cells. The MD trajectories of the NO diffusion in the NiR:NOR complex with the membrane showed that, as a plausible NO transfer mechanism, NO released from NiR rapidly migrates into the membrane, then binds to NOR. These results help us understand the mechanism of the cellular control of the action of cytotoxic NO.

GENeralized-Ensemble SImulation System (GENESIS) is a software package for molecular dynamics (MD) simulation of biological systems. It is designed to extend limitations in system size and accessible time scale by adopting highly parallelized schemes and enhanced conformational sampling algorithms. In this new version, GENESIS 1.1, new functions and advanced algorithms have been added. The all-atom and coarse-grained potential energy functions used in AMBER and GROMACS packages now become available in addition to CHARMM energy functions. The performance of MD simulations has been greatly improved by further optimization, multiple time-step integration, and hybrid (CPU + GPU) computing. The string method and replica-exchange umbrella sampling with flexible collective variable choice are used for finding the minimum free-energy pathway and obtaining free-energy profiles for conformational changes of a macromolecule. These new features increase the usefulness and power of GENESIS for modeling and simulation in biological research.

The effects of crowding in biological environments on biomolecular structure, dynamics, and function remain not well understood. Computer simulations of atomistic models of concentrated peptide and protein systems at different levels of complexity are beginning to provide new insights. Crowding, weak interactions with other macromolecules and metabolites, and altered solvent properties within cellular environments appear to remodel the energy landscape of peptides and proteins in significant ways including the possibility of native state destabilization. Crowding is also seen to affect dynamic properties, both conformational dynamics and diffusional properties of macromolecules. Recent simulations that address these questions are reviewed here and discussed in the context of relevant experiments.

Structures of polyamide 6 are investigated for different hydration levels using molecular dynamics (MD) simulations and quantum vibrational calculations. The MD simulations have shown that hydration leads to an increase in the diffusion coefficient, accompanied by a growth of water clusters in the polymer. The IR difference spectra upon hydration are calculated using a weight-averaged method incorporating anharmonicity of the potential energy surface. The predicted IR difference spectrum for the amide A band is in quantitative agreement with the experiment [Iwamoto, R.; Murase, H. J. Polym. Sci., Part B: Polym. Phys. 2003, 41, 1722−1729]. The proposed method, combined with experimental IR difference spectra, makes it feasible to elucidate the atomistic structure of hydrated polymer materials.

A scoring protocol based on implicit membrane-based scoring functions and a new protocol for optimizing the positioning of proteins inside the membrane was evaluated for its capacity to discriminate native-like states from misfolded decoys. A decoy set previously established by the Baker lab (Proteins: Struct., Funct., Genet. 2006, 62, 1010–1025) was used along with a second set that was generated to cover higher resolution models. The Implicit Membrane Model 1 (IMM1), IMM1 model with CHARMM 36 parameters (IMM1-p36), generalized Born with simple switching (GBSW), and heterogeneous dielectric generalized Born versions 2 (HDGBv2) and 3 (HDGBv3) were tested along with the new HDGB van der Waals (HDGBvdW) model that adds implicit van der Waals contributions to the solvation free energy. For comparison, scores were also calculated with the distance-scaled finite ideal-gas reference (DFIRE) scoring function. Z-scores for native state discrimination, energy vs root-mean-square deviation (RMSD) correlations, and the ability to select the most native-like structures as top-scoring decoys were evaluated to assess the performance of the scoring functions. Ranking of the decoys in the Baker set that were relatively far from the native state was challenging and dominated largely by packing interactions that were captured best by DFIRE with less benefit of the implicit membrane-based models. Accounting for the membrane environment was much more important in the second decoy set where especially the HDGB-based scoring functions performed very well in ranking decoys and providing significant correlations between scores and RMSD, which shows promise for improving membrane protein structure prediction and refinement applications. The new membrane structure scoring protocol was implemented in the MEMScore web server (http://feiglab.org/memscore).

Protein secretion mediated by SecYEG translocon and SecA ATPase is enhanced by membrane-embedded SecDF by using proton motive force. A previous structural study of SecDF indicated that it comprises 12 transmembrane helices that can conduct protons and three periplasmic domains, which form at least two characterized transition states, termed the F and I forms. We report the structures of full-length SecDF in I form at 2.6- to 2.8-Å resolution. The structures revealed that SecDF in I form can generate a tunnel that penetrates the transmembrane region and functions as a proton pathway regulated by a conserved Asp residue of the transmembrane region. In one crystal structure, periplasmic cavity interacts with a molecule, potentially polyethylene glycol, which may mimic a substrate peptide. This study provides structural insights into the Sec protein translocation that allows future analyses to develop a more detailed working model for SecDF.

細胞質の体積の約30%は、生体分子で占められている。このような混み合った環境下では生体分子の立体構造、動態、反応性が希薄溶液環境と大きく異なることは従来の研究から示唆されているが、原子・分子レベルの解像度では十分に理解されていなかった。我々はバクテリア(マイコプラズマ・ジェニタリウム)の細胞質を原子解像度でモデリングし、「京」コンピュータ上で超並列分子動力学シミュレータGENESIS を用いた大規模計算を行い、細胞内の生体分子ダイナミクスを原子レベルで再現した。さらに、得られたデータから蛋白質-代謝物の相互作用や代謝物の細胞内動態を詳細に解析し、新たな知見を得た。

The free energy calculations of complex chemical and biological systems with molecular dynamics (MD) are inefficient due to multiple local minima separated by high-energy barriers. The minima can be escaped using an enhanced sampling method such as metadynamics, which apply bias (i.e., importance sampling) along a set of collective variables (CV), but the maximum number of CVs (or dimensions) is severely limited. We propose a high-dimensional bias potential method (NN2B) based on two machine learning algorithms: the nearest neighbor density estimator (NNDE) and the artificial neural network (ANN) for the bias potential approximation. The bias potential is constructed iteratively from short biased MD simulations accounting for correlation among CVs. Our method is capable of achieving ergodic sampling and calculating free energy of polypeptides with up to 8-dimensional bias potential.

Flexible fitting is a computational algorithm to derive a new conformational model that conforms to low-resolution experimental data by transforming a known structure. A common application is against data from cryo-electron microscopy to obtain conformational models in new functional states. The conventional flexible fitting algorithms cannot derive correct structures in some cases due to the complexity of conformational transitions. In this study, we show the importance of conformational ensemble in the refinement process by performing multiple fittings trials using a variety of different force constants. Application to simulated maps of Ca²⁺ ATPase and diphtheria toxin as well as experimental data of release factor 2 revealed that for these systems, multiple conformations with similar agreement with the density map exist and a large number of fitting trials are necessary to generate good models. Clustering analysis can be an effective approach to avoid over-fitting models. In addition, we show that an automatic adjustment of the biasing force constants during the fitting process, implemented as replica-exchange scheme, can improve the success rate.

Parallelization of molecular dynamics (MD) simulation is essential for investigating conformational dynamics of large biological systems, such as ribosomes, viruses, and multiple proteins in cellular environments. To improve efficiency in the parallel computation, we have to reduce the amount of data transfer between processors by introducing domain decomposition schemes. Also, it is important to optimize the computational balance between real-space non-bonded interactions and reciprocal-space interactions for long-range electrostatic interactions. Here, we introduce a novel parallelization scheme for large-scale MD simulations on massively parallel supercomputers consisting of only CPUs. We make use of a multiple program/multiple data (MPMD) approach for separating the real-space and reciprocal-space computations on different processors. We also utilize the r-RESPA multiple time step integrator on the framework of the MPMD approach in an efficient way: when the reciprocal-space computations are skipped in r-RESPA, processors assigned for them are utilized for half of the real-space computations. The new scheme allows us to use twice as many as processors that are available in the conventional single program approach. The best performances of all-atom MD simulations for 1 million (STMV), 8.5 million (8_STMV), and 28.8 million (27_STMV) atom systems on K computer are 65, 36, and 24 ns/day, respectively. The MPMD scheme can accelerate 23.4, 10.2, and 9.2 ns/day from the maximum performance of single-program approach for STMV, 8_STMV, and 27_STMV systems, respectively, which correspond to 57%, 39%, and 60% speed up. This suggests significant speedups by increasing the number of processors without losing parallel computational efficiency.

Molecular dynamics (MD) simulation of a N-glycan in solution is challenging because of high-energy barriers of the glycosidic linkages, functional group rotational barriers, and numerous intra- and intermolecular hydrogen bonds. In this study, we apply different enhanced conformational sampling approaches, namely, metadynamics (MTD), the replica-exchange MD (REMD), and the recently proposed replica state exchange MTD (RSE-MTD), to a N-glycan in solution and compare the conformational sampling efficiencies of the approaches. MTD helps to cross the high-energy barrier along the ω angle by utilizing a bias potential, but it cannot enhance sampling of the other degrees of freedom. REMD ensures moderate-energy barrier crossings by exchanging temperatures between replicas, while it hardly crosses the barriers along ω. In contrast, RSE-MTD succeeds to cross the high-energy barrier along ω as well as to enhance sampling of the other degrees of freedom. We tested two RSE-MTD schemes: in one scheme, 64 replicas were simulated with the bias potential along ω at different temperatures, while simulations of four replicas were performed with the bias potentials for different CVs at 300 K. In both schemes, one unbiased replica at 300 K was included to compute conformational properties of the glycan. The conformational sampling of the former is better than the other enhanced sampling methods, while the latter shows reasonable performance without spending large computational resources. The latter scheme is likely to be useful when a N-glycan-attached protein is simulated.

The infrared spectrum of H⁺(H₂O)₄ recently observed in a wide spectral range has shown a series of bands in a range of 1700–2500 cm^–1, which can not be understood by the standard harmonic normal mode analysis. Here, we theoretically investigate the origin of these bands with a focus on (1) the possibility of coexistence of multiple isomers in the Eigen [H₃O⁺(H₂O)₃] and Zundel [H₅O₂⁺(H₂O)₂] forms and (2) the effect of anharmonic coupling that gives rise to nonzero intensities for overtones and combination bands. Anharmonic vibrational calculations are carried out for the Eigen and Zundel clusters by the second-order vibrational quasi-degenerate perturbation theory (VQDPT2) based on optimized coordinates. The anharmonic potential energy surface and the dipole moment surfaces are generated by a multiresolution approach combining one-dimensional (1D) grid potential functions derived from CCSD(T)-F12, 2D and 3D grid potential functions derived from B3LYP for important coupling terms, and a quartic force field derived from B3LYP for less important terms. The spectrum calculated for the Eigen cluster is in excellent agreement with the experiment, assigning the bands in the range of 1700–2500 cm^–1 to overtones and combination bands of a H₃O⁺ moiety in line with recent reports [J. Phys. Chem. A 2015, 119, 9425; Science 2016, 354 1131]. On the other hand, characteristic OH stretching bands of the Zundel cluster is found to be absent in the experimental spectrum. We therefore conclude that the experimental spectrum originates solely from the Eigen cluster. Nonetheless, the present calculation for the Eigen cluster poorly reproduces a band observed at 1765 cm^–1. A possible nature of this band is discussed.

We investigate the effects of protein crowder sizes on hydration structure and dynamics in macromolecular crowded systems by all-atom MD simulations. The crowded systems consisting of only small proteins showed larger total surface areas than those of large proteins at the same volume fractions. As a result, more water molecules were trapped within the hydration shells, slowing down water diffusion. The simulation results suggest that the protein crowder size is another factor to determine the effect of macromolecular crowding and to explain the experimental kinetic data of proteins and DNAs in the presence of crowding agents.