Calculation of binding affinities of HIV-1 RT and β-secretase inhibitors using the linear interaction energy method with explicit and continuum solvation approaches

The linear interaction energy (LIE) approach has been applied to estimate the binding free energies of representative sets of HIV-1 RT and β-Secretase inhibitors, using both molecular dynamics (MD) and tethered energy minimization sampling protocols with the OPLS-AA potential, using a range of solvation methodologies. Generalized Born (GB), ‘shell’ and periodic boundary condition (PBC) solvation were used, the latter with reaction field (RF) electrostatics. Poisson-Boltzmann (PB) and GB continuum electrostatics schemes were applied to the simulation trajectories for each solvation type to estimate the electrostatic ligand-water interaction energy in both the free and bound states. Reasonable agreement of the LIE predictions was obtained with respect to experimental binding free energy estimates for both systems: for instance, ‘PB’ fits on MD trajectories carried out with PBC solvation and RF electrostatics led to models with standard errors of 1.11 and 1.03 kcal mol⁻¹ and coefficients of determination, r ² of 0.76 and 0.75 for the HIV-1 RT and β-Secretase sets. However, it was also found that results from MD sampling using PBC solvation provided only slightly better fits than from simulations using shell or Born solvation or tethered energy minimization sampling. Figure Evolution of the running averages for compound H11 (binding to HIV-1RT) of the bound state ligand-water and ligand-protein interaction energies. The ligand-water electrostatic terms are twice the corresponding GB and PB electrostatic solvation free energies. The ligand-receptor van der Waals and Coulombic interaction energies are also shown, in addition to the ligand-water van der Waals interaction term. The terms were calculated (without application of a cut-off) from a trajectory sampled under PBC solvation with reaction field electrostatics Electronic supplementary material The online version of this article (doi:) contains supplementary material, which is available to authorized users.     

High pressure reveals structural determinants for globin hexacoordination: Neuroglobin and myoglobin cases

The influence of pressure on the equilibrium between five‐(5c) and six‐coordination (6c) forms in neuroglobin (Ngb) and myoglobin (Mb) has been examined by means of molecular dynamics (MD) simulations at normal and high pressure. The results show that the main effect of high pressure is to reduce the protein mobility without altering the structure in a significant manner. Moreover, our data suggest that the equilibrium between 5c and 6c states in globins is largely controlled by the structure and dynamics of the C‐D region. Finally, in agreement with the available experimental data, the free energy profiles obtained from steered MD for both proteins indicate that high pressure enhances hexacoordination. In Ngb, the shift in equilibrium is mainly related to an increase in the 6c→5c transition barrier, whereas in Mb such a shift is primarily due to a destabilization of the 5c state. Proteins 2009. © 2008 Wiley‐Liss, Inc.     

Effects of ligand binding on the dynamics of rice nonspecific lipid transfer protein 1: a model from molecular simulations

Plant nonspecific lipid transfer proteins (nsLTPs) are small, basic proteins constituted mainly of alpha-helices and stabilized by four conserved disulfide bridges. They are characterized by the presence of a tunnel-like hydrophobic cavity, capable of transferring various lipid molecules between lipid bilayers in vitro. In this study, molecular dynamics (MD) simulations were performed at room temperature to investigate the effects of lipid binding on the dynamic properties of rice nsLTP1. Rice nsLTP1, either in the free form or complexed with one or two lipids was subjected to MD simulations. The C-terminal loop was very flexible both before and after lipid binding, as revealed by calculating the root-mean-square fluctuation. After lipid binding, the flexibility of some residues that were not in direct contact with lipid molecules increased significantly, indicating an increase of entropy in the region distal from the binding site. Essential dynamics analysis revealed clear differences in motion between unliganded and liganded rice nsLTP1s. In the free form of rice nsLTP1, loop1 exhibited the largest directional motion. This specific essential motion mode diminished after binding one or two lipid molecules. To verify the origin of the essential motion observed in the free form of rice nsLTP1, we performed multiple sequence alignments to probe the intrinsic motion encoded in the primary sequence. We found that the amino acid sequence of loop1 is highly conserved among plant nsLTP1s, thus revealing its functional importance during evolution. Furthermore, the sequence of loop1 is composed mainly of amino acids with short side chains. In this study, we show that MD simulations, together with essential dynamics analysis, can be used to determine structural and dynamic differences of rice nsLTP1 upon lipid binding.     

