Atomistic detailed free‐energy landscape of intrinsically disordered protein studied by multi‐scale divide‐and‐conquer molecular dynamics simulation

Calcineurin (CaN) is a eukaryotic serine/threonine protein phosphatase activated by both Ca²⁺ and calmodulin (CaM), including intrinsically disordered region (IDR). The region undergoes folding into an α‐helix form in the presence Ca²⁺‐loaded CaM. To sample the ordered structure of the IDR by conventional all atom model (AAM) molecular dynamics (MD) simulation, the IDR and Ca²⁺‐loaded CaM must be simultaneously treated. However, it is time‐consuming task because the coupled folding and binding should include repeated binding and dissociation. Then, in this study, we propose novel multi‐scale divide‐and‐conquer MD (MSDC‐MD), which combines AAM‐MD and coarse‐grained model MD (CGM‐MD). To speed up the conformation sampling, MSDC‐MD simulation first treats the IDR by CGM to sample conformations from wide conformation space; then, multiple AAM‐MD in a limited area is initiated using the resultant CGM conformation, which is reconstructed by homology modeling method. To investigate performance, we sampled the ordered conformation of the IDR using MSDC‐MD; the root‐mean‐square distance (RMSD) with respect to the experimental structure was 2.23 Å.     

Characterizing amyloid‐beta protein misfolding from molecular dynamics simulations with explicit water

Extracellular deposition of amyloid‐beta (Aβ) protein, a fragment of membrane glycoprotein called β‐amyloid precursor transmembrane protein (βAPP), is the major characteristic for the Alzheimer's disease (AD). However, the structural and mechanistic information of forming Aβ protein aggregates in a lag phase in cell exterior has been still limited. Here, we have performed multiple all‐atom molecular dynamics simulations for physiological 42‐residue amyloid‐beta protein (Aβ42) in explicit water to characterize most plausible aggregation‐prone structure (APS) for the monomer and the very early conformational transitions for Aβ42 protein misfolding process in a lag phase. Monitoring the early sequential conformational transitions of Aβ42 misfolding in water, the APS for Aβ42 monomer is characterized by the observed correlation between the nonlocal backbone H‐bond formation and the hydrophobic side‐chain exposure. Characteristics on the nature of the APS of Aβ42 allow us to provide new insight into the higher aggregation propensity of Aβ42 over Aβ40, which is in agreement with the experiments. On the basis of the structural features of APS, we propose a plausible aggregation mechanism from APS of Aβ42 to form fibril. The structural and mechanistic observations based on these simulations agree with the recent NMR experiments and provide the driving force and structural origin for the Aβ42 aggregation process to cause AD. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

Confined Water in Amyloid‐Competent Oligomers of the Prion Protein

Conformational switching of the prion protein into the abnormal form involves the formation of (obligatory) molten‐oligomers that mature into ordered amyloid fibrils. The role of water in directing the course of amyloid formation remains poorly understood. Here, we show that the mobility of the water molecules within the on‐pathway oligomers is highly retarded. The water relaxation time within the oligomers was estimated to be ≈1 ns which is about three orders of magnitude slower than the bulk water and resembles the characteristics of (trapped) nano‐confined water. We propose that the coalescence of these obligatory oligomers containing trapped water is entropically favored because of the release of ordered water molecules in the bulk milieu and results in the sequestration of favorable inter‐chain amyloid contacts via nucleated conformational conversion. The dynamic role of water in protein aggregation will have much broader implications in a variety of protein misfolding diseases.     

Fragment‐based quantum mechanical calculation of protein–protein binding affinities

