Free energy calculations and molecular dynamics simulations of wild-type and variants of the DNA-EcoRI complex.

Molecular dynamics simulations and free energy calculations of the wild-type EcoRI-DNA complex and several variants have been performed in aqueous solvent. In general, he theoretical estimations of the free energy differences (DeltaDeltaA) qualitatively agree well with the corresponding experimental data. The modifications which were experimentally found unfavorable compared to the wild-type complex were also found to be so in theoretical estimates. The mutant where the amino group of the base Ade(6) was replaced by a hydrogen atom eliminating one H-bond between the DNA and the protein, was experimentally found to be more stable than the wild-type complex. It was speculated that the modification also caused a structural relaxation in the DNA making DeltaDeltaA favorable. Our theoretical estimate yields a positive DeltaDeltaA in this case, but the difference is small, and no significant local structural relaxation was observed. The major H-bonds between the DNA and the protein in the wild-type complex are found to be maintained in the different mutants although the specific and non-specific interaction energies between the interacting the DNA bases and the protein residues are different in different mutants. The interaction pattern of the other nearby nucleotides are significantly influenced by each modification. Thus, the alteration of the non-specific interactions may also play an indirect role in determining the specificity of the complex. The interaction of the Gua(4) of the DNA with the protein is found to be most sensitive to any alteration in the recognition site. Because Gua(4) is the nucleotide closest to the scissile bond, this extra sensitivity seems to play an important role in altering the functional activity of the complex.     

Free energy surfaces of beta-hairpin and alpha-helical peptides generated by replica exchange molecular dynamics with the AGBNP implicit solvent model

We have studied the potential of mean force of two peptides, one known to adopt a beta-hairpin and the other an alpha-helical conformation in solution. These peptides are, respectively, residues 41-56 of the C-terminus (GEWTYDDATKTFTVTE) of the B1 domain of protein G and the 13 residue C-peptide (KETAAAKFERQHM) of ribonuclease A. Extensive canonical ensemble sampling has been performed using a parallel replica exchange method. The effective potential employed in this work consists of the OPLS all-atom force field (OPLS-AA) and an analytical generalized Born (AGB) implicit solvent model including a novel nonpolar solvation free energy estimator (NP). An additional dielectric screening parameter has been incorporated into the AGBNP model. In the case of the beta-hairpin, the nonpolar solvation free energy estimator provides the necessary effective interactions for the collapse of the hydrophobic core (W43, Y45, F52, and V54), which the more commonly used surface-area-dependent nonpolar model does not provide. For both the beta-hairpin and the alpha-helix, increased dielectric screening reduces the stability of incorrectly formed salt bridges, which tend to disrupt the formation of the hairpin and helix, respectively. The fraction of beta-hairpin and alpha-helix content we obtained using the AGBNP model agrees well with experimental results. The thermodynamic stability of the beta-hairpin from protein G and the alpha-helical C-peptide from ribonuclease A as modeled with the OPLS-AA/AGBNP effective potential reflects the balance between the nonpolar effective potential terms, which drive compaction, and the polar and hydrogen bonding terms, which promote secondary structure formation.     

Efficiency of a second-generation HIV-1 protease inhibitor studied by molecular dynamics and absolute binding free energy calculations

A subnanomolar inhibitor of human immunodeficiency virus type 1 (HIV-1) protease, designated QF34, potently inhibits the wild-type and drug-resistant enzyme. To explain its broad activity, the binding of QF34 to the wild-type HIV-1 protease is investigated by molecular dynamics simulations and compared to the binding of two inhibitors that are used clinically, saquinavir (SQV) and indinavir (IDV). Analysis of the flexibility of protease residues and inhibitor segments in the complex reveals that segments of QF34 were more mobile during the dynamics studies than the segments of SQV and IDV. The dynamics of hydrogen bonding show that QF34 forms a larger number of stable hydrogen bonds than the two inhibitors that are used clinically. Absolute binding free energies were calculated with molecular mechanics-generalized Born surface area (MM-GBSA) methodology using three protocols. The most consistent results were obtained using the single-trajectory approach, due to cancellation of errors and inadequate sampling in the separate-trajectory protocols. For all three inhibitors, energy components in favor of binding include van der Waals and electrostatic terms, whereas polar solvation and entropy terms oppose binding. Decomposition of binding energies reveals that more protease residues contribute significantly to the binding of QF34 than to the binding of SQV and IDV. Moreover, contributions from protease main chains and side chains are balanced in the case of QF34 (52:48 ratio, respectively), whereas side chain contributions prevail in both SQV and IDV (main-chain:side-chain ratios of 41:59 and 45:55, respectively). The presented results help explain the ability of QF34 to inhibit multiple resistant mutants and should be considered in the design of broad-specificity second-generation HIV-1 protease inhibitors.     

