Parametrization of a Generalized Born/Solvent-Accessible Surface Area Model and Applications to the Simulation of Protein Dynamics

ROS1 kinase inhibition reimagined: identifying repurposed drug via virtual screening and molecular dynamics simulations for cancer therapeutics

Precision medicine has revolutionized modern cancer therapeutic management by targeting specific molecular aberrations responsible for the onset and progression of tumorigenesis. ROS proto-oncogene 1 (ROS1) is a receptor tyrosine kinase (RTK) that can induce tumorigenesis through various signaling pathways, such as cell proliferation, survival, migration, and metastasis. It has emerged as a promising therapeutic target in various cancer types. However, there is very limited availability of specific ROS1 inhibitors for therapeutic purposes. Exploring repurposed drugs for rapid and effective treatment is a useful approach. In this study, we utilized an integrated approach of virtual screening and molecular dynamics (MD) simulations of repurposing existing drugs for ROS1 kinase inhibition. Using a curated library of 3648 FDA-approved drugs, virtual screening identified drugs capable of binding to ROS1 kinase domain. The results unveil two hits, Midostaurin and Alectinib with favorable binding profiles and stable interactions with the active site residues of ROS1. These hits were subjected to stability assessment through all-atom MD simulations for 200 ns. MD results showed that Midostaurin and Alectinib were stable with ROS1. Taken together, the study showed a rational framework for the selection of repurposed Midostaurin and Alectinib with ROS1 inhibitory potential for therapeutic management after further validation.

Increased Antiaromaticity through Pentalene Connection in [n]Cyclo-1,5-dibenzopentalenes

Conjugated nanohoops incorporating nonalternant hydrocarbons have altered optoelectronic properties compared to [n]cycloparaphenylenes or other purely aromatic hoops. We synthesized [n]cyclo-1,5-dibenzopentalenes (n = 4, 5), in which nonalternant dibenzo[a,e]pentalenes are connected through their pentalene units. This leads to an increase in antiaromatic character and low-lying LUMO energies. Calculations show puckered or entangled conformations of the precursor macrocyclic Pt-complexes. Our study proves dibenzopentalene as a versatile nonalternant building block for conjugated nanohoops with modifiable antiaromaticity and optoelectronic properties.

Chiral Dibenzopentalene-Based Conjugated Nanohoops through Stereoselective Synthesis

Conjugated nanohoops allow to investigate the effect of radial conjugation and bending on the involved π-systems. They can possess unexpected optoelectronic properties and their radially oriented π-system makes them attractive for host-guest chemistry. Bending the π-subsystems can lead to chiral hoops. Herein, we report the stereoselective synthesis of two enantiomers of chiral conjugated nanohoops by incorporating dibenzo[a,e]pentalenes (DBPs), which are generated in the last synthetic step from enantiomerically pure diketone precursors. Owing to its bent shape, this diketone unit was used as the only bent precursor and novel "corner unit" in the synthesis of the hoops. The [6]DBP[4]Ph-hoops contain six antiaromatic DBP units and four bridging phenylene groups. The small HOMO-LUMO gap and ambipolar electrochemical character of the DBP units is reflected in the optoelectronic properties of the hoop. Electronic circular dichroism spectra and MD simulations showed that the chiral hoop did not racemize even when heated to 110 °C. Due to its large diameter, it was able to accommodate two C60 molecules, as binding studies indicate.

Oxidation-induced destabilization of the fibrinogen αC-domain dimer investigated by molecular dynamics simulations