The electrostatically embedded generalized molecular fractionation with conjugate caps (EE‐GMFCC) method has been successfully utilized for efficient linear‐scaling quantum mechanical (QM) calculation of protein energies. In this work, we applied the EE‐GMFCC method for calculation of binding affinity of Endonuclease colicin–immunity protein complex. The binding free energy changes between the wild‐type and mutants of the complex calculated by EE‐GMFCC are in good agreement with experimental results. The correlation coefficient (R) between the predicted binding energy changes and experimental values is 0.906 at the B3LYP/6‐31G*‐D level, based on the snapshot whose binding affinity is closest to the average result from the molecular mechanics/Poisson–Boltzmann surface area (MM/PBSA) calculation. The inclusion of the QM effects is important for accurate prediction of protein–protein binding affinities. Moreover, the self‐consistent calculation of PB solvation energy is required for accurate calculations of protein–protein binding free energies. This study demonstrates that the EE‐GMFCC method is capable of providing reliable prediction of relative binding affinities for protein–protein complexes. © 2018 Wiley Periodicals, Inc.     

The QM‐MM interface for CHARMM‐deMon

We present a new QM/MM interface for fast and efficient simulations of organic and biological molecules. The CHARMM/deMon interface has been developed and tested to perform minimization and atomistic simulations for multi‐particle systems. The current features of this QM/MM interface include readability for molecular dynamics, tested compatibility with Free Energy Perturbation simulations (FEP) using the dual topology/single coordinate method. The current coupling scheme uses link atoms, but further extensions of the code to incorporate other available schemes are planned. We report the performance of different levels of theory for the treatment of the QM region, while the MM region was represented by a classical force‐field (CHARMM27) or a polarizable force‐field based on a simple Drude model. The current QM/MM implementation can be coupled to the dual‐thermostat method and the VV2 integrator to run molecular dynamics simulations. © 2009 Wiley Periodicals, Inc. J Comput Chem, 2010     

Enhanced stability of the model mini‐protein in amino acid ionic liquids and their aqueous solutions

Using molecular dynamics simulations, the structure of model mini‐protein was thoroughly characterized in the imidazolium‐based amino acid ionic liquids and their aqueous solutions. Complete substitution of water by organic cations and anions further results in hindered conformational flexibility of the mini‐protein. This observation suggests that amino acid‐based ionic liquids are able to defend proteins from thermally induced denaturation. We show by means of radial distributions that the mini‐protein is efficiently solvated by both solvents due to a good mutual miscibility. Amino acid‐based anions prevail in the first coordination sphere of positively charged sites of the mini‐protein whereas water molecules prevail in the first coordination sphere of negatively charged sites of the mini‐protein. © 2015 Wiley Periodicals, Inc.     

Free energies of solvation in the context of protein folding: Implications for implicit and explicit solvent models

Implicit solvent models for biomolecular simulations have been developed to use in place of more expensive explicit models; however, these models make many assumptions and approximations that are likely to affect accuracy. Here, the changes in free energies of solvation upon folding of several fast folding proteins are calculated from previously run μs–ms simulations with a number of implicit solvent models and compared to the values needed to be consistent with the explicit solvent model used in the simulations. In the majority of cases, there is a significant and substantial difference between the values calculated from the two approaches that is robust to the details of the calculations. These differences could only be remedied by selecting values for the model parameters—the internal dielectric constant for the polar term and the surface tension coefficient for the nonpolar term—that were system‐specific or physically unrealistic. We discuss the potential implications of our findings for both implicit and explicit solvent simulations. © 2015 Wiley Periodicals, Inc.     

A noniterative mean-field QM/MM-type approach with a linear response approximation toward an efficient free-energy evaluation