Molecular dynamics guided study of salt bridge length dependence in both fluorinated and non-fluorinated parallel dimeric coiled-coils

The alpha-helical coiled-coil is one of the most common oligomerization motifs found in both native and engineered proteins. To better understand the stability and dynamics of the coiled-coil motifs, including those modified by fluorination, several fluorinated and nonfluorinated parallel dimeric coiled-coil protein structures were designed and modeled. We also attempt to investigate how changing the length and geometry of the important stabilizing salt bridges influences the coiled-coil protein structure. Molecular dynamics (MD) and free energy simulations with AMBER used a particle mesh Ewald treatment of the electrostatics in explicit TIP3P solvent with balanced force field treatments. Preliminary studies with legacy force fields (ff94, ff96, and ff99) show a profound instability of the coiled-coil structures in short MD simulation. Significantly, better behavior is evident with the more balanced ff99SB and ff03 protein force fields. Overall, the results suggest that the coiled-coil structures can readily accommodate the larger acidic arginine or S-2,7-diaminoheptanedoic acid mutants in the salt bridge, whereas substitution of the smaller L-ornithine residue leads to rapid disruption of the coiled-coil structure on the MD simulation time scale. This structural distortion of the secondary structure allows both the formation of large hydration pockets proximal to the charged groups and within the hydrophobic core. Moreover, the increased structural fluctuations and movement lead to a decrease in the water occupancy lifetimes in the hydration pockets. In contrast, analysis of the hydration in the stable dimeric coiled-coils shows high occupancy water sites along the backbone residues with no water occupancy in the hydrophobic core, although transitory water interactions with the salt bridge residues are evident. The simulations of the fluorinated coiled-coils suggest that in some cases fluorination electrostatically stabilizes the intermolecular coiled-coil salt bridges. Structural analyses also reveal different side chain rotamer preferences for leucine when compared with 5,5,5,5',5',5'-hexafluoroleucine mutants. These observed differences in the side chain rotamer populations suggest differential changes in the side chain conformational entropy upon coiled-coil formation when the protein is fluorinated. The free energy of hydration of the isolated 5,5,5,5',5',5'-hexafluoroleucine amino acid is calculated to be 1.1 kcal/mol less stable than leucine; this hydrophobic penalty in the monomer may provide a driving force for coiled-coil dimer formation. Estimation of the ellipticity at 222 nm from a series of snapshots from the MD simulations with DicroCalc shows distinct increases in the ellipticity when the coiled-coil is fluorinated, which suggests that the helicity in the folded coiled-coils is greater when fluorinated.     

Effective energy function for proteins in solution

A Gaussian solvent-exclusion model for the solvation free energy is developed. It is based on theoretical considerations and parametrized with experimental data. When combined with the CHARMM 19 polar hydrogen energy function, it provides an effective energy function (EEF1) for proteins in solution. The solvation model assumes that the solvation free energy of a protein molecule is a sum of group contributions, which are determined from values for small model compounds. For charged groups, the self-energy contribution is accounted for primarily by the exclusion model. Ionic side-chains are neutralized, and a distance-dependent dielectric constant is used to approximate the charge-charge interactions in solution. The resulting EEF1 is subjected to a number of tests. Molecular dynamics simulations at room temperature of several proteins in their native conformation are performed, and stable trajectories are obtained. The deviations from the experimental structures are similar to those observed in explicit water simulations. The calculated enthalpy of unfolding of a polyalanine helix is found to be in good agreement with experimental data. Results reported elsewhere show that EEF1 clearly distinguishes correctly from incorrectly folded proteins, both in static energy evaluations and in molecular dynamics simulations and that unfolding pathways obtained by high-temperature molecular dynamics simulations agree with those obtained by explicit water simulations. Thus, this energy function appears to provide a realistic first approximation to the effective energy hypersurface of proteins.     