Upon activation, fibrinogen is converted to insoluble fibrin, which assembles into long strings called protofibrils. These aggregate laterally to form a fibrin matrix that stabilizes a blood clot. Lateral aggregation of protofibrils is mediated by the αC domain, a partially structured fragment located in a disordered region of fibrinogen. Polymerization of αC domains links multiple fibrin molecules with each other enabling the formation of thick fibrin fibers and a fibrin matrix that is stable but can also be digested by enzymes. However, oxidizing agents produced during the inflammatory response have been shown to cause thinner fibrin fibers resulting in denser clots, which are harder to proteolyze and pose the risk of deep vein thrombosis and lung embolism. Oxidation of Met⁴⁷⁶ located within the αC domain is thought to hinder its ability to polymerize disrupting the lateral aggregation of protofibrils and leading to the observed thinner fibers. How αC domains assemble into polymers is still unclear and yet this knowledge would shed light on the mechanism through which oxidation weakens the lateral aggregation of protofibrils. This study used temperature replica exchange molecular dynamics simulations to investigate the αC-domain dimer and how this is affected by oxidation of Met⁴⁷⁶ . Analysis of the trajectories revealed that multiple stable binding modes were sampled between two αC domains while oxidation decreased the likelihood of dimer formation. Furthermore, the side chain of Met⁴⁷⁶ was observed to act as a docking spot for the binding and this function was impaired by its conversion to methionine sulfoxide.

Molecular dynamics simulation reveals the possible druggable hot-spots of USP7

The plasticity in Ubiquitin Specific Proteases (USP7) inducing conformational changes at important areas has highlighted an intricate mechanism, by which USP7 is regulated. Given the importance of USP7 in oncogenic pathways and immune-oncology, identification of USP7 inhibitors has attracted considerable interest. Despite substantial efforts, the discovery of deubiquitinases (DUBs) inhibitors, knowledge of their binding site and understanding the possible mechanism of action has proven particularly challenging. We disclose the most likely binding site of P5091 (a potent USP7 inhibitor), which reveal a cryptic allosteric site through extensive computational studies in an inhibitor dependent and independent manner. Overall, these findings demonstrate the tractability and druggability of USP7. Through a series of molecular dynamics simulations and detailed quantitative analysis, a dynamically stable allosteric binding site near catalytic center of the inactive state of USP7 (site partially absent in active state), along with two newly identified sites have been revealed, which opens the avenue for rational structure-guided inhibitor designing in USP7 specific-manner.

Electrostatic Interactions in Protein Structure, Folding, Binding, and Condensation

Charged and polar groups, through forming ion pairs, hydrogen bonds, and other less specific electrostatic interactions, impart important properties to proteins. Modulation of the charges on the amino acids, e.g., by pH and by phosphorylation and dephosphorylation, have significant effects such as protein denaturation and switch-like response of signal transduction networks. This review aims to present a unifying theme among the various effects of protein charges and polar groups. Simple models will be used to illustrate basic ideas about electrostatic interactions in proteins, and these ideas in turn will be used to elucidate the roles of electrostatic interactions in protein structure, folding, binding, condensation, and related biological functions. In particular, we will examine how charged side chains are spatially distributed in various types of proteins and how electrostatic interactions affect thermodynamic and kinetic properties of proteins. Our hope is to capture both important historical developments and recent experimental and theoretical advances in quantifying electrostatic contributions of proteins.

Molecular mechanism of tobramycin with human serum albumin for probing binding interactions: multi-spectroscopic and computational approaches

Highlighting novelty: comprehensive in vitro and in silico insights for understanding the novel binding site of TOB with HSA.

Incorporating excluded solvent volume and physical dipoles for computing solvation free energy

The solvation free energy described using the Born equation depends on the solute charge, solute radius, and solvent dielectric constant. However, the dielectric polarization derived from Gauss's law used in the Born equation differs from that obtained from molecular dynamics simulations. Therefore, the adjustment of Born radii is insufficient for fitting the solvation free energy to various solute conformations. In order to mimic the dielectric polarization surrounding a solute in molecular dynamics simulations, the water molecule in the first coordination shell is modeled as a physical dipole in a van der Waals sphere, and the intermediate water is treated as a bulk solvent. The electric dipole of the first-shell water is modeled as positive and negative surface charge layers with fixed charge magnitudes, but with variable separation distance as derived from the distributions of hydrogen and oxygen atoms of water dictated by their orientational distribution functions. An equation that describes the solvation free energy of ions using this solvent scheme with a TIP3P water model is derived, and the values of the solvation free energies of ions estimated from this derived equation are found to be similar to those obtained from the experimental data.

