Refinement of comparative models of protein structure by using multicanonical molecular dynamics simulations

Comparative modelling is a powerful method that easily predicts a considerably accurate structure of a protein by using a template structure having a similar amino-acid sequence to the target protein. However, in the region where the amino-acid sequence is different between the target and the template, the predicted structure remains unreliable. In such a case, the model has to be refined. In the present study, we explored the possibility of a molecular dynamics-based method, using the human SAP Src Homology 2 (SH2) domain as the modelling target. The multicanonical method was used to alleviate the multiple-minima problem and the generalised Born/surface area model was used to reduce the computational cost. In addition, position restraints were imposed on the atoms in the reliable regions to avoid unnecessary conformational sampling. We analyzed the conformational distribution of the ligand-recognition loop of the domain and found that the most populated conformational clusters in the ensemble of the model agreed well with one of the two major clusters in the ensemble of the reference simulation starting from the crystal structure. This demonstrates that the current refinement method can significantly improve the accuracy of an unreliable region in a comparative model.     

Implicit solvent models for flexible protein-protein docking by molecular dynamics simulation

The suitability of three implicit solvent models for flexible protein-protein docking by procedures using molecular dynamics simulation is investigated. The three models are (i) the generalized Born (GB) model implemented in the program AMBER6.0; (ii) a distance-dependent dielectric (DDD) model; and (iii) a surface area-dependent model that we have parameterized and call the NPSA model. This is a distance-dependent dielectric model modified by neutralizing the ionizable side-chains and adding a surface area-dependent solvation term. These solvent models were first tested in molecular dynamics simulations at 300 K of the native structures of barnase, barstar, segment B1 of protein G, and three WW domains. These protein structures display a range of secondary structure contents and stabilities. Then, to investigate the performance of the implicit solvent models in protein docking, molecular dynamics simulations of barnase/barstar complexation, as well as PIN1 WW domain/peptide complexation, were conducted, starting from separated unbound structures. The simulations show that the NPSA model has significant advantages over the DDD and GB models in maintaining the native structures of the proteins and providing more accurate docked complexes.     

Influence of Paleolithic range contraction,admixture and long‐distance dispersal on genetic gradients of modern humans in Asia

Cavalli‐Sforza and coauthors originally explored the genetic variation of modern humans throughout the world and observed an overall east‐west genetic gradient in Asia. However, the specific environmental and population genetics processes causing this gradient were not formally investigated and promoted discussion in recent studies. Here we studied the influence of diverse environmental and population genetics processes on Asian genetic gradients and identified which could have produced the observed gradient. To do so, we performed extensive spatially‐explicit computer simulations of genetic data under the following scenarios: (a) variable levels of admixture between Paleolithic and Neolithic populations, (b) migration through long‐distance dispersal (LDD), (c) Paleolithic range contraction induced by the last glacial maximum (LGM), and (d) Neolithic range expansions from one or two geographic origins (the Fertile Crescent and the Yangzi and Yellow River Basins). Next, we estimated genetic gradients from the simulated data and we found that they were sensible to the analysed processes, especially to the range contraction induced by LGM and to the number of Neolithic expansions. Some scenarios were compatible with the observed east‐west genetic gradient, such as the Paleolithic expansion with a range contraction induced by the LGM or two Neolithic range expansions from both the east and the west. In general, LDD increased the variance of genetic gradients among simulations. We interpreted the obtained gradients as a consequence of both allele surfing caused by range expansions and isolation by distance along the vast east‐west geographic axis of this continent.     

Binding hotspots on K‐ras: Consensus ligand binding sites and other reactive regions from probe‐based molecular dynamics analysis