Validating Solution Ensembles from Molecular Dynamics Simulation by Wide-Angle X-ray Scattering Data

Wide-angle x-ray scattering (WAXS) experiments of biomolecules in solution have become increasingly popular because of technical advances in light sources and detectors. However, the structural interpretation of WAXS profiles is problematic, partly because accurate calculations of WAXS profiles from structural models have remained challenging. In this work, we present the calculation of WAXS profiles from explicit-solvent molecular dynamics (MD) simulations of five different proteins. Using only a single fitting parameter that accounts for experimental uncertainties because of the buffer subtraction and dark currents, we find excellent agreement to experimental profiles both at small and wide angles. Because explicit solvation eliminates free parameters associated with the solvation layer or the excluded solvent, which would require fitting to experimental data, we minimize the risk of overfitting. We further find that the influence from water models and protein force fields on calculated profiles are insignificant up to q ≈ 15 nm^?1. Using a series of simulations that allow increasing flexibility of the proteins, we show that incorporating thermal fluctuations into the calculations significantly improves agreement with experimental data, demonstrating the importance of protein dynamics in the interpretation of WAXS profiles. In addition, free MD simulations up to one microsecond suggest that the calculated profiles are highly sensitive with respect to minor conformational rearrangements of proteins, such as an increased flexibility of a loop or an increase of the radius of gyration by  <  1%. The present study suggests that quantitative comparison between MD simulations and experimental WAXS profiles emerges as an accurate tool to validate solution ensembles of biomolecules.     

Identification of the critical structural determinants of the EF-hand domain arrangements in calcium binding proteins

EF-hand calcium binding proteins (CaBPs) share strong sequence homology, but exhibit great diversity in structure and function. Thus although calmodulin (CaM) and calcineurin B (CNB) both consist of four EF hands, their domain arrangements are quite distinct. CaM and the CaM-like proteins are characterized by an extended architecture, whereas CNB and the CNB-like proteins have a more compact form. In this study, we performed structural alignments and molecular dynamics (MD) simulations on 3 CaM-like proteins and 6 CNB-like proteins, and quantified their distinct structural and dynamical features in an effort to establish how their sequences specify their structures and dynamics. Alignments of the EF2-EF3 region of these proteins revealed that several residues (not restricted to the linker between the EF2 and EF3 motifs) differed between the two groups of proteins. A customized inverse folding approach followed by structural assessments and MD simulations established the critical role of these residues in determining the structure of the proteins. Identification of the critical determinants of the two different EF-hand domain arrangements and the distinct dynamical features relevant to their respective functions provides insight into the relationships between sequence, structure, dynamics and function among these EF-hand CaBPs.     

Subtle functional collective motions in pancreatic-like ribonucleases: from ribonuclease A to angiogenin

The analysis of the dynamic behavior of enzymes is fundamental to structural biology. A direct relationship between protein flexibility and biological function has been shown for bovine pancreatic ribonuclease (RNase A) (Rasmussen et al., Nature 1992;357:423-424). More recently, crystallographic studies have shown that functional motions in RNase A involve the enzyme beta-sheet regions that move concertedly on substrate binding and release (Vitagliano et al., Proteins 2002;46:97-104). These motions have been shown to correspond to intrinsic dynamic properties of the native enzyme by molecular dynamics (MD) simulations. To unveil the occurrence of these collective motions in other members of pancreatic-like superfamily, we carried out MD simulations on human angiogenin (Ang). Essential dynamics (ED) analyses performed on the trajectories reveal that Ang exhibits collective motions similar to RNase A, despite the limited sequence identity (33%) of the two proteins. Furthermore, we show that these collective motions are also present in ensembles of experimentally determined structures of both Ang and RNase A. Finally, these subtle concerted beta-sheet motions were also observed for other two members of the pancreatic-like superfamily by comparing the ligand-bound and ligand-free structures of these enzymes. Taken together, these findings suggest that pancreatic-like ribonucleases share an evolutionary conserved dynamic behavior consisting of subtle beta-sheet motions, which are essential for substrate binding and release.     

