Quantitative Prediction of Dissociation Rates of PYK2 Ligands Using Umbrella Sampling and Milestoning.

We used umbrella sampling and the milestoning simulation method to study the dissociation of multiple ligands from protein kinase PYK2. The activation barriers obtained from the potential of mean force of the umbrella sampling simulations correlated well with the experimental dissociation rates. Using the zero-temperature string method, we obtained optimized paths along the free-energy surfaces for milestoning simulations of three ligands with a similar structure. The milestoning simulations gave an absolute dissociation rate within 2 orders of magnitude of the experimental value for two ligands but at least 3 orders of magnitude too high for the third. Despite the similarity in their structures, the ligands took different pathways to exit from the binding site of PYK2, making contact with different sets of residues. In addition, the protein experienced different conformational changes for dissociation of the three ligands.

Simulation of ligand dissociation kinetics from the protein kinase PYK2.

Early-stage drug discovery projects often focus on equilibrium binding affinity to the target alongside selectivity and other pharmaceutical properties. The kinetics of drug binding are ignored but can have significant influence on drug efficacy. Therefore, increasing attention has been paid on evaluating drug-binding kinetics early in a drug discovery process. Simulating drug-binding kinetics at the atomic level is challenging for the long time scale involved. Here, we used the transition-based reweighting analysis method (TRAM) with the Markov state model to study the dissociation of a ligand from the protein kinase PYK2. TRAM combines biased and unbiased simulations to reduce computational costs. This work used the umbrella sampling technique for the biased simulations. Although using the potential of mean force from umbrella sampling simulations with the transition-state theory over-estimated the dissociation rate by three orders of magnitude, TRAM gave a dissociation rate within an order of magnitude of the experimental value.

Qualitative Prediction of Ligand Dissociation Kinetics from Focal Adhesion Kinase Using Steered Molecular Dynamics.

Most early-stage drug discovery projects focus on equilibrium binding affinity to the target alongside selectivity and other pharmaceutical properties. Since many approved drugs have nonequilibrium binding characteristics, there has been increasing interest in optimizing binding kinetics early in the drug discovery process. As focal adhesion kinase (FAK) is an important drug target, we examine whether steered molecular dynamics (SMD) can be useful for identifying drug candidates with the desired drug-binding kinetics. In simulating the dissociation of 14 ligands from FAK, we find an empirical power-law relationship between the simulated time needed for ligand unbinding and the experimental rate constant for dissociation, with a strong correlation depending on the SMD force used. To improve predictions, we further develop regression models connecting experimental dissociation rate with various structural and energetic quantities derived from the simulations. These models can be used to predict dissociation rates from FAK for related compounds.

Should Virus Capsids Assemble Perfectly? Theory and Observation of Defects.

Although published structural models of viral capsids generally exhibit a high degree of regularity or symmetry, structural defects might be expected because of the fluctuating environment in which capsids assemble and the requirement of some capsids for disassembly before genome delivery. Defective structures are observed in computer simulations, and are evident in single-particle cryoelectron microscopy studies. Here, we quantify the conditions under which defects might be expected, using a statistical mechanics model allowing for ideal, defective, and vacant sites. The model displays a threshold in affinity parameters below which there is an appreciable population of defective capsids. Even when defective sites are not allowed, there is generally some population of vacancies. Analysis of single particles in cryoelectron microscopy micrographs yields a confirmatory ≳15% of defective particles. Our findings suggest structural heterogeneity in virus capsids may be under-appreciated, and also points to a nontraditional strategy for assembly inhibition.

Middle-way flexible docking: Pose prediction using mixed-resolution Monte Carlo in estrogen receptor α.

