Computing Protein-Protein Association Affinity with Hybrid Steered Molecular Dynamics

An insight from computational approach to explore novel, high-affinity phosphodiesterase 10A inhibitors for neurological disorders

The enzyme Phosphodiesterase 10A (PDE10A) plays a regulatory role in the cAMP/protein kinase A (PKA) signaling pathway by means of hydrolyzing cAMP and cGMP. PDE10A emerges as a relevant pharmacological drug target for neurological conditions such as psychosis, schizophrenia, Parkinson's, Huntington's disease, and other memory-related disorders. In the current study, we subjected a set of 1,2,3-triazoles to be explored as PDE10A inhibitors using diverse computational approaches, including molecular docking, classical molecular dynamics (MD) simulations, Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) calculations, steered MD, and umbrella sampling simulations. Molecular docking of cocrystallized ligands papaverine and PFJ, along with a set of in-house synthesized molecules, suggested that molecule 3i haded the highest binding affinity, followed by 3h and 3j. Furthermore, the structural stability studies using MD and MM-PBSA indicated that the 3h and 3j formed stable complexes with PDE10A. The binding free energy of -240.642 kJ/mol and -201.406 kJ/mol was observed for 3h and 3j, respectively. However, the cocrystallized ligands papaverine and PFJ exhibited comparitively higher binding free energy values of -202.030 kJ/mol and -138.764 kJ/mol, respectively. Additionally, steered MD and umbrella sampling simulations provided conclusive evidence that the molecules 3h and 3j could be exploited as promising candidates to target PDE10A.Communicated by Ramaswamy H. Sarma.

The α-Synuclein Monomer May Have Different Misfolding Mechanisms in the Induction of α-Synuclein Fibrils with Different Polymorphs

The aggregation of alpha-synuclein (α-Syn) is closely related to the occurrence of some neurodegenerative diseases such as Parkinson's disease. The misfolding of α-Syn monomer plays a key role in the formation of aggregates and extension of fibril. However, the misfolding mechanism of α-Syn remains elusive. Here, three different α-Syn fibrils (isolated from a diseased human brain, generated by in vitro cofactor-tau induction, and obtained by in vitro cofactor-free induction) were selected for the study. The misfolding mechanisms of α-Syn were uncovered by studying the dissociation of the boundary chains based on the conventional molecular dynamics (MD) and Steered MD simulations. The results showed that the dissociation paths of the boundary chains in the three systems were different. According to the reverse process of dissociation, we concluded that in the human brain system, the binding of the monomer and template starts from the C-terminal and gradually misfolds toward the N-terminal. In the cofactor-tau system, the monomer binding starts from residues 58-66 (contain β3), followed by the C-terminal coil (residues 67-79). Then, the N-terminal coil (residues 36-41) and residues 50-57 (contain β2) bind to the template, followed by residues 42-49 (contain β1). In the cofactor-free system, two misfolding paths were found. One is that the monomer binds to the N/C-terminal (β1/β6) and then binds to the remaining residues. The other one is that the monomer binds sequentially from the C- to N-terminal, similar to the human brain system. Furthermore, in the human brain and cofactor-tau systems, electrostatic interactions (especially from residues 58-66) are the main driving force during the misfolding process, whereas in the cofactor-free system, the contributions of electrostatic and van der Waals interactions are comparable. These results may provide a deeper understanding for the misfolding and aggregation mechanism of α-Syn.

Continuous millisecond conformational cycle of a DEAH box helicase reveals control of domain motions by atomic-scale transitions

Helicases are motor enzymes found in every living organism and viruses, where they maintain the stability of the genome and control against false recombination. The DEAH-box helicase Prp43 plays a crucial role in pre-mRNA splicing in unicellular organisms by translocating single-stranded RNA. The molecular mechanisms and conformational transitions of helicases are not understood at the atomic level. We present a complete conformational cycle of RNA translocation by Prp43 in atomic detail based on molecular dynamics simulations. To enable the sampling of such complex transition on the millisecond timescale, we combined two enhanced sampling techniques, namely simulated tempering and adaptive sampling guided by crystallographic data. During RNA translocation, the center-of-mass motions of the RecA-like domains followed the established inchworm model, whereas the domains crawled along the RNA in a caterpillar-like movement, suggesting an inchworm/caterpillar model. However, this crawling required a complex sequence of atomic-scale transitions involving the release of an arginine finger from the ATP pocket, stepping of the hook-loop and hook-turn motifs along the RNA backbone, and several others. These findings highlight that large-scale domain dynamics may be controlled by complex sequences of atomic-scale transitions.

