Protein surface roughness accounts for binding free energy of Plasmepsin II‐ligand complexes

The calculation of absolute binding affinities for protein‐inhibitor complexes remains as one of the main challenges in computational structure‐based ligand design. The present work explored the calculations of surface fractal dimension (as a measure of surface roughness) and the relationship with experimental binding free energies of Plasmepsin II complexes. Plasmepsin II is an attractive target for novel therapeutic compounds to treat malaria. However, the structural flexibility of this enzyme is a drawback when searching for specific inhibitors. Concerning that, we performed separate explicitly solvated molecular dynamics simulations using the available high‐resolution crystal structures of different Plasmepsin II complexes. Molecular dynamics simulations allowed a better approximation to systems dynamics and, therefore, a more reliable estimation of surface roughness. This constitutes a novel approximation in order to obtain more realistic values of fractal dimension, because previous works considered only x‐ray structures. Binding site fractal dimension was calculated considering the ensemble of structures generated at different simulation times. A linear relationship between binding site fractal dimension and experimental binding free energies of the complexes was observed within 20 ns. Previous studies of the subject did not uncover this relationship. Regression model, coined FD model, was built to estimate binding free energies from binding site fractal dimension values. Leave‐one‐out cross‐validation showed that our model reproduced accurately the absolute binding free energies for our training set (R² = 0.76; <|error|> =0.55 kcal/mol; SD_error = 0.19 kcal/mol). The fact that such a simple model may be applied raises some questions that are addressed in the article.     

GABA Binding to an Insect GABA Receptor: A Molecular Dynamics and Mutagenesis Study

RDL receptors are GABA-activated inhibitory Cys-loop receptors found throughout the insect CNS. They are a key target for insecticides. Here, we characterize the GABA binding site in RDL receptors using computational and electrophysiological techniques. A homology model of the extracellular domain of RDL was generated and GABA docked into the binding site. Molecular dynamics simulations predicted critical GABA binding interactions with aromatic residues F206, Y254, and Y109 and hydrophilic residues E204, S176, R111, R166, S176, and T251. These residues were mutated, expressed in Xenopus oocytes, and their functions assessed using electrophysiology. The data support the binding mechanism provided by the simulations, which predict that GABA forms many interactions with binding site residues, the most significant of which are cation-π interactions with F206 and Y254, H-bonds with E204, S205, R111, S176, T251, and ionic interactions with R111 and E204. These findings clarify the roles of a range of residues in binding GABA in the RDL receptor, and also show that molecular dynamics simulations are a useful tool to identify specific interactions in Cys-loop receptors.Abbreviations used: nACh, nicotinic acetylcholine; AChBP, acetylcholine binding protein; GABA, gamma-aminobutyric acid; MD, molecular dynamics; RDL, resistant to dieldrin; RMSD, root mean-square displacement; RMSF, root mean-square fluctuation     

Cryo-Cooling Effect on DHFR Crystal Studied by Replica-Exchange Molecular Dynamics Simulations

Cryo-cooling is routinely performed before x-ray diffraction image collection to reduce the damage to crystals due to ionizing radiation. It has been suggested that although backbone structures are usually very similar between room temperature and cryo-temperature, cryo-cooling may hamper biologically relevant dynamics. In this study, the crystal of Escherichia coli dihydrofolate reductase is studied with replica-exchange molecular dynamics simulation, and the results are compared with the crystal structure determined at cryo-temperature and room temperature with the time-averaged ensemble method. Although temperature dependence of unit cell compaction and root mean-square fluctuation of Cα is found in accord with experiment, it is found that the protein structure at low temperature can be more heterogeneous than the ensemble of structures reported by using the time-averaged ensemble method, encouraging further development of the time-averaged ensemble method and indicating that data should be examined carefully to avoid overinterpretation of one average structure.     

Structural refinement of membrane proteins by restrained molecular dynamics and solvent accessibility data

