Inclusion of solvation free energy with molecular mechanics energy: alanyl dipeptide as a test case.

A combined force field of molecular mechanics and solvation free energy is tested by carrying out energy minimization and molecular dynamics on several conformations of the alanyl dipeptide. Our results are qualitatively consistent with previous experimental and computational studies, in that the addition of solvation energy stabilizes the C5 conformation of the alanyl dipeptide relative to the C7.     

Free energy surfaces of miniproteins with a betabetaalpha motif: replica exchange molecular dynamics simulation with an implicit solvation model

Designed miniproteins with a betabetaalpha motif, such as BBA5, 1FSD, and 1PSV can serve as a benchmark set to test the validity of all-atom force fields with computer simulation, because they contain all the basic structural elements in protein folding. Unfortunately, it was found that the standard all-atom force fields with the generalized Born (GB) implicit solvation model tend to produce distorted free energy surfaces for the betabetaalpha proteins, not only because energetically those proteins need to be described by more balanced weights of the alpha- and beta-strands, but also because the GB implicit solvation model suffers from overestimated salt bridge effects. In an attempt to resolve these problems, we have modified one of the standard all-atom force fields in conjunction with the GB model, such that each native state of the betabetaalpha proteins is in its free energy minimum state with reasonable energy barriers separating local minima. With this modified energy model, the free energy contour map in each protein was constructed from the replica exchange molecular dynamics REMD simulation. The resulting free energy surfaces are significantly improved in comparison with previous simulation results and consistent with general views on small protein folding behaviors with realistic topology and energetics of all three proteins.     

Complexity of free energy landscapes of peptides revealed by nonlinear principal component analysis

Employing the recently developed hierarchical nonlinear principal component analysis (NLPCA) method of Saegusa et al. (Neurocomputing 2004;61:57-70 and IEICE Trans Inf Syst 2005;E88-D:2242-2248), the complexities of the free energy landscapes of several peptides, including triglycine, hexaalanine, and the C-terminal beta-hairpin of protein G, were studied. First, the performance of this NLPCA method was compared with the standard linear principal component analysis (PCA). In particular, we compared two methods according to (1) the ability of the dimensionality reduction and (2) the efficient representation of peptide conformations in low-dimensional spaces spanned by the first few principal components. The study revealed that NLPCA reduces the dimensionality of the considered systems much better, than did PCA. For example, in order to get the similar error, which is due to representation of the original data of beta-hairpin in low dimensional space, one needs 4 and 21 principal components of NLPCA and PCA, respectively. Second, by representing the free energy landscapes of the considered systems as a function of the first two principal components obtained from PCA, we obtained the relatively well-structured free energy landscapes. In contrast, the free energy landscapes of NLPCA are much more complicated, exhibiting many states which are hidden in the PCA maps, especially in the unfolded regions. Furthermore, the study also showed that many states in the PCA maps are mixed up by several peptide conformations, while those of the NLPCA maps are more pure. This finding suggests that the NLPCA should be used to capture the essential features of the systems.     

Helix-capping interaction in λ cro protein: A free energy simulation analysis

The stability mutant Tyr-26 → Asp was studied in the Cro protein from bacteriophage λ using free energy molecular dynamics simulations. The mutant was calculated to be more stable than the wild type by 3.0 ± 1.7 kcal/mol/monomer, in reasonable agreement with experiment (1.4 kcal/mol/monomer). Moreover, the aspartic acid in the mutant was found to form a capping interation with the amino terminus of the third α-helix of Cro. The simulations were analyzed to understand better the source of the stability of this helix-capping interaction and to examine the results in light of previous explanations of stabilizing helix caps-namely, a model of local unsatisfied hydrogen bonds at the helix termini and the helix macro dipole model. Analysis of the simulations shows that the stabilizing effect of this charged helical cap is due both to favorable hydrogen bonds with backbone NH groups at the helix terminus and to favorable electrostatic interactions (but not hydrogen bonds) with their carbonyls (effectively the next row of local dipoles in the helix). However, electrostatic interactions are weak or negligible with backbone dipolar groups in the helix further away from the terminus. Moreover, the importance of other local electrostatic interactions with polar side chains near the helix terminus, which are neglected in most treatments of this effect, are shown to be important. Thus, the results support a model that is intermediate between the two previous explanations: both unsatisfied hydrogen bonds at the helix terminus and other, local preoriented dipolar groups stabilize the helix cap. These findings suggest that similar interactions with preoriented dipolar groups may be important for cooperativity in other charge–dipole interactions and may be employed to advantage for molecular design. © 1994 Wiley-Liss, Inc.     

