Explicit Solvent Molecular Dynamics Research Articles

A molecular reconstruction approach to site-based 3D-RISM and comparison to GIST hydration thermodynamic maps in an enzyme active site.

Computed, high-resolution, spatial distributions of solvation energy and entropy can provide detailed information about the role of water in molecular recognition. While grid inhomogeneous solvation theory (GIST) provides rigorous, detailed thermodynamic information from explicit solvent molecular dynamics simulations, recent developments in the 3D reference interaction site model (3D-RISM) theory allow many of the same quantities to be calculated in a fraction of the time. However, 3D-RISM produces atomic-site, rather than molecular, density distributions, which are difficult to extract physical meaning from. To overcome this difficulty, we introduce a method to reconstruct molecular density distributions from atomic-site density distributions. Furthermore, we assess the quality of the resulting solvation thermodynamics density distributions by analyzing the binding site of coagulation Factor Xa with both GIST and 3D-RISM. We find good qualitative agreement between the methods for oxygen and hydrogen densities as well as direct solute-solvent energetic interactions. However, 3D-RISM predicts lower energetic and entropic penalties for moving water from the bulk to the binding site.

Edaravone inhibits the conformational transition of amyloid-β42: insights from molecular dynamics simulations

Previous work has shown that edaravone inhibits fibrillogenesis of amyloid-β protein (Aβ). However, the detailed mechanism by which edaravone inhibits the conformational transition of the Aβ42 monomer is not known at the molecular level. Here, explicit-solvent molecular dynamics (MD) simulations were coupled with molecular mechanics–Poisson–Boltzmann surface area (MM–PBSA) method to address the issue. MD simulations confirmed that edaravone inhibits the conformational transition of the Aβ42 monomer in a dose-dependent manner. It was found that the direct interactions between edaravone and Aβ42 are responsible for its inhibiting effects. The analysis of binding free energy using the MM–PBSA method demonstrated that the nonpolar interactions provide favourable contributions (about −71.7 kcal/mol). Conversely, the polar interactions are unfavourable for the binding process. A total of 14 residues were identified as greatly contributing to the binding free energy between edaravone and the Aβ42 monomer. In addition, the intra-peptide hydrophobic interactions were weakened and the salt bridge D23–K28 was interrupted by edaravone. Therefore, the conformational transition was inhibited. Our studies provide molecular-level insights into how edaravone molecules inhibit the conformational transition of the Aβ42 monomer, which may be useful for designing amyloid inhibitors.Communicated by Ramaswamy H. Sarma

Conformational coupling by trans-phosphorylation in calcium calmodulin dependent kinase II.

The calcium calmodulin-dependent protein kinase II (CaMKII) is a dodecameric holoenzyme important for encoding memory. Its activation, triggered by binding of calcium-calmodulin, persists autonomously after calmodulin dissociation. One (receiver) kinase captures and subsequently phosphorylates the regulatory domain peptide of a donor kinase forming a chained dimer as the first stage of autonomous activation. Protein dynamics simulations examined the conformational changes triggered by dimer formation and phosphorylation, aimed to provide a molecular rationale for human mutations that result in learning disabilities. Ensembles generated from X-ray crystal structures were characterized by network centrality and community analysis. Mutual information related collective motions to local fragment dynamics encoded with a structural alphabet. Implicit solvent tCONCOORD conformational ensembles revealed the dynamic architecture of inactive kinase domains was co-opted in the activated dimer but the network hub shifted from the nucleotide binding cleft to the captured peptide. Explicit solvent molecular dynamics (MD) showed nucleotide and substrate binding determinants formed coupled nodes in long-range signal relays between regulatory peptides in the dimer. Strain in the extended captured peptide was balanced by reduced flexibility of the receiver kinase C-lobe core. The relays were organized around a hydrophobic patch between the captured peptide and a key binding helix. The human mutations aligned along the relays. Thus, these mutations could disrupt the allosteric network alternatively, or in addition, to altered binding affinities. Non-binding protein sectors distant from the binding sites mediated the allosteric signalling; providing possible targets for inhibitor design. Phosphorylation of the peptide modulated the dielectric of its binding pocket to strengthen the patch, non-binding sectors, domain interface and temporal correlations between parallel relays. These results provide the molecular details underlying the reported positive kinase cooperativity to enrich the discussion on how autonomous activation by phosphorylation leads to long-term behavioural effects.