We have used probe‐based molecular dynamics (pMD) simulations to search for interaction hotspots on the surface of the therapeutically highly relevant oncogenic K‐Ras G12D. Combining the probe‐based query with an ensemble‐based pocket identification scheme and an analysis of existing Ras‐ligand complexes, we show that (i) pMD is a robust and cost‐effective strategy for binding site identification, (ii) all four of the previously reported ligand binding sites are suitable for structure‐based ligand design, and (iii) in some cases probe binding and expanded sampling of configurational space enable pocket expansion and increase the likelihood of site identification. Furthermore, by comparing the distribution of hotspots in nonpocket‐like regions with known protein‐ and membrane‐interacting interfaces, we propose that pMD has the potential to predict surface patches responsible for protein‐biomolecule interactions. These observations have important implications for future drug design efforts and will facilitate the search for potential interfaces responsible for the proposed transient oligomerization or interaction of Ras with other biomolecules in the cellular milieu. Proteins 2015; 83:898–909. © 2015 Wiley Periodicals, Inc.     

Sex‐specific responses to territorial intrusions in a communication network: Evidence from radio‐tagged great tits

Signals play a key role in the ecology and evolution of animal populations, influencing processes such as sexual selection and conflict resolution. In many species, sexually selected signals have a dual function: attracting mates and repelling rivals. Yet, to what extent males and females under natural conditions differentially respond to such signals remains poorly understood, due to a lack of field studies that simultaneously track both sexes. Using a novel spatial tracking system, we tested whether or not the spatial behavior of male and female great tits (Parus major) changes in relation to the vocal response of a territorial male neighbor to an intruder. We tracked the spatial behavior of male and female great tits (N = 44), 1 hr before and 1 hr after simulating territory intrusions, employing automatized Encounternet radio‐tracking technology. We recorded the spatial and vocal response of the challenged males and quantified attraction and repulsion of neighboring males and females to the intrusion site. We additionally quantified the direct proximity network of the challenged male. The strength of a male's vocal response to an intruder induced sex‐dependent movements in the neighborhood, via female attraction and male repulsion. Stronger vocal responders were older and in better body condition. The proximity networks of the male vocal responders, including the number of sex‐dependent connections and average time spent with connections, however, did not change directly following the intrusion. The effects on neighbor movements suggest that the strength of a male's vocal response can provide relevant social information to both the males and the females in the neighborhood, resulting in both sexes adjusting their spatial behavior in contrasting ways, while the social proximity network remained stable. This study underlines the importance of “silent” eavesdroppers within communication networks for studying the dual functioning and evolution of sexually selected signals.     

Molecular dynamics simulation on the allosteric analysis of the c‐di‐GMP class I riboswitch induced by ligand binding

Riboswitches are RNA molecules that regulate gene expression using conformation change, affected by binding of small molecule ligands. Although a number of ligand‐bound aptamer complex structures have been solved, it is important to know ligand‐free conformations of the aptamers in order to understand the mechanism of specific binding by ligands. In this paper, we use dynamics simulations on a series of models to characterize the ligand‐free and ligand‐bound aptamer domain of the c‐di‐GMP class I (GEMM‐I) riboswitch. The results revealed that the ligand‐free aptamer has a stable state with a folded P2 and P3 helix, an unfolded P1 helix and open binding pocket. The first Mg ions binding to the aptamer is structurally favorable for the successive c‐di‐GMP binding. The P1 helix forms when c‐di‐GMP is successive bound. Three key junctions J1/2, J2/3 and J1/3 in the GEMM‐I riboswitch contributing to the formation of P1 helix have been found. The binding of the c‐di‐GMP ligand to the GEMM‐I riboswitch induces the riboswitch's regulation through the direct allosteric communication network in GEMM‐I riboswitch from the c‐di‐GMP binding sites in the J1/2 and J1/3 junctions to the P1 helix, the indirect ones from those in the J2/3 and P2 communicating to P1 helix via the J1/2 and J1/3 media.     

Effect O6‐guanine alkylation on DNA flexibility studied by comparative molecular dynamics simulations