Structure Based Modeling of Small Molecules Binding to the TLR7 by Atomistic Level Simulations

Toll-Like Receptors (TLR) are a large family of proteins involved in the immune system response. Both the activation and the inhibition of these receptors can have positive effects on several diseases, including viral pathologies and cancer, therefore prompting the development of new compounds. In order to provide new indications for the design of Toll-Like Receptor 7 (TLR7)-targeting drugs, the mechanism of interaction between the TLR7 and two important classes of agonists (imidazoquinoline and adenine derivatives) was investigated through docking and Molecular Dynamics simulations. To perform the computational analysis, a new model for the dimeric form of the receptors was necessary and therefore created. Qualitative and quantitative differences between agonists and inactive compounds were determined. The in silico results were compared with previous experimental observations and employed to define the ligand binding mechanism of TLR7.

Incorporating the excluded solvent volume and surface charges for computing solvation free energy

Gauss's law or Poisson's equation is conventionally used to calculate solvation free energy. However, the near-solute dielectric polarization from Gauss's law or Poisson's equation differs from that obtained from molecular dynamics (MD) simulations. To mimic the near-solute dielectric polarization from MD simulations, the first-shell water was treated as two layers of surface charges, the densities of which are proportional to the electric field at the solvent molecule that is modeled as a hard sphere. The intermediate water was treated as a bulk solvent. An equation describing the solvation free energy of ions using this solvent scheme was derived using the TIP3P water model.

Density functional theory for protein transfer free energy

We cast the problem of protein transfer free energy within the formalism of density functional theory (DFT), treating the protein as a source of external potential that acts upon the solvent. Solvent excluded volume, solvent-accessible surface area, and temperature dependence of the transfer free energy all emerge naturally within this formalism, and may be compared with simplified "back of the envelope" models, which are also developed here. Depletion contributions to osmolyte induced stability range from 5 to 10 kBT for typical protein lengths. The general DFT transfer theory developed here may be simplified to reproduce a Langmuir isotherm condensation mechanism on the protein surface in the limits of short-ranged interactions, and dilute solute. Extending the equation of state to higher solute densities results in non-monotonic behavior of the free energy driving protein or polymer collapse. Effective interaction potentials between protein backbone or side chains and TMAO are obtained, assuming a simple backbone/side chain two-bead model for the protein with an effective 6-12 potential with the osmolyte. The transfer free energy δg shows significant entropy: d(δg)/dT ≈ 20 kB for a 100-residue protein. The application of DFT to effective solvent forces for use in implicit-solvent molecular dynamics is also developed. The simplest DFT expressions for implicit-solvent forces contain both depletion interactions and an "impeded-solvation" repulsive force at larger distances.

Dependence of interaction free energy between solutes on an external electrostatic field

To explore the athermal effect of an external electrostatic field on the stabilities of protein conformations and the binding affinities of protein-protein/ligand interactions, the dependences of the polar and hydrophobic interactions on the external electrostatic field, -Eext, were studied using molecular dynamics (MD) simulations. By decomposing Eext into, along, and perpendicular to the direction formed by the two solutes, the effect of Eext on the interactions between these two solutes can be estimated based on the effects from these two components. Eext was applied along the direction of the electric dipole formed by two solutes with opposite charges. The attractive interaction free energy between these two solutes decreased for solutes treated as point charges. In contrast, the attractive interaction free energy between these two solutes increased, as observed by MD simulations, for Eext = 40 or 60 MV/cm. Eext was applied perpendicular to the direction of the electric dipole formed by these two solutes. The attractive interaction free energy was increased for Eext = 100 MV/cm as a result of dielectric saturation. The force on the solutes along the direction of Eext computed from MD simulations was greater than that estimated from a continuum solvent in which the solutes were treated as point charges. To explore the hydrophobic interactions, Eext was applied to a water cluster containing two neutral solutes. The repulsive force between these solutes was decreased/increased for Eext along/perpendicular to the direction of the electric dipole formed by these two solutes.

GPU/CPU Algorithm for Generalized Born/Solvent-Accessible Surface Area Implicit Solvent Calculations