Conservative and nonconservative mutations in proteins: anomalous mutations in a transport receptor analyzed by free energy and quantum chemical calculations.

Experimental studies on a bacterial sulfate receptor have indicated anomalous relative binding affinities for the mutations Ser130-->Cys,Ser130-->Gly, and Ser130-->Ala. The loss of affinity for sulfate in the former mutation was previously attributed to a greater steric effect on the part of the Cys side chain relative to the Ser side chain, whereas the relatively small loss of binding affinity for the latter two mutations was attributed to the loss of a single hydrogen bond. In this report we present quantum chemical and statistical thermodynamic studies of these mutations. Qualitative results from these studies indicate that for the Ser130-->Cys mutation the large decrease in binding affinity is in part caused by steric effects, but also significantly by the differential work required to polarize the Cys thiol group relative to the Ser hydroxyl group. The Gly mutant cobinds a water molecule in the same location as the Ser side chain resulting in a relatively small decrease in binding affinity. Results for the Ala mutant are in disagreement with experimental results but are likely to be limited by insufficient sampling of configuration space due to physical constraints applied during the simulation.     

Position-resolved free energy of solvation for amino acids in lipid membranes from molecular dynamics simulations

Studies of insertion and interactions of amino acids in lipid membranes are pivotal to our understanding of membrane protein structure and function. Calculating the insertion cost as a function of transmembrane helix sequence is thus an important step towards improved membrane protein prediction and eventually drug design. Here, we present position-dependent free energies of solvation for all amino acid analogs along the membrane normal. The profiles cover the entire region from bulk water to hydrophobic core, and were produced from all-atom molecular dynamics simulations. Experimental differences corresponding to mutations and costs for entire segments match experimental data well, and in addition the profiles provide the spatial resolution currently not available from experiments. Polar side-chains largely maintain their hydration and assume quite ordered conformations, which indicates the solvation cost is mainly entropic. The cost of solvating charged side-chains is not only significantly lower than for implicit solvation models, but also close to experiments, meaning these could well maintain their protonation states inside the membrane. The single notable exception to the experimental agreement is proline, which is quite expensive to introduce in vivo despite its hydrophobicity--a difference possibly explained by kinks making it harder to insert helices in the translocon.     

Enhanced sampling of molecular dynamics simulation of peptides and proteins by double coupling to thermal bath

Low sampling efficiency in conformational space is the well-known problem for conventional molecular dynamics. It greatly increases the difficulty for molecules to find the transition path to native state, and costs amount of CPU time. To accelerate the sampling, in this paper, we re-couple the critical degrees of freedom in the molecule to environment temperature, like dihedrals in generalized coordinates or nonhydrogen atoms in Cartesian coordinate. After applying to ALA dipeptide model, we find that this modified molecular dynamics greatly enhances the sampling behavior in the conformational space and provides more information about the state-to-state transition, while conventional molecular dynamics fails to do so. Moreover, from the results of 16 independent 100?ns simulations by the new method, it shows that trpzip2 has one-half chances to reach the naive state in all the trajectories, which is greatly higher than conventional molecular dynamics. Such an improvement would provide a potential way for searching the conformational space or predicting the most stable states of peptides and proteins.     

Docking study and free energy simulation of the complex between p53 DNA-binding domain and azurin

Molecular interaction between p53 tumor suppressor and the copper protein azurin (AZ) has been demonstrated to enhance p53 stability and hence antitumoral function, opening new perspectives in cancer treatment. While some experimental work has provided evidence for AZ binding to p53, no crystal structure for the p53-AZ complex was solved thus far. In this work the association between AZ and the p53 DNA-binding domain (DBD) was investigated by computational methods. Using a combination of rigid-body protein docking, experimental mutagenesis information, and cluster analysis 10 main p53 DBD-AZ binding modes were generated. The resulting structures were further characterized by molecular dynamics (MD) simulations and free energy calculations. We found that the highest scored docking conformation for the p53 DBD-AZ complex also yielded the most favorable free energy value. This best three-dimensional model for the complex was validated by using a computational mutagenesis strategy. In this structure AZ binds to the flexible L(1) and s(7)-s(8) loops of the p53 DBD and stabilizes them through protein-protein tight packing interactions, resulting in high degree of both surface matching and electrostatic complementarity.     