Mean-field treatment of solvent provides an efficient technique to investigate chemical processes in solution in quantum mechanics/molecular mechanics (QM/MM) framework. In the algorithm, an iterative calculation is required to obtain the self-consistency between QM and MM regions, which is a time-consuming step. In the present study, we have proposed a noniterative approach by introducing a linear response approximation (LRA) into the solvation term in the one-electron part of Fock matrix in a hybrid approach between molecular-orbital calculations and a three-dimensional (3D) integral equation theory for molecular liquids (multicenter molecular Ornstein–Zernike self-consistent field [MC-MOZ-SCF]; Kido et al., J. Chem. Phys. 2015, 143, 014103). To save the computational time, we have also developed a fast method to generate electrostatic potential map near solute and the solvation term in Fock matrix, using Fourier transformation (FT) and real spherical harmonics expansion (RSHE). To numerically validate the LRA and FT-RSHE method, we applied the present approach to water, carbonic acid, and their ionic species in aqueous solution. Molecular properties of the solutes were evaluated by the present approach with four different types of initial wave functions and compared with those by the original (MC-MOZ-SCF). We found that an initial wave function considering solvation effects is needed to appropriately reproduce the properties by MC-MOZ-SCF. Furthermore, a benchmark test for 32 solute molecules was performed to evaluate the accuracy of the present approach for solvation free energy (SFE) and measure the speedup ratio for MC-MOZ-SCF. The error of SFE for MC-MOZ-SCF does not correlate with the SFE but increases in proportion to the electronic reorganization energy. Similar to water and carbonic acid, an initial wave function with solvation effects is also important to make the error small. From the averaged speed up ratio, the present approach is 13.5 times faster than MC-MOZ-SCF. © 2019 Wiley Periodicals, Inc.     

桥函数在统计力学积分方程理论3d-RISM-HNC的应用和溶剂化自由能计算的改进

把氢-桥函数和氧-桥函数应用于统计力学积分方程理论的三维的参考作用点-超链模型(3d-RISM-HNC)中,用以改进极性和非极性溶质的水溶液的热力学性质的计算.用三维和二维图形考察了溶剂水分子的氢原子和氧原子的桥函数在改进溶剂作用点的平均密度分布函数〈ρ_H(r)〉和〈ρ_O(r)〉,和平均超额化学势〈Δμ(r)〉的计算的效果.计算结果表明,氢桥函数和氧桥函数极大地改进了3d-RISM-HNC 方法的精度,把这一方法提高到定量和半定量的水平.研究表明,溶质分子的作用点的超额化学势的径向分布函数〈Δμ(r)〉比平均密度分布函数〈ρ_s(r)〉能够更灵敏地反映桥函数的改进效果.研究表明,为提高3d-RISM-HNC 方法的精度,需要进一步改进桥函数的函数形式和优化其中的参数.     

Transformation between α‐helix and β‐sheet structures of one and two polyglutamine peptides in explicit water molecules by replica‐exchange molecular dynamics simulations

Aggregation of polyglutamine peptides with β‐sheet structures is related to some important neurodegenerative diseases such as Huntington's disease. However, it is not clear how polyglutamine peptides form the β‐sheets and aggregate. To understand this problem, we performed all‐atom replica‐exchange molecular dynamics simulations of one and two polyglutamine peptides with 10 glutamine residues in explicit water molecules. Our results show that two polyglutamine peptides mainly formed helix or coil structures when they are separated, as in the system with one‐polyglutamine peptide. As the interpeptide distance decreases, the intrapeptide β‐sheet structure sometimes appear as an intermediate state, and finally the interpeptide β‐sheets are formed. We also find that the polyglutamine dimer tends to form the antiparallel β‐sheet conformations rather than the parallel β‐sheet, which is consistent with previous experiments and a coarse‐grained molecular dynamics simulation. © 2014 Wiley Periodicals, Inc.     

Spectroscopic Evidence for Polymorphic Aggregates Formed by Amyloid‐β Fragments