There is a vast gulf between the two primary strategies for simulating protein-ligand interactions. Docking methods significantly limit or eliminate protein flexibility to gain great speed at the price of uncontrolled inaccuracy, whereas fully flexible atomistic molecular dynamics simulations are expensive and often suffer from limited sampling. We have developed a flexible docking approach geared especially for highly flexible or poorly resolved targets based on mixed-resolution Monte Carlo (MRMC), which is intended to offer a balance among speed, protein flexibility, and sampling power. The binding region of the protein is treated with a standard atomistic force field, while the remainder of the protein is modeled at the residue level with a Go model that permits protein flexibility while saving computational cost. Implicit solvation is used. Here we assess three facets of the MRMC approach with implications for other docking studies: (i) the role of receptor flexibility in cross-docking pose prediction; (ii) the use of non-equilibrium candidate Monte Carlo (NCMC) and (iii) the use of pose-clustering in scoring. We examine 61 co-crystallized ligands of estrogen receptor α, an important cancer target known for its flexibility. We also compare the performance of the MRMC approach with Autodock smina. Adding protein flexibility, not surprisingly, leads to significantly lower total energies and stronger interactions between protein and ligand, but notably we document the important role of backbone flexibility in the improvement. The improved backbone flexibility also leads to improved performance relative to smina. Somewhat unexpectedly, our implementation of NCMC leads to only modestly improved sampling of ligand poses. Overall, the addition of protein flexibility improves the performance of docking, as measured by energy-ranked poses, but we do not find significant improvements based on cluster information or the use of NCMC. We discuss possible improvements for the model including alternative coarse-grained force fields, improvements to the treatment of solvation, and adding additional types of NCMC moves.

Entire-Dataset Analysis of NMR Fast-Exchange Titration Spectra: A Mg2+ Titration Analysis for HIV-1 Ribonuclease H Domain.

This article communicates our study to elucidate the molecular determinants of weak Mg2+ interaction with the ribonuclease H (RNH) domain of HIV-1 reverse transcriptase in solution. As the interaction is weak (a ligand-dissociation constant >1 mM), nonspecific Mg2+ interaction with the protein or interaction of the protein with other solutes that are present in the buffer solution can confound the observed Mg2+-titration data. To investigate these indirect effects, we monitored changes in the chemical shifts of backbone amides of RNH by recording NMR 1H-15N heteronuclear single-quantum coherence spectra upon titration of Mg2+ into an RNH solution. We performed the titration under three different conditions: (1) in the absence of NaCl, (2) in the presence of 50 mM NaCl, and (3) at a constant 160 mM Cl- concentration. Careful analysis of these three sets of titration data, along with molecular dynamics simulation data of RNH with Na+ and Cl- ions, demonstrates two characteristic phenomena distinct from the specific Mg2+ interaction with the active site: (1) weak interaction of Mg2+, as a salt, with the substrate-handle region of the protein and (2) overall apparent lower Mg2+ affinity in the absence of NaCl compared to that in the presence of 50 mM NaCl. A possible explanation may be that the titrated MgCl2 is consumed as a salt and interacts with RNH in the absence of NaCl. In addition, our data suggest that Na+ increases the kinetic rate of the specific Mg2+ interaction at the active site of RNH. Taken together, our study provides biophysical insight into the mechanism of weak metal interaction on a protein.

Structural integrity of the ribonuclease H domain in HIV-1 reverse transcriptase.

The mature form of reverse transcriptase (RT) is a heterodimer comprising the intact 66-kDa subunit (p66) and a smaller 51-kDa subunit (p51) that is generated by removal of most of the RNase H (RNH) domain from a p66 subunit by proteolytic cleavage between residues 440 and 441. Viral infectivity is eliminated by mutations such as F440A and E438N in the proteolytic cleavage sequence, while normal processing and virus infectivity are restored by a compensatory mutation, T477A, that is located more than 10 Å away from the processing site. The molecular basis for this compensatory effect has remained unclear. We therefore investigated structural characteristics of RNH mutants using computational and experimental approaches. Our Nuclear Magnetic Resonance and Differential Scanning Fluorimetry results show that both F440A and E438N mutations disrupt RNH folding. Addition of the T477A mutation restores correct folding of the RNH domain despite the presence of the F440A or E438N mutations. Molecular dynamics simulations suggest that the T477A mutation affects the processing site by altering relative orientations of secondary structure elements. Predictions of sequence tolerance suggest that phenylalanine and tyrosine are structurally preferred at residues 440 and 441, respectively, which are the P1 and P1' substrate residues known to require bulky side chains for substrate specificity. Interestingly, our study demonstrates that the processing site residues, which are critical for protease substrate specificity and must be exposed to the solvent for efficient processing, also function to maintain proper RNH folding in the p66/p51 heterodimer.