Immune-related proteins induced in the hemolymph after aseptic and septic injury differ in honey bee worker larvae and adults

We have employed the proteomic approach in combination with mass spectrometry to study the immune response of honey bee workers at different developmental stages. Analysis of the hemolymph proteins of noninfected, mock-infected and immune-challenged individuals by polyacrylamide gel electrophoresis showed differences in the protein profiles. We present evidence that in vitro reared honey bee larvae respond with a prominent humoral reaction to aseptic and septic injury as documented by the transient synthesis of the three antimicrobial peptides (AMPs) hymenoptaecin, defensin1, and abaecin. In contrast, young adult worker bees react with a broader spectrum of immune reactions that include the activation of prophenoloxidase and humoral immune responses. At least seven proteins appeared consistently in the hemolymph of immune-challenged bees, three of which are identical to the AMPs induced also in larvae. The other four, i.e., phenoloxidase (PO), peptidoglycan recognition protein-S2, carboxylesterase (CE), and an Apis-specific protein not assigned to any function (HP30), are induced specifically in adult bees and, with the exception of PO, are not expressed after aseptic injury. Structural features of CE and HP30, such as classical leucine zipper motifs, together with their strong simultaneous induction upon challenge with bacteria suggest an important role of the two novel bee-specific immune proteins in response to microbial infections.     

Convergence in peptide folding simulation: multiple trajectories of a potential AIDS pharmacophore

To examine the conformational properties in aqueous solution of a 15-residue peptide that is a potential pharmacophore for AIDS vaccine development, molecular dynamics simulations were performed in water starting from structures determined experimentally in three different organic solvents. Convergence characteristics of the simulation are examined in Cartesian and conformational spaces. In addition, novel analysis tools are employed including a multidimensional scaling method to represent the distance between trajectory frames. As these methods are based on a variety of physical parameters, they provide a useful cross-check on the structural convergence. Theoretical two-dimensional (2D) 1H-NMR spectra are also generated. These are superficially quite different in appearance, demonstrating that backbone similarities difficult to identify by visual inspection of 2D NMR data can be revealed using the methods described here.     

Interfacial water on crystalline silica: a comparative molecular dynamics simulation study

Understanding the properties of interfacial water at solid–liquid interfaces is important in a wide range of applications. Molecular dynamics is becoming a widespread tool for this purpose. Unfortunately, however, the results of such studies are known to strongly depend on the selection of force fields. It is, therefore, of interest to assess the extent by which the implemented force fields can affect the predicted properties of interfacial water. Two silica surfaces, with low and high surface hydroxyl density, respectively, were simulated implementing four force fields. These force fields yield different orientation and flexibility of surface hydrogen atoms, and also different interaction potentials with water molecules. The properties for interfacial water were quantified by calculating contact angles, atomic density profiles, surface density distributions, hydrogen bond density profiles and residence times for water near the solid substrates. We found that at low surface density of hydroxyl groups, the force field strongly affects the predicted contact angle, while at high density of hydroxyl groups, water wets all surfaces considered. From a molecular-level point of view, our results show that the position and intensity of peaks observed from oxygen and hydrogen atomic density profiles are quite different when different force fields are implemented, even when the simulated contact angles are similar. Particularly, the surfaces simulated by the CLAYFF force field appear to attract water more strongly than those simulated by the Bródka and Zerda force field. It was found that the surface density distributions for water strongly depend on the orientation of surface hydrogen atoms. In all cases, we found an elevated number of hydrogen bonds formed between interfacial water molecules. The hydrogen bond density profile does not depend strongly on the force field implemented to simulate the substrate, suggesting that interfacial water assumes the necessary orientation to maximise the number of water–water hydrogen bonds irrespectively of surface properties. Conversely, the residence time for water molecules near the interface strongly depends on the force field and on the flexibility of surface hydroxyl groups. Specifically, water molecules reside for longer times at contact with rigid substrates with high density of hydroxyl groups. These results should be considered when comparisons between simulated and experimental data are attempted.     