Proton transfer dynamics in the propionic acid dimer from path integral molecular dynamics calculations

The double proton transfer process in the cyclic dimer of propionic acid in the gas phase was studied using a path integral molecular dynamics method. Structures, energies and proton trajectories were determined. Very large amplitude motions of the skeleton of a propionic acid molecule were observed during the simulations, and almost free rotation of the C₂H₅ group around the C_α-C bond. A double-well symmetric potential with a very small energy barrier was determined from the free energy profile for the proton motions. Infrared spectra for different isotopomers were calculated, and comparative vibrational analysis was performed. The vibrational results from CPMD appear to be in qualitative agreement with the experimental ones.     

Structure and stability of the insulin dimer investigated by molecular dynamics simulation

Molecular dynamics simulation indicates that the dynamical behaviour of the insulin dimer is asymmetric. Atomic level knowledge of the interaction modes and protein conformation in the solvation state identifies dynamical structures, held by hydrogen bonds that stabilize, mainly in one monomer, the interaction between the chains. Dynamic cross-correlation analysis shows that the two insulin monomers behave asymmetrically and are almost independent. Solvation energy, calculated to evaluate the contribute of each interface residue to the dimer association pattern, well compares with the experimental association state found in protein mutants indicating that this parameter is an important factor to explain the association properties of mutated insulin dimers.     

Absolute binding free energy calculations using molecular dynamics simulations with restraining potentials

The absolute (standard) binding free energy of eight FK506-related ligands to FKBP12 is calculated using free energy perturbation molecular dynamics (FEP/MD) simulations with explicit solvent. A number of features are implemented to improve the accuracy and enhance the convergence of the calculations. First, the absolute binding free energy is decomposed into sequential steps during which the ligand-surrounding interactions as well as various biasing potentials restraining the translation, orientation, and conformation of the ligand are turned "on" and "off." Second, sampling of the ligand conformation is enforced by a restraining potential based on the root mean-square deviation relative to the bound state conformation. The effect of all the restraining potentials is rigorously unbiased, and it is shown explicitly that the final results are independent of all artificial restraints. Third, the repulsive and dispersive free energy contribution arising from the Lennard-Jones interactions of the ligand with its surrounding (protein and solvent) is calculated using the Weeks-Chandler-Andersen separation. This separation also improves convergence of the FEP/MD calculations. Fourth, to decrease the computational cost, only a small number of atoms in the vicinity of the binding site are simulated explicitly, while all the influence of the remaining atoms is incorporated implicitly using the generalized solvent boundary potential (GSBP) method. With GSBP, the size of the simulated FKBP12/ligand systems is significantly reduced, from approximately 25,000 to 2500. The computations are very efficient and the statistical error is small ( approximately 1 kcal/mol). The calculated binding free energies are generally in good agreement with available experimental data and previous calculations (within approximately 2 kcal/mol). The present results indicate that a strategy based on FEP/MD simulations of a reduced GSBP atomic model sampled with conformational, translational, and orientational restraining potentials can be computationally inexpensive and accurate.     

Structural stability of wild type and mutated alpha-keratin fragments: molecular dynamics and free energy calculations

