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.
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 (FeIII) 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 FeIII(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 FeIII(NO) affect the structure around the heme of P450nor leading to an increase in νN–O and νFe–NO.
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 Ca2+-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 G37G38 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 Ca2+ 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+(H2O)4 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 [H3O+(H2O)3] and Zundel [H5O2+(H2O)2] 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 H3O+ 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.
Biological macromolecules function in highly crowded cellular environments. The structure and dynamics of proteins and nucleic acids are well characterized in vitro, but in vivo crowding effects remain unclear. Using molecular dynamics simulations of a comprehensive atomistic model cytoplasm we found that protein-protein interactions may destabilize native protein structures, whereas metabolite interactions may induce more compact states due to electrostatic screening. Protein-protein interactions also resulted in significant variations in reduced macromolecular diffusion under crowded conditions, while metabolites exhibited significant two-dimensional surface diffusion and altered protein-ligand binding that may reduce the effective concentration of metabolites and ligands in vivo. Metabolic enzymes showed weak non-specific association in cellular environments attributed to solvation and entropic effects. These effects are expected to have broad implications for the in vivo functioning of biomolecules. This work is a first step towards physically realistic in silico whole-cell models that connect molecular with cellular biology.
The cytoplasm of a cell is crowded with many different kinds of macromolecules. The macromolecular crowding affects the thermodynamics and kinetics of biological reactions in a living cell, such as protein folding, association, and diffusion. Theoretical and simulation studies using simplified models focus on the essential features of the crowding effects and provide a basis for analyzing experimental data. In most of the previous studies on the crowding effects, a uniform crowder size is assumed, which is in contrast to the inhomogeneous size distribution of macromolecules in a living cell. Here, we evaluate the free energy changes upon macromolecular association in a cell-like inhomogeneous crowding system via a theory of hard-sphere fluids and free energy calculations using Brownian dynamics trajectories. The inhomogeneous crowding model based on 41 different types of macromolecules represented by spheres with different radii mimics the physiological concentrations of macromolecules in the cytoplasm of Mycoplasma genitalium. The free energy changes of macromolecular association evaluated by the theory and simulations were in good agreement with each other. The crowder size distribution affects both specific and nonspecific molecular associations, suggesting that not only the volume fraction but also the size distribution of macromolecules are important factors for evaluating in vivo crowding effects. This study relates in vitro experiments on macromolecular crowding to in vivo crowding effects by using the theory of hard-sphere fluids with crowder-size heterogeneity.
Lipid/water interaction is essential for many biological processes. The water structure at the nonionic lipid interface remains little known, and there is no scope of a priori prediction of water orientation at nonionic interfaces, either. Here, we report our study combining advanced nonlinear spectroscopy and molecular dynamics simulation on the water orientation at the ceramide/water interface. We measured χ(2) spectrum in the OH stretch region of ceramide/isotopically diluted water interface using heterodyne-detected vibrational sum-frequency generation spectroscopy and found that the interfacial water prefers an overall hydrogen-up orientation. Molecular dynamics simulation indicates that this preferred hydrogen-up orientation of water is determined by a delicate balance between hydrogen-up and hydrogen-down orientation induced by lipid–water and intralipid hydrogen bonds. This mechanism also suggests that water orientation at neutral lipid interfaces depends highly on the chemical structure of the lipid headgroup, in contrast to the charged lipid interfaces where the net water orientation is determined solely by the charge of the lipid headgroup.
The graphics processing unit (GPU) has become a popular computational platform for molecular dynamics (MD) simulations of biomolecules. A significant speedup in the simulations of small- or medium-size systems using only a few computer nodes with a single or multiple GPUs has been reported. Because of GPU memory limitation and slow communication between GPUs on different computer nodes, it is not straightforward to accelerate MD simulations of large biological systems that contain a few million or more atoms on massively parallel supercomputers with GPUs. In this study, we develop a new scheme in our MD software, GENESIS, to reduce the total computational time on such computers. Computationally intensive real-space nonbonded interactions are computed mainly on GPUs in the scheme, while less intensive bonded interactions and communication-intensive reciprocal-space interactions are performed on CPUs. On the basis of the midpoint cell method as a domain decomposition scheme, we invent the single particle interaction list for reducing the GPU memory usage. Since total computational time is limited by the reciprocal-space computation, we utilize the RESPA multiple time-step integration and reduce the CPU resting time by assigning a subset of nonbonded interactions on CPUs as well as on GPUs when the reciprocal-space computation is skipped. We validated our GPU implementations in GENESIS on BPTI and a membrane protein, porin, by MD simulations and an alanine-tripeptide by REMD simulations. Benchmark calculations on TSUBAME supercomputer showed that an MD simulation of a million atoms system was scalable up to 256 computer nodes with GPUs.
An accurate theoretical prediction of the vibrational spectrum of polypeptides remains to be a challenge due to (1) their conformational flexibility and (2) non-negligible anharmonic effects. The former makes the search for conformers that contribute to the spectrum difficult, and the latter requires an expensive, quantum mechanical calculation for both electrons and vibrations. Here, we propose a new theoretical approach, which implements an enhanced conformational sampling by the replica-exchange molecular dynamics method, a structural clustering to identify distinct conformations, and a vibrational structure calculation by the second-order vibrational quasi-degenerate perturbation theory (VQDPT2). A systematic mode-selection scheme is developed to reduce the cost of VQDPT2 and the generation of a potential energy surface by the electronic structure calculation. The proposed method is applied to a pentapeptide, SIVSF-NH2, for which the infrared spectrum has recently been measured in the gas phase with high resolution in the OH and NH stretching region. The theoretical spectrum of the lowest energy conformer is obtained with a mean absolute deviation of 11.2 cm-1 from the experimental spectrum. Furthermore, the NH stretching frequencies of the five lowest energy conformers are found to be consistent with the literature values measured for small peptides with a similar secondary structure. Therefore, the proposed method is a promising way to analyze the vibrational spectrum of polypeptides.
Sphingomyelin (SM) is a major sphingolipid in mammalian cells that forms specific lipid domains in combination with cholesterol (Chol). Using molecular-dynamics simulation and density functional theory calculation, we identified a characteristic Raman band of SM at ∼1643 cm-1 as amide I of the SM cluster. Experimental results indicate that this band is sensitive to the hydration of SM and the presence of Chol. We showed that this amide I Raman band can be utilized to examine the membrane distribution of SM. Similarly to SM, ceramide phosphoethanolamine (CerPE) exhibited an amide I Raman band in almost the same region, although CerPE lacks three methyl groups in the phosphocholine moiety of SM. In contrast to SM, the amide I band of CerPE was not affected by Chol, suggesting the importance of the methyl groups of SM in the SM-Chol interaction.
Collective variables (CVs) are often used in molecular dynamics simulations based on enhanced sampling algorithms to investigate large conformational changes of a protein. The choice of CVs in these simulations is essential because it affects simulation results and impacts the free-energy profile, the minimum free-energy pathway (MFEP), and the transition-state structure. Here we examine how many CVs are required to capture the correct transition-state structure during the open-to-close motion of adenylate kinase using a coarse-grained model in the mean forces string method to search the MFEP. Various numbers of large amplitude principal components are tested as CVs in the simulations. The incorporation of local coordinates into CVs, which is possible in higher dimensional CV spaces, is important for capturing a reliable MFEP. The Bayesian measure proposed by Best and Hummer is sensitive to the choice of CVs, showing sharp peaks when the transition-state structure is captured. We thus evaluate the required number of CVs needed in enhanced sampling simulations for describing protein conformational changes.
Bacterial pathogens or cancer cells can acquire multidrug resistance, which causes serious clinical problems. In cells with multidrug resistance, various drugs or antibiotics are extruded across the cell membrane by multidrug transporters. The multidrug and toxic compound extrusion (MATE) transporter is one of the five families of multidrug transporters. MATE from Pyrococcus furiosus uses H+ to transport a substrate from the cytoplasm to the outside of a cell. Crystal structures of MATE from P. furiosus provide essential information on the relevant H+-binding sites (D41 and D184). Hybrid quantum mechanical/molecular mechanical simulations and continuum electrostatic calculations on the crystal structures predict that D41 is protonated in one structure (Straight) and, both D41 and D184 protonated in another (Bent). All-atom molecular dynamics simulations suggest a dynamic equilibrium between the protonation states of the two aspartic acids and that the protonation state affects hydration in the substrate binding cavity and lipid intrusion in the cleft between the N- and C-lobes. This hypothesis is examined in more detail by quantum mechanical/molecular mechanical calculations on snapshots taken from the molecular dynamics trajectories. We find the possibility of two proton transfer (PT) reactions in Straight: the 1st PT takes place between side-chains D41 and D184 through a transient formation of low-barrier hydrogen bonds and the 2nd through another H+ from the headgroup of a lipid that intrudes into the cleft resulting in a doubly protonated (both D41 and D184) state. The 1st PT affects the local hydrogen bond network and hydration in the N-lobe cavity, which would impinge on the substrate-binding affinity. The 2nd PT would drive the conformational change from Straight to Bent. This model may be applicable to several prokaryotic H+-coupled MATE multidrug transporters with the relevant aspartic acids.
Contacts with neighboring molecules in protein crystals inevitably restrict the internal motions of intrinsically flexible proteins. The resultant clear electron densities permit model building, as crystallographic snapshot structures. Although these still images are informative, they could provide biased pictures of the protein motions. If the mobile parts are located at a site lacking direct contacts in rationally designed crystals, then the amplitude of the movements can be experimentally analyzed. We propose a fusion protein method, to create crystal contact-free space (CCFS) in protein crystals and to place the mobile parts in the CCFS. Conventional model building fails when large amplitude motions exist. In this study, the mobile parts appear as smeared electron densities in the CCFS, by suitable processing of the X-ray diffraction data. We applied the CCFS method to a highly mobile presequence peptide bound to the mitochondrial import receptor, Tom20, and a catalytically relevant flexible segment in the oligosaccharyltransferase, AglB. These two examples demonstrated the general applicability of the CCFS method to the analysis of the spatial distribution of motions within protein molecules.
The asymmetric outer membrane of Gram-negative bacteria is formed of the inner leaflet with phospholipids and the outer leaflet with lipopolysaccharides (LPS). Outer membrane protein F (OmpF) is a trimeric porin responsible for the passive transport of small molecules across the outer membrane of Escherichia coli. Here, we report the impact of different levels of heterogeneity in LPS environments on the structure and dynamics of OmpF using all-atom molecular dynamics simulations. The simulations provide insight into the flexibility and dynamics of LPS components that are highly dependent on local environments, with lipid A being the most rigid and O-antigen being the most flexible. Increased flexibility of O-antigen polysaccharides is observed in heterogeneous LPS systems, where the adjacent O-antigen repeating units are weakly interacting and thus more dynamic, compared to homogeneous LPS systems in which LPS interacts strongly with each other with limited overall flexibility due to dense packing. The model systems were validated by comparing molecular-level details of interactions between OmpF surface residues and LPS core sugars with experimental data, establishing the importance of LPS core oligosaccharides in shielding OmpF surface epitopes recognized by monoclonal antibodies. There are LPS environmental influences on the move- ment of bulk ions (K+ and Cl-), but the ion selectivity of OmpF is mainly affected by bulk ion concentration.
This paper reviews various enhanced conformational sampling methods and explicit/implicit solvent/membrane models, as well as their recent applications to the exploration of the structure and dynamics of membranes and membrane proteins. Molecular dynamics simulations have become an essential tool to investigate biological problems, and their success relies on proper molecular models together with efficient conformational sampling methods. The implicit representation of solvent/membrane environments is reasonable approximation to the explicit all-atom models, considering the balance between computational cost and simulation accuracy. Implicit models can be easily combined with replica-exchange molecular dynamics methods to explore a wider conformational space of a protein. Other molecular models and enhanced conformational sampling methods are also briefly discussed. As application examples, we introduce recent simulation studies of glycophorin A, phospholamban, amyloid precursor protein, and mixed lipid bilayers and discuss the accuracy and efficiency of each simulation model and method. This article is part of a Special Issue entitled: Membrane Proteins edited by J.C. Gumbart and Sergei Noskov.
Three-dimensional Fast Fourier Transform (3D FFT) plays an important role in a wide variety of computer simulations and data analyses, including molecular dynamics (MD) simulations. In this study, we develop hybrid (MPI+OpenMP) parallelization schemes of 3D FFT based on two new volumetric decompositions, mainly for the particle mesh Ewald (PME) calculation in MD simulations. In one scheme, (1d_Alltoall), five all-to-all communications in one dimension are carried out, and in the other, (2d_Alltoall), one two-dimensional all-to-all communication is combined with two all-to-all communications in one dimension. 2d_Alltoall is similar to the conventional volumetric decomposition scheme. We performed benchmark tests of 3D FFT for the systems with different grid sizes using a large number of processors on the K computer in RIKEN AICS. The two schemes show comparable performances, and are better than existing 3D FFTs. The performances of 1d_Alltoall and 2d_Alltoall depend on the supercomputer network system and number of processors in each dimension. There is enough leeway for users to optimize performance for their conditions. In the PME method, short-range real-space interactions as well as long-range reciprocal-space interactions are calculated. Our volumetric decomposition schemes are particularly useful when used in conjunction with the recently developed midpoint cell method for short-range interactions, due to the same decompositions of real and reciprocal spaces. The 1d_Alltoall scheme of 3D FFT takes 4.7 ms to simulate one MD cycle for a virus system containing more than 1 million atoms using 32,768 cores on the K computer.
The bacterial SecYEG translocon functions as a conserved protein-conducting channel. Conformational transitions of SecYEG allow protein translocation across the membrane without perturbation of membrane permeability. Here, we report the crystal structures of intact SecYEG at 2.7-Å resolution and of peptide-bound SecYEG at 3.6-Å resolution. The higher-resolution structure revealed that the cytoplasmic loop of SecG covers the hourglass-shaped channel, which was confirmed to also occur in the membrane by disulfide bond formation analysis and molecular dynamics simulation. The cytoplasmic loop may be involved in protein translocation. In addition, the previously unknown peptide-bound crystal structure of SecYEG implies that interactions between the cytoplasmic side of SecY and signal peptides are related to lateral gate opening at the first step of protein translocation. These SecYEG structures therefore provide a number of structural insights into the Sec machinery for further study.
Large conformational changes of multidomain proteins are difficult to simulate using all-atom molecular dynamics (MD) due to the slow time scale. We show that a simple modification of the structure-based coarse-grained (CG) model enables a stable and efficient MD simulation of those proteins. “Motion Tree”, a tree diagram that describes conformational changes between two structures in a protein, provides information on rigid structural units (domains) and the magnitudes of domain motions. In our new CG model, which we call the DoME (domain motion enhanced) model, interdomain interactions are defined as being inversely proportional to the magnitude of the domain motions in the diagram, whereas intradomain interactions are kept constant. We applied the DoME model in combination with the Go model to simulations of adenylate kinase (AdK). The results of the DoME–Go simulation are consistent with an all-atom MD simulation for 10 μs as well as known experimental data. Unlike the conventional Go model, the DoME–Go model yields stable simulation trajectories against temperature changes and conformational transitions are easily sampled despite domain rigidity. Evidently, identification of domains and their interfaces is useful approach for CG modeling of multidomain proteins.
Infrared (IR) and Raman spectra of a sphingomyelin (SM) bilayer have been calculated for the amide I, II and A modes and the double-bonded CC stretching mode by a weight averaged approach, based on an all-atom molecular dynamics (MD) simulation and a vibrational structure calculation. Representative structures and statistical weights of SM clusters connected by hydrogen bonds (HBs) are observed in MD trajectories. After constructing smaller fragments from the SM clusters, the vibrational spectra of the target modes were calculated by normal mode analysis with a correction for anharmonicity, using density functional theory. The final IR and Raman spectra of a SM bilayer were obtained as the weight averages over all SM clusters. The calculated Raman spectrum is in excellent agreement with a recent measurement, providing a clear assignment of the peak in question observed at 1643 cm-1 to the amide I modes of a SM bilayer. The analysis of the IR spectrum has also revealed that the amide bands are sensitive to the water content inside the membrane, since their band positions are strongly modulated by the HB between SM and water molecules. The present study suggests that the amide I band serves as a marker to identify the formation of SM clusters, and opens a new way to detect lipid rafts in the biological membrane.
Data assimilation is a statistical method designed to improve the quality of numerical simulations in combination with real observations. Here, we develop a sequential data assimilation method that incorporates one-dimensional time-series data of smFRET (single-molecule Förster resonance energy transfer) photon-counting into conformational ensembles of biomolecules derived from “replicated” molecular dynamics (MD) simulations. A particle filter using a large number of “replicated” MD simulations with a likelihood function for smFRET photon-counting data is employed to screen the conformational ensembles that match the experimental data. We examine the performance of the method using emulated smFRET data and coarse-grained (CG) MD simulations of a dye-labeled polyproline-20. The method estimates the dynamics of the end-to-end distance from smFRET data as well as revealing that of latent conformational variables. The particle filter is also able to correct model parameter dependence in CG MD simulations. We discuss the applicability of the method to real experimental data for conformational dynamics of biomolecules.
Metadynamics (MTD) is a powerful enhanced sampling method for systems with rugged energy landscapes. It constructs a bias potential in a predefined collective variable (CV) space to overcome barriers between metastable states. In bias-exchange MTD (BE-MTD), multiple replicas approximate the CV space by exchanging bias potentials (replica conditions) with the Metropolis–Hastings (MH) algorithm. We demonstrate that the replica-exchange rates and the convergence of free energy estimates of BE-MTD are improved by introducing the infinite swapping (IS) or the Suwa-Todo (ST) algorithms. Conceptually, IS and ST perform transitions in a replica state space rather than exchanges in a replica condition space. To emphasize this, the proposed scheme is called the replica state exchange MTD (RSE-MTD). Benchmarks were performed with alanine polypeptides in vacuum and water. For the systems tested in this work, there is no significant performance difference between IS and ST.
GENESIS (Generalized‐Ensemble Simulation System) is a new software package for molecular dynamics (MD) simulations of macromolecules. It has two MD simulators, called ATDYN and SPDYN. ATDYN is parallelized based on an atomic decomposition algorithm for the simulations of all‐atom force‐field models as well as coarse‐grained Go‐like models. SPDYN is highly parallelized based on a domain decomposition scheme, allowing large‐scale MD simulations on supercomputers. Hybrid schemes combining OpenMP and MPI are used in both simulators to target modern multicore computer architectures. Key advantages of GENESIS are (1) the highly parallel performance of SPDYN for very large biological systems consisting of more than one million atoms and (2) the availability of various REMD algorithms (T‐REMD, REUS, multi‐dimensional REMD for both all‐atom and Go‐like models under the NVT, NPT, NPAT, and NPγT ensembles). The former is achieved by a combination of the midpoint cell method and the efficient three‐dimensional Fast Fourier Transform algorithm, where the domain decomposition space is shared in real‐space and reciprocal‐space calculations. Other features in SPDYN, such as avoiding concurrent memory access, reducing communication times, and usage of parallel input/output files, also contribute to the performance. We show the REMD simulation results of a mixed (POPC/DMPC) lipid bilayer as a real application using GENESIS. GENESIS is released as free software under the GPLv2 licence and can be easily modified for the development of new algorithms and molecular models.
Sarco(endo)plasmic reticulum Ca2+-ATPase transports two Ca2+ per ATP-hydrolyzed across biological membranes against a large concentration gradient by undergoing large conformational changes. Structural studies with X-ray crystallography revealed functional roles of coupled motions between the cytoplasmic domains and the transmembrane helices in individual reaction steps. Here, we employed "Motion Tree (MT)," a tree diagram that describes a conformational change between two structures, and applied it to representative Ca2+-ATPase structures. MT provides information of coupled rigid-body motions of the ATPase in individual reaction steps. Fourteen rigid structural units, "common rigid domains (CRDs)" are identified from seven MTs throughout the whole enzymatic reaction cycle. CRDs likely act as not only the structural units, but also the functional units. Some of the functional importance has been newly revealed by the analysis. Stability of each CRD is examined on the morphing trajectories that cover seven conformational transitions. We confirmed that the large conformational changes are realized by the motions only in the flexible regions that connect CRDs. The Ca2+-ATPase efficiently utilizes its intrinsic flexibility and rigidity to response different switches like ligand binding/dissociation or ATP hydrolysis. The analysis detects functional motions without extensive biological knowledge of experts, suggesting its general applicability to domain movements in other membrane proteins to deepen the understanding of protein structure and function.
The pathways by which small molecules (substrates or inhibitors) access active sites are a key aspect of the function of enzymes and other proteins. A key problem in designing or altering such proteins is to identify sites for mutation that will have the desired effect on the substrate transport properties. While specific access channels have been invoked in the past, molecular simulations suggest that multiple routes are possible, complicating the analysis. This complexity, however, can be captured by a Markov State Model (MSM) of the ligand diffusion process. We have developed a sensitivity analysis of the resulting rate matrix, which identifies the locations where mutations should have the largest effect on the diffusive on rate. We apply this method to myoglobin, which is the best characterized example both from experiment and simulation. We validate the approach by translating the sensitivity parameter obtained from this method into the CO binding rates in myoglobin upon mutation, resulting in a semi-quantitative correlation with experiments. The model is further validated against an explicit simulation for one of the experimental mutants.
Replica-exchange molecular dynamics (REMD) method is one of the enhanced conformational sampling algorithms widely applied in computational biophysics and biochemistry. In the method, molecular dynamics (MD) simulations for multiple replicas of a target system are performed simultaneously and independently. Every few steps, temperatures or other parameters are exchanged between a pair of replicas according to the modified Metropolis criteria. Replica-Exchange INterface (REIN) is interface software for REMD simulations. It utilizes existing MD software without modification and performs the exchanges of replicas. In this article, we introduce the software design of REIN and demonstrate its applicability through benchmark simulations of alanine pentapeptide in explicit water, as well as the free-energy analysis of N-glycan and Tom20-presequence complex in solution.
A model for the cytoplasm of Mycoplasma genitalium is presented that integrates data from a variety of sources into a physically and biochemically consistent model. Based on gene annotations, core genes expected to be present in the cytoplasm were determined and a metabolic reaction network was reconstructed. The set of cytoplasmic genes and metabolites from the predicted reactions were assembled into a comprehensive atomistic model consisting of proteins with predicted structures, RNA, protein/RNA complexes, metabolites, ions, and solvent. The resulting model bridges between atomistic and cellular scales, between physical and biochemical aspects, and between structural and systems views of cellular systems and is meant as a starting point for a variety of simulation studies.
The dielectric properties of proton(H)–deuteron(D) mixed crystals of the hydrogen-bonded material 5-Bromo-9-hydroxyphenalenone are studied using a novel path integral Monte Carlo (PIMC) method that takes account of the dipole induction effect depending on the relative proton configurations in the surrounding molecules. The induced dipole is evaluated using the fragment molecular orbital method with electron correlation included by second-order Møller–Plesset perturbation theory and long-range corrected density functional theory. The results show a greater influence of CH…O intermolecular weak hydrogen bonding on the induction than for results evaluated with the Hartree–Fock method. The induction correction is incorporated into the PIMC simulations with a model Hamiltonian that consists of long-range dipolar interactions and a transverse term describing proton tunneling. The relationship between the calculated phase transition temperature and H/D mixing ratio is consistent with the experimental phase diagram, indicating that the balance between the proton tunneling and the collective ordering is appropriately described.
Ordering of water structures near the surface of biological membranes has been recently extensively studied using interface-selective techniques like vibrational sum frequency generation (VSFG) spectroscopy. The detailed structures of interface water have emerged for charged lipids, but those for neutral zwitterionic lipids remain obscure. We analyze an all-atom molecular dynamics (MD) trajectory of a hydrated 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine bilayer to characterize the orientation of interface waters in different chemical environments. The structure and dynamics of interfacial waters strongly depend on both their vertical position along the bilayer normal as well as vicinal lipid charged groups. Water orientation in the vicinity of phosphate groups is opposite to that around choline groups. The results are consistent with observed VSFG spectra and demonstrate that a mosaic of water orientation structures exists on the surface of a neutral zwitterionic phospholipid bilayer, reflecting rapid water exchange and the influence of local chemical environments.
A simple methodology for calculating state average energies in the context of vibrational self-consistent field (VSCF) is suggested. The suggested state average energy is employed in the optimization of coordinates for anharmonic vibrational wave function calculations where an orthogonal matrix of transformation between normal and optimized coordinates is variationally optimized. The convergence to the exact limit for approximate vibrational configuration interaction and vibrational coupled cluster wave functions is studied, comparing the performance of state average optimized coordinates, ground state optimized coordinates and standard normal coordinates. Exploratory calculations are presented for water, formaldehyde, the water dimer and trimer, and ethylene.
Adenosine triphosphate (ATP) is an indispensable energy source in cells. In a wide variety of biological phenomena like glycolysis, muscle contraction/relaxation, and active ion transport, chemical energy released from ATP hydrolysis is converted to mechanical forces to bring about large-scale conformational changes in proteins. Investigation of structure–function relationships in these proteins by molecular dynamics (MD) simulations requires modeling of ATP in solution and ATP bound to proteins with accurate force-field parameters. In this study, we derived new force-field parameters for the triphosphate moiety of ATP based on the high-precision quantum calculations of methyl triphosphate. We tested our new parameters on membrane-embedded sarcoplasmic reticulum Ca2+-ATPase and four soluble proteins. The ATP-bound structure of Ca2+-ATPase remains stable during MD simulations, contrary to the outcome in shorter simulations using original parameters. Similar results were obtained with the four ATP-bound soluble proteins. The new force-field parameters were also tested by investigating the range of conformations sampled during replica-exchange MD simulations of ATP in explicit water. Modified parameters allowed a much wider range of conformational sampling compared with the bias toward extended forms with original parameters. A diverse range of structures agrees with the broad distribution of ATP conformations in proteins deposited in the Protein Data Bank. These simulations suggest that the modified parameters will be useful in studies of ATP in solution and of the many ATP-utilizing proteins.
Newly synthesized membrane proteins must be accurately inserted into the membrane, folded and assembled for proper functioning. The protein YidC inserts its substrates into the membrane, thereby facilitating membrane protein assembly in bacteria; the homologous proteins Oxa1 and Alb3 have the same function in mitochondria and chloroplasts, respectively. In the bacterial cytoplasmic membrane, YidC functions as an independent insertase and a membrane chaperone in cooperation with the translocon SecYEG. Here we present the crystal structure of YidC from Bacillus halodurans, at 2.4 Å resolution. The structure reveals a novel fold, in which five conserved transmembrane helices form a positively charged hydrophilic groove that is open towards both the lipid bilayer and the cytoplasm but closed on the extracellular side. Structure-based in vivo analyses reveal that a conserved arginine residue in the groove is important for the insertion of membrane proteins by YidC. We propose an insertion mechanism for single-spanning membrane proteins, in which the hydrophilic environment generated by the groove recruits the extracellular regions of substrates into the low-dielectric environment of the membrane.
The use of variationally optimized coordinates, which minimize the vibrational self-consistent field (VSCF) ground state energy with respect to orthogonal transformations of the coordinates, has recently been shown to improve the convergence of vibrational configuration interaction (VCI) towards the exact full VCI [K. Yagi, M. Keçeli, and S. Hirata, J. Chem. Phys. 137, 204118 (2012)]. The present paper proposes an incorporation of optimized coordinates into the vibrational coupled cluster (VCC), which has in the past been shown to outperform VCI in approximate calculations where similar restricted state spaces are employed in VCI and VCC. An embarrassingly parallel algorithm for variational optimization of coordinates for VSCF is implemented and the resulting coordinates and potentials are introduced into a VCC program. The performance of VCC in optimized coordinates (denoted oc-VCC) is examined through pilot applications to water, formaldehyde, and a series of water clusters (dimer, trimer, and hexamer) by comparing the calculated vibrational energy levels with those of the conventional VCC in normal coordinates and VCI in optimized coordinates. For water clusters, in particular, oc-VCC is found to gain orders of magnitude improvement in the accuracy, exemplifying that the combination of optimized coordinates localized to each monomer with the size-extensive VCC wave function provides a supreme description of systems consisting of weakly interacting sub-systems.
We have developed a new hybrid (MPI+OpenMP) parallelization scheme for molecular dynamics (MD) simulations by combining a cell-wise version of the midpoint method with pair-wise Verlet lists. In this scheme, which we call the midpoint cell method, simulation space is divided into subdomains, each of which is assigned to a MPI processor. Each subdomain is further divided into small cells. The interaction between two particles existing in different cells is computed in the subdomain containing the midpoint cell of the two cells where the particles reside. In each MPI processor, cell pairs are distributed over OpenMP threads for shared memory parallelization. The midpoint cell method keeps the advantages of the original midpoint method, while filtering out unnecessary calculations of midpoint checking for all the particle pairs by single midpoint cell determination prior to MD simulations. Distributing cell pairs over OpenMP threads allows for more efficient shared memory parallelization compared with distributing atom indices over threads. Furthermore, cell grouping of particle data makes better memory access, reducing the number of cache misses. The parallel performance of the midpoint cell method on the K computer showed scalability up to 512 and 32,768 cores for systems of 20,000 and 1 million atoms, respectively. One MD time step for long-range interactions could be calculated within 4.5 ms even for a 1 million atoms system with particle-mesh Ewald electrostatics.
A perturbative extension to optimized coordinate vibrational self-consistent field (oc-VSCF) is proposed based on the quasi-degenerate perturbation theory (QDPT). A scheme to construct the degenerate space (P space) is developed, which incorporates degenerate configurations and alleviates the divergence of perturbative expansion due to localized coordinates in oc-VSCF (e.g., local O–H stretching modes of water). An efficient configuration selection scheme is also implemented, which screens out the Hamiltonian matrix element between the P space configuration (p) and the complementary Q space configuration (q) based on a difference in their quantum numbers (λpq = ∑s |ps − qs |). It is demonstrated that the second-order vibrational QDPT based on optimized coordinates (oc-VQDPT2) smoothly converges with respect to the order of the mode coupling, and outperforms the conventional one based on normal coordinates. Furthermore, an improved, fast algorithm is developed for optimizing the coordinates. First, the minimization of the VSCF energy is conducted in a restricted parameter space, in which only a portion of pairs of coordinates is selectively transformed. A rational index is devised for this purpose, which identifies the important coordinate pairs to mix from others that may remain unchanged based on the magnitude of harmonic coupling induced by the transformation. Second, a cubic force field (CFF) is employed in place of a quartic force field, which bypasses intensive procedures that arise due to the presence of the fourth-order force constants. It is found that oc-VSCF based on CFF together with the pair selection scheme yields the coordinates similar in character to the conventional ones such that the final vibrational energy is affected very little while gaining an order of magnitude acceleration. The proposed method is applied to ethylene and trans-1,3-butadiene. An accurate, multi-resolution potential, which combines the MP2 and coupled-cluster with singles, doubles, and perturbative triples level of electronic structure theory, is generated and employed in the oc-VQDPT2 calculation to obtain the fundamental tones as well as selected overtones/combination tones coupled to the fundamentals through the Fermi resonance. The calculated frequencies of ethylene and trans-1,3-butadiene are found to be in excellent agreement with the experimental values with a mean absolute error of 8 and 9 cm-1, respectively.
Structural information of a transmembrane (TM) helix dimer is useful in understanding molecular mechanisms of important biological phenomena such as signal transduction across the cell membrane. Here, we describe an umbrella sampling (US) scheme for predicting the structure of a TM helix dimer in implicit membrane using the interhelical crossing angle and the TM–TM relative rotation angles as the reaction coordinates. This scheme conducts an efficient conformational search on TM–TM contact interfaces, and its robustness is tested by predicting the structures of glycophorin A (GpA) and receptor tyrosine kinase EphA1 (EphA1) TM dimers. The nuclear magnetic resonance (NMR) structures of both proteins correspond to the global free-energy minimum states in their free-energy landscapes. In addition, using the landscape of GpA as a reference, we also examine the protocols of temperature replica-exchange molecular dynamics (REMD) simulations for structure prediction of TM helix dimers in implicit membrane. A wide temperature range in REMD simulations, for example, 250–1000 K, is required to efficiently obtain a free-energy landscape consistent with the US simulations. The interhelical crossing angle and the TM–TM relative rotation angles can be used as reaction coordinates in multidimensional US and be good measures for conformational sampling of REMD simulations.
We have investigated effects of salt ions on folding events of a helical miniprotein chicken villin headpiece subdomain HP36. Low concentrations of ions alter electrostatic interactions between charged groups of a protein and can change the populations of conformers. Here, we compare two data sets of folding simulations of HP36 in explicit water solvent with or without ions. For efficient sampling of the conformational space of HP36, the multicanonical replica-exchange molecular dynamics method was employed. Our analyses suggest that salt alters salt-bridging nature of the protein at later stages of folding at room temperature. Especially, more nonnative, nonlocal salt bridges are formed at near-native conformations in pure water. Our analyses also show that such salt-bridge formation hinders the fully native hydrophobic-core packing at the final stages of folding.
Conformational sampling is fundamentally important for simulating complex bio-molecular systems. The generalized-ensemble algorithm, especially the temperature replica-exchange molecular dynamics method (T-REMD), is one of the most powerful methods to explore structures of bio-molecules such as proteins, nucleic acids, carbohydrates, and also of lipid membranes. T-REMD simulations have focused on soluble proteins rather than membrane proteins or lipid bilayers, because explicit membranes do not keep their structural integrity at high temperature. Here, we propose a new generalized-ensemble algorithm for membrane systems, which we call the surface-tension REMD method. Each replica is simulated in the NPγT ensemble, and surface tensions in a pair of replicas are exchanged at certain intervals to enhance conformational sampling of the target membrane system. We test the method on two biological membrane systems: a fully hydrated DPPC (1,2-dipalmitoyl-sn-glycero-3-phosphatidylcholine) lipid bilayer and a WALP23–POPC (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine) membrane system. During these simulations, a random walk in surface tension space is realized. Large-scale lateral deformation (shrinking and stretching) of the membranes takes place in all of the replicas without collapse of the lipid bilayer structure. There is accelerated lateral diffusion of DPPC lipid molecules compared with conventional MD simulation, and a much wider range of tilt angle of the WALP23 peptide is sampled due to large deformation of the POPC lipid bilayer and through peptide-lipid interactions. Our method could be applicable to a wide variety of biological membrane systems.
An increasing number of studies are aimed at modeling cellular environments in a comprehensive and realistic fashion. A major challenge in these efforts is how to bridge spatial and temporal scales over many orders of magnitude. Furthermore, there are additional challenges in integrating different aspects ranging from questions about biomolecular stability in crowded environments to the description of reactive processes on cellular scales. In this review, recent studies with models of biomolecules in cellular environments at different levels of detail are discussed in terms of their strengths and weaknesses. In particular, atomistic models, implicit representations of cellular environments, coarse-grained and spheroidal models of biomolecules, as well as the inclusion of reactive processes via reaction–diffusion models are described. Furthermore, strategies for integrating the different models into a comprehensive description of cellular environments are discussed.
The major bottleneck in molecular dynamics (MD) simulations of biomolecules exist in the calculation of pairwise nonbonded interactions like Lennard-Jones and long-range electrostatic interactions. Particle-mesh Ewald (PME) method is able to evaluate long-range electrostatic interactions accurately and quickly during MD simulation. However, the evaluation of energy and gradient includes time-consuming inverse square roots and complementary error functions. To avoid such time-consuming operations while keeping accuracy, we propose a new lookup table for short-range interaction in PME by defining energy and gradient as a linear function of inverse distance squared. In our lookup table approach, densities of table points are inversely proportional to squared pair distances, enabling accurate evaluation of energy and gradient at small pair distances. Regardless of the inverse operation here, the new lookup table scheme allows fast pairwise nonbonded calculations owing to efficient usage of cache memory.
Proton-dependent oligopeptide transporters (POTs) are major facilitator superfamily (MFS) proteins that mediate the uptake of peptides and peptide-like molecules, using the inwardly directed H+ gradient across the membrane. The human POT family transporter peptide transporter 1 is present in the brush border membrane of the small intestine and is involved in the uptake of nutrient peptides and drug molecules such as β-lactam antibiotics. Although previous studies have provided insight into the overall structure of the POT family transporters, the question of how transport is coupled to both peptide and H+ binding remains unanswered. Here we report the high-resolution crystal structures of a bacterial POT family transporter, including its complex with a dipeptide analog, alafosfalin. These structures revealed the key mechanistic and functional roles for a conserved glutamate residue (Glu310) in the peptide binding site. Integrated structural, biochemical, and computational analyses suggested a mechanism for H+-coupled peptide symport in which protonated Glu310 first binds the carboxyl group of the peptide substrate. The deprotonation of Glu310 in the inward open state triggers the release of the bound peptide toward the intracellular space and salt bridge formation between Glu310 and Arg43 to induce the state transition to the occluded conformation.
The formation of like-ion pairs, Na+–Na+ and Cl-–Cl-, in aqueous solution was studied by high-level ab initio methods, classical molecular dynamics (MD), QM/TIP5P, and QM/EFP MD (quantum mechanics/effective fragment potential molecular dynamics). Ab initio calculations on model clusters revealed that the Na+(H2O)nNa+ (n = 2–4) clusters were significantly more stabilized by bridged waters than the corresponding Cl-(H2O)nCl- clusters. QM/EFP MD simulations in solution also predicted a clear local minimum near 3.6 Å only for the Na+–Na+ pair, suggesting that Na+–Na+ pairs may be more likely to form than Cl-–Cl- pairs in solution. Analysis of the hydration structures further showed that two-water bridged Na+–Na+ pairs were dominant at the local minimum. The preferred formation of Na+ like-ion pairs in solution appeared to come from significant short-range effects, in particular, charge delocalization (polarization) between the bridged oxygen p and the vacant valence Na+ orbitals.
An efficient parallel implementation of QM/MM-based replica-exchange molecular dynamics (REMD) as well as umbrella samplings techniques was proposed by adopting the generalized distributed data interface (GDDI). Parallelization speed-up of 40.5 on 48 cores was achieved, making our QM/MM-MD engine a robust tool for studying complex chemical dynamics in solution. They were comparatively used to study the torsional isomerization of hydrogen peroxide in aqueous solution. All results by QM/MM-REMD and QM/MM umbrella sampling techniques yielded nearly identical potentials of mean force (PMFs) regardless of the particular QM theories for solute, showing that the overall dynamics are mainly determined by solvation. Although the entropic penalty of solvent rearrangements exists in cisoid conformers, it was found that both strong intermolecular hydrogen bonding and dipole–dipole interactions preferentially stabilize them in solution, reducing the torsional free-energy barrier at 0° by about 3 kcal/mol as compared to that in gas phase.
Tom20 is located at the outer membrane of mitochondria and functions as a receptor for the N-terminal presequence of mitochondrial-precursor proteins. Recently, three atomic structures of the Tom20-presequence complex were determined using X-ray crystallography and classified into A-, M-, and Y-poses in terms of their presequence-binding modes. Combined with biochemical and NMR data, a dynamic-equilibrium model between the three poses has been proposed. To investigate this mechanism in further detail, we performed all-atom molecular dynamics (MD) simulations and replica-exchange MD (REMD) simulations of the Tom20-presequence complex in explicit water. In the REMD simulations, one major distribution and another minor one were observed in the converged free-energy landscape at 300 K. In the major distribution, structures similar to A- and M-poses exist, whereas those similar to Y-pose are located in the minor one, suggesting that A-pose in solution is more stable than Y-pose. A k-means clustering algorithm revealed a new pose not yet obtained by X-ray crystallography. This structure has double salt bridges between Arg14′ in the presequence and Glu78 or Glu79 in Tom20 and can explain the binding affinity of the complex in previous pull-down assay experiments. Structural clustering and analyses of contacts between Tom20 and the presequence suggest smooth conformational changes from Y- to A-poses through low activation barriers. M-pose lies between Y- and A-poses as a metastable state. The REMD simulations thus provide insights into the energetics of the multiple-binding forms and help to detail the progressive conformational states in the dynamic-equilibrium model based on the experimental data.
The effect of cellular crowding environments on protein structure and stability is a key issue in molecular and cellular biology. The classical view of crowding emphasizes the volume exclusion effect that generally favors compact, native states. Here, results from molecular dynamics simulations and NMR experiments show that protein crowders may destabilize native states via protein–protein interactions. In the model system considered here, mixtures of villin head piece and protein G at high concentrations, villin structures become increasingly destabilized upon increasing crowder concentrations. The denatured states observed in the simulation involve partial unfolding as well as more subtle conformational shifts. The unfolded states remain overall compact and only partially overlap with unfolded ensembles at high temperature and in the presence of urea. NMR measurements on the same systems confirm structural changes upon crowding based on changes of chemical shifts relative to dilute conditions. An analysis of protein–protein interactions and energetic aspects suggests the importance of enthalpic and solvation contributions to the crowding free energies that challenge an entropic-centered view of crowding effects.
The nudged elastic band (NEB) and string methods are widely used to obtain the reaction path of chemical reactions and phase transitions. In these methods, however, it is difficult to define an accurate Lagrangian to generate the conservative forces. On the other hand, the constrained optimization with locally updated planes (CO-LUP) scheme defines target function properly and suitable for micro-iteration optimizations in quantum mechanical/molecular mechanical (QM/MM) systems, which uses the efficient second order QM optimization. However, the method does have problems of inaccurate estimation of reactions and inappropriate accumulation of images around the energy minimum. We introduce three modifications into the CO-LUP scheme to overcome these problems: (1) An improved tangent estimation of the reaction path, which is used in the NEB method, (2) redistribution of images using an energy-weighted interpolation before updating local tangents, and (3) reduction of the number of constraints, in particular translation/rotation constraints, for improved convergence. First, we test the method on the isomerization of alanine dipeptide without QM/MM calculation, showing that the method is comparable to the string method both in accuracy and efficiency. Next, we apply the method for defining the reaction paths of the rearrangement reaction catalyzed by chorismate mutase (CM) and of the phosphoryl transfer reaction catalyzed by cAMP-dependent protein kinase (PKA) using generalized hybrid orbital QM/MM calculations. The reaction energy barrier of CM is in high agreement with the experimental value. The path of PKA reveals that the enzyme reaction is associative and there is a late transfer of the substrate proton to Asp 166, which is in agreement with the recently published result using the NEB method.
The association process of NaCl in aqueous solution was studied by a combination of quantum mechanical calculations on NaCl(H2O)n (n = 1–6) clusters and quantum mechanical/effective fragment potential–molecular dynamics (QM/EFP-MD) simulations for NaCl in 292 EFP waters. The interionic hydration structures (IHSs) were topologically classified as "ring" (R), "half-bridge" (H), and "full-bridge" (F) types on the basis of the quantum mechanical calculations. Subsequent IHS analysis on QM/EFP-MD simulations revealed that the NaCl contact ion pair (CIP) mainly involved R type hydration structures while the solvent-separated ion pair (SSIP) was composed of two different groups of F-type hydration structures. Our IHS analysis also discovered H type hydration even at large separation interionic distances (7 Å), which is denoted as a dissociating ion pair (DIP). The analysis was able to reveal the most complete interionic structures and their reorganizations of the association process. A strong correlation between the IHSs and interionic distance suggests that not only the solvent reorganization but also the local IHS changes are equally important. Mechanistically, it is suggested that the conversion between ring-type and full-bridge hydration structures is the main rate-determining step of ion-pair association.
The effect of protein crowding on the structure and dynamics of water was examined from explicit solvent molecular dynamics simulations of a series of protein G and protein G/villin systems at different protein concentrations. Hydration structure was analyzed in terms of radial distribution functions, three-dimensional hydration sites, and preservation of tetrahedral coordination. Analysis of hydration dynamics focused on self-diffusion rates and dielectric constants as a function of crowding. The results show significant changes in both structure and dynamics of water under highly crowded conditions. The structure of water is altered mostly beyond the first solvation shell. Diffusion rates and dielectric constants are significantly reduced following linear trends as a function of crowding reflecting highly constrained water in crowded environments. The reduced dynamics of diffusion is expected to be strongly related to hydrodynamic properties of crowded cellular environments while the reduced dielectric constant under crowded conditions has implications for the stability of biomolecules in crowded environments. The results from this study suggest a prescription for modeling solvation in simulations of cellular environments.
RATIONALE
A central issue in glycan mass analysis is the ambiguity of structural assignments due to the heterogeneity and complexity of glycan structures. Ion mobility mass spectrometry (IM-MS) has the potential to separate isomeric glycans depending on their unique collisional cross section especially when coupled with hydrophilic interaction liquid chromatography (HILIC).
METHODS
Ten pyridylaminated biantennary N-glycans including isomeric structures were measured by electrospray ionization quadrupole-time-of-flight mass spectrometry with an ion mobility phase. We investigated which adduct ions would be suitable for good separation in the ion mobility phase. The differences in observed drift time of isomeric glycans were assessed by molecular dynamics (MD) simulations in vacuum. Connecting an HILIC system with IM-MS provided another, augmented separation mode.
RESULTS
By selecting doubly protonated precursor ion species, we succeeded in separating a pair of isomeric glycans in the ion mobility phase with reasonable resolution. MD simulations of monogalactosylated glycan isomers indicate that the galactosylated Man α1-3 branch preferentially folds back to the core chitobiose portion to form a compact structure. IM-MS combined with HILIC resulted in even clearer separation of isomeric glycans within 15 min.
CONCLUSIONS
A combination of IM-MS with an HILIC system is eminently suitable for the confident and rapid distinction of glycan structures within a defined mixture.
Replica-exchange molecular dynamics (REMD) method is one of the enhanced conformational sampling techniques in MD simulations of proteins or other systems with rugged-energy landscapes. In REMD method, copies of original simulation system at different temperatures are simulated separately and simultaneously. Every few steps, temperatures between neighboring replicas are exchanged if the Metropolis criteria for their instantaneous potential energies are satisfied. Due to its simplicity and high efficiency in parallel computers, the method has been applied to many biological problems including protein folding, aggregation, receptor-ligand binding, and so on. In the last ten years, continuous effort to improve sampling efficiency of REMD simulations for larger biological systems has been carried out by us and other theoretical scientists. In this review article, we introduce two different approaches in REMD simulations to reduce the computational cost. One is the multicanonical replica-exchange method (MUCAREM) for reducing the number of replicas. In this method, each replica has a different multicanonical weight factor and takes a flat energy distribution to cover a wider potential energy space. Another approach is to employ implicit solvent/membrane models for representing surrounding environments of target proteins in REMD simulations. We show two applications of proteinfolding simulations in explicit solvent using the former approach and a structural prediction of a transmembrane protein dimer using the latter. Finally, we discuss possibilities of REMD method to simulate a large-scale conformational change of protein systems using massively parallel supercomputers.
Nitric oxide reductases (NORs) are membrane proteins that catalyze the reduction of nitric oxide (NO) to nitrous oxide (N2O), which is a critical step of the nitrate respiration process in denitrifying bacteria. Using the recently determined first crystal structure of the cytochrome c-dependent NOR (cNOR) [Hino T, Matsumoto Y, Nagano S, Sugimoto H, Fukumori Y, et al. (2010) Structural basis of biological N2O generation by bacterial nitric oxide reductase. Science 330: 1666–70.], we performed extensive all-atom molecular dynamics (MD) simulations of cNOR within an explicit membrane/solvent environment to fully characterize water distribution and dynamics as well as hydrogen-bonded networks inside the protein, yielding the atomic details of functionally important proton channels. Simulations reveal two possible proton transfer pathways leading from the periplasm to the active site, while no pathways from the cytoplasmic side were found, consistently with the experimental observations that cNOR is not a proton pump. One of the pathways, which was newly identified in the MD simulation, is blocked in the crystal structure and requires small structural rearrangements to allow for water channel formation. That pathway is equivalent to the functional periplasmic cavity postulated in cbb3 oxidase, which illustrates that the two enzymes share some elements of the proton transfer mechanisms and confirms a close evolutionary relation between NORs and C-type oxidases. Several mechanisms of the critical proton transfer steps near the catalytic center are proposed.
Cytochrome c oxidase is a member of the haem copper oxidase superfamily (HCO). HCOs function as the terminal enzymes in the respiratory chain of mitochondria and aerobic prokaryotes, coupling molecular oxygen reduction to transmembrane proton pumping. Integral to the enzyme's function is the transfer of electrons from cytochrome c to the oxidase via a transient association of the two proteins. Electron entry and exit are proposed to occur from the same site on cytochrome c. Here we report the crystal structure of the caa3-type cytochrome oxidase from Thermus thermophilus, which has a covalently tethered cytochrome c domain. Crystals were grown in a bicontinuous mesophase using a synthetic short-chain monoacylglycerol as the hosting lipid. From the electron density map, at 2.36 Å resolution, a novel integral membrane subunit and a native glycoglycerophospholipid embedded in the complex were identified. Contrary to previous electron transfer mechanisms observed for soluble cytochrome c, the structure reveals the architecture of the electron transfer complex for the fused cupredoxin/cytochrome c domain, which implicates different sites on cytochrome c for electron entry and exit. Support for an alternative to the classical proton gate characteristic of this HCO class is presented.
Protein–glycan recognition regulates a wide range of biological and pathogenic processes. Conformational diversity of glycans in solution is apparently incompatible with specific binding to their receptor proteins. One possibility is that among the different conformational states of a glycan, only one conformer is utilized for specific binding to a protein. However, the labile nature of glycans makes characterizing their conformational states a challenging issue. All-atom molecular dynamics (MD) simulations provide the atomic details of glycan structures in solution, but fairly extensive sampling is required for simulating the transitions between rotameric states. This difficulty limits application of conventional MD simulations to small fragments like di- and tri-saccharides. Replica-exchange molecular dynamics (REMD) simulation, with extensive sampling of structures in solution, provides a valuable way to identify a family of glycan conformers. This article reviews recent REMD simulations of glycans carried out by us or other research groups and provides new insights into the conformational equilibria of N-glycans and their alteration by chemical modification. We also emphasize the importance of statistical averaging over the multiple conformers of glycans for comparing simulation results with experimental observables. The results support the concept of "conformer selection" in protein–glycan recognition.
The introduction of bisecting GlcNAc and core fucosylation in N-glycans is essential for fine functional regulation of glycoproteins. In this paper, the effect of these modifications on the conformational properties of N-glycans is examined at the atomic level by performing replica exchange molecular dynamics (REMD) simulations. We simulate four biantennary complex type N-glycans, namely, unmodified, two single substituted with either bisecting GlcNAc or core fucose, and disubstituted forms. By using REMD as an enhanced sampling technique, five distinct conformers in solution, each of which is characterized by its local orientation of the Manα1-6Man glycosidic linkage, are observed for all four N-glycans. The chemical modifications significantly change their conformational equilibria. The number of major conformers is reduced from five to two and from five to four upon the introduction of bisecting GlcNAc and core fucosylation; respectively, The population change is attributed to specific inter residue hydrogen bonds, including water mediated Ones. The experimental NMR data, including nuclear Overhauser enhancement and scalar J-coupling constant's, are well reproduced taking the multiple conformers into account. Our structural model supports the concept of "conformer selection"; which emphasize the, conformational flexibility of N-glycans in protein-glycan interactions.
We show that the hybrid quantum mechanical/effective fragment potential (QM/EFP) can be a very effective and practical quantum mechanical molecular dynamics method, when it is properly combined with well-developed traditional molecular dynamics (MD) techniques. QM/EFP-MD simulations on intra-molecular proton transfer of glycine with 290 EFP waters yielded accurate free energy change and reaction barrier of the zwitterion → neutral form conversion. Water rearrangements turned out to be the main driving force of the proton transfer.
Kaposi's sarcoma-associated herpesvirus (KSHV), a human tumor virus, encodes two homologous membrane-associated E3 ubiquitin ligases, modulator of immune recognition 1 (MIR1) and MIR2, to evade host immunity. Both MIR1 and MIR2 downregulate the surface expression of major histocompatibility complex class I (MHC I) molecules through ubiquitin-mediated endocytosis followed by lysosomal degradation. Since MIR2 additionally downregulates a costimulatory molecule (B7-2) and an integrin ligand (intercellular adhesion molecule 1 [ICAM-1]), MIR2 is thought to be a more important molecule for immune evasion than MIR1; however, the molecular basis of the MIR2 substrate specificity remains unclear. To address this issue, we determined which regions of B7-2 and MIR2 are required for MIR2-mediated B7-2 downregulation. Experiments with chimeras made by swapping domains between human B7-2 and CD8α, a non-MIR2 substrate, and between MIR1 and MIR2 demonstrated a significant contribution of the juxtamembrane (JM) region of B7-2 and the intertransmembrane (ITM) region of MIR2 to MIR2-mediated downregulation. Structure prediction and mutagenesis analyses indicate that Phe119 and Ser120 in the MIR2 ITM region and Asp244 in the B7-2 JM region contribute to the recognition of B7-2 by MIR2. This finding provides new insight into the molecular basis of substrate recognition by MIR family members.
Integrin α5β1 is a major cellular receptor for the extracellular matrix protein fibronectin and plays a fundamental role during mammalian development. A crystal structure of the α5β1 integrin headpiece fragment bound by an allosteric inhibitory antibody was determined at a 2.9-Å resolution both in the absence and presence of a ligand peptide containing the Arg-Gly-Asp (RGD) sequence. The antibody-bound β1 chain accommodated the RGD ligand with very limited structural changes, which may represent the initial step of cell adhesion mediated by nonactivated integrins. Furthermore, a molecular dynamics simulation pointed to an important role for Ca2+ in the conformational coupling between the ligand-binding site and the rest of the molecule. The RGD-binding pocket is situated at the center of a trenchlike exposed surface on the top face of α5β1 devoid of glycosylation sites. The structure also enabled the precise prediction of the acceptor residue for the auxiliary synergy site of fibronectin on the α5 subunit, which was experimentally confirmed by mutagenesis and kinetic binding assays.
The effect of cellular crowding was examined from molecular dynamics simulations of chymotrypsin inhibitor 2 (CI2) in the presence of either lysozyme or bovine serum albumin (BSA) crowder molecules as a complement to recent experimental studies of the same systems (Miklos, A. C.; Sarkar, M.; Wang, Y.; Pielak, G. J. J. Am. Chem. Soc. 2011, 133, 7116). The simulations confirm a destabilization and significantly slowed diffusion of CI2 in the presence of lysozyme and indicate that this observation is a result of extensive, non-specific protein-protein interactions between CI2 and lysozyme. CI2 interacts much less with BSA crowders corresponding to a weak effect of crowding. Energetic analysis suggests an overall favorable crowding free energy in the presence of lysozyme while weaker interactions with BSA appear to be unfavorable.
We have investigated flip-flop motions of diacylglycerol and ceramide in phospholipid bilayers using coarse-grained molecular dynamics simulations. In the simulations, flip-flop motions of diacylglycerol and ceramide in the DAPC membrane are slower than cholesterol. Rates correlate with the number of unsaturated bonds in the membrane phospholipids and hence with fluidity of membranes. These findings qualitatively agree with corresponding experimental data. Statistical analysis of the trajectories suggests that flip-flop can be approximated as a Poisson process. The rate of the transverse movement is influenced by depth of the polar head group in the membrane and extent of interaction with water.
The structure of quinol dependent nitric oxide reductase (qNOR) from Geobacillus stearothermophilus, which catalyzes the reduction of NO to produce a major ozone-depleting gas N2O, has been characterized at 2.5-Å resolution. The overall fold of qNOR is similar to that of cytochrome c dependent NOR (cNOR), and some structural features that are characteristic of cNOR, such as the Ca binding site and hydrophilic cytochrome c domain, are observed in qNOR, even though it harbours no heme c. In contrast to cNOR, structure-based mutagenesis and molecular dynamics simulation studies of qNOR suggest that a water channel from the cytoplasm can serve as a proton transfer pathway for the catalytic reaction. Further structural comparison of qNOR with cNOR and aerobic/micro-aerobic respiratory oxidases elucidates their evolutionary relationship and possible functional conversions.
All-atom molecular dynamics (MD) simulation has become a powerful research tool to investigate structural and dynamical properties of biological membranes and membrane proteins. The lipid structures of simple membrane systems in recent MD simulations are in good agreement with those obtained by experiments. However, for protein-membrane systems, the complexity of protein-lipid interactions makes investigation of lipid structure difficult. Although the area per lipid is one of the essential structural properties in membrane systems, the area in protein-membrane systems cannot be computed easily by conventional approaches like the Voronoi tessellation (VT) method. To overcome this limitation, we propose a new method combining the two-dimensional VT and Monte Carlo integration methods. This approach computes individual surface areas of lipid molecules not only in bulk lipids but also in proximity to membrane proteins. We apply the method to all-atom MD trajectories of the sarcoplasmic reticulum Ca2+-pump and the SecY protein-conducting channel. The calculated lipid surface area is in agreement with experimental values and consistent with other structural parameters of lipid bilayers. We also observe changes in the average area per lipid induced by the conformational transition of the SecY channel. Our method is particularly useful for examining equilibration of lipids around membrane proteins and for analyzing the time course of protein-lipid interactions.