Computer aided drug design in the development of proteolysis targeting chimeras

Proteolysis targeting chimeras represent a class of drug molecules with a number of attractive properties, most notably a potential to work for targets that, so far, have been in-accessible for conventional small molecule inhibitors. Due to their different mechanism of action, and physico-chemical properties, many of the methods that have been designed and applied for computer aided design of traditional small molecule drugs are not applicable for proteolysis targeting chimeras. Here we review recent developments in this field focusing on three aspects: de-novo linker-design, estimation of absorption for beyond-rule-of-5 compounds, and the generation and ranking of ternary complex structures. In spite of this field still being young, we find that a good number of models and algorithms are available, with the potential to assist the design of such compounds in-silico, and accelerate applied pharmaceutical research.

Aquaglyceroporin AQP7's affinity for its substrate glycerol: Have we reached convergence in the computed values of glycerol-aquaglyceroporin affinity?

AQP7 is one of the four human aquaglyceroporins that facilitate glycerol transport across the cell membrane, a biophysical process that is essential in human physiology. Therefore, it is interesting to compute AQP7's affinity for its substrate (glycerol) with reasonable certainty to compare with the experimental data suggesting high affinity in contrast with most computational studies predicting low affinity. In this study aimed at computing the AQP7-glycerol affinity with high confidence, we implemented a direct computation of the affinity from unbiased equilibrium molecular dynamics (MD) simulations of three all-atom systems constituted with 0.16 M, 4.32 M, and 10.23 M atoms, respectively. These three sets of simulations manifested a fundamental physics law that the intrinsic fluctuations of pressure in a system are inversely proportional to the system size (the number of atoms in it). These simulations showed that the computed values of glycerol-AQP7 affinity are dependent upon the system size (the inverse affinity estimations were, respectively, 47.3 mM, 1.6 mM, and 0.92 mM for the three model systems). In this, we obtained a lower bound for the AQP7-glycerol affinity (an upper bound for the dissociation constant). Namely, the AQP7-glycerol affinity is stronger than 1087/M (the dissociation constant is less than 0.92 mM). Additionally, we conducted hyper steered MD (hSMD) simulations to map out the Gibbs free-energy profile. From the free-energy profile, we produced an independent computation of the AQP7-glycerol dissociation constant being approximately 0.18 mM.

AQP7 is one of the four human aquaglyceroporins that facilitate glycerol transport across the cell membrane, a biophysical process that is essential in human physiology.

Serine 477 plays a crucial role in the interaction of the SARS-CoV-2 spike protein with the human receptor ACE2

Since the worldwide outbreak of the infectious disease COVID-19, several studies have been published to understand the structural mechanism of the novel coronavirus SARS-CoV-2. During the infection process, the SARS-CoV-2 spike (S) protein plays a crucial role in the receptor recognition and cell membrane fusion process by interacting with the human angiotensin-converting enzyme 2 (hACE2) receptor. However, new variants of these spike proteins emerge as the virus passes through the disease reservoir. This poses a major challenge for designing a potent antigen for an effective immune response against the spike protein. Through a normal mode analysis (NMA) we identified the highly flexible region in the receptor binding domain (RBD) of SARS-CoV-2, starting from residue 475 up to residue 485. Structurally, the position S477 shows the highest flexibility among them. At the same time, S477 is hitherto the most frequently exchanged amino acid residue in the RBDs of SARS-CoV-2 mutants. Therefore, using MD simulations, we have investigated the role of S477 and its two frequent mutations (S477G and S477N) at the RBD during the binding to hACE2. We found that the amino acid exchanges S477G and S477N strengthen the binding of the SARS-COV-2 spike with the hACE2 receptor.