Molecular dynamics methodologies comprise a vital research tool for structural biology. Molecular dynamics has benefited from technological advances in computing, such as multi-core CPUs and graphics processing units (GPUs), but harnessing the full power of hybrid GPU/CPU computers remains difficult. The generalized Born/solvent-accessible surface area implicit solvent model (GB/SA) stands to benefit from hybrid GPU/CPU computers, employing the GPU for the GB calculation and the CPU for the SA calculation. Here, we explore the computational challenges facing GB/SA calculations on hybrid GPU/CPU computers and demonstrate how NAMD, a parallel molecular dynamics program, is able to efficiently utilize GPUs and CPUs simultaneously for fast GB/SA simulations. The hybrid computation principles demonstrated here are generally applicable to parallel applications employing hybrid GPU/CPU calculations.

Configurational subdiffusion of peptides: a network study

Molecular dynamics (MD) simulation of linear peptides reveals configurational subdiffusion at equilibrium extending from 10⁻¹² to 10⁻⁸ s. Rouse chain and continuous-time random walk models of the subdiffusion are critically discussed. Network approaches to analyzing MD simulations are shown to reproduce the time dependence of the subdiffusive mean squared displacement, which is found to arise from the fractal-like geometry of the accessible volume in the configuration space. Convergence properties of the simulation pertaining to the subdiffusive dynamics are characterized and the effect on the subdiffusive properties of representing the solvent explicitly or implicitly is compared. Non-Markovianity and other factors limiting the range of applicability of the network models are examined.

Energetic decomposition with the generalized-born and Poisson-Boltzmann solvent models: lessons from association of G-protein components

Continuum electrostatic models have been shown to be powerful tools in providing insight into the energetics of biomolecular processes. While the Poisson-Boltzmann (PB) equation provides a theoretically rigorous approach to computing electrostatic free energies of solution in such a model, computational cost makes its use for large ensembles of states impractical. The generalized-Born (GB) approximation provides a much faster alternative, although with a weaker theoretical framework. While much attention has been given to how GB recapitulates PB energetics for the overall stability of a biomolecule or the affinity of a complex, little attention has been given to how the contributions of individual functional groups are captured by the two methods. Accurately capturing these individual electrostatic components is essential both for the development of a mechanistic understanding of biomolecular processes and for the design of variant sequences and structures with desired properties. Here, we present a detailed comparison of the group-wise decomposition of both PB and GB electrostatic free energies of binding, using association of various components of the heterotrimeric-G-protein complex as a model. We find that, while net binding free energies are strongly correlated in the two models, the correlations of individual group contributions are highly variable; in some cases, strong correlation is seen, while in others, there is essentially none. Structurally, the GB model seems to capture the magnitude of direct, short-range electrostatic interactions quite well but performs more poorly with moderate-range "action-at-a-distance" interactions--GB has a tendency to overestimate solvent screening over moderate distances, and to underestimate the costs of desolvating charged groups somewhat removed from the binding interface. Despite this, however, GB does seem to be quite effective as a predictor of those groups that will be computed to be most significant in a PB-based model.

Dehydration-driven solvent exposure of hydrophobic surfaces as a driving force in peptide folding

Recent work has shown that the nature of hydration of pure hydrophobic surfaces changes with the length scale considered: water hydrogen-bonding networks adapt to small exposed hydrophobic species, hydrating or "wetting" them at relatively high densities, whereas larger hydrophobic areas are "dewetted" [Chandler D (2005), Nature 29:640-647]. Here we determine whether this effect is also present in peptides by examining the folding of a beta-hairpin (the 14-residue amyloidogenic prion protein H1 peptide), using microsecond time-scale molecular dynamics simulations. Two simulation models are compared, one explicitly including the water molecules, which may thus adapt locally to peptide configurations, and the other using a popular continuum approximation, the generalized Born/surface area implicit solvent model. The results obtained show that, in explicit solvent, peptide conformers with high solvent-accessible hydrophobic surface area indeed also have low hydration density around hydrophobic residues, whereas a concomitant higher hydration density around hydrophilic residues is observed. This dewetting effect stabilizes the fully folded beta-hairpin state found experimentally. In contrast, the implicit solvent model destabilizes the fully folded hairpin, tending to cluster hydrophobic residues regardless of the size of the exposed hydrophobic surface. Furthermore, the rate of the conformational transitions in the implicit solvent simulation is almost doubled with respect to that of the explicit solvent. The results suggest that dehydration-driven solvent exposure of hydrophobic surfaces may be a significant factor determining peptide conformational equilibria.