Characterization of Hydration Properties in Structural Ensembles of Biomolecules.

Solute-solvent interactions are critical for biomolecular stability and recognition. Explicit solvent molecular dynamics (MD) simulations are routinely used to probe such interactions. However, detailed analyses and interpretation of the hydration patterns seen in MD simulations can be both complex and time-consuming. A variety of approaches/tools to compute and interrogate hydration properties in structural ensembles of proteins, nucleic acids, or in general any molecule are available and are complemented here with a new and free software package ("JAL"). Central to "JAL" is an intuitive atom centric approach of computing hydration properties. In addition to the standard metrics commonly used to understand hydration, "JAL" introduces two nonstandard utilities: a program to rapidly compute buried waters in an MD trajectory and a new method to compute multiwater bridges around a solute. We demonstrate the utility of the package by probing the hydration characteristics of the tumor suppressor protein p53 and the translation initiation factor eif4E. "JAL" is hosted online and can be accessed for free at http://mspc.bii.a-star.edu.sg/minhn/jal.html .

"Mirror"-like Protein Dimers Stabilized by Local Heterogeneity at Protein Surfaces.

Protein aggregation has been observed inside cells and holds true for membraneless organelles. The precise understanding of protein dimerization is a prerequisite for manipulating protein aggregation, which is promising for elevating enzyme concentration to enhance their catalytic performance. Here, the dimerization of two industrially important enzymes of cytochrome P450 (P450) and organophosphorus hydrolase (OPH) is investigated using all-atom explicit solvent molecular dynamics simulations, umbrella sampling, and protein-protein docking calculations. The calculated potentials of mean force of dimer-monomer dissociation demonstrate that the dimeric forms are more stable with the free energy barrier of around 60 kJ/mol for P450 and 101 kJ/mol for OPH. The docking calculations on the OPH dimer evidence the uniqueness of the native orientation. The protein dimers form "mirror"-like orientations with some degree of rotation. Such signature orientations are interpreted based on the predominant polar amino acids in the contact regime. In the dimer conformations, the active sites are exposed. This work highlights the crucial roles of the polar and nonpolar protein surface domains to form enzymatically active protein dimer aggregates. Our work will potentially aid the design of molecules that can deliver and protect native protein function in various environments.

Haptic-Assisted Interactive Molecular Docking Incorporating Receptor Flexibility.

Haptic-assisted interactive docking tools immerse the user in an environment where intuition and knowledge can be used to help guide the docking process. Here we present such a tool where the user "holds" a rigid ligand via a haptic device through which they feel interaction forces with a flexible receptor biomolecule. To ensure forces transmitted through the haptic device are smooth and stable, they must be updated at a rate greater than 500 Hz. Due to this time constraint, the majority of haptic docking tools do not attempt to model the conformational changes that would occur when molecules interact during binding. Our haptic-assisted docking tool, "Haptimol FlexiDock", models a receptor's conformational response to forces of interaction with a ligand while maintaining the required haptic refresh rate. In order to model receptor flexibility we use the method of linear response for which we determine the variance-covariance matrix of atomic fluctuations from the trajectory of an explicit-solvent molecular dynamics simulation of the ligand-free receptor molecule. The key to satisfying the time constraint is an eigenvector decomposition of the variance-covariance matrix which enables a good approximation to the conformational response of the receptor to be calculated rapidly. This exploits a feature of protein dynamics whereby most fluctuation occurs within a relatively small subspace. The method is demonstrated on glutamine binding protein in interaction with glutamine and maltose binding protein in interaction with maltose. For both proteins the movement that occurs when the ligand is docked near to its binding site matches the experimentally determined movement well. It is thought that this tool will be particularly useful for structure-based drug design.

Mean first passage times in variational coarse graining using Markov state models.