Implicit ligand theory for relative binding free energies: II. An estimator based on control variates

Implicit ligand theory describes the relationship between the noncovalent binding free energy and the binding free energy between a ligand and multiple rigid receptor conformations. We have previously shown that if the receptor conformations are sampled from or reweighed to a holo ensemble, the binding free energy relative to the ligand that defines the ensemble can be calculated. Here, we apply a variance reduction technique known as control variates to derive a new statistical estimator for the relative binding free energy. In applications to a data set of 6 reference ligands and 18 test ligands, statistically significant differences between the estimators are not observed for most systems. However, in cases where such differences are observed, the new estimator is more accurate, precise, and converges more quickly. Performance improvements are most consistent where there is a clear correlation, with a correlation coefficient greater than 0.3, between the control variate and the statistic being averaged.

Quantitative study of unsaturated transport of glycerol through aquaglyceroporin that has high affinity for glycerol

The structures of several aquaglyceroporins have been resolved to atomic resolution showing two or more glycerols bound inside a channel and confirming a glycerol-facilitator's affinity for its substrate glycerol. However, the kinetics data of glycerol transport experiments all point to unsaturated transport that is characteristic of low substrate affinity in terms of the Michaelis-Menten kinetics. In this article, we present an in silico-in vitro research focused on AQP3, one of the human aquaglyceroporins that is natively expressed in the abundantly available erythrocytes. We conducted 2.1 μs in silico simulations of AQP3 embedded in a model erythrocyte membrane with intracellular-extracellular asymmetries in leaflet lipid compositions and compartment salt ions. From the equilibrium molecular dynamics (MD) simulations, we elucidated the mechanism of glycerol transport at high substrate concentrations. From the steered MD simulations, we computed the Gibbs free-energy profile throughout the AQP3 channel. From the free-energy profile, we quantified the kinetics of glycerol transport that is unsaturated due to glycerol-glycerol interactions mediated by AQP3 resulting in the concerted movement of two glycerol molecules for the transport of one glycerol molecule across the cell membrane. We conducted in vitro experiments on glycerol uptake into human erythrocytes for a wide range of substrate concentrations and various temperatures. The experimental data quantitatively validated our theoretical-computational conclusions on the unsaturated glycerol transport through AQP3 that has high affinity for glycerol.

Thermodynamic Integration in 3n Dimensions without Biases or Alchemy for Protein Interactions

Thermodynamic integration (TI), a powerful formalism for computing Gibbs free energy, has been implemented for many biophysical processes with alchemical schemes that require delicate human efforts to choose/design biasing potentials for sampling the desired biophysical events and to remove their artifactitious consequences afterwards. Theoretically, an alchemical scheme is exact but practically, an unsophisticated implementation of this exact formula can cause error amplifications. Small relative errors in the input parameters can be amplified many times in their propagation into the computed free energy [due to subtraction of similar numbers such as (105 ± 5)‒(100 ± 5) = 5 ± 7]. In this paper, we present an unsophisticated implementation of TI in 3n dimensions (3nD) (n=1,2,3…) for the potential of mean force along a 3nD path connecting one state in the bound state ensemble to one state in the unbound state ensemble. Fluctuations in these 3nD are integrated in the bound and unbound state ensembles but not along the 3nD path. Using TI3nD, we computed the standard binding free energies of three protein complexes: trometamol in Salmonella effector SpvD (n=1), biotin-avidin (n=2), and Colicin E9 endonuclease with cognate immunity protein Im9 (n=3). We employed three different protocols in three independent computations of E9-Im9 to show TI3nD's robustness. We also computed the hydration energies of ten biologically relevant compounds (n=1 for water, acetamide, urea, glycerol, trometamol, ammonium and n=2 for erythritol, 1,3-propanediol, xylitol, biotin). Each of the 15 computations is accomplishable within one (for hydration) to ten (for E9-Im9) days on an inexpensive GPU workstation. The computed results all agree with the available experimental data.