Virtual substitution scan via single‐step free energy perturbation

With the rapid expansion of our computing power, molecular dynamics (MD) simulations ranging from hundreds of nanoseconds to microseconds or even milliseconds have become increasingly common. The majority of these long trajectories are obtained from plain (vanilla) MD simulations, where no enhanced sampling or free energy calculation method is employed. To promote the “recycling” of these trajectories, we developed the Virtual Substitution Scan (VSS) toolkit as a plugin of the open‐source visualization and analysis software VMD. Based on the single‐step free energy perturbation (sFEP) method, VSS enables the user to post‐process a vanilla MD trajectory for a fast free energy scan of substituting aryl hydrogens by small functional groups. Dihedrals of the functional groups are sampled explicitly in VSS, which improves the performance of the calculation and is found particularly important for certain groups. As a proof‐of‐concept demonstration, we employ VSS to compute the solvation free energy change upon substituting the hydrogen of a benzene molecule by 12 small functional groups frequently considered in lead optimization. Additionally, VSS is used to compute the relative binding free energy of four selected ligands of the T4 lysozyme. Overall, the computational cost of VSS is only a fraction of the corresponding multi‐step FEP (mFEP) calculation, while its results agree reasonably well with those of mFEP, indicating that VSS offers a promising tool for rapid free energy scan of small functional group substitutions. © 2016 Wiley Periodicals, Inc. Biopolymers 105: 324–336, 2016.     

Modeling loop reorganization free energies of acetylcholinesterase: a comparison of explicit and implicit solvent models

The treatment of hydration effects in protein dynamics simulations varies in model complexity and spans the range from the computationally intensive microscopic evaluation to simple dielectric screening of charge-charge interactions. This paper compares different solvent models applied to the problem of estimating the free-energy difference between two loop conformations in acetylcholinesterase. Molecular dynamics (MD) simulations were used to sample potential energy surfaces of the two basins with solvent treated by means of explicit and implicit methods. Implicit solvent methods studied include the generalized Born (GB) model, atomic solvation potential (ASP), and the distance-dependent dieletric constant. By using the linear response approximation (LRA), the explicit solvent calculations determined a free-energy difference that is in excellent agreement with the experimental estimate, while rescoring the protein conformations with GB or the Poisson equation showed inconsistent and inferior results. While the approach of rescoring conformations from explicit water simulations with implicit solvent models is popular among many applications, it perturbs the energy landscape by changing the solvent contribution to microstates without conformational relaxation, thus leading to non-optimal solvation free energies. Calculations applying MD with a GB solvent model produced results of comparable accuracy as observed with LRA, yet the electrostatic free-energy terms were significantly different due to optimization on a potential energy surface favored by an implicit solvent reaction field. The simpler methods of ASP and the distance-dependent scaling of the dielectric constant both produced considerable distortions in the protein internal free-energy terms and are consequently unreliable.     

Effect of the reaction field electrostatic term on the molecular dynamics simulation of the activation domain of procarboxypeptidase B

Molecular dynamics simulations of the activation domain of porcine procarboxypeptidase B (ADBp) were performed in order to examine the effects of the inclusion of a reaction field (RF) term into the calculation of electrostatics forces for highly charged proteins. Two simulations were performed with the GROMOS96 package, studying the influence of counterions on the final results. Comparison with previous results without the inclusion of the RF term (Martí-Renom, M.A., Mas,J.M., Oliva,B., Querol,E. and Avilés,F.X., Protein Engng, 1998, 11, 101-110) shows that the structure is well maintained when the RF term is included. Moreover, the analysis of the trajectories shows that simulations of solvated highly-charged proteins are sensitive to the presence of counterions, the secondary structures being more stable when their charges are neutralized.     