Markov state models (MSMs) provide some of the simplest mathematical and physical descriptions of dynamical and thermodynamical properties of complex systems. However, typically, the large dimensionality of biological systems studied makes them prohibitively expensive to work in fully Markovian regimes. In this case, coarse graining can be introduced to capture the key dynamical processes-slow degrees of the system-and reduce the dimension of the problem. Here, we introduce several possible options for such Markovian coarse graining, including previously commonly used choices: the local equilibrium and the Hummer Szabo approaches. We prove that the coarse grained lower dimensional MSM satisfies a variational principle with respect to its slowest relaxation time scale. This provides an excellent framework for optimal coarse graining, as previously demonstrated. Here, we show that such optimal coarse graining to two or three states has a simple physical interpretation in terms of mean first passage times and fluxes between the coarse grained states. The results are verified numerically using both analytic test potentials and data from explicit solvent molecular dynamics simulations of pentalanine. This approach of optimizing and interpreting clustering protocols has broad applicability and can be used in time series analysis of large data.

Solvent Composition Effects on the Structural Properties of the Aβ42 Monomer from the 3D-RISM-KH Molecular Theory of Solvation.

Structural characterization of amyloid (A)β peptides implicated in Alzheimer's disease is a challenging problem due to their intrinsically disordered nature and their high propensity for aggregation. Only limited information is currently available from experiments on conformational properties and aggregation pathways of the peptides in cellular environments. In silico modeling complements experimental information, providing atomistic insight into structure and dynamics of different Aβ species. All-atom explicit solvent molecular dynamics (MD) simulations with a properly selected force field can deliver reliable structural and dynamic information. In the case of intrinsically disordered Aβ peptides, enhanced sampling simulations beyond the nanosecond time scale are required to obtain statistically meaningful results even for simple solvent conditions. To overcome the challenges of conformational sampling in crowded cellular environments, alternative approaches have to be used, including postprocessing of MD data. In this study, we employ the statistical-mechanical, three-dimensional reference interaction site model with the Kovalenko-Hirata closure integral equation molecular theory of solvation to describe solvent composition effects on the conformational equilibrium in a structural ensemble of the Aβ42 (covering residues 1-42) monomer based on a statistical reweighting technique. The methodology enables a computationally efficient prediction on how different factors in the cellular environment, such as solvent composition, nonpolar solvation, and macromolecular crowding, affect the structural properties of the monomer. Similarities have been identified between changes in the structural ensemble caused by nonpolar solvation and crowded environments modeled by ionic solution with large negative ions. In particular, both solvent conditions reduce the random coil content and enhance the helical structure content of the monomer. In contrast to the previous studies, which reported increased α-helical content of peptides in crowded environments, this work attributes these structural features to the difference in solvent exposure of hydrophilic residues of the monomer for different secondary structure elements, rather than to (entropic) excluded volume effects.

Heterogeneous Solvation in Distinctive Protein-Protein Interfaces Revealed by Molecular Dynamics Simulations

Water, despite being a driving force in biochemical processes, has an elusively complex microscopic behavior. While water can increase its local density near amphiphilic protein surfaces, water is also thought to evaporate from hydrophobic surfaces and cavities, an effect known as “dewetting”. The existence and extent of dewetting effects remains elusive due to the difficulty in observing clear “drying” transitions in experiments or simulations. Here, we use explicit solvent molecular dynamics (MD) simulations to study the molecular solvation at the binding interfaces of two distinctive molecular complexes: the highly hydrophilic barnase-barstar and the highly hydrophobic MDM2-p53. Our simulations, in conjunction with simple volumetric analyses, reveal a strikingly different water behavior at the binding interfaces of these two molecular complexes. In both complexes, we observe significant changes in the water local density as the two proteins approach, supporting the existence of a clear dewetting transition in the case of MDM2-p53, with an onset distance of 5.6-7.6Å. Furthermore, the solvation analysis reported herein is a valuable tool to capture and quantify persistent or transient dewetting events in future explicit solvent MD simulations.

Structure and Dynamics of Alzheimer's Associated Amyloid-Beta Peptide