We present an approach for incorporating solvent accessibility data from electron paramagnetic resonance experiments in the structural refinement of membrane proteins through restrained molecular dynamics simulations. The restraints have been parameterized from oxygen (ΠO₂) and nickel-ethylenediaminediacetic acid (ΠNiEdda) collision frequencies, as indicators of lipid or aqueous exposed spin-label sites. These are enforced through interactions between a pseudoatom representation of the covalently attached Nitroxide spin-label and virtual “solvent” particles corresponding to O₂ and NiEdda in the surrounding environment. Interactions were computed using an empirical potential function, where the parameters have been optimized to account for the different accessibilities of the spin-label pseudoatoms to the surrounding environment. This approach, “pseudoatom-driven solvent accessibility refinement”, was validated by refolding distorted conformations of the Streptomyces lividans potassium channel (KcsA), corresponding to a range of 2-30 Å root mean-square deviations away from the native structure. Molecular dynamics simulations based on up to 58 electron paramagnetic resonance restraints derived from spin-label mutants were able to converge toward the native structure within 1-3 Å root mean-square deviations with minimal computational cost. The use of energy-based ranking and structure similarity clustering as selection criteria helped in the convergence and identification of correctly folded structures from a large number of simulations. This approach can be applied to a variety of integral membrane protein systems, regardless of oligomeric state, and should be particularly useful in calculating conformational changes from a known reference crystal structure.     

Light Activation of Rhodopsin: Insights from Molecular Dynamics Simulations Guided by Solid-State NMR Distance Restraints

Structural restraints provided by solid-state NMR measurements of the metarhodopsin II intermediate are combined with molecular dynamics simulations to help visualize structural changes in the light activation of rhodopsin. Since the timescale for the formation of the metarhodopsin II intermediate (> 1 ms) is beyond that readily accessible by molecular dynamics, we use NMR distance restraints derived from ¹³C dipolar recoupling measurements to guide the simulations. The simulations yield a working model for how photoisomerization of the 11-cis retinylidene chromophore bound within the interior of rhodopsin is coupled to transmembrane helix motion and receptor activation. The mechanism of activation that emerges is that multiple switches on the extracellular (or intradiscal) side of rhodopsin trigger structural changes that converge to disrupt the ionic lock between helices H3 and H6 on the intracellular side of the receptor.     

TASSER_WT: a protein structure prediction algorithm with accurate predicted contact restraints for difficult protein targets

To improve the prediction accuracy in the regime where template alignment quality is poor, an updated version of TASSER_2.0, namely TASSER_WT, was developed. TASSER_WT incorporates more accurate contact restraints from a new method, COMBCON. COMBCON uses confidence-weighted contacts from PROSPECTOR_3.5, the latest version, PROSPECTOR_4, and a new local structural fragment-based threading algorithm, STITCH, implemented in two variants depending on expected fragment prediction accuracy. TASSER_WT is tested on 622 Hard proteins, the most difficult targets (incorrect alignments and/or templates and incorrect side-chain contact restraints) in a comprehensive benchmark of 2591 nonhomologous, single domain proteins ≤200 residues that cover the PDB at 35% pairwise sequence identity. For 454 of 622 Hard targets, COMBCON provides contact restraints with higher accuracy and number of contacts per residue. As contact coverage with confidence weight ≥3 (F^wt≥3_cov) increases, the more improved are TASSER_WT models. When F^wt≥3_cov > 1.0 and > 0.4, the average root mean-square deviation of TASSER_WT (TASSER_2.0) models is 4.11 Å (6.72 Å) and 5.03 Å (6.40 Å), respectively. Regarding a structure prediction as successful when a model has a TM-score to the native structure ≥0.4, when F^wt≥3_cov > 1.0 and > 0.4, the success rate of TASSER_WT (TASSER_2.0) is 98.8% (76.2%) and 93.7% (81.1%), respectively.     

An Unusual Hydrophobic Core Confers Extreme Flexibility to HEAT Repeat Proteins

Alpha-solenoid proteins are suggested to constitute highly flexible macromolecules, whose structural variability and large surface area is instrumental in many important protein-protein binding processes. By equilibrium and nonequilibrium molecular dynamics simulations, we show that importin-β, an archetypical α-solenoid, displays unprecedentedly large and fully reversible elasticity. Our stretching molecular dynamics simulations reveal full elasticity over up to twofold end-to-end extensions compared to its bound state. Despite the absence of any long-range intramolecular contacts, the protein can return to its equilibrium structure to within 3 Å backbone RMSD after the release of mechanical stress. We find that this extreme degree of flexibility is based on an unusually flexible hydrophobic core that differs substantially from that of structurally similar but more rigid globular proteins. In that respect, the core of importin-β resembles molten globules. The elastic behavior is dominated by nonpolar interactions between HEAT repeats, combined with conformational entropic effects. Our results suggest that α-solenoid structures such as importin-β may bridge the molecular gap between completely structured and intrinsically disordered proteins.     