Alkylation of guanine at the O6 atom is a highly mutagenic DNA lesion because it alters the coding specificity of the base causing G:C to A:T transversion mutations. Specific DNA repair enzymes, e.g. O⁶‐alkylguanin‐DNA‐Transferases (AGT), recognize and repair such damage after looping out the damaged base to transfer it into the enzyme active site. The exact mechanism how the repair enzyme identifies a damaged site within a large surplus of undamaged DNA is not fully understood. The O⁶‐alkylation of guanine may change the deformability of DNA which may facilitate the initial binding of a repair enzyme at the damaged site. In order to characterize the effect of O⁶‐methyl‐guanine (O⁶‐MeG) containing base pairs on the DNA deformability extensive comparative molecular dynamics (MD) simulations on duplex DNA with central G:C, O⁶‐MeG:C or O⁶‐MeG:T base pairs were performed. The simulations indicate significant differences in the helical deformability due to the presence of O⁶‐MeG compared to regular undamaged DNA. This includes enhanced base pair opening, shear and stagger motions and alterations in the backbone fine structure caused in part by transient rupture of the base pairing at the damaged site and transient insertion of water molecules. It is likely that the increased opening motions of O⁶‐MeG:C or O⁶‐MeG:T base pairs play a decisive role for the induced fit recognition or for the looping out of the damaged base by repair enzymes. © 2014 Wiley Periodicals, Inc. Biopolymers 103: 23–32, 2015.     

Thermal stability of single‐domain antibodies estimated by molecular dynamics simulations

Single‐domain antibodies (sdAbs) function like regular antibodies, however, consist of only one domain. Because of their low molecular weight, sdAbs have advantages with respect to production and delivery to their targets and for applications such as antibody drugs and biosensors. Thus, sdAbs with high thermal stability are required. In this work, we chose seven sdAbs, which have a wide range of melting temperature (T_m) values and known structures. We applied molecular dynamics (MD) simulations to estimate their relative stability and compared them with the experimental data. High‐temperature MD simulations at 400 K and 500 K were executed with simulations at 300 K as a control. The fraction of native atomic contacts, Q, measured for the 400 K simulations showed a fairly good correlation with the T_m values. Interestingly, when the residues were classified by their hydrophobicity and size, the Q values of hydrophilic residues exhibited an even better correlation, suggesting that stabilization is correlated with favorable interactions of hydrophilic residues. Measuring the Q value on a per‐residue level enabled us to identify residues that contribute significantly to the instability and thus demonstrating how our analysis can be used in a mutant case study.     

Enantioseparation of Four Organophosphonate Derivatives on N‐(3,5‐Dinitrobenzoyl)‐ l‐leucine–n‐Propylamide Stationary Phase by Molecular Modeling

Four groups of organophosphonate derivatives enantiomers were separated on N‐(3,5‐dinitrobenzoyl)‐S‐leucine chiral stationary phase. The three‐dimensional structures of the complexes between the single enantiotopic chiral compounds and chiral stationary phase have been studied using molecular model and molecular dynamics simulation. Detailed results regarding the conformation, auto‐docking, and thermodynamic estimation are presented. The elution order of the enantiomer could be determined from the energy. The predicted chiral discrimination was obtained by computational results. Chirality 25:101–106, 2013. © 2012 Wiley Periodicals, Inc.     

One short cysteine‐rich sequence pattern – two different disulfide‐bonded structures – a molecular dynamics simulation study