Understanding the structure of amyloid‐β (Aβ) aggregates is a key step towards elucidating the pathology of Alzheimer’s disease. In this work, three fragments of the Aβ_1–42 protein, Aβ_1–25 (DAEFRHDSGYEVHHQKLVFFAEDVG), Aβ_25–35 (GSNKGAIIGLM), and Aβ_33–42 (GLMVGGVVIA), were synthesized, and their aggregated structures were examined by linear infrared spectroscopy in the amide‐I (mainly the C?O stretching) region. The structures of the formed aggregates were found to be both sequence and pH dependent. The results suggest that instead of forming matured fibrils, as in the case of full‐length Aβ_1–42, both Aβ_1–25 and Aβ_33–42 form a mixture of threadlike β‐sheet fibril, soluble β‐sheet oligomer, and random coil structures. The β‐sheet conformations were found to be mainly antiparallel for the former and both parallel and antiparallel for the latter. However, the Aβ_25–35 fragment was found to form assembled fibrils containing predominantly parallel β‐sheets. The conformation and morphology of the aggregates were also confirmed by circular dichroism measurements and transmission electron microscopy. Factors influencing the structures of the aggregates formed by the Aβ fragments were discussed.     

Folding Dynamics of 10‐Residue β‐Hairpin Peptide Chignolin

Short peptides that fold into β‐hairpins are ideal model systems for investigating the mechanism of protein folding because their folding process shows dynamics typical of proteins. We performed folding, unfolding, and refolding molecular dynamics simulations (total of 2.7 μs) of the 10‐residue β‐hairpin peptide chignolin, which is the smallest β‐hairpin structure known to be stable in solution. Our results revealed the folding mechanism of chignolin, which comprises three steps. First, the folding begins with hydrophobic assembly. It brings the main chain together; subsequently, a nascent turn structure is formed. The second step is the conversion of the nascent turn into a tight turn structure along with interconversion of the hydrophobic packing and interstrand hydrogen bonds. Finally, the formation of the hydrogen‐bond network and the complete hydrophobic core as well as the arrangement of side‐chain–side‐chain interactions occur at approximately the same time. This three‐step mechanism appropriately interprets the folding process as involving a combination of previous inconsistent explanations of the folding mechanism of the β‐hairpin, that the first event of the folding is formation of hydrogen bonds and the second is that of the hydrophobic core, or vice versa.     

Parameterization of the Hamiltonian Dielectric Solvent (HADES) Reaction‐Field Method for the Solvation Free Energies of Amino Acid Side‐Chain Analogs

Optimization of the Hamiltonian dielectric solvent (HADES) method for biomolecular simulations in a dielectric continuum is presented with the goal of calculating accurate absolute solvation free energies while retaining the model’s accuracy in predicting conformational free‐energy differences. The solvation free energies of neutral and polar amino acid side‐chain analogs calculated by using HADES, which may optionally include nonpolar contributions, were optimized against experimental data to reach a chemical accuracy of about 0.5 kcal mol^?1. The new parameters were evaluated for charged side‐chain analogs. The HADES results were compared with explicit‐solvent, generalized Born, Poisson–Boltzmann, and QM‐based methods. The potentials of mean force (PMFs) between pairs of side‐chain analogs obtained by using HADES and explicit‐solvent simulations were used to evaluate the effects of the improved parameters optimized for solvation free energies on intermolecular potentials.     

Interaction identification of Zif268 and TATAZF proteins with GC‐/AT‐rich DNA sequence: A theoretical study

