Incorporating replacement free energy of binding‐site waters in molecular docking

Binding‐site water molecules play a crucial role in protein‐ligand recognition, either being displaced upon ligand binding or forming water bridges to stabilize the complex. However, rigorously treating explicit binding‐site waters is challenging in molecular docking, which requires to fully sample ensembles of waters and to consider the free energy cost of replacing waters. Here, we describe a method to incorporate structural and energetic properties of binding‐site waters into molecular docking. We first developed a solvent property analysis (SPA) program to compute the replacement free energies of binding‐site water molecules by post‐processing molecular dynamics trajectories obtained from ligand‐free protein structure simulation in explicit water. Next, we implemented a distance‐dependent scoring term into DOCK scoring function to take account of the water replacement free energy cost upon ligand binding. We assessed this approach in protein targets containing important binding‐site waters, and we demonstrated that our approach is reliable in reproducing the crystal binding geometries of protein‐ligand‐water complexes, as well as moderately improving the ligand docking enrichment performance. In addition, SPA program (free available to academic users upon request) may be applied in identifying hot‐spot binding‐site residues and structure‐based lead optimization. Proteins 2014; 82:1765–1776. © 2014 Wiley Periodicals, Inc.     

Conformational fluctuations of the AXH monomer of Ataxin‐1

In this paper, we report the results of molecular dynamics simulations of AXH monomer of Ataxin‐1. The AXH domain plays a crucial role in Ataxin‐1 aggregation, which accompanies the initiation and progression of Spinocerebellar ataxia type 1. Our simulations involving both classical and replica exchange molecular dynamics, followed by principal component analysis of the trajectories obtained, reveal substantial conformational fluctuations of the protein structure, especially in the N‐terminal region. We show that these fluctuations can be generated by thermal noise since the free energy barriers between conformations are small enough for thermally stimulated transitions. In agreement with the previous experimental findings, our results can be considered as a basis for a future design of ataxin aggregation inhibitors that will require several key conformations identified in the present study as molecular targets for ligand binding. Proteins 2016; 84:52–59. © 2015 Wiley Periodicals, Inc.     

Mechanisms for allosteric activation of protease DegS by ligand binding and oligomerization as revealed from molecular dynamics simulations

A local perturbation of a protein may lead to functional changes at some distal site, a phenomenon denoted as allostery. Here, we study the allosteric control of a protease using molecular dynamics simulations. The system considered is the bacterial protein DegS which includes a protease domain activated on ligand binding to an adjacent PDZ domain. Starting from crystallographic structures of DegS homo‐trimers, we perform simulations of the ligand‐free and ‐bound state of DegS at equilibrium. Considering a single protomer only, the trimeric state was mimicked by applying restraints on the residues in contact with other protomers in the DegS trimer. In addition, the bound state was also simulated without any restraints to mimic the monomer. Our results suggest that not only ligand release but also disassembly of a DegS trimer inhibits proteolytic activity. Considering various observables for structural changes, we infer allosteric pathways from the interface with other protomers to the active site. Moreover, we study how ligand release leads to (i) catalytically relevant changes involving residues 199–201 and (ii) a transition from a stretched to a bent conformation for residues 217–219 (which prohibits proper substrate binding). Finally, based on ligand‐induced C_α shifts we identify residues in contact with other protomers in the DegS trimer that likely transduce the perturbation from ligand release from a given protomer to adjacent protomers. These residues likely play a key role in the experimentally known effect of ligand release from a protomer on the proteolytic activity of the other protomers. Proteins 2016; 84:1690–1705. © 2016 Wiley Periodicals, Inc.     

Structure‐based engineering of streptavidin monomer with a reduced biotin dissociation rate