Tabulation as a high-resolution alternative to coarse-graining protein interactions: Initial application to virus capsid subunits.

Traditional coarse-graining based on a reduced number of interaction sites often entails a significant sacrifice of chemical accuracy. As an alternative, we present a method for simulating large systems composed of interacting macromolecules using an energy tabulation strategy previously devised for small rigid molecules or molecular fragments [S. Lettieri and D. M. Zuckerman, J. Comput. Chem. 33, 268-275 (2012); J. Spiriti and D. M. Zuckerman, J. Chem. Theory Comput. 10, 5161-5177 (2014)]. We treat proteins as rigid and construct distance and orientation-dependent tables of the interaction energy between them. Arbitrarily detailed interactions may be incorporated into the tables, but as a proof-of-principle, we tabulate a simple α-carbon Go-like model for interactions between dimeric subunits of the hepatitis B viral capsid. This model is significantly more structurally realistic than previous models used in capsid assembly studies. We are able to increase the speed of Monte Carlo simulations by a factor of up to 6700 compared to simulations without tables, with only minimal further loss in accuracy. To obtain further enhancement of sampling, we combine tabulation with the weighted ensemble (WE) method, in which multiple parallel simulations are occasionally replicated or pruned in order to sample targeted regions of a reaction coordinate space. In the initial study reported here, WE is able to yield pathways of the final ∼25% of the assembly process.

Tunable Coarse Graining for Monte Carlo Simulations of Proteins via Smoothed Energy Tables: Direct and Exchange Simulations.

Many commonly used coarse-grained models for proteins are based on simplified interaction sites and consequently may suffer from significant limitations, such as the inability to properly model protein secondary structure without the addition of restraints. Recent work on a benzene fluid (Lettieri S.; Zuckerman D. M.J. Comput. Chem.2012, 33, 268-275) suggested an alternative strategy of tabulating and smoothing fully atomistic orientation-dependent interactions among rigid molecules or fragments. Here we report our initial efforts to apply this approach to the polar and covalent interactions intrinsic to polypeptides. We divide proteins into nearly rigid fragments, construct distance and orientation-dependent tables of the atomistic interaction energies between those fragments, and apply potential energy smoothing techniques to those tables. The amount of smoothing can be adjusted to give coarse-grained models that range from the underlying atomistic force field all the way to a bead-like coarse-grained model. For a moderate amount of smoothing, the method is able to preserve about 70-90% of the α-helical structure while providing a factor of 3-10 improvement in sampling per unit computation time (depending on how sampling is measured). For a greater amount of smoothing, multiple folding-unfolding transitions of the peptide were observed, along with a factor of 10-100 improvement in sampling per unit computation time, although the time spent in the unfolded state was increased compared with less smoothed simulations. For a ß hairpin, secondary structure is also preserved, albeit for a narrower range of the smoothing parameter and, consequently, for a more modest improvement in sampling. We have also applied the new method in a "resolution exchange" setting, in which each replica runs a Monte Carlo simulation with a different degree of smoothing. We obtain exchange rates that compare favorably to our previous efforts at resolution exchange (Lyman E.; Zuckerman D. M.J. Chem. Theory Comput.2006, 2, 656-666).

DNA binding and bending by Sac7d is stepwise.

Which came first? Using thermodynamic cycles constructed from experimental measures of the binding free energy and free energy simulations, we have shown that binding of Sulfolobus acidocaldarius protein Sac7d to DNA occurs before DNA bending, thus indicating that a conformational selection mechanism is unlikely to be operative in Sac7d binding.

DNA Bending through Roll Angles Is Independent of Adjacent Base Pairs.

We have studied DNA bending for a wide range of DNA sequences by two-dimensional adaptive umbrella sampling simulations on adjacent roll angles. Calculated free energy surfaces are largely additive and can be well approximated by the sum of the one-dimensional free energy surfaces. Cooperativity between adjacent roll angles was found to be negligible: less than 1.0 kcal/mol and a small fraction of the overall bending energy. Our calculations validate the assumptions underlying many popular coarse-grained models for DNA bending, and demonstrate their theoretical validity for investigating DNA bending.

DNA Bending through Large Angles Is Aided by Ionic Screening.