The objective of this study is to investigate the influence of point mutations on the structural stability of coiled coil fragments of the human hair intermediate filament by molecular dynamics simulations and free energy calculations. Mutations in the helix termination motif of human hair keratin gene hHb6 seem to be connected to the hereditary hair dystrophy Monilethrix. The most common mutations reported are Glu413Lys and Glu413Asp, located at the C-terminal end of the coiled coil 2B rod domain of the IF. According to our simulations, significant conformational changes of the side chains at the mutation and neighboring sites occur due to the Glu413Lys mutation. Furthermore, the differences in electrostatic interactions cause a large change in free energy during transformation of Glu413 to Lys calculated by the thermodynamic integration approach. It is speculated that the structural rearrangement necessary to adapt the interactions in the mutated coiled coil leads to changes in the IF assembly or its stability. The second mutation, Glu413Asp, only leads to a small value of the calculated free energy difference that is within the error limits of the simulations. Thus, it has to be concluded that this mutation does not affect the coiled coil stability.     

Free energy surface for rotamers of cis-enol malonaldehyde in aqueous solution studied by molecular dynamics calculations

Two-dimensional free energy surfaces for four rotamers of cis-enol malonaldehyde in water have been investigated by umbrella sampling molecular dynamics (MD) calculations. Biasing potential used in the umbrella sampling calculation was adopted to be the minus of conformational free energy preliminary obtained by the thermodynamic integration MD calculations for the rigid malonaldehyde whose stretching and bending were all fixed. The calculated free energy surface shows that, in water, a rotamer that has an intramolecular hydrogen bond is most stable among the rotamers. This is the same as that in vacuum, while order of relative stability of the other three rotamers is different in water and in vacuum. Inclusion of intramolecular vibrations changed the free energy surface little, i.e. at most 2.6 kJ/mol, which is much smaller than the solvation free energy. Free energy barriers from the most stable intramolecular hydrogen bonded rotamer to the others are lowered by hydration but they are still very high, >50 kJ/mol, such that the malonaldehyde molecule spends most of its time in water taking this conformation. Thus, reaction coordinate for intramolecular proton transfer reaction in water may be constructed assuming this rotamer.     

Computational investigation of the HIV-1 Rev multimerization using molecular dynamics simulations and binding free energy calculations

The HIV Rev protein mediates the nuclear export of viral mRNA, and is thereby essential for the production of late viral proteins in the replication cycle. Rev forms a large organized multimeric protein-protein complex for proper functioning. Recently, the three-dimensional structures of a Rev dimer and tetramer have been resolved and provide the basis for a thorough structural analysis of the binding interaction. Here, molecular dynamics (MD) and binding free energy calculations were performed to elucidate the forces thriving dimerization and higher order multimerization of the Rev protein. It is found that despite the structural differences between each crystal structure, both display a similar behavior according to our calculations. Our analysis based on a molecular mechanics-generalized Born surface area (MM/GBSA) and a configurational entropy approach demonstrates that the higher order multimerization site is much weaker than the dimerization site. In addition, a quantitative hot spot analysis combined with a mutational analysis reveals the most contributing amino acid residues for protein interactions in agreement with experimental results. Additional residues were found in each interface, which are important for the protein interaction. The investigation of the thermodynamics of the Rev multimerization interactions performed here could be a further step in the development of novel antiretrovirals using structure based drug design. Moreover, the variability of the angle between each Rev monomer as measured during the MD simulations suggests a role of the Rev protein in allowing flexibility of the arginine rich domain (ARM) to accommodate RNA binding.     

Ion transport through membrane-spanning nanopores studied by molecular dynamics simulations and continuum electrostatics calculations

Narrow hydrophobic regions are a common feature of biological channels, with possible roles in ion-channel gating. We study the principles that govern ion transport through narrow hydrophobic membrane pores by molecular dynamics simulation of model membranes formed of hexagonally packed carbon nanotubes. We focus on the factors that determine the energetics of ion translocation through such nonpolar nanopores and compare the resulting free-energy barriers for pores with different diameters corresponding to the gating regions in closed and open forms of potassium channels. Our model system also allows us to compare the results from molecular dynamics simulations directly to continuum electrostatics calculations. Both simulations and continuum calculations show that subnanometer wide pores pose a huge free-energy barrier for ions, but a small increase in the pore diameter to approximately 1 nm nearly eliminates that barrier. We also find that in those wider channels the ion mobility is comparable to that in the bulk phase. By calculating local electrostatic potentials, we show that the long range Coulomb interactions of ions are strongly screened in the wide water-filled channels. Whereas continuum calculations capture the overall energetics reasonably well, the local water structure, which is not accounted for in this model, leads to interesting effects such as the preference of hydrated ions to move along the pore wall rather than through the center of the pore.     