Calculating the absolute binding free energy of the insulin dimer in an explicit solvent

Insulin is a significant hormone in the regulation of glucose level in the blood. Its monomers bind to each other to form dimers or hexamers through a complex process. To study the binding of the insulin dimer, we first calculate its absolute binding free energy by the steered molecular dynamics method and the confinement method based on a fictitious thermodynamic cycle. After considering some special correction terms, the final calculated binding free energy at 298 K is −8.97 ± 1.41 kcal mol⁻¹, which is close to the experimental value of −7.2 ± 0.8 kcal mol⁻¹. Furthermore, we discuss the important residue–residue interactions between the insulin monomers, including hydrophobic interactions, π–π interactions and hydrogen bond interactions. The analysis reveals five key residues, Vla^B12, Tyr^B16, Phe^B24, Phe^B25, and Tyr^B26, for the dimerization of the insulin. We also perform MM-PBSA calculations for the wild-type dimer and some mutants and study the roles of the key residues by the change of the binding energy of the insulin dimer.

In this paper, we calculate the absolute binding free energy of an insulin dimer by steered MD method. The result of −8.97 kcal mol⁻¹ is close to the experimental value −7.2 kcal mol⁻¹. We also analyze the residue–residue interactions.

Disclosing the Template-Induced Misfolding Mechanism of Tau Protein by Studying the Dissociation of the Boundary Chain from the Formed Tau Fibril Based on a Steered Molecular Dynamics Simulation

The level of tau aggregation into neurofibrillary tangles, including paired helical filament (PHF) and straight filament (SF), is closely associated with Alzheimer's disease. Despite the pathological importance of misfolding and aggregation of tau, the corresponding mechanism remains unclear. Therefore, to uncover the misfolding mechanism of the tau monomer upon induction of formed PHF and SF, in this study, a conventional molecular dynamics simulation combined with a steered molecular dynamics simulation was performed to study the dissociation of the boundary chain. Interestingly, our results show that the dissociation mechanisms of the boundary chain in PHF and SF are different. In PHF, the boundary chain begins to dissociate from regions β2 and β3 and ends at β8. However, in SF, it is simultaneously dissociated from β1 and β8 and ends at β5. The dissociation of the boundary chain is the reverse of template-induced misfolding of the monomer. Therefore, we can deduce the misfolding mechanism of the monomer upon induction of the template. For PHF, β8 first interacts with the template by hydrophobic interaction. Then β7, β6, β5, β4, and β1 sequentially bind to the template by electrostatic and hydrophobic interactions. After β1 binds to the template, β2 and β3 very quickly bind to the template through hydrophobic interaction. For SF, β5 of the monomer first interacts with the template by electrostatic attraction. Then β4 and β6, β3 and β7, and β2 and β8 bind to the template in turn. Finally, β1 and β8 are fully bound to the template by hydrophobic interaction. The obtained results will be vital for understanding the earlier events during misfolding and aggregation of tau.

On the mechanisms of protein interactions: predicting their affinity from unbound tertiary structures

Motivation

The characterization of the protein–protein association mechanisms is crucial to understanding how biological processes occur. It has been previously shown that the early formation of non-specific encounters enhances the realization of the stereospecific (i.e. native) complex by reducing the dimensionality of the search process. The association rate for the formation of such complex plays a crucial role in the cell biology and depends on how the partners diffuse to be close to each other. Predicting the binding free energy of proteins provides new opportunities to modulate and control protein–protein interactions. However, existing methods require the 3D structure of the complex to predict its affinity, severely limiting their application to interactions with known structures.

Results

We present a new approach that relies on the unbound protein structures and protein docking to predict protein–protein binding affinities. Through the study of the docking space (i.e. decoys), the method predicts the binding affinity of the query proteins when the actual structure of the complex itself is unknown. We tested our approach on a set of globular and soluble proteins of the newest affinity benchmark, obtaining accuracy values comparable to other state-of-art methods: a 0.4 correlation coefficient between the experimental and predicted values of ΔG and an error < 3 Kcal/mol.

Availability and implementation