Alzheimer's disease and other neurodegenerative illnesses are among the leading causes of death in the developed world, yet there are currently no treatments available. One of the characteristics of Alzheimer's disease is the formation of neurotoxic fibrils and oligomers through the aggregation of the intrinsically disordered protein (IDP) amyloid-β. While understanding this process has seen widespread experimental interest, it is becoming evident that molecular simulations can provide a unique, in-depth description of these events. A major challenge in this context has been in modelling not only the thermodynamic and structural, but also kinetic properties of the system. The latter requires a discretisation into individual states, which is typically difficult to achieve for IDPs. However, many recent developments in the analysis of dynamical systems as well as the growing accessibility of computational resources have made addressing the problem of generating high quality ensembles possible. We here present an in-depth description of an amyloid-β 42 ensemble acquired from all-atom, explicit solvent molecular dynamics simulations totaling 240 μs using a Markov state model and validate our results using data from NMR measurements. We then correlate our findings with information on the subsequences known to be important for aggregation. Taken together, these results could aid in the understanding of the oligomerisation and fibrillisation processes, and subsequently inform drug developments targeting specific epitopes along the sequence of this aggregation-prone disordered protein.

A Dynamic Hydrophobic Core and Surface Salt Bridges Thermostabilize a Designed Three-Helix Bundle

Thermostable proteins are advantageous in industrial applications, as pharmaceuticals or biosensors, and as templates for directed evolution. As protein-design methodologies improve, bioengineers are able to design proteins to perform a desired function. Although many rationally designed proteins end up being thermostable, how to intentionally design de novo, thermostable proteins is less clear. UVF is a de novo-designed protein based on the backbone structure of the Engrailed homeodomain (EnHD) and is highly thermostable (Tm > 99°C vs. 52°C for EnHD). Although most proteins generally have polar amino acids on their surfaces and hydrophobic amino acids buried in their cores, protein engineers followed this rule exactly when designing UVF. To investigate the contributions of the fully hydrophobic core versus the fully polar surface to UVF’s thermostability, we built two hybrid, chimeric proteins combining the sets of buried and surface residues from UVF and EnHD. Here, we determined a structural, dynamic, and thermodynamic explanation for UVF’s thermostability by performing 4 μs of all-atom, explicit-solvent molecular dynamics simulations at 25 and 100°C, Tanford-Kirkwood solvent accessibility Monte Carlo electrostatic calculations, and a thermodynamic analysis of 40 temperature runs by the weighted-histogram analysis method of heavy-atom, structure-based models of UVF, EnHD, and both chimeric proteins. Our models showed that UVF was highly dynamic because of its fully hydrophobic core, leading to a smaller loss of entropy upon folding. The charged residues on its surface made favorable electrostatic interactions that contributed enthalpically to its thermostability. In the chimeric proteins, both the hydrophobic core and charged surface independently imparted thermostability.

β-Branched Amino Acids Stabilize Specific Conformations of Cyclic Hexapeptides

Cyclic peptides (CPs) are a promising class of molecules for drug development, particularly as inhibitors of protein-protein interactions. Predicting low-energy structures and global structural ensembles of individual CPs is critical for the design of bioactive molecules, but these are challenging to predict and difficult to verify experimentally. In our previous work, we used explicit-solvent molecular dynamics simulations with enhanced sampling methods to predict the global structural ensembles of cyclic hexapeptides containing different permutations of glycine, alanine, and valine. One peptide, cyclo-(VVGGVG) or P7, was predicted to be unusually well structured. In this work, we synthesized P7, along with a less well-structured control peptide, cyclo-(VVGVGG) or P6, and characterized their global structural ensembles in water using NMR spectroscopy. The NMR data revealed a structural ensemble similar to the prediction for P7 and showed that P6 was indeed much less well-structured than P7. We then simulated and experimentally characterized the global structural ensembles of several P7 analogs and discovered that β-branching at one critical position within P7 is important for overall structural stability. The simulations allowed deconvolution of thermodynamic factors that underlie this structural stabilization. Overall, the excellent correlation between simulation and experimental data indicates that our simulation platform will be a promising approach for designing well-structured CPs and also for understanding the complex interactions that control the conformations of constrained peptides and other macrocycles.

Aflatoxin B₁⁻Formamidopyrimidine DNA Adducts: Relationships between Structures, Free Energies, and Melting Temperatures.