We recently reported the engineering of monomeric streptavidin, mSA, corresponding to one subunit of wild type (wt) streptavidin tetramer. The monomer was designed by homology modeling, in which the streptavidin and rhizavidin sequences were combined to engineer a high affinity binding pocket containing residues from a single subunit only. Although mSA is stable and binds biotin with nanomolar affinity, its fast off rate (k_off) creates practical challenges during applications. We obtained a 1.9 Å crystal structure of mSA bound to biotin to understand their interaction in detail, and used the structure to introduce targeted mutations to improve its binding kinetics. To this end, we compared mSA to shwanavidin, which contains a hydrophobic lid containing F43 in the binding pocket and binds biotin tightly. However, the T48F mutation in mSA, which introduces a comparable hydrophobic lid, only resulted in a modest 20–40% improvement in the measured k_off. On the other hand, introducing the S25H mutation near the bicyclic ring of bound biotin increased the dissociation half life (t_½) from 11 to 83 min at 20°C. Molecular dynamics (MD) simulations suggest that H25 stabilizes the binding loop L3,4 by interacting with A47, and protects key intermolecular hydrogen bonds by limiting solvent entry into the binding pocket. Concurrent T48F or T48W mutation clashes with H25 and partially abrogates the beneficial effects of H25. Taken together, this study suggests that stabilization of the binding loop and solvation of the binding pocket are important determinants of the dissociation kinetics in mSA. Proteins 2013. © 2013 Wiley Periodicals, Inc.     

Lipid‐associated aggregate formation of superoxide dismutase‐1 is initiated by membrane‐targeting loops

Copper‐Zinc superoxide dismutase 1 (SOD1) is a homodimeric enzyme that protects cells from oxidative damage. Hereditary and sporadic amyotrophic lateral sclerosis may be linked to SOD1 when the enzyme is destabilized through mutation or environmental stress. The cytotoxicity of demetallated or apo‐SOD1 aggregates may be due to their ability to cause defects within cell membranes by co‐aggregating with phospholipids. SOD1 monomers may associate with the inner cell membrane to receive copper ions from membrane‐bound copper chaperones. But how apo‐SOD1 interacts with lipids is unclear. We have used atomistic molecular dynamics simulations to reveal that flexible electrostatic and zinc‐binding loops in apo‐SOD1 dimers play a critical role in the binding of 1‐octanol clusters and phospholipid bilayer, without any significant unfolding of the protein. The apo‐SOD1 monomer also associates with phospholipid bilayer via its zinc‐binding loop rather than its exposed hydrophobic dimerization interface. Our observed orientation of the monomer on the bilayer would facilitate its association with a membrane‐bound copper chaperone. The orientation also suggests how membrane‐bound monomers could act as seeds for membrane‐associated SOD1 aggregation. Proteins 2014; 82:3194–3209. © 2014 Wiley Periodicals, Inc.     

Identification of a determinant of epidermal growth factor receptor ligand-binding specificity using a truncated, high-affinity form of the ectodomain

Murine and human epidermal growth factor receptors (EGFRs) bind human EGF (hEGF), mouse EGF (mEGF), and human transforming growth factor alpha (hTGF-alpha) with high affinity despite the significant differences in the amino acid sequences of the ligands and the receptors. In contrast, the chicken EGFR can discriminate between mEGF (and hEGF) and hTGF-alpha and binds the EGFs with approximately 100-fold lower affinity. The regions responsible for this poor binding are known to be Arg(45) in hEGF and the L2 domain in the chicken EGFR. In this study we have produced a truncated form of the hEGFR ectodomain comprising residues 1-501 (sEGFR501), which, unlike the full-length hEGFR ectodomain (residues 1-621, sEGFR621), binds hEGF and hTGF-alpha with high affinity (K(D) = 13-21 and 35-40 nM, respectively). sEGFR501 was a competitive inhibitor of EGF-stimulated mitogenesis, being almost 10-fold more effective than the full-length EGFR ectodomain and three times more potent than the neutralizing anti-EGFR monoclonal antibody Mab528. Analytical ultracentrifugation showed that the primary EGF binding sites on sEGFR501 were saturated at an equimolar ratio of ligand and receptor, leading to the formation of a 2:2 EGF:sEGFR501 dimer complex. We have used sEGFR501 to generate three mutants with single position substitutions at Glu(367), Gly(441), or Glu(472) to Lys, the residue found in the corresponding positions in the chicken EGFR. All three mutants bound hTGF-alpha and were recognized by Mab528. However, mutant Gly(441)Lys showed markedly reduced binding to hEGF, implicating Gly(441), in the L2 domain, as part of the binding site that recognizes Arg(45) of hEGF.     