The binding affinity predictor is implemented and available at http://sbi.upf.edu/BADock and https://github.com/badocksbi/BADock.

Supplementary information

Supplementary data are available at Bioinformatics online.

Gibbs Free-Energy Gradient along the Path of Glucose Transport through Human Glucose Transporter 3

Fourteen glucose transporters (GLUTs) play essential roles in human physiology by facilitating glucose diffusion across the cell membrane. Due to its central role in the energy metabolism of the central nervous system, GLUT3 has been thoroughly investigated. However, the Gibbs free-energy gradient (what drives the facilitated diffusion of glucose) has not been mapped out along the transport path. Some fundamental questions remain. Here we present a molecular dynamics study of GLUT3 embedded in a lipid bilayer to quantify the free-energy profile along the entire transport path of attracting a β-d-glucose from the interstitium to the inside of GLUT3 and, from there, releasing it to the cytoplasm by Arrhenius thermal activation. From the free-energy profile, we elucidate the unique Michaelis–Menten characteristics of GLUT3, low K_M and high V_MAX, specifically suitable for neurons’ high and constant demand of energy from their low-glucose environments. We compute GLUT3’s binding free energy for β-d-glucose to be −4.6 kcal/mol in agreement with the experimental value of −4.4 kcal/mol (K_M = 1.4 mM). We also compute the hydration energy of β-d-glucose, −18.0 kcal/mol vs the experimental data, −17.8 kcal/mol. In this, we establish a dynamics-based connection from GLUT3’s crystal structure to its cellular thermodynamics with quantitative accuracy. We predict equal Arrhenius barriers for glucose uptake and efflux through GLUT3 to be tested in future experiments.

Atomistic Analysis of ToxN and ToxI Complex Unbinding Mechanism

ToxIN is a triangular structure formed by three protein toxins (ToxNs) and three specific noncoding RNA antitoxins (ToxIs). To respond to stimuli, ToxI is preferentially degraded, releasing the ToxN. Thus, the dynamic character is essential in the normal function interactions between ToxN and ToxI. Here, equilibrated molecular dynamics (MD) simulations were performed to study the stability of ToxN and ToxI. The results indicate that ToxI adjusts the conformation of 3' and 5' termini to bind to ToxN. Steered molecular dynamics (SMD) simulations combined with the recently developed thermodynamic integration in 3nD (TI3nD) method were carried out to investigate ToxN unbinding from the ToxIN complex. The potentials of mean force (PMFs) and atomistic pictures suggest the unbinding mechanism as follows: (1) dissociation of the 5' terminus from ToxN, (2) missing the interactions involved in the 3' terminus of ToxI without three nucleotides (G31, A32, and A33), (3) starting to unfold for ToxI, (4) leaving the binding package of ToxN for three nucleotides of ToxI, (5) unfolding of ToxI. This work provides information on the structure-function relationship at the atomistic level, which is helpful for designing new potent antibacterial drugs in the future.

Simple Entropy Terms for End-Point Binding Free Energy Calculations

We introduce a number of computationally inexpensive modifications to the MM/PBSA and MM/GBSA estimators for binding free energies, which are based on average receptor-ligand interaction energies in simulations of a noncovalent complex, to improve the treatment of entropy: second- and higher-order terms in a cumulant expansion and a confining potential on ligand external degrees of freedom. We also consider a filter for snapshots where ligands have drifted from the initial binding pose. The variations were tested on six sets of systems for which binding modes and free energies have previously been experimentally determined. For some data sets, none of the tested estimators led to results significantly correlated with measured free energies. In data sets with nontrivial correlation, a ligand RMSD cutoff of 3 Å and a second-order truncation of the cumulant expansion was found to be comparable or better than the average interaction energy by several statistical metrics.

Affinity and path of binding xylopyranose unto E. coli xylose permease