Molecular dynamics simulation of interleukin-2 and its complex and determination of the binding free energy

Interleukin-2 (IL-2) protein belongs to the signal modulator cytokine's family and therefore it is prevalent for immunological responses. It has been identified as a centrally important potential drug target for the inhibition of protein-protein interactions; so as to suppress the immunological responses associated with autoimmune, inflammatory and immunological diseases, and cancer. In the present work, we have performed two independent 100?ns of molecular dynamics (MD) simulations on the apo IL-2 protein and its ligand-bound complex (with a potent inhibitor FRG), to study the effect of inhibitor binding on the dynamics and stability of the protein. The calculation of binding free energy via post-processing end state method of Molecular Mechanics Poisson Boltzmann Surface Area (MM-PBSA) and Molecular Mechanics Generalised Born Surface Area (MM-GBSA) has inferred a good correlation in accordance with the already reported experimental data, demonstrating that the free energy of binding calculated by the two methods has no significant difference. The investigation of individual components of free energy revealed that the association of IL-2 protein with FRG ligand is primarily driven by the van der Waals energy contribution that represents the non-polar/hydrophobic energy contribution as dominant in this case of ligand binding.     

Crabrolin,a natural antimicrobial peptide: structural properties

A joint application of experimental and computational approaches has revealed the exceptionally high attitude of crabrolin, a 13‐residue peptide with sequence FLPLILRKIVTAL‐NH₂, to adopt alpha‐helix conformation not only in membrane‐mimicking solvents but also in the presence of a not negligible amount of water. Our study shows that this propensity essentially resides in the intrinsic thermodynamic stability of alpha‐helix conformation whose kinetic stability is drastically reduced in water solvent. Our analysis suggests that this is due to two effects enhanced by water: a more local effect consisting of the demolition of intra‐peptide H‐bonds, essential for the alpha‐helix formation, and a bulk – electrostatic – effect favoring conformational states more polar than alpha‐helix. Copyright © 2017 European Peptide Society and John Wiley & Sons, Ltd.     

Deciphering the crucial residues involved in heterodimerization of Bak peptide and anti-apoptotic proteins for apoptosis

B-cell lymphoma 2 (Bcl-2) family proteins are the central regulators of apoptosis, functioning via mitochondrial outer membrane permeabilization. The family members are involved in several stages of apoptosis regulation. The overexpression of the anti-apoptotic proteins leads to several cancer pathological conditions. This overexpression is modulated or inhibited by heterodimerization of pro-apoptotic BH3 domain or BH3-only peptides to the hydrophobic groove present at the surface of anti-apoptotic proteins. Additionally, the heterodimerization displayed differences in binding affinity profile among the pro-apoptotic peptides binding to anti-apoptotic proteins. In light of discovering the novel peptide/drug molecules that contain the potential to inhibit specific anti-apoptotic protein, it is necessary to understand the molecular basis of recognition between the protein and its binding partner (peptide or ligand) along with its binding energies. Therefore, the present work focused on deciphering the molecular basis of recognition between pro-apoptotic Bak peptide binding to different anti-apoptotic (Bcl-xL, Bfl-1, Bcl-W, Mcl-1, and Bcl-2) proteins using advanced Molecular Dynamics (MD) approach such as Molecular Mechanics-Generalized Born Solvent Accessible. The results from our investigation revealed that the predicted binding free energies showed excellent correlation with the experimental values (r² = .95). The electrostatic (ΔG_ele) contributions are the major component that drives the interaction between Bak peptides and different anti-apoptotic peptides. Additionally, van der Waals (ΔG_vdw) energies also play an indispensible role in determining the binding free energy. Furthermore, the decomposition analysis highlighted the comprehensive information about the energy contributions of hotspot residues involved in stabilizing the interaction between Bak peptide and different anti-apoptotic proteins.     

Molecular dynamics simulation of the phosphorylation‐induced conformational changes of a tau peptide fragment