Thermal stabilities of DNA duplexes containing Gua (g), α- (a) or β-anomer of formamidopyrimidine-N7-9-hydroxy-aflatoxin B1 (b) differ markedly (Tm: ), but the underlying molecular origin of this experimentally observed phenomenon is yet to be identified and determined. Here, by employing explicit-solvent molecular dynamics simulations coupled with free-energy calculations using a combined linear-interaction-energy/linear-response-approximation approach, we explain the quantitative differences in T in terms of three structural features (bulkiness, order, and compactness) and three energetical contributions (non-polar, electrostatic, and preorganized-electrostatic), and thus advance the current understanding of the relationships between structures, free energies, and thermal stabilities of DNA double helices.

A consensus protocol for the in silico optimisation of antibody fragments.

We present an in silico mutagenetic protocol for improving the binding affinity of single domain antibodies (or nanobodies, VHHs). The method iteratively attempts random mutations in the interacting region of the protein and evaluates the resulting binding affinity towards the target by scoring, with a collection of scoring functions, short explicit solvent molecular dynamics trajectories of the binder-target complexes. The acceptance/rejection of each attempted mutation is carried out by a consensus decision-making algorithm, which considers all individual assessments derived from each scoring function. The method was benchmarked by evolving a single complementary determining region (CDR) of an anti-HER2 VHH hit obtained by direct panning of a phage display library. The optimised VHH mutant showed significantly enhanced experimental affinity with respect to the original VHH it matured from. The protocol can be employed as it is for the optimization of peptides, antibody fragments, and (given enough computational power) larger antibodies.

Temperature-induced structural changes of apo-lactoferrin and their functional implications: a molecular dynamics simulation study

ABSTRACTLactoferrin (Lf) is an iron-binding glycoprotein present in secretory fluids from human and bovine sources. Sequence alignment was employed to identify a region on the C-lobe of Lf capable of binding to bacterial cell surfaces, followed by all-atom explicit solvent molecular dynamics simulations to study the conformational changes of Lf after exposure to three processing temperatures: pasteurisation (72°C), spray drying (90°C) and ultra-high temperature (UHT) (127°C). Below 90°C, the simulations indicate relatively minor changes in overall protein structure. At UHT conditions (127°C), however, marked disruptions to protein structure were found as demonstrated by a substantial decrease in protein dimensions due to collapse in the inter-lobe region. There was also a marked increase in residue fluctuations in several regions of known functional importance, including antibacterial, iron-binding, and putative membrane binding regions, the latter of which is stabilised by a triplet of hydrophobic residues comprised of Leu446, Trp448 and Leu451 at low temperature, but which are disrupted under UHT conditions. A unique network analysis confirmed these results as demonstrated by large clusters of residues with increased dynamical correlation in the N-terminal lobe.

Counterion Effects on Aggregate Structure of 12-Hydroxystearate Salts in Hexane: A Quantum Mechanical and Molecular Dynamics Simulation Study.

Salts of 12-hydroxystearate are important organogelators and grease thickeners, but a structural rationale for their rheological properties remains elusive. We use quantum mechanical calculations and molecular dynamics (MD) simulations to analyze aggregate structures for (1) ( R)-12-hydroxystearic acid (( R)-12HSA), (2) lithium ( R)-12-hydroxystearate (( R)-Li12HS), and (3) sodium ( R)-12-hydroxystearate (( R)-Na12HS). First, quantum mechanical calculations were used to establish the structure and complexation energies of dimers of acetic acid, lithium acetate, and sodium acetate. The expected acetic acid dimer is predicted, and both the lithium acetate and sodium acetate dimer formed a C2 h-symmetric structure. All dimers were sufficiently stable to allow modeling them as pseudocovalent complexes in all-atom, explicit solvent MD. After microsecond-long MD, all systems produced strong ringlike ordered nuclei. The C2 h lithium salt molecules produced aggregates that had the most efficient packing at the head group and a higher frequency of hydroxyl hydrogen bonding compared to the sodium salt. This ordering propensity explains the high melting temperature of ( R)-Li12HS. Also, the higher frequency of hydrogen bonding leads to fewer solvent-exposed hydrogen bond partners. This explains why lithium is a common counterion in high-temperature and water-resistant greases.

Multi-strand β-sheet of Alzheimer Aβ(1–40) folds to β-strip helix: implication for protofilament formation