Glucose transporters (GLUTs), expressed in all types of human cells, are responsible for the uptake of sugars as the primary energy source for the normal functions of good cells and for the abnormal growth of cancer cells. The E. coli xylose permease (XylE), a homologue of human GLUTs, has been investigated more thoroughly than other major facilitator proteins in the current literature. In this paper, we present a molecular dynamics (MD) study of an all-atom model system to elucidate the atomistic details and the free-energy landscape along the path of binding a xylopyranose (XYP) from the extracellular space to the inside of the transporter protein XylE. From the MD simulations, the Gibbs free energy of binding was found to be -4.4kcal/mol in agreement with the experimental value of -4.7kcal/mol. The accuracy of our study is further shown in the computed hydration energy of XYP of -14.6kcal/mol in comparison with the experimental data of -15.0kcal/mol. Along the binding path, the Gibbs free energy of the XYP-XylE complex first rises from zero in the dissociated state to approximately 4 kcal/mol in the transition state (when XylE slightly increases its opening toward the extracellular side to accommodate XYP) before dropping down to -9.0 kcal/mol in the bound state. These quantitative insights indicate the fast equilibration between the bound and the unbound states of XylE and XYP. They also serve as an atomistic-dynamic corroboration of the experimental conclusion that XylE is a high-affinity sugar transporter.

Simulation of Reversible Protein-Protein Binding and Calculation of Binding Free Energies Using Perturbed Distance Restraints

Virtually all biological processes depend on the interaction between proteins at some point. The correct prediction of biomolecular binding free-energies has many interesting applications in both basic and applied pharmaceutical research. While recent advances in the field of molecular dynamics (MD) simulations have proven the feasibility of the calculation of protein–protein binding free energies, the large conformational freedom of proteins and complex free energy landscapes of binding processes make such calculations a difficult task. Moreover, convergence and reversibility of resulting free-energy values remain poorly described. In this work, an easy-to-use, yet robust approach for the calculation of standard-state protein–protein binding free energies using perturbed distance restraints is described. In the binding process the conformations of the proteins were restrained, as suggested earlier. Two approaches to avoid end-state problems upon release of the conformational restraints were compared. The method was evaluated by practical application to a small model complex of ubiquitin and the very flexible ubiquitin-binding domain of human DNA polymerase ι (UBM2). All computed free energy differences were closely monitored for convergence, and the calculated binding free energies had a mean unsigned deviation of only 1.4 or 2.5 kJ·mol^–1 from experimental values. Statistical error estimates were in the order of thermal noise. We conclude that the presented method has promising potential for broad applicability to quantitatively describe protein–protein and various other kinds of complex formation.

Elongation affinity, activation barrier, and stability of Aβ42 oligomers/fibrils in physiological saline

Amyloid-beta (Aβ) peptides, Aβ40 and the more neurotoxic Aβ42, have been the subject of many research efforts for Alzheimer's disease. In two recent independent investigations, the atomistic structure of Aβ42 fibril has been clearly established in the S-shaped conformation consisting of three β-sheets stabilized by salt bridges formed between the Lys28 sidechain and the C-terminus of Ala42. This structure distinctively differs from the long-known structure of Aβ40 in the β-hairpin shaped conformation consisting of two β-sheets. Recent in silico investigations based on all-atom models have reached closer agreement with the in vitro measurements of Aβ40 thermodynamics. In this study, we present an in silico investigation of Aβ42 thermodynamics. Using the established force field parameters in seven sets of all-atom simulations, we examined the stability of small Aβ42 oligomers in physiological saline. We computed the elongation affinity of the S-shaped Aβ42 fibril, reaching agreement with the experimental data. We also estimated the Arrhenius activation barrier along the elongation pathway (from the disordered conformation of a free Aβ42 peptide to its S-shaped conformation on a fibril) that amounts to about 16 kcal/mol, which is consistent with the experimental data. Based on these quantitative agreements, we conclude that aggregation of Aβ42 peptides into fibrils is thermodynamically slow without precipitation by extrinsic factors such as heparan sulfate proteoglycan and highlight the possibility to prevent Aβ42 aggregation by eliminating some precipitation factors or by increasing competitive agents to capture and transport free Aβ42 peptides from the cerebrospinal fluid.

Computing the binding affinity of a ligand buried deep inside a protein with the hybrid steered molecular dynamics