Insights into mechanism of pyrido[2,3‐d]pyrimidines as DYRK1A inhibitors based on molecular dynamic simulations

DYRK1A is characterized by the early development and regulation of neuronal proliferation, and its over expression gives rise to neurological abnormalities. As the promising DYRK1A inhibitors, the binding mechanism between DYRK1A and pyrido[2,3‐d]pyrimidines derivatives at molecular level are still veiled. In this article, it was achieved to get the structural insights into pyrido[2,3‐d]pyrimidines derivatives as DYRK1A inhibitors by means of comprehensive computational approaches involving molecular docking, molecular dynamics simulation, free energy calculation, and energy decomposition analysis. The calculated energy values were highly consistent with the experimental activities. Based on the individual energy terms analysis, the van der Waals interaction was the major leading force in the DYRK1A–ligand interaction. Lys188 was the important residue that formed the hydrogen bond, which improved the inhibitory activity. Furthermore, four novel inhibitors with higher predicted activity were designed based on the obtained findings and confirmed by molecular simulations. Our study is expected to provide significant drug design strategy for the development of more promising DYRK1A inhibitors. Proteins 2016; 84:1108–1123. © 2016 Wiley Periodicals, Inc.     

Hydration energy landscape of the active site cavity in cytochrome P450cam

Hydration of protein cavities influences protein stability, dynamics, and function. Protein active sites usually contain water molecules that, upon ligand binding, are either displaced into bulk solvent or retained to mediate protein–ligand interactions. The contribution of water molecules to ligand binding must be accounted for to compute accurate values of binding affinities. This requires estimation of the extent of hydration of the binding site. However, it is often difficult to identify the water molecules involved in the binding process when ligands bind on the surface of a protein. Cytochrome P450cam is, therefore, an ideal model system because its substrate binds in a buried active site, displacing partially disordered solvent, and the protein is well characterized experimentally. We calculated the free energy differences for having five to eight water molecules in the active site cavity of the unliganded enzyme from molecular dynamics simulations by thermodynamic integration employing a three-stage perturbation scheme. The computed free energy differences between the hydration states are small (within 12 kJ mol⁻¹) but distinct. Consistent with the crystallographic determination and studies employing hydrostatic pressure, we calculated that, although ten water molecules could in principle occupy the volume of the active site, occupation by five to six water molecules is thermodynamically most favorable. Proteins 32:381–396, 1998. © 1998 Wiley-Liss, Inc.     

Specificity of a protein–protein interface: Local dynamics direct substrate recognition of effector caspases

Proteases are prototypes of multispecific protein–protein interfaces. Proteases recognize and cleave protein and peptide substrates at a well‐defined position in a substrate binding groove and a plethora of experimental techniques provide insights into their substrate recognition. We investigate the caspase family of cysteine proteases playing a key role in programmed cell death and inflammation, turning caspases into interesting drug targets. Specific ligand binding to one particular caspase is difficult to achieve, as substrate specificities of caspase isoforms are highly similar. In an effort to rationalize substrate specificity of two closely related caspases, we investigate the substrate promiscuity of the effector Caspases 3 and 7 by data mining (cleavage entropy) and by molecular dynamics simulations. We find a strong correlation between binding site rigidity and substrate readout for individual caspase subpockets explaining more stringent substrate readout of Caspase 7 via its narrower conformational space. Caspase 3 subpockets S3 and S4 show elevated local flexibility explaining the more unspecific substrate readout of that isoform in comparison to Caspase 7. We show by in silico exchange mutations in the S3 pocket of the proteases that a proline residue in Caspase 7 contributes to the narrowed conformational space of the binding site. These findings explain the substrate specificities of caspases via a mechanism of conformational selection and highlight the crucial importance of binding site local dynamics in substrate recognition of proteases. Proteins 2014; 82:546–555. © 2013 Wiley Periodicals, Inc.     