Comparative study of generalized born models: Born radii and peptide folding

In this study, we have implemented four analytical generalized Born (GB) models and investigated their performance in conjunction with the GROMOS96 force field. The four models include that of Still and co-workers, the HCT model of Cramer, Truhlar, and co-workers, a modified form of the AGB model of Levy and co-workers, and the GBMV2 model of Brooks and co-workers. The models were coded independently and implemented in the GROMOS software package and in TINKER. They were compared in terms of their ability to reproduce the results of Poisson-Boltzmann (PB) calculations and in their performance in the ab initio peptide folding of two peptides, one that forms a beta-hairpin in solution and one that forms an alpha-helix. In agreement with previous work, the GBMV2 model is most successful in reproducing PB results while the other models tend to underestimate the effective Born radii of buried atoms. In contrast, stochastic dynamics simulations on the folding of the two peptides, the C-terminus beta-hairpin of the B1 domain of protein G and the alanine-based alpha-helical peptide 3K(I), suggest that the simpler GB models are more effective in sampling conformational space. Indeed, the Still model used in conjunction with the GROMOS96 force field is able to fold the hairpin peptide to a native-like structure without the benefit of enhanced sampling techniques. This is due in part to the properties of the united-atom GROMOS96 force field which appears to be more flexible, and hence to sample more efficiently, than force fields such as OPLSAA. Our results suggest a general strategy which involves using different combinations of force fields and solvent models in different applications, for example, using GROMOS96 and a simple GB model in sampling and OPLSAA and a more accurate GB model in refinement. The fact that various methods have been implemented in a unified way should facilitate the testing and subsequent use of different methods to evaluate conformational free energies in different applications. Our results also bear on some general issues involved in peptide folding and structure prediction which are addressed in the Discussion.

A generalized Born formalism for heterogeneous dielectric environments: application to the implicit modeling of biological membranes

Reliable computer simulations of complex biological environments such as integral membrane proteins with explicit water and lipid molecules remain a challenging task. We propose a modification of the standard generalized Born theory of homogeneous solvent for modeling the heterogeneous dielectric environments such as lipid/water interfaces. Our model allows the representation of biological membranes in the form of multiple layered dielectric regions with dielectric constants that are different from the solute cavity. The proposed new formalism is shown to predict the electrostatic component of solvation free energy with a relative error of 0.17% compared to exact finite-difference solutions of the Poisson equation for a transmembrane helix test system. Molecular dynamics simulations of melittin and bacteriorhodopsin are carried out and performed over 10 ns and 7 ns of simulation time, respectively. The center of melittin along the membrane normal in these stable simulations is in excellent agreement with the relevant experimental data. Simulations of bacteriorhodopsin started from the experimental structure remained stable and in close agreement with experiment. We also examined the free energy profiles of water and amino acid side chain analogs upon membrane insertion. The results with our implicit membrane model agree well with the experimental transfer free energy data from cyclohexane to water as well as explicit solvent simulations of water and selected side chain analogs.

Implicit solvation based on generalized Born theory in different dielectric environments

In this paper we are investigating the effect of the dielectric environment on atomic Born radii used in generalized Born (GB) methods. Motivated by the Kirkwood expression for the reaction field of a single off-center charge in a spherical cavity, we are proposing extended formalisms for the calculation of Born radii as a function of external and internal dielectric constants. We demonstrate that reaction field energies calculated from environmentally dependent Born radii lead to much improved agreement with Poisson-Boltzmann solutions for low dielectric external environments, such as biological membranes or organic solvent, compared to previous methods where the calculation of Born radii does not depend on the environment. We also examine how this new approach can be applied for the calculation of transfer free energies from vacuum to a given external dielectric for a system with an internal dielectric larger than one. This has not been possible with standard GB theory but is relevant when scoring minimized or average structures with implicit solvent.