Computing the ligand-protein binding affinity (or the Gibbs free energy) with chemical accuracy has long been a challenge for which many methods/approaches have been developed and refined with various successful applications. False positives and, even more harmful, false negatives have been and still are a common occurrence in practical applications. Inevitable in all approaches are the errors in the force field parameters we obtain from quantum mechanical computation and/or empirical fittings for the intra- and inter-molecular interactions. These errors propagate to the final results of the computed binding affinities even if we were able to perfectly implement the statistical mechanics of all the processes relevant to a given problem. And they are actually amplified to various degrees even in the mature, sophisticated computational approaches. In particular, the free energy perturbation (alchemical) approaches amplify the errors in the force field parameters because they rely on extracting the small differences between similarly large numbers. In this paper, we develop a hybrid steered molecular dynamics (hSMD) approach to the difficult binding problems of a ligand buried deep inside a protein. Sampling the transition along a physical (not alchemical) dissociation path of opening up the binding cavity---pulling out the ligand---closing back the cavity, we can avoid the problem of error amplifications by not relying on small differences between similar numbers. We tested this new form of hSMD on retinol inside cellular retinol-binding protein 1 and three cases of a ligand (a benzylacetate, a 2-nitrothiophene, and a benzene) inside a T4 lysozyme L99A/M102Q(H) double mutant. In all cases, we obtained binding free energies in close agreement with the experimentally measured values. This indicates that the force field parameters we employed are accurate and that hSMD (a brute force, unsophisticated approach) is free from the problem of error amplification suffered by many sophisticated approaches in the literature.

1,3-propanediol binds deep inside the channel to inhibit water permeation through aquaporins

Aquaporins and aquaglyceroporins (AQPs) are membrane channel proteins responsible for transport of water and for transport of glycerol in addition to water across the cell membrane, respectively. They are expressed throughout the human body and also in other forms of life. Inhibitors of human AQPs have been sought for therapeutic treatment for various medical conditions including hypertension, refractory edema, neurotoxic brain edema, and so forth. Conducting all-atom molecular dynamics simulations, we computed the binding affinity of acetazolamide to human AQP4 that agrees closely with in vitro experiments. Using this validated computational method, we found that 1,3-propanediol (PDO) binds deep inside the AQP4 channel to inhibit that particular aquaporin efficaciously. Furthermore, we used the same method to compute the affinities of PDO binding to four other AQPs and one aquaglyceroporin whose atomic coordinates are available from the protein data bank (PDB). For bovine AQP1, human AQP2, AQP4, AQP5, and Plasmodium falciparum PfAQP whose structures were resolved with high resolution, we obtained definitive predictions on the PDO dissociation constant. For human AQP1 whose PDB coordinates are less accurate, we estimated the dissociation constant with a rather large error bar. Taking into account the fact that PDO is generally recognized as safe by the US FDA, we predict that PDO can be an effective diuretic which directly modulates water flow through the protein channels. It should be free from the serious side effects associated with other diuretics that change the hydro-homeostasis indirectly by altering the osmotic gradients.

Molecular dynamics simulations on the effect of size and shape on the interactions between negative Au18(SR)14, Au102(SR)44 and Au144(SR)60 nanoparticles in physiological saline

Molecular dynamics simulations employing all-atom force fields have become a reliable way to study binding interactions quantitatively for a wide range of systems. In this work, we employ two recently developed methods for the calculation of dissociation constants KD between gold nanoparticles (AuNPs) of different sizes in a near-physiological environment through the potential of mean force (PMF) formalism: the method of geometrical restraints developed by Woo et al. and formalized by Gumbart et al. and the method of hybrid Steered Molecular Dynamics (hSMD). Obtaining identical results (within the margin of error) from both approaches on the negatively charged Au18(SR)14 NP, functionalized by the negatively charged 4-mercapto-benzoate (pMBA) ligand, we draw parallels between their energetic and entropic interactions. By applying the hSMD method on Au102(SR)44 and Au144(SR)60, both of them near-spherical in shape and functionalized by pMBA, we study the effects of size and shape on the binding interactions. Au18 binds weakly with KD = 13mM as a result of two opposing effects: its large surface curvature hindering the formation of salt bridges, and its large ligand density on preferential orientations favoring their formation. On the other hand, Au102 binds more strongly with KD = 30μM and Au144 binds the strongest with KD = 3.2nM.