X-ray fiber diffraction experiments on Alzheimer Aβ(1–40) fibrils formed in an assembly process thought to simulate a portion of the pathophysiological process in Alzheimer's disease, indicated protofilaments with tilted β-strands rather than strands oriented perpendicular to the fibril axis as is usually interpreted from cross-β patterns. The protofilament width and tilt angle determined by these experiments were used to predict a β-strip helix model – a β-helix-like structure in which multiple identical polypeptide molecules assemble in-register to form a helical sheet structure such that the outer strands 1 and m join with a register shift t – with m = 11 and t = 22. Starting from untwisted β-sheets comprising 10, 11, and 12 strands, multiple explicit solvent molecular dynamics (MD) simulations were performed to determine whether the sheets form β-strip helices matching the dimensions of the experimentally measured protofilament. In the simulations, the predicted 11-strand sheets curled up to form a closed β-strip helix-like structure with dimensions matching experimental values, whereas the 10- and 12-strand sheets did not form a closed helical structure. The 12-strand structure did, however, show similarity to a cross-β structure determined by a solid-state NMR experiment. The 11-strand β-strip helix resembles a trans-membrane β-barrel which could explain the ability of small oligomers of Aβ(1–40) to form toxic ion channels. A further consequence of opposite sides of the 11-strand strip coming together at a register shift of 22 is end-to-end joins between neighboring β-strip helices, resulting in a protofilament that keeps growing in both directions.Communicated by Ramaswamy H. Sarma

New Binding Mode of SLURP Protein to a7 Nicotinic Acetylcholine Receptor Revealed by Computer Simulations

SLURP-1 is a member of three-finger toxin-like proteins. Their characteristic feature is a set of three beta strands extruding from hydrophobic core stabilized by disulfide bonds. Each beta-strand carries a flexible loop, which is responsible for recognition. SLURP-1 was recently shown to act as an endogenous growth regulator of keratinocytes and tumor suppressor by reducing cell migration and invasion by antagonizing the pro-malignant effects of nicotine. This effect is achieved through allosteric interaction with alpha7 nicotinic acetylcholine receptors (alpha-7 nAChRs) in an antagonist-like manner. Moreover, this interaction is unaffected by several well-known agents specifically alpha-bungarotoxin. In this work, we carry out the conformational analysis of the SLURP-1 by a microsecond-long full-atom explicit solvent molecular dynamics simulations followed by clustering, to identify representative states. To achieve this timescale we employed a GPU-accelerated version of GROMACS modeling package. To avoid human bias in clustering we used a non-parametric clustering algorithm Affinity Propagation adapted for biomolecules and HPC environments. Then, we applied protein-protein molecular docking of the ten most massive clusters to alpha7-nAChRs in order to test if structural variability can affect binding. Docking simulations revealed the unusual binding mode of one of the minor SLURP-1 conformations.

Correction to Large-Scale Analysis of 48 DNA and 48 RNA Tetranucleotides Studied by 1 μs Explicit-Solvent Molecular Dynamics Simulations

Correction to Large-Scale Analysis of 48 DNA and 48 RNA Tetranucleotides Studied by 1 μs Explicit-Solvent Molecular Dynamics Simulations

Correction to Direct Comparison of Amino Acid and Salt Interactions with Double-Stranded and Single-Stranded DNA from Explicit-Solvent Molecular Dynamics Simulations.

Given the ubiquitous nature of protein–DNA interactions, it is important to understand the interaction thermodynamics of individual amino acid side chains for DNA. One way to assess these preferences is to perform molecular dynamics (MD) simulations. Here we report MD simulations of 20 amino acid side chain analogs interacting simultaneously with both a 70-base-pair double-stranded DNA and with a 70-nucleotide single-stranded DNA. The relative preferences of the amino acid side chains for dsDNA and ssDNA match well with values deduced from crystallographic analyses of protein–DNA complexes. The estimated apparent free energies of interaction for ssDNA, on the other hand, correlate well with previous simulation values reported for interactions with isolated nucleobases, and with experimental values reported for interactions with guanosine. Comparisons of the interactions with dsDNA and ssDNA indicate that, with the exception of the positively charged side chains, all types of amino acid side chain interact...