Binding interactions between the core central domain of 16S rRNA and the ribosomal protein S15 determined by molecular dynamics simulations

The goal of the current study is to utilize molecular dynamic (MD) simulations to investigate the dynamic behavior of 16S rRNA in the presence and absence of S15 and to identify the binding interactions between these two molecules. The simulations show that: (i) 16S rRNA remains in a highly folded structure when it is bound to S15; (ii) in the absence of S15, 16S rRNA significantly alters its conformation and transiently forms conformations that are similar to the bound structure that make it available for binding with S15; (iii) the unbound rRNA spends the majority of its time in extended conformations. The formation of the extended conformations is a result of the molecule reaching a lower electrostatic energy and the formation of the highly folded, crystal-like conformation is a result of achieving a lower solvation energy. In addition, our MD simulations show that 16S rRNA and S15 bind across the major groove of helix 22 (H22) via electrostatic interactions. The negatively charged phosphate groups of G658, U740, G741 and G742 bind to the positively charged S15 residues Lys7, Arg34 and Arg37. The current study provides a dynamic view of the binding of 16S rRNA with S15.     

Long time dynamics of Met-enkephalin: comparison of explicit and implicit solvent models

Met-enkephalin is one of the smallest opiate peptides. Yet, its dynamical structure and receptor docking mechanism are still not well understood. The conformational dynamics of this neuron peptide in liquid water are studied here by using all-atom molecular dynamics (MD) and implicit water Langevin dynamics (LD) simulations with AMBER potential functions and the three-site transferable intermolecular potential (TIP3P) model for water. To achieve the same simulation length in physical time, the full MD simulations require 200 times as much CPU time as the implicit water LD simulations. The solvent hydrophobicity and dielectric behavior are treated in the implicit solvent LD simulations by using a macroscopic solvation potential, a single dielectric constant, and atomic friction coefficients computed using the accessible surface area method with the TIP3P model water viscosity as determined here from MD simulations for pure TIP3P water. Both the local and the global dynamics obtained from the implicit solvent LD simulations agree very well with those from the explicit solvent MD simulations. The simulations provide insights into the conformational restrictions that are associated with the bioactivity of the opiate peptide dermorphin for the delta-receptor.     

Molecular dynamics simulations of a highly charged peptide from an SH3 domain: possible sequence-function relationship

A seven-residue peptide that is highly conserved in SH3 domains despite being far from the active site has been shown by NMR to be stable in solution. This peptide, biologically important because it is a likely folding nucleus for SH3 domains, provides a challenging subject for molecular dynamics because it is highly charged. We present stable, 10-ns simulations of both the native-like diverging turn structure and a helical model. Free energies of these two conformations, estimated through MM-PBSA analysis using several force fields, suggest a comparable free energy (DeltaDeltaG < or =6 kcal/mol) for native and helix conformations. NOE intensities calculated from the native trajectory reproduce experimental data quite well, suggesting that the conformations sampled by the trajectory reasonably represent those observed in the NMR experiment. The molecular dynamics results, as well as sequence analysis of a diverse 690-member family of SH3 domain proteins, suggest that the presence of two elements is essential for formation of the diverging turn structure: a pair of residues with low helical propensity in the turn region and, as previously recognized, two hydrophobic residues to close the end of the diverging turn. Thus, these two sequence features may form the structural basis for the function of this peptide as a folding nucleus in this family of proteins.     

An evaluation of Poisson-Boltzmann electrostatic free energy calculations through comparison with experimental mutagenesis data