Role of tyrosine hot‐spot residues at the interface of colicin E9 and immunity protein 9: A comparative free energy simulation study

The endonuclease activity of the bacterial colicin 9 enzyme is controlled by the specific and high‐affinity binding of immunity protein 9 (Im9). Molecular dynamics simulation studies in explicit solvent were used to investigate the free energy change associated with the mutation of two hot‐spot interface residues [tyrosine (Tyr): Tyr54 and Tyr55] of Im9 to Ala. In addition, the effect of several other mutations (Leu33Ala, Leu52Ala, Val34Ala, Val37Ala, Ser48Ala, and Ile53Ala) with smaller influence on binding affinity was also studied. Good qualitative agreement of calculated free energy changes and experimental data on binding affinity of the mutations was observed. The simulation studies can help to elucidate the molecular details on how the mutations influence protein–protein binding affinity. The role of solvent and conformational flexibility of the partner proteins was studied by comparing the results in the presence or absence of solvent and with or without positional restraints. Restriction of the conformational mobility of protein partners resulted in significant changes of the calculated free energies but of similar magnitude for isolated Im9 and for the complex and therefore in only modest changes of binding free energy differences. Although the overall binding free energy change was similar for the two Tyr–Ala mutations, the physical origin appeared to be different with solvation changes contributing significantly to the Tyr55Ala mutation and to a loss of direct protein–protein interactions dominating the free energy change due to the Tyr54Ala mutation. Proteins 2013. © 2012 Wiley Periodicals, Inc.     

Molecular dynamics simulations reveal structural insights into inhibitor binding modes and functionality in human Group IIA phospholipase A2

Human Group IIA phospholipase A₂ (hGIIA) promotes inflammation in immune‐mediated pathologies by regulating the arachidonic acid pathway through both catalysis‐dependent and ‐independent mechanisms. The hGIIA crystal structure, both alone and inhibitor‐bound, together with structures of closely related snake‐venom‐derived secreted phospholipase enzymes has been well described. However, differentiation of biological and nonbiological contacts and the relevance of structures determined from snake venom enzymes to human enzymes are not clear. We employed molecular dynamics (MD) and docking approaches to understand the binding of inhibitors that selectively or nonselectively block the catalysis‐independent mechanism of hGIIA. Our results indicate that hGIIA behaves as a monomer in the solution environment rather than a dimer arrangement that is in the asymmetric unit of some crystal structures. The binding mode of a nonselective inhibitor, KH064, was validated by a combination of the experimental electron density and MD simulations. The binding mode of the selective pentapeptide inhibitor FLSYK to hGIIA was stipulated to be different to that of the snake venom phospholipases A₂ of Daboia russelli pulchella (svPLA₂). Our data suggest that the application of MD approaches to crystal structure data is beneficial in evaluating the robustness of conclusions drawn based on crystal structure data alone. Proteins 2017; 85:827–842. © 2016 Wiley Periodicals, Inc.     

Structural insight into the UNC‐45–myosin complex

The UNC‐45 chaperone protein interacts with and affects the folding, stability, and the ATPase activity of myosins. It plays a critical role in the cardiomyopathy development and in the breast cancer tumor growth. Here we propose the first structural model of the UNC‐45–myosin complex using various in silico methods. Initially, the human UNC‐45B binding epitope was identified and the protein was docked to the cardiac myosin (MYH7) motor domain. The final UNC45B–MYH7 structure was obtained by performing of total 630 ns molecular dynamics simulations. The results indicate a complex formation, which is mainly stabilized by electrostatic interactions. Remarkably, the contact surface area is similar to that of the myosin‐actin complex. A significant interspecies difference in the myosin binding epitope is observed. Our results reveal the structural basis of MYH7 exons 15–16 hypertrophic cardiomyopathy mutations and provide directions for drug targeting. Proteins 2013; 81:1212–1221. © 2013 Wiley Periodicals, Inc.     