Subdiffraction-limit study of Kaede diffusion and spatial distribution in live Escherichia coli

Photoactivation localization microscopy (PALM) is used to study the spatial distribution and diffusion of single copies of the protein Kaede in the cytoplasm of live Escherichia coli under moderate growth conditions (67 min doubling time). The spatial distribution of Kaede is uniform within the cytoplasm. The cytoplasmic radius of 380 ± 30 nm varies little from cell to cell. Single-particle tracking using 4 ms exposure times reveals negatively curved plots of mean-square displacement versus time. A detailed comparison with Monte Carlo simulations in a spherocylindrical volume shows that the curvature can be quantitatively understood in terms of free diffusion within a confining volume. The mean diffusion coefficient across cells is <D_Kaede> = 7.3 ± 1.1 μm²·s⁻¹, consistent with a homotetrameric form of Kaede. The distribution of squared displacements along the long axis for individual Kaede molecules is consistent with homogeneous diffusion. However, for longer cells, a spatial map of one-step estimates of the diffusion coefficient along x suggests that diffusion is ∼20–40% faster within nucleoids than in the ribosome-rich region lying between nucleoid lobes at the cell mid-plane. Fluorescence recovery after photobleaching yielded <D_FRAP> = 8.3 ± 1.6 μm²·s⁻¹, in agreement with the single-particle tracking results.     

A one-dimensional reaction coordinate for identification of transition states from explicit solvent P(fold)-like calculations

A properly identified transition state ensemble (TSE) in a molecular dynamics (MD) simulation can reveal a tremendous amount about how a protein folds and offer a point of comparison to experimentally derived Phi(F) values, which reflect the degree of structure in these transient states. In one such method of TSE identification, dubbed P(fold), MD simulations of individual protein structures taken from an unfolding trajectory are used to directly assess an input structure's probability of folding before unfolding, and P(fold) is, by definition, 0.5 for the TSE. Other, less computationally intensive methods, such as multidimensional scaling (MDS) of the pairwise root mean-squared deviation (RMSD) matrix of the conformations sampled in a thermal unfolding trajectory, have also been used to identify the TSE. Identification of the TSE is made from the original MD simulation without the need to run further simulations. Here we present a P(fold)-like study and describe methods for identification of the TSE through the derivation of a high fidelity, bounded, one-dimensional reaction coordinate for protein folding. These methods are applied to the engrailed homeodomain. The TSE identified by this approach is essentially identical to the TSE identified previously by MDS of the pairwise RMSD matrix. However, the cost of performing P(fold), or even our reduced P(fold)-like calculations, is at least 36,000 times greater than the MDS method.     

Molecular dynamics simulations of a protein-protein dimer: particle-mesh Ewald electrostatic model yields far superior results to standard cutoff model.

In this article we present two 1000 ps molecular dynamics simulations on the rat micro-glutathione S-transferase dimeric enzyme in complex with the product 1-(S-glutathionyl)-2,4-dinitrobenzene, in a periodic box with explicit solvent molecules, and investigate the effect of long-range electrostatics models on the structure and dynamics of the dimer and its components. One simulation used the standard cutoff method (10A), whilst the other used the particle-mesh Ewald (PME) method. We monitored the root mean-square atomic deviation (RMSD) from the initial crystal structure to examine the convergence of both simulations, as well as several other structural parameters such as the distance between active sites, rigid body rotation between domains in subunits, radius of gyration, B-factors, number of hydrogen bonds and salt bridges and solvent-accessible surface area. For example, with the PME method, the dimer structure remains much closer to the initial crystallographic structure with an average RMSD of 1.3A +/- 0.1A and 1.0A +/- 0.1A for all heavy and backbone atoms, respectively, in the last 200 ps; the respective values for the cutoff simulation are 4.7A +/- 0.3A and 4.2A +/- 0.3A. The large deviations observed in the cutoff simulation severely affected the stability of the enzyme dimer and its complex with the bound product. This finding is contrary to that found in a similar study of the monomeric protein ubiquitin [Fox, T. & Kollman, P. A. Proteins Struct. Func. Genet. 25, 315-334 (1996)]. Unlike the earlier published work, the present study provides evidence that the standard cutoff method is not generally valid for the study of protein complexes, or their subunits.     