We used adaptive umbrella sampling on a modified version of the roll angle to simulate the bending of DNA dodecamers. Simulations were carried out with the AMBER and CHARMM force fields for 10 sequences in which the central base pair step was varied. On long length scales, the DNA behavior was found to be consistent with the worm-like chain model. Persistence lengths calculated directly from the simulated structures and indirectly through the use of sequence-independent coarse-grained models based on simulation data were similar to literature values. On short length scales, the free energy cost of bending DNA was found to be consistent with the worm-like chain model for small and intermediate bending angles. At large angles, the bending free energy as a function of the roll angle became linear, suggesting a relative increase in flexibility at larger roll angles. Counterions congregated on the concave side of the highly bent DNA and screened the repulsion of the phosphate groups, facilitating the bending.

Cy3-DNA stacking interactions strongly depend on the identity of the terminal basepair.

We characterized the effect of the first basepair on the conformational dynamics of the fluorescent dye Cy3 attached to the 5' end of double-stranded DNA using gaussian-mixture adaptive umbrella sampling simulations. In the simulations, the sampling of all five dihedral angles along the linker was enhanced, so that both stacked and unstacked states were sampled. The affinity of Cy3 for a T·A basepair (with the dye attached to T) was found to be significantly less than for the other basepairs. This was verified experimentally by measuring the activation energies for cis-trans isomerization of the dye. The simulation and experimental results indicate the existence of partially unstacked conformations amenable to photoisomerization. The simulations also showed that stacking of Cy3 straightens the DNA while stabilizing the first basepair. Our findings indicate that fluorescence is modulated by Cy3-DNA interactions in a sequence-dependent manner.

Mechanism of the calcium-induced trans-cis isomerization of a non-prolyl peptide bond in Clostridium histolyticum collagenase.

cis peptide bonds in proteins are often rate-limiting steps in protein folding or conformational change and are frequently stabilized by metal ions. In the collagen-binding domain of Clostridium histolyticum collagenase, the binding of calcium ions triggers the formation of a cis peptide bond. We present free energy simulations of the formation of this cis peptide bond using a combined quantum mechanics/molecular mechanics approach together with adaptive umbrella sampling. From these simulations, we have determined that the calcium ions not only stabilize the cis peptide bond thermodynamically but also catalyze its formation; the free energy barrier to the formation of the cis peptide bond decreases from 21.4 kcal/mol in the absence of calcium ions to 10.3 kcal/mol in their presence. Two principal factors contribute to this reduction in the energy barrier. The calcium ions electrostatically stabilize the lone pair on the nitrogen atom that forms during the isomerization. In addition, their attraction to acidic amino acid side chains and formation of a hydrogen bond network constrain the peptide backbone in a way that makes it easier for the nitrogen to pyramidalize. Factors that explain the observed cooperativity of calcium binding are discussed.

Modulation of protein stability by O-glycosylation in a designed Gc-MAF analog.

The post-translational modification of proteins by the covalent attachment of carbohydrates to specific side chains, or glycosylation, is emerging as a crucial process in modulating the function of proteins. In particular, the dynamic processing of the oligosaccharide can correlate with a change in function. For example, a potent macrophage-activating factor, Gc-MAF, is obtained from serum vitamin D binding protein (VDBP) by stepwise processing of the oligosaccharide attached to Thr 420 to the core alpha-GalNAc moiety. In previous work we designed a miniprotein analog of Gc-MAF, MM1, by grafting the glycosylated loop of Gc-MAF on a stable scaffold. GalNAc-MM1 showed native-like activity on macrophages (Bogani 2006, J. Am. Chem. Soc. 128 7142-43). Here, we present data on the thermodynamic stability and conformational dynamics of the mono- and diglycosylated forms. We observed an unusual trend: each glycosylation event destabilized the protein by about 1 kcal/mol. This effect is matched by an increase in the mobility of the glycosylated forms, as evaluated by molecular dynamics simulations. An analysis of the solvent-accessible surface area shows that glycosylation causes the three-helix bundle to adopt conformations in which the hydrophobic residues are more solvent exposed. The number of hydrophobic contacts is also affected. These two factors, which are ultimately explained with a change in occupancy for conformers of specific side chains, may contribute to the observed destabilization.