For systems involving highly and oppositely charged proteins, electrostatic forces dominate association and contribute to biomolecular complex stability. Using experimental or theoretical alanine-scanning mutagenesis, it is possible to elucidate the contribution of individual ionizable amino acids to protein association. We evaluated our electrostatic free energy calculations by comparing calculated and experimental data for alanine mutants of five protein complexes. We calculated Poisson-Boltzmann electrostatic free energies based on a thermodynamic cycle, which incorporates association in a reference (Coulombic) and solvated (solution) state, as well as solvation effects. We observe that Coulombic and solvation free energy values correlate with experimental data in highly and oppositely charged systems, but not in systems comprised of similarly charged proteins. We also observe that correlation between solution and experimental free energies is dependent on dielectric coefficient selection for the protein interior. Free energy correlations improve as protein dielectric coefficient increases, suggesting that the protein interior experiences moderate dielectric screening, despite being shielded from solvent. We propose that higher dielectric coefficients may be necessary to more accurately predict protein-protein association. Additionally, our data suggest that Coulombic potential calculations alone may be sufficient to predict relative binding of protein mutants.     

Multi‐scale characterization of the energy landscape of proteins with application to the C3D/Efb‐C complex

We present a novel multi‐level methodology to explore and characterize the low energy landscape and the thermodynamics of proteins. Traditional conformational search methods typically explore only a small portion of the conformational space of proteins and are hard to apply to large proteins due to the large amount of calculations required. In our multi‐scale approach, we first provide an initial characterization of the equilibrium state ensemble of a protein using an efficient computational conformational sampling method. We then enrich the obtained ensemble by performing short Molecular Dynamics (MD) simulations on selected conformations from the ensembles as starting points. To facilitate the analysis of the results, we project the resulting conformations on a low‐dimensional landscape to efficiently focus on important interactions and examine low energy regions. This methodology provides a more extensive sampling of the low energy landscape than an MD simulation starting from a single crystal structure as it explores multiple trajectories of the protein. This enables us to obtain a broader view of the dynamics of proteins and it can help in understanding complex binding, improving docking results and more. In this work, we apply the methodology to provide an extensive characterization of the bound complexes of the C3d fragment of human Complement component C3 and one of its powerful bacterial inhibitors, the inhibitory domain of Staphylococcus aureus extra‐cellular fibrinogen‐binding domain (Efb‐C) and two of its mutants. We characterize several important interactions along the binding interface and define low free energy regions in the three complexes. Proteins 2010. © 2009 Wiley‐Liss, Inc.     

MM/PBSA analysis of molecular dynamics simulations of bovine β‐lactoglobulin: Free energy gradients in conformational transitions?

The pH-driven opening and closure of beta-lactoglobulin EF loop, acting as a lid and closing the internal cavity of the protein, has been studied by molecular dynamics (MD) simulations and free energy calculations based on molecular mechanics/Poisson-Boltzmann (PB) solvent-accessible surface area (MM/PBSA) methodology. The forms above and below the transition pH differ presumably only in the protonation state of residue Glu89. MM/PBSA calculations are able to reproduce qualitatively the thermodynamics of the transition. The analysis of MD simulations using a combination of MM/PBSA methodology and the colony energy approach is able to highlight the driving forces implied in the transition. The analysis suggests that global rearrangements take place before the equilibrium local conformation is reached. This conclusion may bear general relevance to conformational transitions in all lipocalins and proteins in general.     

Linear response theory: An alternative to PB and GB methods for the analysis of molecular dynamics trajectories?

We explore the use of classical Linear Response Theory (LRT) as an alternative strategy to the use of Molecular Mechanics/Poisson-Boltzmann strategies to compute the solvation free energy of macromolecules from molecular dynamics simulations using an explicit representation of solvent. The method reproduces well the free energy of solvation of standard amino acid side chains, small peptides, and proteins. The use of a fully discrete representation of solvent avoids the possible problems of continuum models to represent the solvation of systems containing tightly bound water molecules.