Molecular dynamics (MD) simulations for Zif268 (a zinc‐finger‐protein binding specifically to the GC‐rich DNA)‐d(A₁G₂C₃G₄T₅G₆G₇G₈C₉A₁₀C₁₁)₂ and TATA_ZF (a zinc‐finger‐protein recognizing the AT‐rich DNA)‐d(A₁C₂G₃C₄T₅A₆T₇A₈A₉A₁₀A₁₁G₁₂G₁₃)₂ complexes have been performed for investigating the DNA binding affinities and specific recognitions of zinc fingers to GC‐rich and AT‐rich DNA sequences. The binding free energies for the two systems have been further analyzed by using the molecular mechanics Poisson‐Boltzmann surface area (MM‐PBSA) method. The calculations of the binding free energies reveal that the affinity energy of Zif268‐DNA complex is larger than that of TATA_ZF‐DNA one. The affinity between the zinc‐finger‐protein and DNA is mainly driven by more favorable van‐der‐Waals and nonpolar/solvation interactions in both complexes. However, the affinity energy difference of the two binding systems is mainly caused by the difference of van‐der‐Waals interactions and entropy components. The decomposition analysis of MM‐PBSA free energies on each residue of the proteins predicts that the interactions between the residues with the positive charges and DNA favor the binding process; while the interactions between the residues with the negative charges and DNA behave in the opposite way. The interhydrogen‐bonds at the protein‐DNA interface and the induced intrafinger hydrogen bonds between the residues of protein for the Zif268‐DNA complex have been identified at some key contact sites. However, only the interhydrogen‐bonds between the residues of protein and DNA for TATA_ZF‐DNA complex have been found. The interactions of hydrogen‐bonds, electrostatistics and van‐der‐Waals type at some new contact sites have been identified. Moreover, the recognition characteristics of the two studied zinc‐finger‐proteins have also been discussed. © 2010 Wiley Periodicals, Inc. J Comput Chem, 2011     

A highly efficient hybrid method for calculating the hydration free energy of a protein

We develop a new method for calculating the hydration free energy (HFE) of a protein with any net charge. The polar part of the energetic component in the HFE is expressed as a linear combination of four geometric measures (GMs) of the protein structure and the generalized Born (GB) energy plus a constant. The other constituents in the HFE are expressed as linear combinations of the four GMs. The coefficients (including the constant) in the linear combinations are determined using the three‐dimensional reference interaction site model (3D‐RISM) theory applied to sufficiently many protein structures. Once the coefficients are determined, the HFE and its constituents of any other protein structure are obtained simply by calculating the four GMs and GB energy. Our method and the 3D‐RISM theory give perfectly correlated results. Nevertheless, the computation time required in our method is over four orders of magnitude shorter.     

The role of protein “Stability patches” in molecular recognition: A case study of the human growth hormone‐receptor complex

Dynamic characteristics of protein surfaces are among the factors determining their functional properties, including their potential participation in protein‐protein interactions. The presence of clusters of static residues—“stability patches” (SPs)—is a characteristic of protein surfaces involved in intermolecular recognition. The mechanism, by with SPs facilitate molecular recognition, however, remains unclear. Analyzing the surface dynamic properties of the growth hormone and of its high‐affinity variant we demonstrated that reshaping of the SPs landscape may be among the factors accountable for the improved affinity of this variant to the receptor. We hypothesized that SPs facilitate molecular recognition by moderating the conformational entropy of the unbound state, diminishing enthalpy–entropy compensation upon binding, and by augmenting the favorable entropy of desolvation. SPs mapping emerges as a valuable tool for investigating the structural basis of the stability of protein complexes and for rationalizing experimental approaches, such as affinity maturation, aimed at improving it. © 2015 Wiley Periodicals, Inc.     

Net charge changes in the calculation of relative ligand‐binding free energies via classical atomistic molecular dynamics simulation

The calculation of binding free energies of charged species to a target molecule is a frequently encountered problem in molecular dynamics studies of (bio‐)chemical thermodynamics. Many important endogenous receptor‐binding molecules, enzyme substrates, or drug molecules have a nonzero net charge. Absolute binding free energies, as well as binding free energies relative to another molecule with a different net charge will be affected by artifacts due to the used effective electrostatic interaction function and associated parameters (e.g., size of the computational box). In the present study, charging contributions to binding free energies of small oligoatomic ions to a series of model host cavities functionalized with different chemical groups are calculated with classical atomistic molecular dynamics simulation. Electrostatic interactions are treated using a lattice‐summation scheme or a cutoff‐truncation scheme with Barker–Watts reaction‐field correction, and the simulations are conducted in boxes of different edge lengths. It is illustrated that the charging free energies of the guest molecules in water and in the host strongly depend on the applied methodology and that neglect of correction terms for the artifacts introduced by the finite size of the simulated system and the use of an effective electrostatic interaction function considerably impairs the thermodynamic interpretation of guest‐host interactions. Application of correction terms for the various artifacts yields consistent results for the charging contribution to binding free energies and is thus a prerequisite for the valid interpretation or prediction of experimental data via molecular dynamics simulation. Analysis and correction of electrostatic artifacts according to the scheme proposed in the present study should therefore be considered an integral part of careful free‐energy calculation studies if changes in the net charge are involved. © 2013 The Authors Journal of Computational Chemistry Published by Wiley Periodicals, Inc.     