Ranking protein-protein docking results using steered molecular dynamics and potential of mean force calculations

Crystallization of protein-protein complexes can often be problematic and therefore computational structural models are often relied on. Such models are often generated using protein-protein docking algorithms, where one of the main challenges is selecting which of several thousand potential predictions represents the most near-native complex. We have developed a novel technique that involves the use of steered molecular dynamics (sMD) and umbrella sampling to identify near-native complexes among protein-protein docking predictions. Using this technique, we have found a strong correlation between our predictions and the interface RMSD (iRMSD) in ten diverse test systems. On two of the systems, we investigated if the prediction results could be further improved using potential of mean force calculations. We demonstrated that a near-native (<2.0 Å iRMSD) structure could be identified in the top-1 ranked position for both systems. © 2016 Wiley Periodicals, Inc.

Molecular dynamics study of human carbonic anhydrase II in complex with Zn(2+) and acetazolamide on the basis of all-atom force field simulations

Human carbonic anhydrase II (hCAII) represents an ultimate example of the perfectly efficient metalloenzymes, which is capable of catalyzing the hydration of carbon dioxide with a rate approaching the diffusion controlled limit. Extensive experimental studies of this physiologically important metalloprotein have been done to elucidate the fundamentals of its enzymatic actions: what residues anchor the Zn(2+) (or another divalent cation) at the bottom of the binding pocket; how the relevant residues work concertedly with the divalent cation in the reversible conversions between CO2 and HCO3(-); what are the protonation states of the relevant residues and acetazolamide, an inhibitor complexed with hCAII, etc. In this article, we present a detailed computational study on the basis of the all-atom CHARMM force field where Zn(2+) is represented with a simple model of divalent cation using the transferrable parameters available from the current literature. We compute the hydration free energy of Zn(2+), the characteristics of hCAII-Zn(2+) complexation, and the absolute free energy of binding acetazolamide to the hCAII-Zn(2+) complex. In each of these three problems, our computed results agree with the experimental data within the known margin of error without making any case-by-case adjustments to the parameters. The quantitatively accurate insights we gain in this all-atom molecular dynamics study should be helpful in the search and design of more specific inhibitors of this and other carbonic anhydrases.

1,3-Propanediol binds inside the water-conducting pore of aquaporin 4: Does this efficacious inhibitor have sufficient potency?

Among the thirteen types of water channel proteins, aquaporins (AQPs), which play various essential roles in human physiology, AQP4 is richly expressed in cells of the central nervous system and implicated in pathological conditions such as brain edema. Therefore, researchers have been looking for ways to inhibit AQP4's water-conducting function. Many small molecules have been investigated for their interactions with the residues that form the AQP4 channel entry vestibule on the extracellular side and their interruption of waters entering into the conducting pore. Conducting all-atom simulations on the basis of CHARMM 36 force field, we study one such inhibitor, 5-acetamido-1,3,4-thiadiazole-2-sulfonamide (AZM), to achieve quantitative agreement between the computed and the experimentally measured values of AZM-AQP4 binding affinity. Using the same method, we examine the possibility of plugging up the AQP4 channel around the Asn-Pro-Ala motifs located near the channel center because a small molecule bound there would totally occlude water conduction through AQP4. We compute the binding affinities of 1,2-ethanediol (EDO) and 1,3-propanediol (PDO) inside the AQP4 conducting pore and identify the specificities of the interactions. The EDO-AQP4 interaction is weak with a dissociation constant of 80 mM. The PDO-AQP4 interaction is rather strong with a dissociation constant of 328 μM, which indicates that PDO is an efficacious AQP4 inhibitor with sufficiently high potency. Considering the fact that PDO is classified by the US Food and Drug Administration as generally safe, we predict that 1,3-propanediol could be an effective drug for brain edema and other AQP4-correlated neurological conditions.