Atomic-Level Mechanisms for Phospholamban Regulation of the Calcium Pump

We performed protein pKa calculations and molecular dynamics (MD) simulations of the calcium pump (sarcoplasmic reticulum Ca²⁺-ATPase (SERCA)) in complex with phospholamban (PLB). X-ray crystallography studies have suggested that PLB locks SERCA in a low-Ca²⁺-affinity E2 state that is incompatible with metal-ion binding, thereby blocking the conversion toward a high-Ca²⁺-affinity E1 state. Estimation of pKa values of the acidic residues in the transport sites indicates that at normal intracellular pH (7.1–7.2), PLB-bound SERCA populates an E1 state that is deprotonated at residues E309 and D800 yet protonated at residue E771. We performed three independent microsecond-long MD simulations to evaluate the structural dynamics of SERCA-PLB in a solution containing 100 mM K⁺ and 3 mM Mg²⁺. Principal component analysis showed that PLB-bound SERCA lies exclusively along the structural ensemble of the E1 state. We found that the transport sites of PLB-bound SERCA are completely exposed to the cytosol and that K⁺ ions bind transiently (≤5 ns) and nonspecifically (nine different positions) to the two transport sites, with a total occupancy time of K⁺ in the transport sites of 80%. We propose that PLB binding to SERCA populates a novel (to our knowledge) E1 intermediate, E1⋅H⁺₇₇₁. This intermediate serves as a kinetic trap that controls headpiece dynamics and depresses the structural transitions necessary for Ca²⁺-dependent activation of SERCA. We conclude that PLB-mediated regulation of SERCA activity in the heart results from biochemical and structural transitions that occur primarily in the E1 state of the pump.     

Characterizing a partially ordered miniprotein through folding molecular dynamics simulations: Comparison with the experimental data

The villin headpiece helical subdomain (HP36) is one of the best known model systems for computational studies of fast‐folding all‐α miniproteins. HP21 is a peptide fragment—derived from HP36—comprising only the first and second helices of the full domain. Experimental studies showed that although HP21 is mostly unfolded in solution, it does maintain some persistent native‐like structure as indicated by the analysis of NMR‐derived chemical shifts. Here we compare the experimental data for HP21 with the results obtained from a 15‐μs long folding molecular dynamics simulation performed in explicit water and with full electrostatics. We find that the simulation is in good agreement with the experiment and faithfully reproduces the major experimental findings, namely that (a) HP21 is disordered in solution with <10% of the trajectory corresponding to transiently stable structures, (b) the most highly populated conformer is a native‐like structure with an RMSD from the corresponding portion of the HP36 crystal structure of <1 Å, (c) the simulation‐derived chemical shifts—over the whole length of the trajectory—are in reasonable agreement with the experiment giving reduced χ² values of 1.6, 1.4, and 0.8 for the Δδ¹³C^α, Δδ¹³CO, and Δδ¹³C^β secondary shifts, respectively (becoming 0.8, 0.7, and 0.3 when only the major peptide conformer is considered), and finally, (d) the secondary structure propensity scores are in very good agreement with the experiment and clearly indicate the higher stability of the first helix. We conclude that folding molecular dynamics simulations can be a useful tool for the structural characterization of even marginally stable peptides.     

Coiled-Coil Response to Mechanical Force: Global Stability and Local Cracking

Coiled coils are important structural motifs formed by two or more amphipathic α-helices that twist into a supercoil. These motifs are found in a wide range of proteins, including motor proteins and structural proteins, that are known to transmit mechanical loads. We analyze atomically detailed simulations of coiled-coil cracking under load with Milestoning. Milestoning is an approach that captures the main features of the process in a network, quantifying kinetics and thermodynamics. A 112-residue segment of the β-myosin S2 domain was subjected to constant-magnitude (0–200 pN) and constant-direction tensile forces in molecular dynamics simulations. Twenty 20 ns straightforward simulations at several load levels revealed that initial single-residue cracking events (Ψ > 90°) at loads <100 pN were accompanied by rapid refolding without either intra- or interhelix unfolding propagation. Only initial unfolding events at the highest load (200 pN) regularly propagated along and between helices. Analysis of hydrophobic interactions and of interhelix hydrogen bonds did not show significant variation as a function of load. Unfolding events were overwhelmingly located in the vicinity of E929, a charged residue in a hydrophobic position of the heptad repeat. Milestoning network analysis of E929 cracking determined that the mean first-passage time ranges from 20 ns (200 pN) to 80 ns (50 pN), which is ∼20 times the mean first-passage time of an isolated helix with the same sequence.     