Analysis of the stability of looped-out and stacked-in conformations of an adenine bulge in DNA using a continuum model for solvent and ions

A combination of conformational search, energy minimization, and energetic evaluation using a continuum solvent treatment has been employed to study the stability of various conformations of the DNA fragment d(CGCAGAA)/d(TTCGCG) containing a single adenine bulge. The extra-helical (looped-out) bulge conformation derived from a published x-ray structure and intra-helical (stacked bulge base) model structures partially based on nuclear magnetic resonance (NMR) data were used as start structures for the conformational search. Solvent-dependent contributions to the stability of the conformations were calculated from the solvent exposed molecular surface area and by using the finite difference Poisson-Boltzmann approach. Three classes (I-III) of bulge conformations with calculated low energies can be distinguished. The lowest-energy conformations were found in class I, corresponding to structures with the bulge base stacked between flanking helices, and class II, composed of structures forming a triplet of the bulge base and a flanking base pair. All extra-helical bulge structures, forming class III, were found to be less stable compared with the lowest energy structures of class I and II. The results are consistent with NMR data on an adenine bulge in the same sequence context indicating an intra-helical or triplet bulge conformation in solution. Although the total energies and total electrostatic energies of the low-energy conformations show only relatively modest variations, the energetic contributions to the stability were found to vary significantly among the classes of bulge structures. All intra-helical bulge structures are stabilized by a more favorable Coulomb charge-charge interaction but destabilized by a larger electrostatic reaction field contribution compared with all extra-helical and most triplet bulge structures. Van der Waals packing interactions and nonpolar surface-area-dependent contributions appear to favor triplet class II structures and to a lesser degree also the intra-helical stacked bulge conformations. The large conformational variation found for class III conformers might add a favorable entropic contribution to the stability of the extra-helical bulge form.     

Insight into binding mechanisms of EGFR allosteric inhibitors using molecular dynamics simulations and free energy calculations

Abstract

Lung cancer is the leading cause of cancer death, and epidermal growth factor receptor (EGFR) kinase domain mutations are a common cause of non-small-cell lung cancer (NSCLC), a major subtype of lung cancers. Patients harboring most of these mutations respond well to the EGFR inhibitors Gefitinib and Erlotinib initially, but soon develop resistance to them due to the emergence of the gatekeeper mutation T790M. The new-generation inhibitors such as AZD9291, HM61713, CO-1686 and WZ4002 can overcome T790M through covalent binding to Cys 797, but ultimately lose their efficacy upon the emergence of the C797S mutation that abolishes the covalent bonding. Allosteric inhibitors EAI001 and EAI045 are a new type of EGFR inhibitors that bind to EGFR away from the ATP-binding site and not relying on Cys 797. In this study, molecular dynamics simulations and free energy calculations were carried out on EAI001 and EAI045 in complex with EGFR, revealing the detailed inhibitory mechanism of EAI001 and EAI045 as EGFR allosteric inhibitor, which was expected to provide a basis for rational drug design of the EGFR allosteric inhibitors.

Communicated by Ramaswamy H. Sarma     

Ion solvation and structural stability in a sodium channel investigated by molecular dynamics calculations

The stability and ion binding properties of the homo-tetrameric pore domain of a prokaryotic, voltage-gated sodium channel are studied by extensive all-atom molecular dynamics simulations, with the channel protein being embedded in a fully hydrated lipid bilayer. It is found that Na(+) ion presents in a mostly hydrated state inside the wide pore of the selectivity filter of the sodium channel, in sharp contrast to the nearly fully dehydrated state for K(+) ions in potassium channels. Our results also indicate that Na(+) ions make contact with only one or two out of the four polypeptide chains forming the selectivity filter, and surprisingly, the selectivity filter exhibits robust stability for various initial ion configurations even in the absence of ions. These findings are quite different from those in potassium channels. Furthermore, an electric field above 0.5V/nm is suggested to be able to induce Na(+) permeation through the selectivity filter.     