The nematocyst walls of Hydra are formed by proteins containing small cysteine‐rich domains (CRDs) of ~25 amino acids. The first CRD of nematocyst outer all antigen (NW1) and the C‐terminal CRD of minicollagen‐1 (Mcol1C) contain six cysteines at identical sequence positions, however adopt different disulfide bonded structures. NW1 shows the disulfide connectivities C2‐C14/C6‐C19/C10‐C18 and Mcol1C C2‐C18/C6‐C14/C10‐C19. To analyze if both show structural preferences in the open, non‐disulfide bonded form, which explain the formation of either disulfide connectivity pattern, molecular dynamics (MD) simulations at different temperatures were performed. NW1 maintained in the 100‐ns MD simulations at 283 K a rather compact fold that is stabilized by specific hydrogen bonds. The Mcol1C structure fluctuated overall more, however stayed most of the time also rather compact. The analysis of the backbone Φ/ψ angles indicated different turn propensities for NW1 and Mcol1C, which mostly can be explained based on published data about the influence of different amino acid side chains on the local backbone conformation. Whereas a folded precursor mechanism may be considered for NW1, Mcol1C may fold according to the quasi‐stochastic folding model involving disulfide bond reshuffling and conformational changes, locking the native disulfide conformations. The study further demonstrates the power of MD simulations to detect local structural preferences in rather dynamic systems such as the open, non‐disulfide bonded forms of NW1 and Mcol1C, which complement published information from NMR backbone residual dipolar couplings. Because the backbone structural preferences encoded by the amino acid sequence embedding the cysteines influence which disulfide connectivities are formed, the data are generally interesting for a better understanding of oxidative folding and the design of disulfide stabilized therapeutics. Copyright © 2015 European Peptide Society and John Wiley & Sons, Ltd.     

Folding simulations of Trp‐cage mini protein in explicit solvent using biasing potential replica‐exchange molecular dynamics simulations

Replica exchange molecular dynamics (RexMD) simulations are frequently used for studying structure formation and dynamics of peptides and proteins. A significant drawback of standard temperature RexMD is, however, the rapid increase of the replica number with increasing system size to cover a desired temperature range. A recently developed Hamiltonian RexMD method has been used to study folding of the Trp‐cage protein. It employs a biasing potential that lowers the backbone dihedral barriers and promotes peptide backbone transitions along the replica coordinate. In two independent applications of the biasing potential RexMD method including explicit solvent and starting from a completely unfolded structure the formation of near‐native conformations was observed after 30–40 ns simulation time. The conformation representing the most populated cluster at the final simulation stage had a backbone root mean square deviation of ～1.3 Å from the experimental structure. This was achieved with a very modest number of five replicas making it well suited for peptide and protein folding and refinement studies including explicit solvent. In contrast, during five independent continuous 70 ns molecular dynamics simulations formation of collapsed states but no near native structure formation was observed. The simulations predict a largely collapsed state with a significant helical propensity for the helical domain of the Trp‐cage protein already in the unfolded state. Hydrogen bonded bridging water molecules were identified that could play an active role by stabilizing the arrangement of the helical domain with respect to the rest of the chain already in intermediate states of the protein. Proteins 2009. © 2008 Wiley‐Liss, Inc.     

Identification of potent inhibitors against snake venom metalloproteinase (SVMP) using molecular docking and molecular dynamics studies