Importance of the CMAP correction to the CHARMM22 protein force field: dynamics of hen lysozyme

The recently developed CMAP correction to the CHARMM22 force field (C22) is evaluated from 25 ns molecular dynamics simulations on hen lysozyme. Substantial deviations from experimental backbone root mean-square fluctuations and N-H NMR order parameters obtained in the C22 trajectories (especially in the loops) are eliminated by the CMAP correction. Thus, the C22/CMAP force field yields improved dynamical and structural properties of proteins in molecular dynamics simulations.     

Lipid bilayer structure determined by the simultaneous analysis of neutron and X-ray scattering data

Quantitative structures were obtained for the fully hydrated fluid phases of dioleoylphosphatidylcholine (DOPC) and dipalmitoylphosphatidylcholine (DPPC) bilayers by simultaneously analyzing x-ray and neutron scattering data. The neutron data for DOPC included two solvent contrasts, 50% and 100% D₂O. For DPPC, additional contrast data were obtained with deuterated analogs DPPC_d62, DPPC_d13, and DPPC_d9. For the analysis, we developed a model that is based on volume probability distributions and their spatial conservation. The model's design was guided and tested by a DOPC molecular dynamics simulation. The model consistently captures the salient features found in both electron and neutron scattering density profiles. A key result of the analysis is the molecular surface area, A. For DPPC at 50°C A = 63.0 Å², whereas for DOPC at 30°C A = 67.4 Å², with estimated uncertainties of 1 Å². Although A for DPPC agrees with a recently reported value obtained solely from the analysis of x-ray scattering data, A for DOPC is almost 10% smaller. This improved method for determining lipid areas helps to reconcile long-standing differences in the values of lipid areas obtained from stand-alone x-ray and neutron scattering experiments and poses new challenges for molecular dynamics simulations.     

Robust probabilistic superposition and comparison of protein structures

Background  

Protein structure comparison is a central issue in structural bioinformatics. The standard dissimilarity measure for protein structures is the root mean square deviation (RMSD) of representative atom positions such as α-carbons. To evaluate the RMSD the structures under comparison must be superimposed optimally so as to minimize the RMSD. How to evaluate optimal fits becomes a matter of debate, if the structures contain regions which differ largely - a situation encountered in NMR ensembles and proteins undergoing large-scale conformational transitions.     

Molecular Dynamics Study of High Temperature Phase-Separation in a H2O/N2 Mixture with Exp-6 Interactions

Abstract

Equimolar H₂O/N₂ fluid mixture was studied by molecular dynamics simulations for NVT ensemble. Calculations were performed with the modified Buckingham (exp-6) potentials at T = 2000 K. Particular attention was given to the phase separation at very high pressures relevant to a detonation environment. Calculations of pair correlation functions and local mole fractions clearly indicated the occurrence of the fluid separation into N₂-rich and H₂O-rich phase. The density at the phase boundary between homogeneous and inhomogeneous phase-separated state was determined to be p = 1.35 g/cm³ on the basis of the static cross correlation factor which is defined by the sum of the local mole fractions. The ratio of the self-diffusion coefficients of N₂ and H₂O at p < 1.35 g/cm³ was found to be approximately equal to the value predicted by the kinetic theory of the ideal gas, whereas the ratio was close to unity at the phase-separated state (p > 1.35 g/cm³). In addition, two distinctive behaviors of the system could be observed for the relaxation from the initial uniform mixture to the phase-separated fluid: at lower densities (1.35 < p < 2.0 g/cm³) the fluid mixture began to relax into the phase-separated system without obvious incubation time, while clear incubation period was associated for the separation at higher densities. During this incubation period, discontinuous jumps in the mean square displacements were found.