Transition State Ensemble for the Folding of B Domain of Protein A: A Comparison of Distributed Molecular Dynamics Simulations with Experiments

Folding pathways of the B domain of staphylococcal protein A have been sampled with a distributed computing approach. Starting from an extended structure, the method employs an index measuring topological similarity to the native structure to selectively sample trajectory branches leading to the native fold. Unperturbed and continuous folding trajectories are drawn on a physics-based atomic potential energy surface with an implicit solvent. The sampled folding trajectories demonstrate a similar sequence of events: the earlier stage involves a partial formation of helix 2 and to a less extent of helix 1 at their N terminals, followed by the hydrophobic collapse between residues F14, I17, and L18 on helix 1 and residues R28, F31, and I32 on helix 2, which results in the rigidification of the helix turn from R28 to I32. Helix 2 is then able to extend, allowing for the formation to turn 2. The above description explains one experimental result why a G30A mutant of the protein was observed to be the fastest folder among proteins of its size. And the ensemble of structures right before the final collapse is in good agreement with the transition state ensemble mapped by another recent experiment with Fersht Phi values. We emphasize that because the approach here does not provide quantifications of the free energy landscape, our model of the transition state ensemble emerges from comparisons of simulations and previous experimental results rather from the simulation results alone. On the other hand, as our approach does not rely on a low-dimensional free energy surface, it can complement methods based on the construction of free energy surfaces.

Recent advances in the development and application of implicit solvent models in biomolecule simulations

Advances have recently been made in the development of implicit solvent methodologies and their application to the modeling of biomolecules, particularly with regard to generalized Born approaches, dielectric screening function formulations and models based on solvent-accessible surface areas. Interesting new developments include more refined non-polar solvation energy estimators, and implicit methods for modeling low-dielectric and heterogeneous environments such as membrane systems. These have been successfully applied to molecular dynamics simulations, the scoring of protein conformations, and the calculation of binding affinities and folding free energy landscapes.

Comparative study of generalized Born models: protein dynamics

In this work, we compare the results of molecular dynamics simulations involving the application of three generalized Born (GB) models to 10 different proteins. The three GB models, the Still, HCT, and modified analytical generalized Born models, were implemented in the computationally efficient gromacs package. The performance of each model was assessed from the backbone rms deviation from the native structure, the number of native hydrogen bonds retained in the simulation, and the experimental and calculated radius of gyration. Analysis of variance (ANOVA) was used to analyze the results of the simulations. The rms deviation measure was found to be unable to distinguish the quality of the results obtained with the three different GB models, whereas the number of native hydrogen bonds and radius of gyration yielded a statistically meaningful discrimination among models. Our results suggest that, of the three, modified analytical generalized Born yields the best agreement between calculated and experimentally derived structures. More generally, our study highlights the need both to evaluate the effects of different variables on the results of simulations and to verify that the results of molecular dynamics simulations are statistically meaningful.

Solution structure of AF-6 PDZ domain and its interaction with the C-terminal peptides from Neurexin and Bcr