A molecular dynamics study of the binary complexes of APP,JIP1, and the cargo binding domain of KLC

Mutations in the amyloid precursor protein (APP) are responsible for the formation of amyloid‐β peptides. These peptides play a role in Alzheimer's and other dementia‐related diseases. The cargo binding domain of the kinesin‐1 light chain motor protein (KLC1) may be responsible for transporting APP either directly or via interaction with C‐jun N‐terminal kinase‐interacting protein 1 (JIP1). However, to date there has been no direct experimental or computational assessment of such binding at the atomistic level. We used molecular dynamics and free energy estimations to gauge the affinity for the binary complexes of KLC1, APP, and JIP1. We find that all binary complexes (KLC1:APP, KLC1:JIP1, and APP:JIP1) contain conformations with favorable binding free energies. For KLC1:APP the inclusion of approximate entropies reduces the favorability. This is likely due to the flexibility of the 42‐residue APP protein. In all cases we analyze atomistic/residue driving forces for favorable interactions. Proteins 2017; 85:221–234. © 2016 Wiley Periodicals, Inc.     

Folding and stability of helical bundle proteins from coarse‐grained models

We develop a coarse‐grained model where solvent is considered implicitly, electrostatics are included as short‐range interactions, and side‐chains are coarse‐grained to a single bead. The model depends on three main parameters: hydrophobic, electrostatic, and side‐chain hydrogen bond strength. The parameters are determined by considering three level of approximations and characterizing the folding for three selected proteins (training set). Nine additional proteins (containing up to 126 residues) as well as mutated versions (test set) are folded with the given parameters. In all folding simulations, the initial state is a random coil configuration. Besides the native state, some proteins fold into an additional state differing in the topology (structure of the helical bundle). We discuss the stability of the native states, and compare the dynamics of our model to all atom molecular dynamics simulations as well as some general properties on the interactions governing folding dynamics. Proteins 2013; 81:1200–1211. © 2013 Wiley Periodicals, Inc.     

Recruiting monomer for dimer formation: resolving the antagonistic mechanisms of novel immune check point inhibitors against Programmed Death Ligand-1 in cancer immunotherapy

The design of small molecule antagonists against Programmed Death Ligand-1 (PD-L1) has been the recent highlight of the immune checkpoint blockade therapy. This interventive approach has been potentiated by the development of BMS compounds; BMS-1001 and BMS-1166, which exert their therapeutic activities by inducing dimerisation of PD-L1; a molecular mechanism that has remained unclear. For the first time, we resolve the dynamical events that underlie the antagonistic mechanisms of BMS-1001 and BMS-1166 when bound to PD-L1 using an all-atom molecular dynamics (MD) simulations approach and free binding energy Molecular Mechanics/Poisson-Boltzmann Surface Area (MM/PBSA) calculations. Time-scale dynamical findings revealed that upon binding a PD-L1 monomer, the BMS-compounds gradually facilitated the ‘inbound’ motion of another PD-L1 monomer in the same conformational phase space up till dimer formation. Moreover, the non-liganded PD-L1 monomer exhibited the highest structural flexibility and atomistic motions relative to the BMS-liganded monomer as revealed by post-MD trajectory analyses using root mean square deviation (RMSD) and root mean square fluctuations (RMSF) parameters. Trajectory investigations into ligand motions also revealed that the BMS compounds exhibited mechanistic transitions from the monomeric binding site (monomer A) where they were initially bound, to the second monomeric site (monomer B) where they were strongly bound, followed by eventual high-affinity interactions at the tunnel-like binding cleft formed upon the dimerisation of both PD-L1 monomers. These findings present a model that describes the mechanism by which the BMS compounds induce PD-L1 dimerisation and could further enhance the design of highly selective and novel monomeric recruiters of PD-L1 in cancer immunotherapy.