Update on phosphate and charged post‐translationally modified amino acid parameters in the GROMOS force field

In this study, we propose newly derived parameters for phosphate ions in the context of the GROMOS force field parameter sets. The non‐bonded parameters used up to now lead to a hydration free energy, which renders the dihydrogen phosphate ion too hydrophobic when compared to experimentally derived values, making a reparametrization of the phosphate moiety necessary. Phosphate species are of great importance in biomolecular simulations not only because of their crucial role in the backbone of nucleic acids but also as they represent one of the most important types of post‐translational modifications to protein side‐chains and are an integral part in many lipids. Our re‐parametrization of the free dihydrogen phosphate (H PO ) and three derivatives (methyl phosphate, dimethyl phosphate, and phenyl phosphate) leads, in conjunction with the previously updated charged side‐chains in the GROMOS parameter set 54A8, to new nucleic acid backbone parameters and a 54A8 version of the widely used GROMOS protein post‐translational modification parameter set. © 2017 Wiley Periodicals, Inc.     

Identification of thermostabilizing mutations for a membrane protein whose three‐dimensional structure is unknown

We recently developed a physics‐based method for identifying thermostabilizing mutations of a membrane protein. The method uses a free‐energy function F where the importance of translational entropy of hydrocarbon groups within the lipid bilayer is emphasized. All of the possible mutations can rapidly be examined. The method was illustrated for the adenosine A_2a receptor (A_2aR) whose three‐dimensional (3D) structure experimentally determined was utilized as the wild‐type structure. Nine single mutations and a double mutation predicted to be stabilizing or destabilizing were checked by referring to the experimental results: The success rate was remarkably high. In this work, we postulate that the 3D structure of A_2aR is unknown. We construct candidate models for the 3D structure using the homology modeling and select the model giving the lowest value to the change in F on protein folding. The performance achieved is only slightly lower than that in the recent work. © 2016 Wiley Periodicals, Inc.     

Localization and quantification of hydrophobicity: The molecular free energy density (MolFESD) concept and its application to sweetness recognition

A method for the localization, the quantification, and the analysis of hydrophobicity of a molecule or a molecular fragment is presented. It is shown that the free energy of solvation for a molecule or the transfer free energy from one solvent to another can be represented by a surface integral of a scalar quantity, the molecular free energy surface density (MolFESD), over the solvent accessible surface of that molecule. This MolFESD concept is based on a model approach where the solvent molecules are considered to be small in comparison to the solute molecule, and the solvent can be represented by a continuous medium with a given dielectric constant. The transfer energy surface density for a 1-octanol/water system is empirically determined employing a set of atomic increment contributions and distance dependent membership functions measuring the contribution of the increments to the surface value of the MolFESD. The MolFESD concept can be well used for the quantification of the purely hydrophobic contribution to the binding constants of molecule-receptor complexes. This is demonstrated with the sweeteners sucrose and sucralose and various halogen derivatives. Therein the relative sweetness, which is assumed to be proportional to the binding constant, nicely correlates to the surface integral over the positive, hydrophobic part of the MolFESD, indicating that the sweetness receptor can be characterized by a highly flexible hydrophobic pocket instead of a localized binding site.