AF-6 is a key molecule essential for structure organization of cell-cell junction of polarized epithelia. It belongs to a novel cell-cell adhesion system. The AF-6 PDZ domain mediates interactions by binding to a specific amino acid sequence in target proteins. Here we report the solution structure of the AF-6 PDZ domain determined by NMR. Previously, the AF-6 PDZ domain was considered to be a class II PDZ domain. However we found that a unique hydrophilic amino acid, Gln70, at position alphaB1 makes the alphaB/betaB groove of the AF-6 PDZ domain significantly different from that of the canonical class II PDZ domain. The AF-6 PDZ domain does not have the second hydrophobic binding pocket, and the N-terminal end of alphaB is closer to betaB. Using BIACORE and NMR chemical shift perturbation experiments, we have studied the binding characteristics of the PDZ domain to the C-terminal peptide of Neurexin, KKNKDKEYYV, and that of Bcr, KRQSILFSTEV. The C-terminal peptide of Neurexin is a class II ligand, whereas that of Bcr is a class I ligand. The dissociation constants of these ligands were 4.08 x 10(-7) and 2.23 x 10(-6) m, respectively. Each of the four C-terminal positions in Neurexin and Bcr may contribute to the interaction. The three-dimensional models of the AF-6 PDZ-Neurexin C-terminal peptide complex and the AF-6 PDZ-Bcr C-terminal peptide complex were built up by molecular dynamics simulations. Unlike the canonical class II PDZ domain, Ala74 at alphaB5 rather than the residue at alphaB1 makes direct hydrophobic contact with the side chain of Tyr at the -2 position of the ligand.

Free energy surfaces of β-hairpin and α-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.

Performance comparison of generalized born and Poisson methods in the calculation of electrostatic solvation energies for protein structures

This study compares generalized Born (GB) and Poisson (PB) methods for calculating electrostatic solvation energies of proteins. A large set of GB and PB implementations from our own laboratories as well as others is applied to a series of protein structure test sets for evaluating the performance of these methods. The test sets cover a significant range of native protein structures of varying size, fold topology, and amino acid composition as well as nonnative extended and misfolded structures that may be found during structure prediction and folding/unfolding studies. We find that the methods tested here span a wide range from highly accurate and computationally demanding PB-based methods to somewhat less accurate but more affordable GB-based approaches and a few fast, approximate PB solvers. Compared with PB solvation energies, the latest, most accurate GB implementations were found to achieve errors of 1% for relative solvation energies between different proteins and 0.4% between different conformations of the same protein. This compares to accurate PB solvers that produce results with deviations of less than 0.25% between each other for both native and nonnative structures. The performance of the best GB methods is discussed in more detail for the application for force field-based minimizations or molecular dynamics simulations.

Continuum description of ionic and dielectric shielding for molecular-dynamics simulations of proteins in solution

We extend our continuum description of solvent dielectrics in molecular-dynamics (MD) simulations, which has provided an efficient and accurate solution of the Poisson equation, to ionic solvents as described by the linearized Poisson-Boltzmann (LPB) equation. We start with the formulation of a general theory for the electrostatics of an arbitrarily shaped molecular system, which consists of partially charged atoms and is embedded in a LPB continuum. This theory represents the reaction field induced by the continuum in terms of charge and dipole densities localized within the molecular system. Because these densities cannot be calculated analytically for systems of arbitrary shape, we introduce an atom-based discretization and a set of carefully designed approximations. This allows us to represent the densities by charges and dipoles located at the atoms. Coupled systems of linear equations determine these multipoles and can be rapidly solved by iteration during a MD simulation. The multipoles yield the reaction field forces and energies. Finally, we scrutinize the quality of our approach by comparisons with an analytical solution restricted to perfectly spherical systems and with results of a finite difference method.

How well can we predict native contacts in proteins based on decoy structures and their energies?

One strategy for ab initio protein structure prediction is to generate a large number of possible structures (decoys) and select the most fitting ones based on a scoring or free energy function. The conformational space of a protein is huge, and chances are rare that any heuristically generated structure will directly fall in the neighborhood of the native structure. It is desirable that, instead of being thrown away, the unfitting decoy structures can provide insights into native structures so prediction can be made progressively. First, we demonstrate that a recently parameterized physics-based effective free energy function based on the GROMOS96 force field and a generalized Born/surface area solvent model is, as several other physics-based and knowledge-based models, capable of distinguishing native structures from decoy structures for a number of widely used decoy databases. Second, we observe a substantial increase in correlations of the effective free energies with the degree of similarity between the decoys and the native structure, if the similarity is measured by the content of native inter-residue contacts in a decoy structure rather than its root-mean-square deviation from the native structure. Finally, we investigate the possibility of predicting native contacts based on the frequency of occurrence of contacts in decoy structures. For most proteins contained in the decoy databases, a meaningful amount of native contacts can be predicted based on plain frequencies of occurrence at a relatively high level of accuracy. Relative to using plain frequencies, overwhelming improvements in sensitivity of the predictions are observed for the 4_state_reduced decoy sets by applying energy-dependent weighting of decoy structures in determining the frequency. There, approximately 80% native contacts can be predicted at an accuracy of approximately 80% using energy-weighted frequencies. The sensitivity of the plain frequency approach is much lower (20% to 40%). Such improvements are, however, not observed for the other decoy databases. The rationalization and implications of the results are discussed.