Snake venom metalloproteinase (SVMP) (Echis coloratus (Carpet viper) is a multifunctional enzyme that is involved in producing several symptoms that follow a snakebite, such as severe local hemorrhage, nervous system effects and tissue necrosis. Because the three-dimensional (3D) structure of SVMP is not known, models were constructed, and the best model was selected based on its stereo-chemical quality. The stability of the modeled protein was analyzed through molecular dynamics (MD) simulation studies. Structure-based virtual screening was performed, and 15 potential molecules with the highest binding energies were selected. Further analysis was carried out with induced fit docking, Prime/MM–GBSA (ΔG_Bind calculations), quantum-polarized ligand docking, and density functional theory calculations. Further, the stability of the lead molecules in the SVMP-active site was examined using MD simulation. The results showed that the selected lead molecules were highly stable in the active site of SVMP. Hence, these molecules could potentially be selective inhibitors of SVMP. These lead molecules can be experimentally validated, and their backbone structural scaffold could serve as building blocks in designing drug-like molecules for snake antivenom.     

Strong “bottom‐up” influences on small mammal populations: State‐space model analyses from long‐term studies

“Bottom‐up” influences, that is, masting, plus population density and climate, commonly influence woodland rodent demography. However, “top‐down” influences (predation) also intervene. Here, we assess the impacts of masting, climate, and density on rodent populations placed in the context of what is known about “top‐down” influences. To explain between‐year variations in bank vole Myodes glareolus and wood mouse Apodemus sylvaticus population demography, we applied a state‐space model to 33 years of catch‐mark‐release live‐trapping, winter temperature, and precise mast‐collection data. Experimental mast additions aided interpretation. Rodent numbers in European ash Fraxinus excelsior woodland were estimated (May/June, November/December). December–March mean minimum daily temperature represented winter severity. Total marked adult mice/voles (and juveniles in May/June) provided density indices validated against a model‐generated population estimate; this allowed estimation of the structure of a time‐series model and the demographic impacts of the climatic/biological variables. During two winters of insignificant fruit‐fall, 6.79 g/m² sterilized ash seed (as fruit) was distributed over an equivalent woodland similarly live‐trapped. September–March fruit‐fall strongly increased bank vole spring reproductive rate and winter and summer population growth rates; colder winters weakly reduced winter population growth. September–March fruit‐fall and warmer winters marginally increased wood mouse spring reproductive rate and September–December fruit‐fall weakly elevated summer population growth. Density dependence significantly reduced both species' population growth. Fruit‐fall impacts on demography still appeared after a year. Experimental ash fruit addition confirmed its positive influence on bank vole winter population growth with probable moderation by colder temperatures. The models show the strong impact of masting as a “bottom‐up” influence on rodent demography, emphasizing independent masting and weather influences; delayed effects of masting; and the importance of density dependence and its interaction with masting. We conclude that these rodents show strong “bottom‐up” and density‐dependent influences on demography moderated by winter temperature. “Top‐down” influences appear weak and need further investigation.     

Effect of DNA on the conformational dynamics of the endonucleases I‐DmoI as provided by molecular dynamics simulations

The conformational behavior of the wild‐type endonucleases I‐DmoI and two of its mutants has been studied in the presence and in the absence of DNA target sequences by means of extended molecular dynamics simulations. Our results show that in the absence of DNA, the three protein forms explore a similar essential conformational space, whereas when bound to the same DNA target sequence of 25 base pairs, they diversify and restrain the subspace explored. In addition, the differences in the essential subspaces explored by the residues near the catalytic site for both the bound and unbound forms are discussed in background of the experimental protein activity.     

Insights into the anthrax lethal factor–substrate interaction and selectivity using docking and molecular dynamics simulations

The anthrax toxin of the bacterium Bacillus anthracis consists of three distinct proteins, one of which is the anthrax lethal factor (LF). LF is a gluzincin Zn‐dependent, highly specific metalloprotease with a molecular mass of ～90 kDa that cleaves most isoforms of the family of mitogen‐activated protein kinase kinases (MEKs/MKKs) close to their amino termini, resulting in the inhibition of one or more signaling pathways. Previous studies on the crystal structures of uncomplexed LF and LF complexed with the substrate MEK2 or a MKK‐based synthetic peptide provided structure‐activity correlations and the basis for the rational design of efficient inhibitors. However, in the crystallographic structures, the substrate peptide was not properly oriented in the active site because of the absence of the catalytic zinc atom. In the current study, docking and molecular dynamics calculations were employed to examine the LF‐MEK/MKK interaction along the catalytic channel up to a distance of 20 Å from the zinc atom. This residue‐specific view of the enzyme‐substrate interaction provides valuable information about: (i) the substrate selectivity of LF and its inactivation of MEKs/MKKs (an issue highly important not only to anthrax infection but also to the pathogenesis of cancer), and (ii) the discovery of new, previously unexploited, hot‐spots of the LF catalytic channel that are important in the enzyme/substrate binding and interaction.     

Global fold of human cannabinoid type 2 receptor probed by solid‐state 13C‐, 15N‐MAS NMR and molecular dynamics simulations

The global fold of human cannabinoid type 2 (CB₂) receptor in the agonist‐bound active state in lipid bilayers was investigated by solid‐state ¹³C‐ and ¹⁵N magic‐angle spinning (MAS) NMR, in combination with chemical‐shift prediction from a structural model of the receptor obtained by microsecond‐long molecular dynamics (MD) simulations. Uniformly ¹³C‐ and ¹⁵N‐labeled CB₂ receptor was expressed in milligram quantities by bacterial fermentation, purified, and functionally reconstituted into liposomes. ¹³C MAS NMR spectra were recorded without sensitivity enhancement for direct comparison of C_α, C_β, and C?O bands of superimposed resonances with predictions from protein structures generated by MD. The experimental NMR spectra matched the calculated spectra reasonably well indicating agreement of the global fold of the protein between experiment and simulations. In particular, the ¹³C chemical shift distribution of C_α resonances was shown to be very sensitive to both the primary amino acid sequence and the secondary structure of CB₂. Thus the shape of the C_α band can be used as an indicator of CB₂ global fold. The prediction from MD simulations indicated that upon receptor activation a rather limited number of amino acid residues, mainly located in the extracellular Loop 2 and the second half of intracellular Loop 3, change their chemical shifts significantly (≥1.5 ppm for carbons and ≥5.0 ppm for nitrogens). Simulated two‐dimensional ¹³C_α(i)? ¹³C?O(i) and ¹³C?O(i)? ¹⁵NH(i + 1) dipolar‐interaction correlation spectra provide guidance for selective amino acid labeling and signal assignment schemes to study the molecular mechanism of activation of CB₂ by solid‐state MAS NMR. Proteins 2014; 82:452–465. © 2013 Wiley Periodicals, Inc.     

Theoretical investigation on the glycan‐binding specificity of Agrocybe cylindracea galectin using molecular modeling and molecular dynamics simulation studies

Galectins are β‐galactoside binding proteins which have the ability to serve as potent antitumor, cancer biomarker, and induce tumor cell apoptosis. Agrocybe cylindracea galectin (ACG) is a fungal galectin which specifically recognizes α(2,3)‐linked sialyllactose at the cell surface that plays extensive roles in the biological recognition processes. To investigate the change in glycan‐binding specificity upon mutations, single point and double point site‐directed in silico mutations are performed at the binding pocket of ACG. Molecular dynamics (MD) simulation studies are carried out for the wild‐type (ACG) and single point (ACG1) and double point (ACG2) mutated ACGs to investigate the dynamics of substituted mutants and their interactions with the receptor sialyllactose. Plausible binding modes are proposed for galectin–sialylglycan complexes based on the analysis of hydrogen bonding interactions, total pair‐wise interaction energy between the interacting binding site residues and sialyllactose and binding free energy of the complexes using molecular mechanics–Poisson–Boltzmann surface area. Our result shows that high contribution to the binding in different modes is due to the direct and water‐mediated hydrogen bonds. The binding specificity of double point mutant Y59R/N140Q of ACG2 is found to be high, and it has 26 direct and water‐mediated hydrogen bonds with a relatively low‐binding free energy of −47.52 ± 5.2 kcal/mol. We also observe that the substituted mutant Arg59 is crucial for glycan‐binding and for the preference of α(2,3)‐linked sialyllactose at the binding pocket of ACG2 galectin. When compared with the wild‐type and single point mutant, the double point mutant exhibits enhanced affinity towards α(2,3)‐linked sialyllactose, which can be effectively used as a model for biological cell marker in cancer therapeutics. Copyright © 2015 John Wiley & Sons, Ltd.