Aggregation of the microtubule associated protein tau (MAPT) within neurons of the brain is the leading cause of tauopathies such as Alzheimer's disease. MAPT is a phospho‐protein that is selectively phosphorylated by a number of kinases in vivo to perform its biological function. However, it may become pathogenically hyperphosphorylated, causing aggregation into paired helical filaments and neurofibrillary tangles. The phosphorylation induced conformational change on a peptide of MAPT (htau_225?250) was investigated by performing molecular dynamics simulations with different phosphorylation patterns of the peptide (pThr231 and/or pSer235) in different simulation conditions to determine the effect of ionic strength and phosphate charge. All phosphorylation patterns were found to disrupt a nascent terminal β‐sheet pattern (²²⁶VAVVR²³⁰ and ²⁴⁴QTAPVP²⁴⁹), replacing it with a range of structures. The double pThr231/pSer235 phosphorylation pattern at experimental ionic strength resulted in the best agreement with NMR structural characterization, with the observation of a transient α‐helix (²³⁹AKSRLQT²⁴⁵). PPII helical conformations were only found sporadically throughout the simulations. Proteins 2014; 82:1907–1923. © 2014 Wiley Periodicals, Inc.     

Anisotropy and anharmonicity of atomic fluctuations in proteins: analysis of a molecular dynamics simulation

Positional probability density functions (pdf) for the atomic fluctuations are determined from a molecular dynamics simulation for hen egg-white lysozyme. Most atoms are found to have motions that are highly anisotropic but only slightly anharmonic. The largest deviations from harmonic motion are in the direction of the largest rms fluctuations in the local principal axis frame. Backbone atoms tend to be more nearly harmonic than sidechain atoms. The atoms with the largest anharmonicities tend to have pdfs with multiple peaks, each of which is close to harmonic. Several model pdfs are evaluated on the basis of how well they fit probability densities from the dynamics simulations when parameterized in terms of the moments of the distribution. Gram-Charlier and Edgeworth perturbation expansions, which have been successful in describing the motions of small molecules in crystals, are shown to be inadequate for the distributions found in the dynamics of proteins. Multipeaked distribution functions are found to be more appropriate.     

Glycogen phosphorylase inhibitors: a free energy perturbation analysis of glucopyranose spirohydantoin analogues

GP catalyzes the phosphorylation of glycogen to Glc-1-P. Because of its fundamental role in the metabolism of glycogen, GP has been the target for a systematic structure-assisted design of inhibitory compounds, which could be of value in the therapeutic treatment of type 2 diabetes mellitus. The most potent catalytic-site inhibitor of GP identified to date is spirohydantoin of glucopyranose (hydan). In this work, we employ MD free energy simulations to calculate the relative binding affinities for GP of hydan and two spirohydantoin analogues, methyl-hydan and n-hydan, in which a hydrogen atom is replaced by a methyl- or amino group, respectively. The results are compared with the experimental relative affinities of these ligands, estimated by kinetic measurements of the ligand inhibition constants. The calculated binding affinity for methyl-hydan (relative to hydan) is 3.75 +/- 1.4 kcal/mol, in excellent agreement with the experimental value (3.6 +/- 0.2 kcal/mol). For n-hydan, the calculated value is 1.0 +/- 1.1 kcal/mol, somewhat smaller than the experimental result (2.3 +/- 0.1 kcal/mol). A free energy decomposition analysis shows that hydan makes optimum interactions with protein residues and specific water molecules in the catalytic site. In the other two ligands, structural perturbations of the active site by the additional methyl- or amino group reduce the corresponding binding affinities. The computed binding free energies are sensitive to the preference of a specific water molecule for two well-defined positions in the catalytic site. The behavior of this water is analyzed in detail, and the free energy profile for the translocation of the water between the two positions is evaluated. The results provide insights into the role of water molecules in modulating ligand binding affinities. A comparison of the interactions between a set of ligands and their surrounding groups in X-ray structures is often used in the interpretation of binding free energy differences and in guiding the design of new ligands. For the systems in this work, such an approach fails to estimate the order of relative binding strengths, in contrast to the rigorous free energy treatment.     

Protein-inhibitor flexible docking by a multicanonical sampling: native complex structure with the lowest free energy and a free-energy barrier distinguishing the native complex from the others