Evaluation of a fast implicit solvent model for molecular dynamics simulations.

A solvation term based on the solvent accessible surface area (SASA) is combined with the CHARMM polar hydrogen force field for the efficient simulation of peptides and small proteins in aqueous solution. Only two atomic solvation parameters are used: one is negative for favoring the direct solvation of polar groups and the other positive for taking into account the hydrophobic effect on apolar groups. To approximate the water screening effects on the intrasolute electrostatic interactions, a distance-dependent dielectric function is used and ionic side chains are neutralized. The use of an analytical approximation of the SASA renders the model extremely efficient (i.e., only about 50% slower than in vacuo simulations). The limitations and range of applicability of the SASA model are assessed by simulations of proteins and structured peptides. For the latter, the present study and results reported elsewhere show that with the SASA model it is possible to sample a significant amount of folding/unfolding transitions, which permit the study of the thermodynamics and kinetics of folding at an atomic level of detail.     

Free energy calculations of ALS‐causing SOD1 mutants reveal common perturbations to stability and dynamics along the maturation pathway

With over 150 heritable mutations identified as disease‐causative, superoxide dismutase 1 (SOD1) has been a main target of amyotrophic lateral sclerosis (ALS) research and therapeutic efforts. However, recent evidence has suggested that neither loss of function nor protein aggregation is responsible for promoting neurotoxicity. Furthermore, there is no clear pattern to the nature or the location of these mutations that could suggest a molecular mechanism behind SOD1‐linked ALS. Here, we utilize reliable and accurate computational techniques to predict the perturbations of 10 such mutations to the free energy changes of SOD1 as it matures from apo monomer to metallated dimer. We find that the free energy perturbations caused by these mutations strongly depend on maturational progress, indicating the need for state‐specific therapeutic targeting. We also find that many mutations exhibit similar patterns of perturbation to native and non‐native maturation, indicating strong thermodynamic coupling between the dynamics at various sites of maturation within SOD1. These results suggest the presence of an allosteric network in SOD1 which is vulnerable to disruption by these mutations. Analysis of these perturbations may contribute to uncovering a unifying molecular mechanism which explains SOD1‐linked ALS and help to guide future therapeutic efforts.     

Protein molecular dynamics with the generalized Born/ACE solvent model

Implicit solvent models are increasingly important for the study of proteins in aqueous solution. Here, the generalized Born (GB) solvent polarization model as implemented in the analytical ACE potential [Schaefer and Karplus (1996) J Phys Chem 100:1578] is used to perform molecular dynamics simulations of two small, homologous proteins: the immunoglobulin-binding domain of streptococcal protein G and the Ras binding domain of Raf. Several model parameterizations are compared through more than 60 ns of simulation. Results are compared with two simpler solvent models-an accessible surface area model and a distant-dependent dielectric model, with finite-difference Poisson calculations, with existing explicit solvent simulations, and with experimental data. The simpler models yield stable but distorted structures. The best GB/ACE implementation uses a set of atomic Voronoi volumes reported recently, obtained by averaging over a large database of crystallographic protein structures. A 20% reduction is applied to the volumes, compensating in an average sense for an excessive de-screening of individual charges inherent in the ACE self-energy and for an undersolvation of dipolar groups inherent in the GB screening function. This GB/ACE parameterization yields stable trajectories on the 0.5-1-ns time scale that deviate moderately (approximately 1.5-2.5 A) from the X-ray structure, reproduce approximately the surface distribution of charged, polar, and hydrophobic groups, and reproduce accurately backbone flexibility as measured by amide NMR-order parameters. Over longer time scales (1.5-3 ns), some of the protein G runs escape from the native energy basin and deviate strongly (3 A) from the native structure. The conformations sampled during the transition out of the native energy basin are overstabilized by the GB/ACE solvation model, as compared with a numerical treatment of the full dielectric continuum model.