Molecular dynamics simulations of peptides and proteins with amplified collective motions

We present a novel method that uses the collective modes obtained with a coarse-grained model/anisotropic network model to guide the atomic-level simulations. Based on this model, local collective modes can be calculated according to a single configuration in the conformational space of the protein. In the molecular dynamics simulations, the motions along the slowest few modes are coupled to a higher temperature by the weak coupling method to amplify the collective motions. This amplified-collective-motion (ACM) method is applied to two test systems. One is an S-peptide analog. We realized the refolding of the denatured peptide in eight simulations out of 10 using the method. The other system is bacteriophage T4 lysozyme. Much more extensive domain motions between the N-terminal and C-terminal domain of T4 lysozyme are observed in the ACM simulation compared to a conventional simulation. The ACM method allows for extensive sampling in conformational space while still restricting the sampled configurations within low energy areas. The method can be applied in both explicit and implicit solvent simulations, and may be further applied to important biological problems, such as long timescale functional motions, protein folding/unfolding, and structure prediction.

Molecular dynamics simulations of peptides and proteins with a continuum electrostatic model based on screened Coulomb potentials

A continuum electrostatics approach for molecular dynamics (MD) simulations of macromolecules is presented and analyzed for its performance on a peptide and a globular protein. The approach incorporates the screened Coulomb potential (SCP) continuum model of electrostatics, which was reported earlier. The model was validated in a broad set of tests some of which were based on Monte Carlo simulations that included single amino acids, peptides, and proteins. The implementation for large-scale MD simulations presented in this article is based on a pairwise potential that makes the electrostatic model suitable for fast analytical calculation of forces. To assess the suitability of the approach, a preliminary validation is conducted, which consists of (i) a 3-ns MD simulation of the immunoglobulin-binding domain of streptococcal protein G, a 56-residue globular protein and (ii) a 3-ns simulation of Dynorphin, a biological peptide of 17 amino acids. In both cases, the results are compared with those obtained from MD simulations using explicit water (EW) molecules in an all-atom representation. The initial structure of Dynorphin was assumed to be an alpha-helix between residues 1 and 9 as suggested from NMR measurements in micelles. The results obtained in the MD simulations show that the helical structure collapses early in the simulation, a behavior observed in the EW simulation and consistent with spectroscopic data that suggest that the peptide may adopt mainly an extended conformation in water. The dynamics of protein G calculated with the SCP implicit solvent model (SCP-ISM) reveals a stable structure that conserves all the elements of secondary structure throughout the entire simulation time. The average structures calculated from the trajectories with the implicit and explicit solvent models had a cRMSD of 1.1 A, whereas each average structure had a cRMSD of about 0.8A with respect to the X-ray structure. The main conformational differences of the average structures with respect to the crystal structure occur in the loop involving residues 8-14. Despite the overall similarity of the simulated dynamics with EW and SCP models, fluctuations of side-chains are larger when the implicit solvent is used, especially in solvent exposed side-chains. The MD simulation of Dynorphin was extended to 40 ns to study its behavior in an aqueous environment. This long simulation showed that the peptide has a tendency to form an alpha-helical structure in water, but the stabilization free energy is too weak, resulting in frequent interconversions between random and helical conformations during the simulation time. The results reported here suggest that the SCP implicit solvent model is adequate to describe electrostatic effects in MD simulation of both peptides and proteins using the same set of parameters. It is suggested that the present approach could form the basis for the development of a reliable and general continuum approach for use in molecular biology, and directions are outlined for attaining this long-term goal.