Flexible docking between a protein (lysozyme) and an inhibitor (tri-N-acetyl-D-glucosamine, tri-NAG) was carried out by an enhanced conformational sampling method, multicanonical molecular dynamics simulation. We used a flexible all-atom model to express lysozyme, tri-NAG, and water molecules surrounding the two bio-molecules. The advantages of this sampling method are as follows: the conformation of system is widely sampled without trapping at energy minima, a thermally equilibrated conformational ensemble at an arbitrary temperature can be reconstructed from the simulation trajectory, and the thermodynamic weight can be assigned to each sampled conformation. During the simulation, exchanges between the binding and free (i.e., unbinding) states of the protein and the inhibitor were repeatedly observed. The conformational ensemble reconstructed at 300 K involved various conformational clusters. The main outcome of the current study is that the most populated conformational cluster (i.e., the cluster of the lowest free energy) was assigned to the native complex structure (i.e., the X-ray complex structure). The simulation also produced non-native complex structures, where the protein and the inhibitor bound with different modes from that of the native complex structure, as well as the unbinding structures. A free-energy barrier (i.e., activation free energy) was clearly detected between the native complex structures and the other structures. The thermal fluctuations of tri-NAG in the lowest free-energy complex correlated well with the X-ray B-factors of tri-NAG in the X-ray complex structure. The existence of the free-energy barrier ensures that the lowest free-energy structure can be discriminated naturally from the other structures. In other words, the multicanonical molecular dynamics simulation can predict the native complex structure without any empirical objective function. The current study also manifested that the flexible all-atom model and the physico-chemically defined atomic-level force field can reproduce the native complex structure. A drawback of the current method is that it requires a time consuming computation due to the exhaustive conformational sampling. We discussed a possibility for combining the current method with conventional docking methods.     

The binding mechanism of a novel nicotinamide isostere inhibiting with TNKSs: a molecular dynamic simulation and binding free energy calculation

Tankyrases (TNKSs), a member of human poly (ADP-ribose) polymerase (PARP) protein superfamily, plays a key role in regulation of cell proliferation. Among the representative proteins of the PARPs family, it is found that the inhibitors have high selectivity for Tankyrase1 (TNKS1). The specific binding modes are investigated between the TNKS1 protein and nicotinamide isostere (ISX) which functions as an inhibitor of TNKS1. The stabilities of ISX-TNKS1 and AVA939-TNKS1 complexes are estimated by molecular dynamics (MD) simulations and free energy calculations; a good agreement with experimental results is reached. On the basis of the calculated results of MD simulations, we found that the inhibitors influence the conformational flexibility of TNKS1 and the XAV939 binding drive the peptide Ile1228-Gly1229-Gly1230 to form a helical structure while the ISX binding drive the peptide to form a turn structure. Moreover, the formed important hydrogen bonds of Tyr1203 residue with XVA939 and WAT1551 with ISX enhance stabilities of the complexes, and the electrostatic interactions in XAV939-TNKS1 and van der Waals interactions in ISX-TNKS1 system are main driving forces for affinity. According to the results of the decomposition of binding free energy, it is obvious that the residues Try1224 and Lys1220 make the most favorable contributions to the binding in, respectively, ISX and XAV939 complexes. Taken together, the obtained results are useful for studying the binding mechanisms of TNKSs and inhibitors and for designing potent inhibitors.     

The severity of osteogenesis imperfecta: a comparison to the relative free energy differences of collagen model peptides

Molecular dynamics simulations were carried out to calculate free energy differences between the folded and unfolded states of wild type and mutant collagen model peptides. The calculated stability of the collagen models was compared with the severity of osteogenesis imperfecta. Free energy differences of Gly → Xaa (Xaa: Ser, Cys, Glu, and Asp) mutations between the wild type and the mutants at position 15 of the model peptide were 3.8, 4.2, 5.6, and 8.8 kcal/mol, respectively. The corresponding free energy differences of a second Gly mutation at the same position in different chains were, on average, 1.3, 1.5, 2.9, and 5.4 kcal/mol, respectively. Free energy simulations were also performed to estimate the relative stability between an oxidized form and a reduced form of the mutants containing two Cys residues, which indicated that the mutant of the collagen-like peptide containing an intramolecular disulfide bond was more stable than the mutant containing one Cys residue but less stable than the wild type. The calculated free energy differences between an oxidized and a reduced form of the mutants containing two Cys residues are 0.8 and 2.6 kcal/mol for the disulfide bonds between Chains A and B and between Chains A and C, respectively.